From 5431451db37b6c1f73ee5b816fad0775927d11f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Thu, 9 Apr 2020 17:44:17 +0200 Subject: [PATCH] Block decomposition: refactor algorithms on the variable dependency graph In particular, move them into a separate class rather than a namespace. --- src/Makefile.am | 4 +- src/MinimumFeedbackSet.cc | 455 -------------------------------- src/MinimumFeedbackSet.hh | 60 ----- src/ModelTree.cc | 148 ++++------- src/VariableDependencyGraph.cc | 468 +++++++++++++++++++++++++++++++++ src/VariableDependencyGraph.hh | 90 +++++++ 6 files changed, 615 insertions(+), 610 deletions(-) delete mode 100644 src/MinimumFeedbackSet.cc delete mode 100644 src/MinimumFeedbackSet.hh create mode 100644 src/VariableDependencyGraph.cc create mode 100644 src/VariableDependencyGraph.hh diff --git a/src/Makefile.am b/src/Makefile.am index e880905a..5583a2ae 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -42,8 +42,8 @@ dynare_m_SOURCES = \ Statement.hh \ ExprNode.cc \ ExprNode.hh \ - MinimumFeedbackSet.cc \ - MinimumFeedbackSet.hh \ + VariableDependencyGraph.cc \ + VariableDependencyGraph.hh \ DynareMain.cc \ MacroExpandModFile.cc \ CodeInterpreter.hh \ diff --git a/src/MinimumFeedbackSet.cc b/src/MinimumFeedbackSet.cc deleted file mode 100644 index a162ce33..00000000 --- a/src/MinimumFeedbackSet.cc +++ /dev/null @@ -1,455 +0,0 @@ -/* - * Copyright © 2009-2020 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 . - */ - -#include - -#include "MinimumFeedbackSet.hh" - -using namespace boost; - -namespace MFS -{ - using color_t = map::vertex_descriptor, boost::default_color_type>; - using vector_vertex_descriptor_t = vector; - - 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 &circuit_stack) - { - property_map::type v_index = get(vertex_index, g); - color[u] = gray_color; - graph_traits::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 &circuit_stack, AdjacencyList_t &g) - { - // Initialize color map to white - color_t color; - graph_traits::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; - auto v_index = get(vertex_index, G); - cout << "Graph\n" - << "-----\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, int n) - { - AdjacencyList_t G(n); - auto v_index = get(vertex_index, G); - auto v_index1 = get(vertex_index1, G); - for (int i = 0; i < n; i++) - { - put(v_index, vertex(i, G), i); - put(v_index1, vertex(i, G), i); - } - for (int i = 0; i < n; i++) - for (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 select_index) - { - int n = select_index.size(); - AdjacencyList_t G(n); - auto v_index = get(vertex_index, G); - auto v_index1 = get(vertex_index1, G); - auto v1_index = get(vertex_index, G1); - map reverse_index; - set::iterator it; - 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; - auto 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 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 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; - int i = 1; - while (i < static_cast(liste.size()) && agree) - { - int j = i + 1; - while (j < static_cast(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 - auto 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 - auto 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 &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 - auto v_index = get(vertex_index, G); - cout << "store v[*it] = " << v_index[*it]+1 << "\n"; -#endif - auto v_index1 = get(vertex_index1, G); - feed_back_vertices.insert(v_index1[*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 &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 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*/ - 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 (static_cast(in_degree(*it, G) + out_degree(*it, G)) > max_degree) - { - max_degree = in_degree(*it, G) + out_degree(*it, G); - max_degree_index = it; - } - - auto v_index1 = get(vertex_index1, G); - feed_back_vertices.insert(v_index1[*max_degree_index]); - cut_++; -#ifdef verbose - auto 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 &feedback_vertices, vector< int> &Reordered_Vertices) - { - AdjacencyList_t G(G1); - auto v_index = get(vertex_index, G); - set fv; - for (auto its = feedback_vertices.begin(); its != feedback_vertices.end(); ++its) - fv.insert(*its); - int i = 0; - for (auto 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"; - } -} diff --git a/src/MinimumFeedbackSet.hh b/src/MinimumFeedbackSet.hh deleted file mode 100644 index 9e10ab86..00000000 --- a/src/MinimumFeedbackSet.hh +++ /dev/null @@ -1,60 +0,0 @@ -/* - * Copyright © 2009-2020 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 . - */ - -#ifndef _MINIMUMFEEDBACKSET_HH -#define _MINIMUMFEEDBACKSET_HH - -#include -#include - -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wold-style-cast" -#include -#pragma GCC diagnostic pop - -using namespace std; - -namespace MFS -{ - using VertexProperty_t = boost::property>>>>; - using AdjacencyList_t = boost::adjacency_list; - - //! Extracts a subgraph - /*! - \param[in] G1 The original graph - \param[in] select_index The vertex indices to select - \return The subgraph - - The property vertex_index of the subgraph contains indices of the original - graph, the property vertex_index1 contains new contiguous indices specific - to the subgraph. - */ - AdjacencyList_t extract_subgraph(AdjacencyList_t &G1, set select_index); - //! Return the feedback set - AdjacencyList_t Minimal_set_of_feedback_vertex(set &feed_back_vertices, const AdjacencyList_t &G); - //! Reorder the recursive variables - /*! They appear first in a quasi triangular form and they are followed by the feedback variables */ - void Reorder_the_recursive_variables(const AdjacencyList_t &G1, set &feedback_vertices, vector &Reordered_Vertices); -}; - -#endif // _MINIMUMFEEDBACKSET_HH diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 769665f9..fc8f8996 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -18,14 +18,13 @@ */ #include "ModelTree.hh" -#include "MinimumFeedbackSet.hh" +#include "VariableDependencyGraph.hh" #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wold-style-cast" #pragma GCC diagnostic ignored "-Wsign-compare" #pragma GCC diagnostic ignored "-Wmaybe-uninitialized" #include #include -#include #include #pragma GCC diagnostic pop @@ -740,93 +739,56 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob int nb_var = variable_reordered.size(); int n = nb_var - prologue - epilogue; - MFS::AdjacencyList_t G2(n); - - // It is necessary to manually initialize vertex_index property since this graph uses listS and not vecS as underlying vertex container - auto v_index = get(boost::vertex_index, G2); - for (int i = 0; i < n; i++) - put(v_index, vertex(i, G2), i); + /* Construct the graph representing the dependencies between all + variables that do not belong to the prologue or the epilogue. */ + VariableDependencyGraph G(n); vector reverse_equation_reordered(nb_var), reverse_variable_reordered(nb_var); - for (int i = 0; i < nb_var; i++) { reverse_equation_reordered[equation_reordered[i]] = i; reverse_variable_reordered[variable_reordered[i]] = i; } + jacob_map_t tmp_normalized_contemporaneous_jacobian; if (cutoff == 0) - { - set> endo; - for (int i = 0; i < nb_var; i++) - { - endo.clear(); - equations[i]->collectEndogenous(endo); - for (const auto &it : endo) - tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1; - } - } + for (int i = 0; i < nb_var; i++) + { + set> endo; + equations[i]->collectEndogenous(endo); + for (const auto &it : endo) + tmp_normalized_contemporaneous_jacobian[{ i, it.first }] = 1; + } else tmp_normalized_contemporaneous_jacobian = static_jacobian; for (const auto &[key, value] : tmp_normalized_contemporaneous_jacobian) if (reverse_equation_reordered[key.first] >= prologue && reverse_equation_reordered[key.first] < nb_var - epilogue && reverse_variable_reordered[key.second] >= prologue && reverse_variable_reordered[key.second] < nb_var - epilogue && key.first != endo2eq[key.second]) - add_edge(vertex(reverse_equation_reordered[endo2eq[key.second]]-prologue, G2), - vertex(reverse_equation_reordered[key.first]-prologue, G2), - G2); + add_edge(vertex(reverse_equation_reordered[endo2eq[key.second]]-prologue, G), + vertex(reverse_equation_reordered[key.first]-prologue, G), + G); - vector endo2block(num_vertices(G2)), discover_time(num_vertices(G2)); - boost::iterator_property_map::type, int, int &> endo2block_map(&endo2block[0], get(boost::vertex_index, G2)); + /* Compute the mapping between endogenous and blocks, using a strongly + connected components (SCC) decomposition */ + auto [num_scc, endo2block] = G.sortedStronglyConnectedComponents(); - // Compute strongly connected components - int num = strong_components(G2, endo2block_map); + vector> blocks(num_scc, { 0, 0 }); - vector> blocks(num, { 0, 0 }); - - // Create directed acyclic graph associated to the strongly connected components - using DirectedGraph = boost::adjacency_list; - DirectedGraph dag(num); - - for (int i = 0; i < static_cast(num_vertices(G2)); i++) - { - MFS::AdjacencyList_t::out_edge_iterator it_out, out_end; - MFS::AdjacencyList_t::vertex_descriptor vi = vertex(i, G2); - for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out) - { - int t_b = endo2block_map[target(*it_out, G2)]; - int s_b = endo2block_map[source(*it_out, G2)]; - if (s_b != t_b) - add_edge(s_b, t_b, dag); - } - } - - // Compute topological sort of DAG (ordered list of unordered SCC) - deque ordered2unordered; - topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order - - // Construct mapping from unordered SCC to ordered SCC - vector unordered2ordered(num); - for (int i = 0; i < num; i++) - unordered2ordered[ordered2unordered[i]] = i; - - //This vector contains for each block: - // - first set = equations belonging to the block, - // - second set = the feeback variables, - // - third vector = the reordered non-feedback variables. - vector, set, vector>> components_set(num); + /* This vector contains for each block: + – first set = equations belonging to the block, + — second set = the feeback variables, + — third vector = the reordered non-feedback variables. */ + vector, set, vector>> components_set(num_scc); for (int i = 0; i < static_cast(endo2block.size()); i++) { - endo2block[i] = unordered2ordered[endo2block[i]]; blocks[endo2block[i]].first++; get<0>(components_set[endo2block[i]]).insert(i); } - auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block, num); + auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock(endo2block, num_scc); - vector tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered); - int order = prologue; - //Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set + // Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set if (select_feedback_variable) { for (int i = 0; i < n; i++) @@ -836,17 +798,18 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob || equation_lag_lead[equation_reordered[i+prologue]].second > 0 || equation_lag_lead[equation_reordered[i+prologue]].first > 0 || mfs == 0) - add_edge(vertex(i, G2), vertex(i, G2), G2); + add_edge(vertex(i, G), vertex(i, G), G); } else for (int i = 0; i < n; i++) if (Equation_Type[equation_reordered[i+prologue]].first == EquationType::solve || mfs == 0) - add_edge(vertex(i, G2), vertex(i, G2), G2); + add_edge(vertex(i, G), vertex(i, G), G); - //Determines the dynamic structure of each equation - vector n_static(prologue+num+epilogue, 0), n_forward(prologue+num+epilogue, 0), - n_backward(prologue+num+epilogue, 0), n_mixed(prologue+num+epilogue, 0); + // Determines the dynamic structure of each equation + vector n_static(prologue+num_scc+epilogue, 0), n_forward(prologue+num_scc+epilogue, 0), + n_backward(prologue+num_scc+epilogue, 0), n_mixed(prologue+num_scc+epilogue, 0); + vector tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered); for (int i = 0; i < prologue; i++) if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0) n_mixed[i]++; @@ -857,24 +820,23 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0) n_static[i]++; - //For each block, the minimum set of feedback variable is computed - // and the non-feedback variables are reordered to get - // a sub-recursive block without feedback variables + /* For each block, the minimum set of feedback variable is computed and the + non-feedback variables are reordered to get a sub-recursive block without + feedback variables. */ - for (int i = 0; i < num; i++) + int order = prologue; + for (int i = 0; i < num_scc; i++) { - MFS::AdjacencyList_t G = MFS::extract_subgraph(G2, get<0>(components_set[i])); - set feed_back_vertices; - MFS::AdjacencyList_t G1 = MFS::Minimal_set_of_feedback_vertex(feed_back_vertices, G); - auto v_index = get(boost::vertex_index, G); + auto subG = G.extractSubgraph(get<0>(components_set[i])); + auto [G1, feed_back_vertices] = subG.minimalSetOfFeedbackVertices(); + auto v_index1 = get(boost::vertex_index1, subG); get<1>(components_set[i]) = feed_back_vertices; blocks[i].second = feed_back_vertices.size(); - vector Reordered_Vertice; - MFS::Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice); + auto reordered_vertices = subG.reorderRecursiveVariables(feed_back_vertices); - //First we have the recursive equations conditional on feedback variables + // First we have the recursive equations conditional on feedback variables for (int j = 0; j < 4; j++) - for (int its : Reordered_Vertice) + for (int its : reordered_vertices) { bool something_done = false; if (j == 2 && variable_lag_lead[tmp_variable_reordered[its+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its+prologue]].second != 0) @@ -905,36 +867,36 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob } } - get<2>(components_set[i]) = Reordered_Vertice; - //Second we have the equations related to the feedback variables + get<2>(components_set[i]) = reordered_vertices; + // Second we have the equations related to the feedback variables for (int j = 0; j < 4; j++) for (int feed_back_vertice : feed_back_vertices) { bool something_done = false; - if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0) + if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second != 0) { n_mixed[prologue+i]++; something_done = true; } - else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second != 0) + else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second != 0) { n_forward[prologue+i]++; something_done = true; } - else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0) + else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second == 0) { n_backward[prologue+i]++; something_done = true; } - else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]].second == 0) + else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]].second == 0) { n_static[prologue+i]++; something_done = true; } if (something_done) { - equation_reordered[order] = tmp_equation_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]; - variable_reordered[order] = tmp_variable_reordered[v_index[vertex(feed_back_vertice, G)]+prologue]; + equation_reordered[order] = tmp_equation_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]; + variable_reordered[order] = tmp_variable_reordered[v_index1[vertex(feed_back_vertice, subG)]+prologue]; order++; } } @@ -942,13 +904,13 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob for (int i = 0; i < epilogue; i++) if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0) - n_mixed[prologue+num+i]++; + n_mixed[prologue+num_scc+i]++; else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0) - n_forward[prologue+num+i]++; + n_forward[prologue+num_scc+i]++; else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0) - n_backward[prologue+num+i]++; + n_backward[prologue+num_scc+i]++; else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0) - n_static[prologue+num+i]++; + n_static[prologue+num_scc+i]++; inv_equation_reordered.resize(nb_var); inv_variable_reordered.resize(nb_var); diff --git a/src/VariableDependencyGraph.cc b/src/VariableDependencyGraph.cc new file mode 100644 index 00000000..f0187704 --- /dev/null +++ b/src/VariableDependencyGraph.cc @@ -0,0 +1,468 @@ +/* + * Copyright © 2009-2020 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 . + */ + +#include +#include + +#include "VariableDependencyGraph.hh" + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wold-style-cast" +#pragma GCC diagnostic ignored "-Wsign-compare" +#pragma GCC diagnostic ignored "-Wmaybe-uninitialized" +#include +#include +#pragma GCC diagnostic pop + +using namespace boost; + +VariableDependencyGraph::VariableDependencyGraph(int n) : base(n) +{ + /* It is necessary to manually initialize the vertex_index property since + this graph uses listS and not vecS as underlying vertex container */ + auto v_index = get(vertex_index, *this); + for (int i = 0; i < n; i++) + put(v_index, vertex(i, *this), i); +} + +void +VariableDependencyGraph::suppress(vertex_descriptor vertex_to_eliminate) +{ + clear_vertex(vertex_to_eliminate, *this); + remove_vertex(vertex_to_eliminate, *this); +} + +void +VariableDependencyGraph::suppress(int vertex_num) +{ + suppress(vertex(vertex_num, *this)); +} + +void +VariableDependencyGraph::eliminate(vertex_descriptor vertex_to_eliminate) +{ + if (in_degree(vertex_to_eliminate, *this) > 0 && out_degree(vertex_to_eliminate, *this) > 0) + { + in_edge_iterator it_in, in_end; + out_edge_iterator it_out, out_end; + for (tie(it_in, in_end) = in_edges(vertex_to_eliminate, *this); it_in != in_end; ++it_in) + for (tie(it_out, out_end) = out_edges(vertex_to_eliminate, *this); it_out != out_end; ++it_out) + { + edge_descriptor ed; + bool exist; + tie(ed, exist) = edge(source(*it_in, *this), target(*it_out, *this), *this); + if (!exist) + add_edge(source(*it_in, *this), target(*it_out, *this), *this); + } + } + suppress(vertex_to_eliminate); +} + +bool +VariableDependencyGraph::hasCycleDFS(vertex_descriptor u, color_t &color, vector &circuit_stack) const +{ + auto v_index = get(vertex_index, *this); + color[u] = gray_color; + out_edge_iterator vi, vi_end; + for (tie(vi, vi_end) = out_edges(u, *this); vi != vi_end; ++vi) + if (color[target(*vi, *this)] == white_color && hasCycleDFS(target(*vi, *this), color, circuit_stack)) + { + // cycle detected, return immediately + circuit_stack.push_back(v_index[target(*vi, *this)]); + return true; + } + else if (color[target(*vi, *this)] == gray_color) + { + // *vi is an ancestor! + circuit_stack.push_back(v_index[target(*vi, *this)]); + return true; + } + color[u] = black_color; + return false; +} + +bool +VariableDependencyGraph::hasCycle() const +{ + // Initialize color map to white + color_t color; + vector circuit_stack; + vertex_iterator vi, vi_end; + for (tie(vi, vi_end) = vertices(*this); vi != vi_end; ++vi) + color[*vi] = white_color; + + // Perform depth-first search + for (tie(vi, vi_end) = vertices(*this); vi != vi_end; ++vi) + if (color[*vi] == white_color && hasCycleDFS(*vi, color, circuit_stack)) + return true; + + return false; +} + +void +VariableDependencyGraph::print() const +{ + vertex_iterator it, it_end; + auto v_index = get(vertex_index, *this); + cout << "Graph\n" + << "-----\n"; + for (tie(it, it_end) = vertices(*this); it != it_end; ++it) + { + cout << "vertex[" << v_index[*it] + 1 << "] <-"; + in_edge_iterator it_in, in_end; + for (tie(it_in, in_end) = in_edges(*it, *this); it_in != in_end; ++it_in) + cout << v_index[source(*it_in, *this)] + 1 << " "; + cout << "\n ->"; + out_edge_iterator it_out, out_end; + for (tie(it_out, out_end) = out_edges(*it, *this); it_out != out_end; ++it_out) + cout << v_index[target(*it_out, *this)] + 1 << " "; + cout << "\n"; + } +} + +VariableDependencyGraph +VariableDependencyGraph::extractSubgraph(const set &select_index) const +{ + int n = select_index.size(); + VariableDependencyGraph G(n); + auto v_index = get(vertex_index, *this); + auto v_index1_G = get(vertex_index1, G); // Maps new vertices to original indices + map reverse_index; // Maps orig indices to new ones + set::iterator it; + int i; + for (it = select_index.begin(), i = 0; i < n; i++, ++it) + { + reverse_index[get(v_index, vertex(*it, *this))] = i; + put(v_index1_G, vertex(i, G), get(v_index, vertex(*it, *this))); + } + for (it = select_index.begin(), i = 0; i < n; i++, ++it) + { + out_edge_iterator it_out, out_end; + auto vi = vertex(*it, *this); + for (tie(it_out, out_end) = out_edges(vi, *this); it_out != out_end; ++it_out) + { + int ii = v_index[target(*it_out, *this)]; + if (select_index.find(ii) != select_index.end()) + add_edge(vertex(reverse_index[get(v_index, source(*it_out, *this))], G), + vertex(reverse_index[get(v_index, target(*it_out, *this))], G), G); + } + } + return G; +} + +bool +VariableDependencyGraph::vertexBelongsToAClique(vertex_descriptor vertex) const +{ + vector liste; + bool agree = true; + in_edge_iterator it_in, in_end; + out_edge_iterator it_out, out_end; + tie(it_in, in_end) = in_edges(vertex, *this); + tie(it_out, out_end) = out_edges(vertex, *this); + while (it_in != in_end && it_out != out_end && agree) + { + agree = (source(*it_in, *this) == target(*it_out, *this) + && source(*it_in, *this) != target(*it_in, *this)); //not a loop + liste.push_back(source(*it_in, *this)); + ++it_in; + ++it_out; + } + if (agree) + { + if (it_in != in_end || it_out != out_end) + agree = false; + int i = 1; + while (i < static_cast(liste.size()) && agree) + { + int j = i + 1; + while (j < static_cast(liste.size()) && agree) + { + edge_descriptor ed; + bool exist1, exist2; + tie(ed, exist1) = edge(liste[i], liste[j], *this); + tie(ed, exist2) = edge(liste[j], liste[i], *this); + agree = (exist1 && exist2); + j++; + } + i++; + } + } + return agree; +} + +bool +VariableDependencyGraph::eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree() +{ + bool something_has_been_done = false; + bool not_a_loop; + int i; + vertex_iterator it, it1, ita, it_end; + for (tie(it, it_end) = vertices(*this), i = 0; it != it_end; ++it, i++) + { + int in_degree_n = in_degree(*it, *this); + int out_degree_n = out_degree(*it, *this); + 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! + { + in_edge_iterator it_in, in_end; + for (tie(it_in, in_end) = in_edges(*it, *this); it_in != in_end; ++it_in) + if (source(*it_in, *this) == target(*it_in, *this)) + { +#ifdef verbose + cout << v_index[source(*it_in, *this)] << " == " << v_index[target(*it_in, *this)] << "\n"; +#endif + not_a_loop = false; + } + } + if (not_a_loop) + { +#ifdef verbose + auto v_index1 = get(vertex_index1, *this); + cout << "->eliminate vertex[" << v_index1[*it] + 1 << "]\n"; +#endif + eliminate(*it); +#ifdef verbose + print(); +#endif + something_has_been_done = true; + if (i > 0) + it = ita; + else + { + tie(it, it_end) = vertices(*this); + i--; + } + } + } + ita = it; + } + return something_has_been_done; +} + +bool +VariableDependencyGraph::eliminationOfVerticesBelongingToAClique() +{ + vertex_iterator it, it1, ita, it_end; + bool something_has_been_done = false; + int i; + for (tie(it, it_end) = vertices(*this), i = 0; it != it_end; ++it, i++) + { + if (vertexBelongsToAClique(*it)) + { +#ifdef verbose + auto v_index1 = get(vertex_index1, *this); + cout << "eliminate vertex[" << v_index1[*it] + 1 << "]\n"; +#endif + eliminate(*it); + something_has_been_done = true; + if (i > 0) + it = ita; + else + { + tie(it, it_end) = vertices(*this); + i--; + } + } + ita = it; + } + return something_has_been_done; +} + +bool +VariableDependencyGraph::suppressionOfVerticesWithLoop(set &feed_back_vertices) +{ + bool something_has_been_done = false; + vertex_iterator it, it_end, ita; + int i = 0; + for (tie(it, it_end) = vertices(*this); it != it_end; ++it, i++) + { + edge_descriptor ed; + bool exist; + tie(ed, exist) = edge(*it, *it, *this); + if (exist) + { +#ifdef verbose + auto v_index1 = get(vertex_index1, *this); + cout << "store v[*it] = " << v_index1[*it]+1 << "\n"; +#endif + auto v_index = get(vertex_index, *this); + feed_back_vertices.insert(v_index[*it]); + suppress(*it); + something_has_been_done = true; + if (i > 0) + it = ita; + else + { + tie(it, it_end) = vertices(*this); + i--; + } + } + ita = it; + } + return something_has_been_done; +} + +pair> +VariableDependencyGraph::minimalSetOfFeedbackVertices() const +{ + bool something_has_been_done = true; + int cut_ = 0; + set feed_back_vertices; + VariableDependencyGraph G(*this); + while (num_vertices(G) > 0) + { + while (something_has_been_done && num_vertices(G) > 0) + { + //Rule 1 + something_has_been_done = G.eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree() /*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 = (G.eliminationOfVerticesBelongingToAClique() || something_has_been_done); +#ifdef verbose + cout << "2 something_has_been_done=" << something_has_been_done << "\n"; +#endif + + //Rule 3 + something_has_been_done = (G.suppressionOfVerticesWithLoop(feed_back_vertices) || something_has_been_done); +#ifdef verbose + cout << "3 something_has_been_done=" << something_has_been_done << "\n"; +#endif + } + if (!G.hasCycle()) + { +#ifdef verbose + cout << "has_cycle=false\n"; +#endif + return { G, feed_back_vertices }; + } + 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 */ + int max_degree = 0, num = 0; + vertex_iterator it, it_end, max_degree_index; + for (tie(it, it_end) = vertices(G); it != it_end; ++it, num++) + if (static_cast(in_degree(*it, G) + out_degree(*it, G)) > max_degree) + { + max_degree = in_degree(*it, G) + out_degree(*it, G); + max_degree_index = it; + } + + auto v_index = get(vertex_index, G); + feed_back_vertices.insert(v_index[*max_degree_index]); + cut_++; +#ifdef verbose + auto v_index1 = get(vertex_index1, G); + cout << "--> cut vertex " << v_index1[*max_degree_index] + 1 << "\n"; +#endif + G.suppress(*max_degree_index); + something_has_been_done = true; + } + } +#ifdef verbose + cout << "cut_=" << cut_ << "\n"; +#endif + return { G, feed_back_vertices }; +} + +vector +VariableDependencyGraph::reorderRecursiveVariables(const set &feedback_vertices) const +{ + vector reordered_vertices; + VariableDependencyGraph G(*this); + auto v_index = get(vertex_index, G); + set> fv; + for (auto its = feedback_vertices.begin(); its != feedback_vertices.end(); ++its) + fv.insert(*its); + int i = 0; + for (auto its = fv.begin(); its != fv.end(); ++its, i++) + G.suppress(*its); + bool something_has_been_done = true; + while (something_has_been_done) + { + something_has_been_done = false; + 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]); + G.suppress(*it); + 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"; + return reordered_vertices; +} + +pair> +VariableDependencyGraph::sortedStronglyConnectedComponents() const +{ + vector vertex2scc(num_vertices(*this)); + auto v_index = get(vertex_index, *this); + + // Compute SCCs and create mapping from vertices to unordered SCCs + int num_scc = strong_components(static_cast(*this), make_iterator_property_map(vertex2scc.begin(), v_index)); + + // Create directed acyclic graph (DAG) associated to the SCCs + adjacency_list dag(num_scc); + for (int i = 0; i < static_cast(num_vertices(*this)); i++) + { + out_edge_iterator it_out, out_end; + auto vi = vertex(i, *this); + for (tie(it_out, out_end) = out_edges(vi, *this); it_out != out_end; ++it_out) + { + int t_b = vertex2scc[v_index[target(*it_out, *this)]]; + int s_b = vertex2scc[v_index[source(*it_out, *this)]]; + if (s_b != t_b) + add_edge(s_b, t_b, dag); + } + } + + /* Compute topological sort of DAG (ordered list of unordered SCC) + Note: the order is reversed. */ + vector reverseOrdered2unordered; + topological_sort(dag, back_inserter(reverseOrdered2unordered)); + + // Construct mapping from unordered SCC to ordered SCC + vector unordered2ordered(num_scc); + for (int j = 0; j < num_scc; j++) + unordered2ordered[reverseOrdered2unordered[num_scc-j-1]] = j; + + // Update the mapping of vertices to (now sorted) SCCs + for (int i = 0; i < static_cast(num_vertices(*this)); i++) + vertex2scc[i] = unordered2ordered[vertex2scc[i]]; + + return { num_scc, vertex2scc }; +} diff --git a/src/VariableDependencyGraph.hh b/src/VariableDependencyGraph.hh new file mode 100644 index 00000000..6974c020 --- /dev/null +++ b/src/VariableDependencyGraph.hh @@ -0,0 +1,90 @@ +/* + * Copyright © 2009-2020 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 . + */ + +#ifndef _VARIABLEDEPENDENCYGRAPH_HH +#define _VARIABLEDEPENDENCYGRAPH_HH + +#include +#include + +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wold-style-cast" +#include +#pragma GCC diagnostic pop + +using namespace std; + +using VertexProperty_t = boost::property>>>>; + +/* Class used to store a graph representing dependencies between variables. + Used in the block decomposition. */ +class VariableDependencyGraph + : public boost::adjacency_list +{ +public: + using color_t = map::vertex_descriptor, boost::default_color_type>; + using base = boost::adjacency_list; + + VariableDependencyGraph(int n); + //! Extracts a subgraph + /*! + \param[in] select_index The vertex indices to select + \return The subgraph + + The property vertex_index1 of the subgraph contains indices of the original + graph. + */ + VariableDependencyGraph extractSubgraph(const set &select_index) const; + //! Return the feedback set + pair> minimalSetOfFeedbackVertices() const; + //! Reorder the recursive variables + /*! They appear first in a quasi triangular form and they are followed by the feedback variables */ + vector reorderRecursiveVariables(const set &feedback_vertices) const; + /* Computes the strongly connected components (SCCs) of the graph, and sort them + topologically. + Returns the number of SCCs, and a mapping of vertex indices to sorted SCC + indices. */ + pair> sortedStronglyConnectedComponents() const; +private: + // Remove a vertex (including all edges to and from it); takes a vertex descriptor + void suppress(vertex_descriptor vertex_to_eliminate); + // Remove a vertex (including all edges to and from it); takes a vertex index + void suppress(int vertex_num); + /* Remove a vertex, but keeping the paths that go through it (i.e. by adding + edges that directly connect vertices that would otherwise be connected + through the vertex to be removed) */ + void eliminate(vertex_descriptor vertex_to_eliminate); + // Internal helper for hasCycle() + bool hasCycleDFS(vertex_descriptor u, color_t &color, vector &circuit_stack) const; + // Determine whether the graph has a cycle + bool hasCycle() const; + // Print on stdout a description of the graph + void print() const; + bool vertexBelongsToAClique(vertex_descriptor vertex) const; + bool eliminationOfVerticesWithOneOrLessIndegreeOrOutdegree(); + bool eliminationOfVerticesBelongingToAClique(); + // The suppressed vertices are stored in feedback set + bool suppressionOfVerticesWithLoop(set &feed_back_vertices); +}; + +#endif // _VARIABLEDEPENDENCYGRAPH_HH