Skip to content
Snippets Groups Projects
Select Git revision
  • bcad4f31f24098bc875acac00ac8b07888d8c08f
  • master default protected
  • 6.x protected
  • madysson
  • 5.x protected
  • asm
  • time-varying-information-set
  • 4.6 protected
  • dynare_minreal
  • dragonfly
  • various_fixes
  • 4.5 protected
  • clang+openmp
  • exo_steady_state
  • declare_vars_in_model_block
  • julia
  • error_msg_undeclared_model_vars
  • static_aux_vars
  • slice
  • aux_func
  • penalty
  • 6.4 protected
  • 6.3 protected
  • 6.2 protected
  • 6.1 protected
  • 6.0 protected
  • 6-beta2 protected
  • 6-beta1 protected
  • 5.5 protected
  • 5.4 protected
  • 5.3 protected
  • 5.2 protected
  • 5.1 protected
  • 5.0 protected
  • 5.0-rc1 protected
  • 4.7-beta3 protected
  • 4.7-beta2 protected
  • 4.7-beta1 protected
  • 4.6.4 protected
  • 4.6.3 protected
  • 4.6.2 protected
41 results

MinimumFeedbackSet.cc

Blame
  • MinimumFeedbackSet.cc 16.13 KiB
    /*
     * Copyright (C) 2009-2012 Dynare Team
     *
     * This file is part of Dynare.
     *
     * Dynare is free software: you can redistribute it and/or modify
     * it under the terms of the GNU General Public License as published by
     * the Free Software Foundation, either version 3 of the License, or
     * (at your option) any later version.
     *
     * Dynare is distributed in the hope that it will be useful,
     * but WITHOUT ANY WARRANTY; without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     * GNU General Public License for more details.
     *
     * You should have received a copy of the GNU General Public License
     * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
     */
    
    #include <iostream>
    
    #include "MinimumFeedbackSet.hh"
    
    namespace MFS
    {
      void
      Suppress(AdjacencyList_t::vertex_descriptor vertex_to_eliminate, AdjacencyList_t &G)
      {
        clear_vertex(vertex_to_eliminate, G);
        remove_vertex(vertex_to_eliminate, G);
      }
    
      void
      Suppress(int vertex_num, AdjacencyList_t &G)
      {
        Suppress(vertex(vertex_num, G), G);
      }
    
      void
      Eliminate(AdjacencyList_t::vertex_descriptor vertex_to_eliminate, AdjacencyList_t &G)
      {
        if (in_degree(vertex_to_eliminate, G) > 0 && out_degree(vertex_to_eliminate, G) > 0)
          {
            AdjacencyList_t::in_edge_iterator it_in, in_end;
            AdjacencyList_t::out_edge_iterator it_out, out_end;
            for (tie(it_in, in_end) = in_edges(vertex_to_eliminate, G); it_in != in_end; ++it_in)
              for (tie(it_out, out_end) = out_edges(vertex_to_eliminate, G); it_out != out_end; ++it_out)
                {
                  AdjacencyList_t::edge_descriptor ed;
                  bool exist;
                  tie(ed, exist) = edge(source(*it_in, G), target(*it_out, G), G);
                  if (!exist)
                    add_edge(source(*it_in, G), target(*it_out, G), G);
                }
          }
        Suppress(vertex_to_eliminate, G);
      }
    
      bool
      has_cycle_dfs(AdjacencyList_t &g, AdjacencyList_t::vertex_descriptor u, color_t &color, vector<int> &circuit_stack)
      {
        property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, g);
        color[u] = gray_color;
        graph_traits<AdjacencyList_t>::out_edge_iterator vi, vi_end;
        for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi)
          if (color[target(*vi, g)] == white_color && has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
            {
              // cycle detected, return immediately
              circuit_stack.push_back(v_index[target(*vi, g)]);
              return true;
            }
          else if (color[target(*vi, g)] == gray_color)
            {
              // *vi is an ancestor!
              circuit_stack.push_back(v_index[target(*vi, g)]);
              return true;
            }
        color[u] = black_color;
        return false;
      }
    
      bool
      has_cycle(vector<int> &circuit_stack, AdjacencyList_t &g)
      {
        // Initialize color map to white
        color_t color;
        graph_traits<AdjacencyList_t>::vertex_iterator vi, vi_end;
        for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
          color[*vi] = white_color;
    
        // Perform depth-first search
        for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
          if (color[*vi] == white_color && has_cycle_dfs(g, *vi, color, circuit_stack))
            return true;
    
        return false;
      }
    
      void
      Print(AdjacencyList_t &G)
      {
        AdjacencyList_t::vertex_iterator  it, it_end;
        property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
        cout << "Graph\n";
        cout << "-----\n";
        for (tie(it, it_end) = vertices(G); it != it_end; ++it)
          {
            cout << "vertex[" << v_index[*it] + 1 << "] <-";
            AdjacencyList_t::in_edge_iterator it_in, in_end;
            for (tie(it_in, in_end) = in_edges(*it, G); it_in != in_end; ++it_in)
              cout << v_index[source(*it_in, G)] + 1 << " ";
            cout << "\n       ->";
            AdjacencyList_t::out_edge_iterator it_out, out_end;
            for (tie(it_out, out_end) = out_edges(*it, G); it_out != out_end; ++it_out)
              cout << v_index[target(*it_out, G)] + 1 << " ";
            cout << "\n";
          }
      }
    
      AdjacencyList_t
      AM_2_AdjacencyList(bool *AM, unsigned int n)
      {
        AdjacencyList_t G(n);
        property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
        property_map<AdjacencyList_t, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
        for (unsigned int i = 0; i < n; i++)
          {
            put(v_index, vertex(i, G), i);
            put(v_index1, vertex(i, G), i);
          }
        for (unsigned int i = 0; i < n; i++)
          for (unsigned int j = 0; j < n; j++)
            if (AM[i*n+j])
              add_edge(vertex(j, G), vertex(i, G), G);
        return G;
      }
    
      AdjacencyList_t
      extract_subgraph(AdjacencyList_t &G1, set<int> select_index)
      {
        unsigned int n = select_index.size();
        AdjacencyList_t G(n);
        property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
        property_map<AdjacencyList_t, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
        property_map<AdjacencyList_t, vertex_index_t>::type v1_index = get(vertex_index, G1);
        map<int, int> reverse_index;
        set<int>::iterator it;
        unsigned int i;
        for (it = select_index.begin(), i = 0; i < n; i++, ++it)
          {
            reverse_index[get(v1_index, vertex(*it, G1))] = i;
            put(v_index, vertex(i, G), get(v1_index, vertex(*it, G1)));
            put(v_index1, vertex(i, G), i);
          }
        for (it = select_index.begin(), i = 0; i < n; i++, ++it)
          {
            AdjacencyList_t::out_edge_iterator it_out, out_end;
            AdjacencyList_t::vertex_descriptor vi = vertex(*it, G1);
            for (tie(it_out, out_end) = out_edges(vi, G1); it_out != out_end; ++it_out)
              {
                int ii = v1_index[target(*it_out, G1)];
                if (select_index.find(ii) != select_index.end())
                  add_edge(vertex(reverse_index[get(v1_index, source(*it_out, G1))], G), vertex(reverse_index[get(v1_index, target(*it_out, G1))], G), G);
              }
          }
        return G;
      }
    
      vector_vertex_descriptor_t
      Collect_Doublet(AdjacencyList_t::vertex_descriptor vertex, AdjacencyList_t &G)
      {
        AdjacencyList_t::in_edge_iterator it_in, in_end;
        AdjacencyList_t::out_edge_iterator it_out, out_end;
        vector<AdjacencyList_t::vertex_descriptor> Doublet;
        if (in_degree(vertex, G) > 0 && out_degree(vertex, G) > 0)
          for (tie(it_in, in_end) = in_edges(vertex, G); it_in != in_end; ++it_in)
            for (tie(it_out, out_end) = out_edges(vertex, G); it_out != out_end; ++it_out)
              if (source(*it_in, G) == target(*it_out, G) && source(*it_in, G) != target(*it_in, G))                                                                                                                                                                                                                                                                                                       // not a loop
                Doublet.push_back(source(*it_in, G));
        return Doublet;
      }
    
      bool
      Vertex_Belong_to_a_Clique(AdjacencyList_t::vertex_descriptor vertex, AdjacencyList_t &G)
      {
        vector<AdjacencyList_t::vertex_descriptor> liste;
        bool agree = true;
        AdjacencyList_t::in_edge_iterator it_in, in_end;
        AdjacencyList_t::out_edge_iterator it_out, out_end;
        tie(it_in, in_end) = in_edges(vertex, G);
        tie(it_out, out_end) = out_edges(vertex, G);
        while (it_in != in_end && it_out != out_end && agree)
          {
            agree = (source(*it_in, G) == target(*it_out, G) && source(*it_in, G) != target(*it_in, G));  //not a loop
            liste.push_back(source(*it_in, G));
            ++it_in;
            ++it_out;
          }
        if (agree)
          {
            if (it_in != in_end || it_out != out_end)
              agree = false;
            unsigned int i = 1;
            while (i < liste.size() && agree)
              {
                unsigned int j = i + 1;
                while (j < liste.size() && agree)
                  {
                    AdjacencyList_t::edge_descriptor ed;
                    bool exist1, exist2;
                    tie(ed, exist1) = edge(liste[i], liste[j], G);
                    tie(ed, exist2) = edge(liste[j], liste[i], G);
                    agree = (exist1 && exist2);
                    j++;
                  }
                i++;
              }
          }
        return agree;
      }
    
      bool
      Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_t &G)
      {
        bool something_has_been_done = false;
        bool not_a_loop;
        int i;
        AdjacencyList_t::vertex_iterator  it, it1, ita, it_end;
        for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
          {
            int in_degree_n = in_degree(*it, G);
            int out_degree_n = out_degree(*it, G);
            if (in_degree_n <= 1 || out_degree_n <= 1)
              {
                not_a_loop = true;
                if (in_degree_n >= 1 && out_degree_n >= 1) // Do not eliminate a vertex if it loops on itself!
                  {
                    AdjacencyList_t::in_edge_iterator it_in, in_end;
                    for (tie(it_in, in_end) = in_edges(*it, G); it_in != in_end; ++it_in)
                      if (source(*it_in, G) == target(*it_in, G))
                        {
    #ifdef verbose
                          cout << v_index[source(*it_in, G)] << " == " << v_index[target(*it_in, G)] << "\n";
    #endif
                          not_a_loop = false;
                        }
                  }
                if (not_a_loop)
                  {
    #ifdef verbose
                    property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
                    cout << "->eliminate vertex[" << v_index[*it] + 1 << "]\n";
    #endif
                    Eliminate(*it, G);
    #ifdef verbose
                    Print(G);
    #endif
                    something_has_been_done = true;
                    if (i > 0)
                      it = ita;
                    else
                      {
                        tie(it, it_end) = vertices(G);
                        i--;
                      }
                  }
              }
            ita = it;
          }
        return something_has_been_done;
      }
    
      bool
      Elimination_of_Vertex_belonging_to_a_clique_Step(AdjacencyList_t &G)
      {
        AdjacencyList_t::vertex_iterator  it, it1, ita, it_end;
        bool something_has_been_done = false;
        int i;
        for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
          {
            if (Vertex_Belong_to_a_Clique(*it, G))
              {
    #ifdef verbose
                property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
                cout << "eliminate vertex[" << v_index[*it] + 1 << "]\n";
    #endif
                Eliminate(*it, G);
                something_has_been_done = true;
                if (i > 0)
                  it = ita;
                else
                  {
                    tie(it, it_end) = vertices(G);
                    i--;
                  }
              }
            ita = it;
          }
        return something_has_been_done;
      }
    
      bool
      Suppression_of_Vertex_X_if_it_loops_store_in_set_of_feedback_vertex_Step(set<int> &feed_back_vertices, AdjacencyList_t &G)
      {
        bool something_has_been_done = false;
        AdjacencyList_t::vertex_iterator  it, it_end, ita;
        int i = 0;
        for (tie(it, it_end) = vertices(G); it != it_end; ++it, i++)
          {
            AdjacencyList_t::edge_descriptor ed;
            bool exist;
            tie(ed, exist) = edge(*it, *it, G);
            if (exist)
              {
    #ifdef verbose
                property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
                cout << "store v[*it] = " << v_index[*it]+1 << "\n";
    #endif
                property_map<AdjacencyList_t, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
                feed_back_vertices.insert(v_index1[*it]);
                /*property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
                  feed_back_vertices.insert(v_index[*it] );*/
                Suppress(*it, G);
                something_has_been_done = true;
                if (i > 0)
                  it = ita;
                else
                  {
                    tie(it, it_end) = vertices(G);
                    i--;
                  }
              }
            ita = it;
          }
        return something_has_been_done;
      }
    
      AdjacencyList_t
      Minimal_set_of_feedback_vertex(set<int> &feed_back_vertices, const AdjacencyList_t &G1)
      {
        bool something_has_been_done = true;
        int cut_ = 0;
        feed_back_vertices.clear();
        AdjacencyList_t G(G1);
        while (num_vertices(G) > 0)
          {
            while (something_has_been_done && num_vertices(G) > 0)
              {
                //Rule 1
                something_has_been_done = (Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(G) /*or something_has_been_done*/);
    #ifdef verbose
                cout << "1 something_has_been_done=" << something_has_been_done << "\n";
    #endif
    
                //Rule 2
                something_has_been_done = (Elimination_of_Vertex_belonging_to_a_clique_Step(G) || something_has_been_done);
    #ifdef verbose
                cout << "2 something_has_been_done=" << something_has_been_done << "\n";
    #endif
    
                //Rule 3
                something_has_been_done = (Suppression_of_Vertex_X_if_it_loops_store_in_set_of_feedback_vertex_Step(feed_back_vertices, G) || something_has_been_done);
    #ifdef verbose
                cout << "3 something_has_been_done=" << something_has_been_done << "\n";
    #endif
              }
            vector<int> circuit;
            if (!has_cycle(circuit, G))
              {
    #ifdef verbose
                cout << "has_cycle=false\n";
    #endif
                //sort(feed_back_vertices.begin(), feed_back_vertices.end());
                return G;
              }
            if (num_vertices(G) > 0)
              {
                /*if nothing has been done in the five previous rule then cut the vertex with the maximum in_degree+out_degree*/
                unsigned int max_degree = 0, num = 0;
                AdjacencyList_t::vertex_iterator  it, it_end, max_degree_index;
                for (tie(it, it_end) = vertices(G); it != it_end; ++it, num++)
                  {
                    if (in_degree(*it, G) + out_degree(*it, G) > max_degree)
                      {
                        max_degree = in_degree(*it, G) + out_degree(*it, G);
                        max_degree_index = it;
                      }
                  }
                property_map<AdjacencyList_t, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
                feed_back_vertices.insert(v_index1[*max_degree_index]);
                /*property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
                  feed_back_vertices.insert(v_index[*max_degree_index]);*/
                //cout << "v_index1[*max_degree_index] = " << v_index1[*max_degree_index] << "\n";
                cut_++;
    #ifdef verbose
                property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
                cout << "--> cut vertex " << v_index[*max_degree_index] + 1 << "\n";
    #endif
                Suppress(*max_degree_index, G);
                something_has_been_done = true;
              }
          }
    #ifdef verbose
        cout << "cut_=" << cut_ << "\n";
    #endif
        //sort(feed_back_vertices.begin(), feed_back_vertices.end());
        return G;
      }
    
      struct rev
      {
        bool
        operator()(const int a, const int b) const
        {
          return (a > b);
        }
      };
    
      void
      Reorder_the_recursive_variables(const AdjacencyList_t &G1, set<int> &feedback_vertices, vector< int> &Reordered_Vertices)
      {
        AdjacencyList_t G(G1);
        property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
        set<int>::iterator its, ita;
        set<int, rev> fv;
        for (its = feedback_vertices.begin(); its != feedback_vertices.end(); its++)
          fv.insert(*its);
        int i = 0;
        for (its = fv.begin(); its != fv.end(); ++its, i++)
          Suppress(*its, G);
        bool something_has_been_done = true;
        while (something_has_been_done)
          {
            something_has_been_done = false;
            AdjacencyList_t::vertex_iterator  it, it_end, ita;
            for (tie(it, it_end) = vertices(G), i = 0; it != it_end; ++it, i++)
              {
                if (in_degree(*it, G) == 0)
                  {
                    Reordered_Vertices.push_back(v_index[*it]);
                    Suppress(*it, G);
                    something_has_been_done = true;
                    if (i > 0)
                      it = ita;
                    else
                      {
                        tie(it, it_end) = vertices(G);
                        i--;
                      }
                  }
                ita = it;
              }
          }
        if (num_vertices(G))
          cout << "Error in the computation of feedback vertex set\n";
      }
    }