stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Block decomposition: refactor algorithms on the variable dependency graph 

In particular, move them into a separate class rather than a namespace.
issue#70
Sébastien Villemot 2020-04-09 17:44:17 +02:00
parent c38f9aa7c3
commit 5431451db3
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 615 additions and 610 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/Makefile.am
View File

@ -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 \

455
src/MinimumFeedbackSet.cc
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include "MinimumFeedbackSet.hh"
using namespace boost;
namespace MFS
{
using color_t = map<boost::graph_traits<AdjacencyList_t>::vertex_descriptor, boost::default_color_type>;
using vector_vertex_descriptor_t = vector<AdjacencyList_t::vertex_descriptor>;
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;
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<int> 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<int, int> reverse_index;
set<int>::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<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;
int i = 1;
while (i < static_cast<int>(liste.size()) && agree)
{
int j = i + 1;
while (j < static_cast<int>(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<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
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<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*/
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<int>(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<int> &feedback_vertices, vector< int> &Reordered_Vertices)
{
AdjacencyList_t G(G1);
auto v_index = get(vertex_index, G);
set<int, rev> 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";
}
}

60
src/MinimumFeedbackSet.hh
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef _MINIMUMFEEDBACKSET_HH
#define _MINIMUMFEEDBACKSET_HH
#include <map>
#include <vector>
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wold-style-cast"
#include <boost/graph/adjacency_list.hpp>
#pragma GCC diagnostic pop
using namespace std;
namespace MFS
{
using VertexProperty_t = boost::property<boost::vertex_index_t, int,
boost::property<boost::vertex_index1_t, int,
boost::property<boost::vertex_degree_t, int,
boost::property<boost::vertex_in_degree_t, int,
boost::property<boost::vertex_out_degree_t, int >>>>>;
using AdjacencyList_t = boost::adjacency_list<boost::listS, boost::listS, boost::bidirectionalS, VertexProperty_t>;
//! 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<int> select_index);
//! Return the feedback set
AdjacencyList_t Minimal_set_of_feedback_vertex(set<int> &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<int> &feedback_vertices, vector<int> &Reordered_Vertices);
};
#endif // _MINIMUMFEEDBACKSET_HH

148
src/ModelTree.cc
View File

@ -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 <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>
#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<int> 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<pair<int, int>> 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<pair<int, int>> 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<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
boost::iterator_property_map<int *, boost::property_map<MFS::AdjacencyList_t, boost::vertex_index_t>::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<pair<int, int>> blocks(num_scc, { 0, 0 });
vector<pair<int, int>> blocks(num, { 0, 0 });
// Create directed acyclic graph associated to the strongly connected components
using DirectedGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS>;
DirectedGraph dag(num);
for (int i = 0; i < static_cast<int>(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<int> 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<int> 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<tuple<set<int>, set<int>, vector<int>>> 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<tuple<set<int>, set<int>, vector<int>>> components_set(num_scc);
for (int i = 0; i < static_cast<int>(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<int> 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<int> 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<int> 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<int> 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<int> 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<int> 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);

468
src/VariableDependencyGraph.cc Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <algorithm>
#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 <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>
#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<int> &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<int> 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<int> &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<int, int> reverse_index; // Maps orig indices to new ones
set<int>::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<vertex_descriptor> 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<int>(liste.size()) && agree)
{
int j = i + 1;
while (j < static_cast<int>(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<int> &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, set<int>>
VariableDependencyGraph::minimalSetOfFeedbackVertices() const
{
bool something_has_been_done = true;
int cut_ = 0;
set<int> 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<int>(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<int>
VariableDependencyGraph::reorderRecursiveVariables(const set<int> &feedback_vertices) const
{
vector<int> reordered_vertices;
VariableDependencyGraph G(*this);
auto v_index = get(vertex_index, G);
set<int, greater<int>> 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<int, vector<int>>
VariableDependencyGraph::sortedStronglyConnectedComponents() const
{
vector<int> 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<base>(*this), make_iterator_property_map(vertex2scc.begin(), v_index));
// Create directed acyclic graph (DAG) associated to the SCCs
adjacency_list<vecS, vecS, directedS> dag(num_scc);
for (int i = 0; i < static_cast<int>(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<int> reverseOrdered2unordered;
topological_sort(dag, back_inserter(reverseOrdered2unordered));
// Construct mapping from unordered SCC to ordered SCC
vector<int> 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<int>(num_vertices(*this)); i++)
vertex2scc[i] = unordered2ordered[vertex2scc[i]];
return { num_scc, vertex2scc };
}

90
src/VariableDependencyGraph.hh Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef _VARIABLEDEPENDENCYGRAPH_HH
#define _VARIABLEDEPENDENCYGRAPH_HH
#include <map>
#include <vector>
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wold-style-cast"
#include <boost/graph/adjacency_list.hpp>
#pragma GCC diagnostic pop
using namespace std;
using VertexProperty_t = boost::property<boost::vertex_index_t, int,
boost::property<boost::vertex_index1_t, int,
boost::property<boost::vertex_degree_t, int,
boost::property<boost::vertex_in_degree_t, int,
boost::property<boost::vertex_out_degree_t, int>>>>>;
/* Class used to store a graph representing dependencies between variables.
Used in the block decomposition. */
class VariableDependencyGraph
: public boost::adjacency_list<boost::listS, boost::listS, boost::bidirectionalS, VertexProperty_t>
{
public:
using color_t = map<boost::graph_traits<VariableDependencyGraph>::vertex_descriptor, boost::default_color_type>;
using base = boost::adjacency_list<boost::listS, boost::listS, boost::bidirectionalS, VertexProperty_t>;
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<int> &select_index) const;
//! Return the feedback set
pair<VariableDependencyGraph, set<int>> minimalSetOfFeedbackVertices() const;
//! Reorder the recursive variables
/*! They appear first in a quasi triangular form and they are followed by the feedback variables */
vector<int> reorderRecursiveVariables(const set<int> &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<int, vector<int>> 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<int> &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<int> &feed_back_vertices);
};
#endif // _VARIABLEDEPENDENCYGRAPH_HH