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

New implementation of block decomposition & feedback variables using Boost for DynamicModel 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2671 ac1d8469-bf42-47a9-8791-bf33cf982152
issue#70
ferhat 2009-05-12 23:03:40 +00:00
parent f01306a29b
commit 916e950df9
13 changed files with 1052 additions and 276 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

528
BlockTriangular.cc
View File

@ -25,15 +25,26 @@
#include <cstdlib>
#include <cstring>
#include <cmath>
using namespace std;
#include <algorithm>
#include <functional>
#include "MinimumFeedbackSet.hh"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
//------------------------------------------------------------------------------
#include "BlockTriangular.hh"
//------------------------------------------------------------------------------
using namespace std;
using namespace boost;
using namespace MFS;
BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) :
symbol_table(symbol_table_arg),
normalization(symbol_table_arg),
incidencematrix(symbol_table_arg)
symbol_table(symbol_table_arg),
normalization(symbol_table_arg),
incidencematrix(symbol_table_arg)
{
bt_verbose = 0;
ModelBlock = NULL;
@ -45,19 +56,19 @@ BlockTriangular::BlockTriangular(const SymbolTable &symbol_table_arg) :
//------------------------------------------------------------------------------
// Find the prologue and the epilogue of the model
void
BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0)
BlockTriangular::Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM0)
{
bool modifie = 1;
int i, j, k, l = 0;
/*Looking for a prologue */
*prologue = 0;
prologue = 0;
while (modifie)
{
modifie = 0;
for (i = *prologue;i < n;i++)
for (i = prologue;i < n;i++)
{
k = 0;
for (j = *prologue;j < n;j++)
for (j = prologue;j < n;j++)
{
if (IM[i*n + j])
{
@ -65,23 +76,23 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
l = j;
}
}
if ((k == 1) && IM0[Index_Equ_IM[i].index*n + Index_Var_IM[l].index])
if ((k == 1) && IM0[Index_Equ_IM[i]*n + Index_Var_IM[l]])
{
modifie = 1;
incidencematrix.swap_IM_c(IM, *prologue, i, l, Index_Var_IM, Index_Equ_IM, n);
(*prologue)++;
incidencematrix.swap_IM_c(IM, prologue, i, l, Index_Var_IM, Index_Equ_IM, n);
prologue++;
}
}
}
*epilogue = 0;
epilogue = 0;
modifie = 1;
while (modifie)
{
modifie = 0;
for (i = *prologue;i < n - *epilogue;i++)
for (i = prologue;i < n - epilogue;i++)
{
k = 0;
for (j = *prologue;j < n - *epilogue;j++)
for (j = prologue;j < n - epilogue;j++)
{
if (IM[j*n + i])
{
@ -89,17 +100,148 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n
l = j;
}
}
if ((k == 1) && IM0[Index_Equ_IM[l].index*n + Index_Var_IM[i].index])
if ((k == 1) && IM0[Index_Equ_IM[l]*n + Index_Var_IM[i]])
{
modifie = 1;
incidencematrix.swap_IM_c(IM, n - (1 + *epilogue), l, i, Index_Var_IM, Index_Equ_IM, n);
(*epilogue)++;
incidencematrix.swap_IM_c(IM, n - (1 + epilogue), l, i, Index_Var_IM, Index_Equ_IM, n);
epilogue++;
}
}
}
}
//------------------------------------------------------------------------------
// Find a matching between equations and endogenous variables
bool
BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector<int> &Index_Equ_IM) const
{
int n = equation_number - prologue - epilogue;
typedef adjacency_list<vecS, vecS, undirectedS> BipartiteGraph;
/*
Vertices 0 to n-1 are for endogenous (using type specific ID)
Vertices n to 2*n-1 are for equations (using equation no.)
*/
BipartiteGraph g(2 * n);
// Fill in the graph
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
if (IM0[(i+prologue) * equation_number+j+prologue])
add_edge(i + n, j, g);
// Compute maximum cardinality matching
typedef vector<graph_traits<BipartiteGraph>::vertex_descriptor> mate_map_t;
mate_map_t mate_map(2*n);
bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
//cout << "check = " << check << "\n";
if (check)
{
// Check if all variables are normalized
mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
if (it != mate_map.begin() + n)
{
if (verbose)
{
cerr << "ERROR: Could not normalize static model. Variable "
<< symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
<< " is not in the maximum cardinality matching." << endl;
exit(EXIT_FAILURE);
}
return false;
}
for (int i = 0; i < n; i++)
{
int j = Index_Equ_IM[mate_map[i]-n+prologue];
Index_Equ_IM[mate_map[i]-n+prologue] = Index_Equ_IM[i+prologue];
Index_Equ_IM[i+prologue] = j;
for (int k=0; k<n; k++)
{
j = IM[(mate_map[i]-n+prologue)*equation_number+k +prologue];
IM[(mate_map[i]-n+prologue)*equation_number+k +prologue] = IM[(i+prologue)*equation_number+k +prologue];
IM[(i+prologue)*equation_number+k +prologue] = j;
}
}
}
return check;
}
void
BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, bool verbose_) const
{
int n = nb_var - prologue - epilogue;
bool *AMp;
AMp = (bool*) malloc(n*n*sizeof(bool));
//transforms the incidence matrix of the complet model into an adjancent matrix of the non-recursive part of the model
for (int i=prologue; i<nb_var - epilogue; i++)
for (int j=prologue; j<nb_var - epilogue; j++)
if (Index_Equ_IM[i]!=Index_Var_IM[j])
AMp[(i-prologue)*n+j-prologue] = IM[i*nb_var + j];
else
AMp[(i-prologue)*n+j-prologue] = 0;
//In a first step we compute the strong components of the graph representation of the static model.
// This insures that block are dynamically recursives.
GraphvizDigraph G2 = AM_2_GraphvizDigraph(AMp, n);
vector<int> component(num_vertices(G2)), discover_time(num_vertices(G2));
//write_graphviz(cout, G2);
int num = strong_components(G2, &component[0]);
blocks = vector<pair<int, int> >(num, make_pair(0,0));
//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<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
for (unsigned int i=0; i<component.size(); i++)
{
blocks[component[i]].first++;
components_set[component[i]].first.insert(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
vector<int> tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
int order = prologue;
bool *SIM;
SIM = (bool*)malloc(nb_var*nb_var*sizeof(bool));
memcpy(SIM, IM, nb_var*nb_var*sizeof(bool));
for (int i=0; i<num; i++)
{
AdjacencyList_type G = GraphvizDigraph_2_AdjacencyList(G2, components_set[i].first);
set<int> feed_back_vertices;
AdjacencyList_type G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
components_set[i].second.first = feed_back_vertices;
blocks[i].second = feed_back_vertices.size();
vector<int> Reordered_Vertice;
Reordered_Vertice = Reorder_the_recursive_variables(G, feed_back_vertices);
for (vector<int>::iterator its=Reordered_Vertice.begin();its != Reordered_Vertice.end(); its++)
{
Index_Equ_IM[order] = tmp_Index_Equ_IM[*its+prologue];
Index_Var_IM[order] = tmp_Index_Var_IM[*its+prologue];
order++;
}
components_set[i].second.second = Reordered_Vertice;
for (set<int>::iterator its=feed_back_vertices.begin();its != feed_back_vertices.end(); its++)
{
Index_Equ_IM[order] = tmp_Index_Equ_IM[v_index[vertex(*its, G)]+prologue];
Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue];
order++;
}
}
}
void
BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock)
{
@ -145,9 +287,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
{
ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]=new temporary_terms_type ();
ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]->clear();
ModelBlock->Block_List[count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index;
ModelBlock->Block_List[count_Block].Variable[i] = Index_Var_IM[*count_Equ].index;
i_1 = Index_Var_IM[*count_Equ].index;
ModelBlock->Block_List[count_Block].Equation[i] = Index_Equ_IM[*count_Equ];
ModelBlock->Block_List[count_Block].Variable[i] = Index_Var_IM[*count_Equ];
i_1 = Index_Var_IM[*count_Equ];
for (k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
@ -158,7 +300,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
{
for (j = 0;j < size;j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr()])
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j]*symbol_table.endo_nbr()])
{
tmp_variable_evaluated[i_1] = true;
tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
@ -177,7 +319,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
{
for (j = 0;j < size;j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr()])
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j]*symbol_table.endo_nbr()])
{
tmp_variable_evaluated[i_1] = true;
tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
@ -208,7 +350,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if (Cur_IM)
{
i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.endo_nbr();
i_1 = Index_Equ_IM[first_count_equ+i] * symbol_table.endo_nbr();
for (j = 0;j < symbol_table.endo_nbr();j++)
if (Cur_IM[i_1 + j])
{
@ -244,7 +386,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
Cur_IM = incidencematrix.Get_IM(k, eExogenous);
if (Cur_IM)
{
i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.exo_nbr();
i_1 = Index_Equ_IM[first_count_equ+i] * symbol_table.exo_nbr();
for (j=0;j<symbol_table.exo_nbr();j++)
if (Cur_IM[i_1 + j])
{
@ -331,10 +473,10 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
{
for (j = first_count_equ;j < size + first_count_equ;j++)
{
i_1 = Index_Var_IM[j].index;
i_1 = Index_Var_IM[j];
m = 0;
for (k = first_count_equ;k < size + first_count_equ;k++)
if (IM[i_1 + Index_Equ_IM[k].index*symbol_table.endo_nbr()])
if (IM[i_1 + Index_Equ_IM[k] * symbol_table.endo_nbr()])
m++;
if (m > 0)
{
@ -345,9 +487,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
m = 0;
for (j = first_count_equ;j < size + first_count_equ;j++)
{
i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr();
i_1 = Index_Equ_IM[j] * symbol_table.endo_nbr();
for (k = first_count_equ;k < size + first_count_equ;k++)
if (IM[Index_Var_IM[k].index + i_1])
if (IM[Index_Var_IM[k] + i_1])
{
if (i == Lag)
{
@ -357,9 +499,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].IM_lead_lag[i].u[m] = li;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index;
tmp_variable_evaluated[Index_Var_IM[k].index] = true;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j];
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k];
tmp_variable_evaluated[Index_Var_IM[k]] = true;
l++;
m++;
li++;
@ -369,15 +511,15 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
m = 0;
for (j = first_count_equ;j < size + first_count_equ;j++)
{
i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr();
i_1 = Index_Equ_IM[j] * symbol_table.endo_nbr();
for (k = 0;k < symbol_table.endo_nbr();k++)
if ((!tmp_variable_evaluated[Index_Var_IM[k].index]) && IM[Index_Var_IM[k].index + i_1])
if ((!tmp_variable_evaluated[Index_Var_IM[k]]) && IM[Index_Var_IM[k] + i_1])
{
ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_other_endo[m] = l;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_other_endo[m] = j - first_count_equ;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_other_endo[m] = k;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index_other_endo[m] = Index_Equ_IM[j].index;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index_other_endo[m] = Index_Var_IM[k].index;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index_other_endo[m] = Index_Equ_IM[j];
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index_other_endo[m] = Index_Var_IM[k];
l++;
m++;
}
@ -390,7 +532,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
m = 0;
for (j = first_count_equ;j < size + first_count_equ;j++)
{
i_1 = Index_Equ_IM[j].index * symbol_table.exo_nbr();
i_1 = Index_Equ_IM[j] * symbol_table.exo_nbr();
for (k = 0; k<tmp_nb_exo; k++)
{
if (IM[ModelBlock->Block_List[count_Block].Exogenous[k]+i_1])
@ -398,7 +540,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Exogenous[m] = k;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Exogenous_Index[m] = ModelBlock->Block_List[count_Block].Exogenous[k];
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X[m] = j - first_count_equ;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X_Index[m] = Index_Equ_IM[j].index;
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X_Index[m] = Index_Equ_IM[j];
m++;
}
}
@ -418,56 +560,54 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
void
BlockTriangular::Free_Block(Model_Block* ModelBlock) const
{
int blk, i;
for (blk = 0;blk < ModelBlock->Size;blk++)
{
{
int blk, i;
for (blk = 0;blk < ModelBlock->Size;blk++)
{
free(ModelBlock->Block_List[blk].Equation);
free(ModelBlock->Block_List[blk].Variable);
free(ModelBlock->Block_List[blk].Exogenous);
free(ModelBlock->Block_List[blk].Own_Derivative);
free(ModelBlock->Block_List[blk].Other_Endogenous);
for (i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++)
{
if (incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size*/)
{
free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].u_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index_other_endo);
}
if (incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo*/)
{
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X);
}
}
free(ModelBlock->Block_List[blk].IM_lead_lag);
for(i=0; i<ModelBlock->Block_List[blk].Size; i++)
delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
delete(ModelBlock->Block_List[blk].Temporary_InUse);
}
free(ModelBlock->Block_List);
free(ModelBlock);
free(Index_Equ_IM);
free(Index_Var_IM);
}
free(ModelBlock->Block_List[blk].Equation);
free(ModelBlock->Block_List[blk].Variable);
free(ModelBlock->Block_List[blk].Exogenous);
free(ModelBlock->Block_List[blk].Own_Derivative);
free(ModelBlock->Block_List[blk].Other_Endogenous);
for (i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++)
{
if (incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size*/)
{
free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].u_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index_other_endo);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index_other_endo);
}
if (incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size_exo*/)
{
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index);
free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X);
}
}
free(ModelBlock->Block_List[blk].IM_lead_lag);
for (i=0; i<ModelBlock->Block_List[blk].Size; i++)
delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
delete(ModelBlock->Block_List[blk].Temporary_InUse);
}
free(ModelBlock->Block_List);
free(ModelBlock);
}
t_type
BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, block_result_t* res, vector<BinaryOpNode *> equations )
BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> equations )
{
int i=0;
NodeID lhs, rhs;
@ -481,23 +621,23 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
t_type Type;
bool *Cur_IM;
BlockSimulationType Simulation_Type , prev_Type=UNKNOWN;
for ( i=0; i<prologue+res->n_sets+epilogue; i++)
for ( i=0; i<prologue+(int)blocks.size()+epilogue; i++)
{
int first_count_equ = count_equ;
if (i < prologue)
Blck_Size = 1;
else if (i < prologue+res->n_sets)
else if (i < prologue+(int)blocks.size())
{
Blck_Size = res->sets_f[res->ordered[blck_count_simult]] - res->sets_s[res->ordered[blck_count_simult]] + 1;
Blck_Size = blocks[blck_count_simult].first;
blck_count_simult++;
}
else if (i < prologue+res->n_sets+epilogue)
else if (i < prologue+(int)blocks.size()+epilogue)
Blck_Size = 1;
Lag = Lead = 0;
for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
{
int i_1 = Index_Var_IM[count_equ].index;
int i_1 = Index_Var_IM[count_equ];
for (int k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
@ -505,7 +645,7 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
{
for (int j = 0;j < Blck_Size;j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr()])
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j] * symbol_table.endo_nbr()])
{
if (k > Lead)
Lead = k;
@ -540,13 +680,13 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
if (Blck_Size == 1)
{
temporary_terms.clear();
eq_node = equations[Index_Equ_IM[first_count_equ].index];
eq_node = equations[Index_Equ_IM[first_count_equ]];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
tmp_s.str("");
tmp_output.str("");
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
tmp_s << "y(it_, " << Index_Var_IM[first_count_equ].index+1 << ")";
tmp_s << "y(it_, " << Index_Var_IM[first_count_equ]+1 << ")";
//Determine whether the equation could be evaluated rather than to be solved
if (tmp_output.str()==tmp_s.str())
{
@ -571,7 +711,7 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
{
if ( ((prev_Type == EVALUATE_FORWARD_R || prev_Type == EVALUATE_FORWARD) && (Simulation_Type == EVALUATE_FORWARD_R || Simulation_Type == EVALUATE_FORWARD))
|| ((prev_Type == EVALUATE_BACKWARD_R || prev_Type == EVALUATE_BACKWARD) && (Simulation_Type == EVALUATE_BACKWARD_R || Simulation_Type == EVALUATE_BACKWARD))
)
)
{
BlockSimulationType c_Type = (Type[Type.size()-1]).first;
int c_Size = (Type[Type.size()-1]).second;
@ -593,150 +733,107 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
}
//------------------------------------------------------------------------------
// Normalize each equation of the model (endgenous_i = f_i(endogenous_1, ..., endogenous_n) - in order to apply strong connex components search algorithm -
// and find the optimal blocks triangular decomposition
bool
BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0, jacob_map j_m, vector<BinaryOpNode *> equations )
void
BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM_0, jacob_map j_m, vector<BinaryOpNode *> equations )
{
int i, j, Nb_TotalBlocks, Nb_RecursBlocks, Nb_SimulBlocks;
BlockType Btype;
int count_Block, count_Equ;
block_result_t* res;
Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set));
//block_result_t* res;
//Equation_set* Equation_gr = (Equation_set*) malloc(sizeof(Equation_set));
bool* SIM0, *SIM00;
SIM0 = (bool*)malloc(n * n * sizeof(bool));
memcpy(SIM0,IM_0,n*n*sizeof(bool));
Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0);
free(SIM0);
if (bt_verbose)
{
cout << "prologue : " << *prologue << " epilogue : " << *epilogue << "\n";
cout << "IM_0\n";
incidencematrix.Print_SIM(IM_0, eEndogenous);
cout << "IM\n";
incidencematrix.Print_SIM(IM, eEndogenous);
for (i = 0;i < n;i++)
cout << "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index << " Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i].index << "\n";
}
int counted=0;
if (*prologue+*epilogue<n)
if (prologue+epilogue<n)
{
if (Do_Normalization)
cout << "Normalizing the model ...\n";
double* max_val=(double*)malloc(n*sizeof(double));
memset(max_val,0,n*sizeof(double));
for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
{
cout << "Normalizing the model ...\n";
if (mixing)
{
double* max_val=(double*)malloc(n*sizeof(double));
memset(max_val,0,n*sizeof(double));
for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
{
if (fabs(iter->second)>max_val[iter->first.first])
max_val[iter->first.first]=fabs(iter->second);
}
for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
iter->second/=max_val[iter->first.first];
free(max_val);
bool OK=false;
double bi=0.99999999;
int suppressed=0;
while (!OK && bi>1e-14)
{
int suppress=0;
SIM0 = (bool*)malloc(n * n * sizeof(bool));
memset(SIM0,0,n*n*sizeof(bool));
SIM00 = (bool*)malloc(n * n * sizeof(bool));
memset(SIM00,0,n*n*sizeof(bool));
for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
{
if (fabs(iter->second)>bi)
{
SIM0[iter->first.first*n+iter->first.second]=1;
if (!IM_0[iter->first.first*n+iter->first.second])
{
cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << " " << iter->second << "\n";
}
}
else
suppress++;
}
for (i = 0;i < n;i++)
for (j = 0;j < n;j++)
{
SIM00[i*n + j] = SIM0[Index_Equ_IM[i].index * n + Index_Var_IM[j].index];
}
free(SIM0);
if (suppress!=suppressed)
{
OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
if (!OK)
normalization.Free_Equation(n-*prologue-*epilogue,Equation_gr);
}
suppressed=suppress;
if (!OK)
//bi/=1.07;
bi/=3;
counted++;
if (bi>1e-14)
free(SIM00);
}
if (!OK)
{
normalization.Set_fp_verbose(true);
OK=normalization.Normalize(n, *prologue, *epilogue, SIM00, Index_Equ_IM, Equation_gr, 1, IM);
cout << "Error\n";
exit(EXIT_FAILURE);
}
}
else
normalization.Normalize(n, *prologue, *epilogue, IM, Index_Equ_IM, Equation_gr, 0, 0);
if (fabs(iter->second)>max_val[iter->first.first])
max_val[iter->first.first]=fabs(iter->second);
}
else
normalization.Gr_to_IM_basic(n, *prologue, *epilogue, IM, Equation_gr, false);
for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
iter->second/=max_val[iter->first.first];
free(max_val);
bool OK=false;
double bi=0.99999999;
int suppressed=0;
while (!OK && bi>1e-14)
{
int suppress=0;
SIM0 = (bool*)malloc(n * n * sizeof(bool));
memset(SIM0,0,n*n*sizeof(bool));
SIM00 = (bool*)malloc(n * n * sizeof(bool));
memset(SIM00,0,n*n*sizeof(bool));
for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
{
if (fabs(iter->second)>bi)
{
SIM0[iter->first.first*n+iter->first.second]=1;
if (!IM_0[iter->first.first*n+iter->first.second])
{
cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << " " << iter->second << "\n";
}
}
else
suppress++;
}
for (i = 0;i < n;i++)
for (j = 0;j < n;j++)
{
SIM00[i*n + j] = SIM0[Index_Equ_IM[i] * n + Index_Var_IM[j]];
}
free(SIM0);
if (suppress!=suppressed)
OK = Compute_Normalization(IM, n, prologue, epilogue, false, SIM00, Index_Equ_IM);
suppressed=suppress;
if (!OK)
//bi/=1.07;
bi/=3;
counted++;
if (bi>1e-14)
free(SIM00);
}
if (!OK)
Compute_Normalization(IM, n, prologue, epilogue, true, SIM00, Index_Equ_IM);
}
//cout << "Finding the optimal block decomposition of the model ..." << counted << "\n";
cout << "Finding the optimal block decomposition of the model ...\n";
if (*prologue+*epilogue<n)
{
if (bt_verbose)
blocks.Print_Equation_gr(Equation_gr);
res = blocks.sc(Equation_gr);
normalization.Free_Equation(n-*prologue-*epilogue,Equation_gr);
if (bt_verbose)
blocks.block_result_print(res);
}
else
{
res = (block_result_t*)malloc(sizeof(*res));
res->n_sets=0;
}
free(Equation_gr);
blocks.block_result_to_IM(res, IM, *prologue, symbol_table.endo_nbr(), Index_Equ_IM, Index_Var_IM);
t_type Type = Reduce_Blocks_and_type_determination(*prologue, *epilogue, res, equations);
vector<pair<int, int> > blocks;
if (prologue+epilogue<n)
Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, false);
t_type Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations);
i = 0;
j = 0;
Nb_SimulBlocks = 0;
int Nb_fv = 0;
for (t_type::const_iterator it = Type.begin(); it!=Type.end(); it++)
{
if (it->first==SOLVE_FORWARD_COMPLETE || it->first==SOLVE_BACKWARD_COMPLETE || it->first==SOLVE_TWO_BOUNDARIES_COMPLETE)
{
Nb_SimulBlocks++;
if (it->second>j)
j=it->second;
{
j=it->second;
Nb_fv = blocks[Nb_SimulBlocks-1].second;
}
}
}
Nb_TotalBlocks = Type.size();
Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
cout << Nb_TotalBlocks << " block(s) found:\n";
cout << " " << Nb_RecursBlocks << " recursive block(s) and " << res->n_sets << " simultaneous block(s). \n";
cout << " the largest simultaneous block has " << j << " equation(s). \n";
cout << " " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n";
cout << " the largest simultaneous block has " << j << " equation(s)\n" <<
" and " << Nb_fv << " feedback variable(s).\n";
ModelBlock->Size = Nb_TotalBlocks;
ModelBlock->Periods = periods;
@ -745,9 +842,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
count_Equ = count_Block = 0;
for (t_type::const_iterator it = Type.begin(); it!=Type.end(); it++)
{
if (count_Equ<*prologue)
if (count_Equ<prologue)
Btype = PROLOGUE;
else if (count_Equ<n-*epilogue)
else if (count_Equ<n-epilogue)
if (it->second==1)
Btype = PROLOGUE;
else
@ -756,17 +853,9 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
Btype = EPILOGUE;
Allocate_Block(it->second, &count_Equ, count_Block++, Btype, it->first, ModelBlock);
}
//exit(-1);
if (res->n_sets)
blocks.block_result_free(res);
else
free(res);
return 0;
}
//------------------------------------------------------------------------------
// normalize each equation of the dynamic model
// and find the optimal block triangular decomposition of the static model
@ -796,15 +885,16 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_
cout << "incidence matrix for the static model (unsorted) \n";
incidencematrix.Print_SIM(SIM, eEndogenous);
}
Index_Equ_IM = (simple*)malloc(symbol_table.endo_nbr() * sizeof(*Index_Equ_IM));
//Index_Equ_IM = (simple*)malloc(symbol_table.endo_nbr() * sizeof(*Index_Equ_IM));
Index_Equ_IM = vector<int>(symbol_table.endo_nbr());
for (i = 0;i < symbol_table.endo_nbr();i++)
{
Index_Equ_IM[i].index = i;
Index_Equ_IM[i] = i;
}
Index_Var_IM = (simple*)malloc(symbol_table.endo_nbr() * sizeof(*Index_Var_IM));
Index_Var_IM = vector<int>(symbol_table.endo_nbr());
for (i = 0;i < symbol_table.endo_nbr();i++)
{
Index_Var_IM[i].index = i;
Index_Var_IM[i] = i;
}
if (ModelBlock != NULL)
Free_Block(ModelBlock);
@ -813,7 +903,7 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_
SIM_0 = (bool*)malloc(symbol_table.endo_nbr() * symbol_table.endo_nbr() * sizeof(*SIM_0));
for (i = 0;i < symbol_table.endo_nbr()*symbol_table.endo_nbr();i++)
SIM_0[i] = Cur_IM[i];
Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr(), &prologue, &epilogue, Index_Var_IM, Index_Equ_IM, 1, 1, SIM_0, j_m, equations);
Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr(), prologue, epilogue, Index_Var_IM, Index_Equ_IM, SIM_0, j_m, equations);
free(SIM_0);
free(SIM);
}

28
BlockTriangular.hh
View File

@ -34,31 +34,39 @@
//! Matrix of doubles for representing jacobian
//! Sparse matrix of double to store the values of the Jacobian
typedef map<pair<int ,int >,double> jacob_map;
typedef vector<pair<BlockSimulationType, int> > t_type;
//! Create the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
//! Creates the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
class BlockTriangular
{
//friend class IncidenceMatrix;
private:
//! Find equations and endogenous variables belonging to the prologue and epilogue of the model
void Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM0);
//! Allocates and fills the Model structure describing the content of each block
void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock);
//! Finds a matching between equations and endogenous variables
bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector<int> &Index_Var_IM) const;
//! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables
void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, bool verbose_) const;
//! Tries to merge the consecutive blocks in a single block and determine the type of each block: recursive, simultaneous, ...
t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> equations );
public:
const SymbolTable &symbol_table;
BlockTriangular(const SymbolTable &symbol_table_arg);
//! Frees the Model structure describing the content of each block
void Free_Block(Model_Block* ModelBlock) const;
//BlockTriangular(const IncidenceMatrix &incidence_matrix_arg);
//const SymbolTable &symbol_table;
Blocks blocks;
Normalization normalization;
IncidenceMatrix incidencematrix;
void Normalize_and_BlockDecompose_Static_0_Model(const jacob_map &j_m, vector<BinaryOpNode *> equations);
bool Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int* prologue, int* epilogue, simple* Index_Var_IM, simple* Index_Equ_IM, bool Do_Normalization, bool mixing, bool* IM_0 , jacob_map j_m, vector<BinaryOpNode *> equations);
void Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n, simple* Index_Var_IM, simple* Index_Equ_IM, bool* IM0);
void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock);
void Free_Block(Model_Block* ModelBlock) const;
t_type Reduce_Blocks_and_type_determination(int prologue, int epilogue, block_result_t* res, vector<BinaryOpNode *> equations );
simple *Index_Equ_IM;
simple *Index_Var_IM;
void Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM_0 , jacob_map j_m, vector<BinaryOpNode *> equations);
vector<int> Index_Equ_IM;
vector<int> Index_Var_IM;
int prologue, epilogue;
bool bt_verbose;
//int endo_nbr, exo_nbr;

20
DynamicModel.cc
View File

@ -21,6 +21,7 @@
#include <cstdlib>
#include <cassert>
#include "DynamicModel.hh"
// For mkdir() and chdir()
@ -124,7 +125,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++)
@ -134,7 +135,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eExogenous, var), lag)));
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
}*/
//jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -172,7 +173,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;i++)
@ -183,7 +184,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
//it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
}*/
//jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -1736,7 +1737,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
<< endl
<< jacobian_output.str()
<< "end" << endl;
if (second_derivatives.size())
{
// Writing initialization instruction for matrix g2
@ -1780,7 +1781,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput) const
<< " {" << endl
<< jacobian_output.str()
<< " }" << endl;
if (second_derivatives.size())
{
DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl
@ -2150,6 +2151,9 @@ DynamicModel::BlockLinear(Model_Block *ModelBlock)
}
}
void
DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives,
const eval_context_type &eval_context, bool no_tmp_terms)
@ -2210,6 +2214,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
block_triangular.incidencematrix.Print_IM(eEndogenous);
}
block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations);
BlockLinear(block_triangular.ModelBlock);
if (!no_tmp_terms)
computeTemporaryTermsOrdered(block_triangular.ModelBlock);
@ -2372,7 +2378,7 @@ DynamicModel::computeDynJacobianCols(bool jacobianExo)
const int &deriv_id = it->second;
SymbolType type = symbol_table.getType(symb_id);
int tsid = symbol_table.getTypeSpecificID(symb_id);
switch(type)
{
case eEndogenous:

1
DynamicModel.hh
View File

@ -110,6 +110,7 @@ private:
//! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms();
public:
DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants);
//! Adds a variable node

14
IncidenceMatrix.cc
View File

@ -211,16 +211,16 @@ IncidenceMatrix::Print_IM(SymbolType type) const
//------------------------------------------------------------------------------
// Swap rows and columns of the incidence matrix
void
IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const
IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int n) const
{
int tmp_i, j;
bool tmp_b;
/* We exchange equation (row)...*/
if(pos1 != pos2)
{
tmp_i = Index_Equ_IM[pos1].index;
Index_Equ_IM[pos1].index = Index_Equ_IM[pos2].index;
Index_Equ_IM[pos2].index = tmp_i;
tmp_i = Index_Equ_IM[pos1];
Index_Equ_IM[pos1] = Index_Equ_IM[pos2];
Index_Equ_IM[pos2] = tmp_i;
for(j = 0;j < n;j++)
{
tmp_b = SIM[pos1 * n + j];
@ -231,9 +231,9 @@ IncidenceMatrix::swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Inde
/* ...and variables (column)*/
if(pos1 != pos3)
{
tmp_i = Index_Var_IM[pos1].index;
Index_Var_IM[pos1].index = Index_Var_IM[pos3].index;
Index_Var_IM[pos3].index = tmp_i;
tmp_i = Index_Var_IM[pos1];
Index_Var_IM[pos1] = Index_Var_IM[pos3];
Index_Var_IM[pos3] = tmp_i;
for(j = 0;j < n;j++)
{
tmp_b = SIM[j * n + pos1];

2
IncidenceMatrix.hh
View File

@ -46,7 +46,7 @@ public:
void Free_IM() const;
void Print_IM(SymbolType type) const;
void Print_SIM(bool* IM, SymbolType type) const;
void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, simple* Index_Var_IM, simple* Index_Equ_IM, int n) const;
void swap_IM_c(bool *SIM, int pos1, int pos2, int pos3, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int n) const;
int Model_Max_Lead, Model_Max_Lag;
int Model_Max_Lead_Endo, Model_Max_Lag_Endo, Model_Max_Lead_Exo, Model_Max_Lag_Exo;
private:

1
Makefile.in
View File

@ -28,6 +28,7 @@ MAIN_OBJS = \
ExprNode.o \
ModelNormalization.o \
ModelBlocks.o \
MinimumFeedbackSet.o \
IncidenceMatrix.o \
BlockTriangular.o \
ModelGraph.o \

637
MinimumFeedbackSet.cc Normal file
View File

@ -0,0 +1,637 @@
#include "MinimumFeedbackSet.hh"
namespace MFS
{
void
Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G)
{
/*clear all in and out edges of vertex_to_eliminate
and remove vertex_to_eliminate from the graph*/
clear_vertex(vertex_to_eliminate, G);
remove_vertex(vertex_to_eliminate, G);
}
void
Suppress(int vertex_num, AdjacencyList_type& G)
{
Suppress(vertex(vertex_num, G), G);
}
void
Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G)
{
/*before the vertex i suppression replace all edges e_k_i and e_i_j by e_k_j*/
if (in_degree (vertex_to_eliminate, G) > 0 and out_degree (vertex_to_eliminate, G) > 0)
{
AdjacencyList_type::in_edge_iterator it_in, in_end;
AdjacencyList_type::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_type::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_type& g, AdjacencyList_type::vertex_descriptor u, color_type& color, vector<int> &circuit_stack)
{
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, g);
color[u] = gray_color;
graph_traits<AdjacencyList_type>::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)
{
if (has_cycle_dfs(g, target(*vi, g), color, circuit_stack))
{
circuit_stack.push_back(v_index[target(*vi, g)]);
return true; // cycle detected, return immediately
}
}
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_cylce(AdjacencyList_type& g, vector<int> &circuit_stack)
{
color_type color;
graph_traits<AdjacencyList_type>::vertex_iterator vi, vi_end;
for (tie(vi, vi_end) = vertices(g); vi != vi_end; vi++)
color[*vi] = white_color;
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, g);
for (tie(vi, vi_end) = vertices(g); vi != vi_end; vi++)
if (color[*vi] == white_color)
if (has_cycle_dfs(g, *vi, color, circuit_stack))
return true;
return false;
}
bool
has_cycle(vector<int> &circuit_stack, AdjacencyList_type& G)
{
return has_cylce(G, circuit_stack);
}
void
Print(AdjacencyList_type& G)
{
AdjacencyList_type::vertex_iterator it, it_end, it_begin;
property_map<AdjacencyList_type, 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_type::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_type::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_type
AM_2_AdjacencyList(bool* AM, unsigned int n)
{
AdjacencyList_type G(n);
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
property_map<AdjacencyList_type, 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;
}
void
Print(GraphvizDigraph& G)
{
GraphvizDigraph::vertex_iterator it, it_end, it_begin;
property_map<GraphvizDigraph, 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 << "] ->";
GraphvizDigraph::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";
}
}
GraphvizDigraph
AM_2_GraphvizDigraph(bool* AM, unsigned int n)
{
GraphvizDigraph G(n);
property_map<GraphvizDigraph, vertex_index_t>::type v_index = get(vertex_index, G);
/*for (unsigned int i = 0;i < n;i++)
cout << "v_index[" << i << "] = " << v_index[i] << "\n";*/
//put(v_index, vertex(i, G), i);
//v_index[/*vertex(i,G)*/i]["v_index"]=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_type
GraphvizDigraph_2_AdjacencyList(GraphvizDigraph& G1, set<int> select_index)
{
unsigned int n = select_index.size();
//cout << "n=" << n << "\n";
AdjacencyList_type G(n);
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
property_map<AdjacencyList_type, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
property_map<GraphvizDigraph, vertex_index_t>::type v1_index = get(vertex_index, G1);
set<int>::iterator it = select_index.begin();
map<int,int> reverse_index;
for (unsigned int i = 0;i < n;i++, it++)
{
reverse_index[v1_index[*it]]=i;
put(v_index, vertex(i, G), v1_index[*it]);
put(v_index1, vertex(i, G), i);
}
unsigned int i;
for (it = select_index.begin(), i = 0;i < n;i++, it++)
{
GraphvizDigraph::out_edge_iterator it_out, out_end;
GraphvizDigraph::vertex_descriptor vi = vertex(*it, G1);
for (tie(it_out, out_end) = out_edges(vi, G1); it_out != out_end; ++it_out)
{
int ii = target(*it_out, G1);
if (select_index.find(ii) != select_index.end())
{
/*cout << "*it=" << *it << " i = " << i << " ii=" << ii << " n=" << n << " *it_out=" << *it_out << "\n";
cout << "source(*it_out, G1) = " << source(*it_out, G1) << " target(*it_out, G1) = " << target(*it_out, G1) << "\n";
cout << "vertex(source(*it_out, G1), G) = " << vertex(source(*it_out, G1), G) << " vertex(target(*it_out, G1), G) = " << vertex(target(*it_out, G1), G) << "\n";*/
add_edge( vertex(reverse_index[source(*it_out, G1)],G), vertex(reverse_index[target(*it_out, G1)], G), G);
//add_edge(vertex(source(*it_out, G1), G) , vertex(target(*it_out, G1), G), G);
}
}
}
return G;
}
vector_vertex_descriptor
Collect_Doublet(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G)
{
/*collect all doublet (for each edge e_i_k there is an edge e_k_i with k!=i) in the graph
and return the vector of doublet*/
AdjacencyList_type::in_edge_iterator it_in, in_end;
AdjacencyList_type::out_edge_iterator it_out, out_end;
vector<AdjacencyList_type::vertex_descriptor> Doublet;
if (in_degree(vertex, G) > 0 and 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) and 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_type::vertex_descriptor vertex, AdjacencyList_type& G)
{
/*Detect all the clique (all vertex in a clique are related to each other) in the graph*/
vector<AdjacencyList_type::vertex_descriptor> liste;
bool agree = true;
AdjacencyList_type::in_edge_iterator it_in, in_end;
AdjacencyList_type::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 and it_out != out_end and agree)
{
agree = (source(*it_in, G) == target(*it_out, G) and 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 or it_out != out_end)
agree = false;
unsigned int i = 1;
while (i < liste.size() and agree)
{
unsigned int j = i + 1;
while (j < liste.size() and agree)
{
AdjacencyList_type::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 and exist2);
j++;
}
i++;
}
}
return agree;
}
bool
Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G)
{
/*Graph reduction: eliminating purely intermediate variables or variables outside of any circuit*/
bool something_has_been_done = false;
bool not_a_loop;
int i;
AdjacencyList_type::vertex_iterator it, it1, ita, it_end, it_begin;
tie(it, it_end) = vertices(G);
it_begin = it;
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
for ( 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 or out_degree_n <= 1)
{
not_a_loop = true;
if (in_degree_n >= 1 and out_degree_n >= 1) //do not eliminate a vertex if it loops on its self!
{
AdjacencyList_type::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_type, 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_type& G)
{
/*Graphe reduction: elimination of a vertex inside a clique*/
AdjacencyList_type::vertex_iterator it, it1, ita, it_end, it_begin;
bool something_has_been_done = false;
int i;
tie(it, it_end) = vertices(G);
it_begin = it;
for (i = 0;it != it_end; ++it, i++)
{
if (Vertex_Belong_to_a_Clique(*it, G))
{
#ifdef verbose
property_map<AdjacencyList_type, 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_Edge_i_j_if_not_a_loop_and_if_for_all_i_k_edge_we_have_a_k_j_edge_Step(AdjacencyList_type& G) //Suppression
{
bool something_has_been_done = false;
AdjacencyList_type::vertex_iterator it, it_end;
int i = 0;
bool agree;
for (tie(it, it_end) = vertices(G);it != it_end; ++it, i++)
{
AdjacencyList_type::in_edge_iterator it_in, in_end;
AdjacencyList_type::out_edge_iterator it_out, out_end, it_out1, ita_out;
int j = 0;
for (tie(ita_out = it_out, out_end) = out_edges(*it, G); it_out != out_end; ++it_out, j++)
{
AdjacencyList_type::edge_descriptor ed;
bool exist;
tie(ed, exist) = edge(target(*it_out, G), source(*it_out, G) , G);
if (!exist)
{
agree = true;
for (tie(it_out1, out_end) = out_edges(*it, G); it_out1 != out_end; ++it_out1)
{
bool exist;
tie(ed, exist) = edge(target(*it_out1, G), target(*it_out, G) , G);
if (target(*it_out1, G) != target(*it_out, G) and !exist)
agree = false;
}
if (agree)
{
something_has_been_done = true;
remove_edge(*it_out, G);
if (out_degree(*it, G) == 0)
break;
if (j > 0)
{
it_out = ita_out;
tie(it_out1, out_end) = out_edges(*it, G);
}
else
{
tie(it_out, out_end) = out_edges(*it, G);
j--;
}
}
}
ita_out = it_out;
}
}
return something_has_been_done;
}
bool
Suppression_of_all_in_Edge_in_i_if_not_a_loop_and_if_all_doublet_i_eq_Min_inDegree_outDegree_Step(AdjacencyList_type& G)
{
bool something_has_been_done = false;
AdjacencyList_type::vertex_iterator it, it_end;
int i = 0;
for (tie(it, it_end) = vertices(G);it != it_end; ++it, i++)
{
AdjacencyList_type::in_edge_iterator it_in, in_end, it_in1, ita_in;
vector<AdjacencyList_type::vertex_descriptor> doublet = Collect_Doublet(*it, G);
if (doublet.size() == (unsigned int) min(in_degree(*it, G), out_degree(*it, G)))
{
int j = 0;
if (in_degree(*it, G))
for (tie(ita_in = it_in, in_end) = in_edges(*it, G); it_in != in_end; ++it_in, j++)
{
vector<AdjacencyList_type::vertex_descriptor>::iterator it1 = doublet.begin();
bool not_a_doublet = true;
while (it1 != doublet.end() and not_a_doublet)
{
if (target(*it_in, G) == *it1)
not_a_doublet = false;
it1++;
}
if (not_a_doublet and source(*it_in, G) != target(*it_in, G))
{
#ifdef verbose
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
cout << "remove_edge(" << v_index[source(*it_in, G)] << ", " << v_index[target(*it_in, G)] << ", G) j=" << j << " it_in == in_end : " << (it_in == in_end) << " in_degree(*it, G)=" << in_degree(*it, G) << ";\n";
#endif
something_has_been_done = true;
remove_edge(source(*it_in, G), target(*it_in, G), G);
cout << " in_degree(*it, G)=" << in_degree(*it, G) << ";\n";
if (in_degree(*it, G) == 0)
break;
if (j > 0)
{
it_in = ita_in;
}
else
{
tie(it_in, in_end) = in_edges(*it, G);
j--;
}
}
ita_in = it_in;
}
}
}
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_type& G)
{
/*If a vertex loop on itself it's a feedback variable
we eliminate it from the graph and store the vertex
in the minimum feedback set*/
bool something_has_been_done = false;
AdjacencyList_type::vertex_iterator it, it_end, it_begin, ita;
int i = 0;
tie(it, it_end) = vertices(G);
it_begin = it;
for (;it != it_end; ++it, i++)
{
AdjacencyList_type::edge_descriptor ed;
bool exist;
tie(ed, exist) = edge(*it, *it , G);
if (exist)
{
#ifdef verbose
property_map<AdjacencyList_type, vertex_index_t>::type v_index = get(vertex_index, G);
cout << "store v[*it] = " << v_index[*it]+1 << "\n";
#endif
property_map<AdjacencyList_type, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
feed_back_vertices.insert(v_index1[*it] );
/*property_map<AdjacencyList_type, 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_type
Minimal_set_of_feedback_vertex(set<int> &feed_back_vertices,const AdjacencyList_type& G1)
{
bool something_has_been_done = true;
int cut_ = 0;
AdjacencyList_type G(G1);
while (num_vertices(G) > 0)
{
while (something_has_been_done and 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) or 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_Edge_i_j_if_not_a_loop_and_if_for_all_i_k_edge_we_have_a_k_j_edge_Step(G) or something_has_been_done);
#ifdef verbose
cout << "3 something_has_been_done=" << something_has_been_done << "\n";
#endif
//Rule 4
//something_has_been_done=(Suppression_of_all_in_Edge_in_i_if_not_a_loop_and_if_all_doublet_i_eq_Min_inDegree_outDegree_Step(G) or something_has_been_done);
#ifdef verbose
cout << "4 something_has_been_done=" << something_has_been_done << "\n";
#endif
//Rule 5
something_has_been_done = (Suppression_of_Vertex_X_if_it_loops_store_in_set_of_feedback_vertex_Step(feed_back_vertices, G) or something_has_been_done);
#ifdef verbose
cout << "5 something_has_been_done=" << something_has_been_done << "\n";
#endif
}
vector<int> circuit;
if (!has_cycle(circuit, G))
{
#ifdef verobse
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_type::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_type, vertex_index1_t>::type v_index1 = get(vertex_index1, G);
feed_back_vertices.insert(v_index1[*max_degree_index]);
/*property_map<AdjacencyList_type, 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_type, 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);
}
};
vector<int>
Reorder_the_recursive_variables(const AdjacencyList_type& G1, set<int> &feedback_vertices)
{
AdjacencyList_type G(G1);
property_map<AdjacencyList_type, 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++)
{
//cout << "supress " << v_index[vertex(*its, G)]+1 << " " << *its << "\n";
Suppress(*its, G);
}
vector< int> Reordered_Vertices;
bool something_has_been_done = true;
while (something_has_been_done)
{
something_has_been_done = false;
AdjacencyList_type::vertex_iterator it, it_end, it_begin, ita;
tie(it, it_end) = vertices(G);
int i = 0;
for (it_begin = it;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";
return Reordered_Vertices;
}
}

50
MinimumFeedbackSet.hh Normal file
View File

@ -0,0 +1,50 @@
#ifndef MFE_BOOST
#define MFE_BOOST
#include <iostream>
#include <map>
#include <vector>
#include <boost/graph/graphviz.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <stdio.h>
#include <stdlib.h>
using namespace std;
using namespace boost;
//#define verbose
typedef property<vertex_index_t, int,
property<vertex_index1_t, int,
property<vertex_degree_t, int,
property<vertex_in_degree_t, int,
property<vertex_out_degree_t, int > > > > > VertexProperty;
typedef adjacency_list<listS, listS, bidirectionalS, VertexProperty> AdjacencyList_type;
typedef map<graph_traits<AdjacencyList_type>::vertex_descriptor,default_color_type> color_type;
typedef vector<AdjacencyList_type::vertex_descriptor> vector_vertex_descriptor;
namespace MFS
{
void Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G);
vector_vertex_descriptor Collect_Doublet(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G);
bool Vertex_Belong_to_a_Clique(AdjacencyList_type::vertex_descriptor vertex, AdjacencyList_type& G);
bool Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G); //Graph reduction: eliminating purely intermediate variables or variables outside of any circuit
bool Elimination_of_Vertex_belonging_to_a_clique_Step(AdjacencyList_type& G); //Graphe reduction: eliminaion of Cliques
bool Suppression_of_Edge_i_j_if_not_a_loop_and_if_for_all_i_k_edge_we_have_a_k_j_edge_Step(AdjacencyList_type& G); //Suppression
bool Suppression_of_all_in_Edge_in_i_if_not_a_loop_and_if_all_doublet_i_eq_Min_inDegree_outDegree_Step(AdjacencyList_type& G);
bool Suppression_of_Vertex_X_if_it_loops_store_in_set_of_feedback_vertex_Step(vector<pair<int, AdjacencyList_type::vertex_descriptor> > &looping_variable, AdjacencyList_type& G);
void Print(AdjacencyList_type& G);
AdjacencyList_type AM_2_AdjacencyList(bool* AMp,unsigned int n);
void Print(GraphvizDigraph& G);
GraphvizDigraph AM_2_GraphvizDigraph(bool* AM, unsigned int n);
AdjacencyList_type GraphvizDigraph_2_AdjacencyList(GraphvizDigraph& G1, set<int> select_index);
bool has_cycle_dfs(AdjacencyList_type& g, AdjacencyList_type::vertex_descriptor u, color_type& color, vector<int> &circuit_stack);
bool has_cylce(AdjacencyList_type& g, vector<int> &circuit_stack, int size);
bool has_cycle(vector<int> &circuit_stack, AdjacencyList_type& G);
AdjacencyList_type Minimal_set_of_feedback_vertex(set<int> &feed_back_vertices, const AdjacencyList_type& G);
void Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G);
void Suppress(int vertex_num, AdjacencyList_type& G);
vector<int> Reorder_the_recursive_variables(const AdjacencyList_type& G1, set<int> &feed_back_vertices);
};
#endif

16
ModelBlocks.cc
View File

@ -249,25 +249,19 @@ Blocks::block_result_print(block_result_t *r)
void
Blocks::block_result_to_IM(block_result_t *r,bool* IM,int prologue, int n,simple* Index_Equ_IM,simple* Index_Var_IM)
Blocks::block_result_to_IM(block_result_t *r,bool* IM,int prologue, int n,vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM)
{
int i, j, k, l;
bool* SIM=(bool*)malloc(n*n*sizeof(*SIM));
simple* Index_Equ_IM_tmp=(simple*)malloc(n*sizeof(*Index_Equ_IM_tmp));
simple* Index_Var_IM_tmp=(simple*)malloc(n*sizeof(*Index_Var_IM_tmp));
vector<int> Index_Equ_IM_tmp(Index_Equ_IM), Index_Var_IM_tmp(Index_Var_IM);
for(i=0;i<n*n;i++)
SIM[i]=IM[i];
for(i=0;i<n;i++)
{
Index_Equ_IM_tmp[i].index=Index_Equ_IM[i].index;
Index_Var_IM_tmp[i].index=Index_Var_IM[i].index;
}
l=prologue;
for(i = 0; i < r->n_sets; i++)
{
for(j = r->sets_s[r->ordered[i]]; j <= r->sets_f[r->ordered[i]]; j++)
{
Index_Equ_IM[l].index=Index_Equ_IM_tmp[r->vertices[j]+prologue].index;
Index_Equ_IM[l]=Index_Equ_IM_tmp[r->vertices[j]+prologue];
for(k=0;k<n;k++)
SIM[l*n+k]=IM[(r->vertices[j]+prologue)*n+k];
l++;
@ -280,14 +274,12 @@ Blocks::block_result_to_IM(block_result_t *r,bool* IM,int prologue, int n,simple
{
for(j = r->sets_s[r->ordered[i]]; j <= r->sets_f[r->ordered[i]]; j++)
{
Index_Var_IM[l].index=Index_Var_IM_tmp[r->vertices[j]+prologue].index;
Index_Var_IM[l]=Index_Var_IM_tmp[r->vertices[j]+prologue];
for(k=0;k<n;k++)
IM[k*n+l]=SIM[(k*n+r->vertices[j]+prologue)];
l++;
}
}
free(Index_Equ_IM_tmp);
free(Index_Var_IM_tmp);
free(SIM);
}

2
ModelBlocks.hh
View File

@ -43,7 +43,7 @@ public:
void block_result_print(block_result_t *r);
void Print_Equation_gr(Equation_set* Equation);
//! Converts the output of Tarjan algorithm into reordered incidence matrices
void block_result_to_IM(block_result_t *r,bool* IM,int prologue, int n,simple* Index_Equ_IM,simple* Index_Var_IM);
void block_result_to_IM(block_result_t *r,bool* IM,int prologue, int n,vector<int> &Index_Equ_IM,vector<int> &Index_Var_IM);
Equation_vertex *vertices;
int *block_vertices, *sets_s, *sets_f;
int *block_stack, *sp, tos;

19
ModelNormalization.cc
View File

@ -348,12 +348,12 @@ Normalization::Gr_to_IM_basic(int n0, int prologue, int epilogue, bool* IM, Equa
}
void
Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* Index_Equ_IM, Equation_set *Equation, bool mixing, bool* IM_s)
Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, vector<int> &Index_Equ_IM, Equation_set *Equation, bool mixing, bool* IM_s)
{
int i, j, n, l;
Edge *e1, *e2;
Equation_set* Equation_p;
simple* Index_Equ_IM_tmp = (simple*)malloc(n0 * sizeof(*Index_Equ_IM_tmp));
vector<int> Index_Equ_IM_tmp(Index_Equ_IM);
bool* SIM = (bool*)malloc(n0 * n0 * sizeof(bool));
#ifdef DEBUG
cout << "in Gr_to_IM\n";
@ -363,14 +363,12 @@ Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* In
{
for(i = 0;i < n0*n0;i++)
SIM[i] = IM_s[i];
for(i = 0;i < n0;i++)
Index_Equ_IM_tmp[i].index = Index_Equ_IM[i].index;
for(i = 0;i < n;i++)
{
/*Index_Var_IM[j+prologue].index=Index_Var_IM_tmp[Equation->Number[j].matched+prologue].index;*/
if(fp_verbose)
cout << "Equation->Number[" << i << "].matched=" << Equation->Number[i].matched << "\n";
Index_Equ_IM[i + prologue].index = Index_Equ_IM_tmp[Equation->Number[i].matched + prologue].index;
Index_Equ_IM[i + prologue] = Index_Equ_IM_tmp[Equation->Number[i].matched + prologue];
for(j = 0;j < n0;j++)
SIM[(i + prologue)*n0 + j] = IM_s[(Equation->Number[i].matched + prologue) * n0 + j];
}
@ -382,12 +380,12 @@ Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* In
for(i = 0;i < n0*n0;i++)
SIM[i] = IM[i];
for(i = 0;i < n0;i++)
Index_Equ_IM_tmp[i].index = Index_Equ_IM[i].index;
Index_Equ_IM_tmp[i] = Index_Equ_IM[i];
for(j = 0;j < n;j++)
{
if(fp_verbose)
cout << "Equation->Number[" << j << "].matched=" << Equation->Number[j].matched << "\n";
Index_Equ_IM[j + prologue].index = Index_Equ_IM_tmp[Equation->Number[j].matched + prologue].index;
Index_Equ_IM[j + prologue] = Index_Equ_IM_tmp[Equation->Number[j].matched + prologue];
for(i = 0;i < n0;i++)
SIM[(i)*n0 + j + prologue] = IM[(i) * n0 + Equation->Number[j].matched + prologue];
}
@ -395,7 +393,6 @@ Normalization::Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* In
IM[i] = SIM[i];
}
free(SIM);
free(Index_Equ_IM_tmp);
//cout << "mixing=" << mixing << "\n";
if(mixing)
{
@ -506,7 +503,7 @@ Normalization::Set_fp_verbose(bool ok)
}
bool
Normalization::Normalize(int n, int prologue, int epilogue, bool* IM, simple* Index_Equ_IM, Equation_set* Equation, bool mixing, bool* IM_s)
Normalization::Normalize(int n, int prologue, int epilogue, bool* IM, vector<int> &Index_Equ_IM, Equation_set* Equation, bool mixing, bool* IM_s)
{
int matchingSize, effective_n;
int save_fp_verbose=fp_verbose;
@ -531,7 +528,7 @@ Normalization::Normalize(int n, int prologue, int epilogue, bool* IM, simple* In
cout << "\n and the following variables:\n - ";
for(i = 0; i < Variable->size; i++)
if(Variable->Number[i].matched == -1)
cout << symbol_table.getName(Index_Equ_IM[i].index) << " ";
cout << symbol_table.getName(Index_Equ_IM[i]) << " ";
cout << "\n could not be normalized\n";
//ErrorHandling(n, IM, Index_Equ_IM);
//system("PAUSE");
@ -544,7 +541,7 @@ Normalization::Normalize(int n, int prologue, int epilogue, bool* IM, simple* In
{
OutputMatching(Equation);
for(int i = 0;i < n;i++)
cout << "Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i].index /*<< " == " "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index*/ << "\n";
cout << "Index_Equ_IM[" << i << "]=" << Index_Equ_IM[i]/*<< " == " "Index_Var_IM[" << i << "]=" << Index_Var_IM[i].index*/ << "\n";
}
}
Free_Other(Variable);

10
ModelNormalization.hh
View File

@ -46,12 +46,6 @@ struct Equation_set
int edges;
};
//! Stores result of block decomposition for ONE equation or ONE variable
struct simple
{
//! New {variable, equation} index after reordering
int index;
};
//! Computes the model normalization
class Normalization
@ -76,7 +70,7 @@ private:
};
public:
Normalization(const SymbolTable &symbol_table_arg);
bool Normalize(int n, int prologue, int epilogue, bool* IM, simple* Index_Var_IM, Equation_set* Equation,bool mixing, bool* IM_s);
bool Normalize(int n, int prologue, int epilogue, bool* IM, vector<int> &Index_Var_IM, Equation_set* Equation,bool mixing, bool* IM_s);
void Gr_to_IM_basic(int n0, int prologue, int epilogue, bool* IM, Equation_set *Equation,bool transpose);
const SymbolTable &symbol_table;
void Set_fp_verbose(bool ok);
@ -90,7 +84,7 @@ private:
void MaximumMatching(Equation_set *Equation, Variable_set *Variable);
int MeasureMatching(Equation_set *Equation);
void OutputMatching(Equation_set* Equation);
void Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, simple* Index_Var_IM, Equation_set *Equation,bool mixing, bool* IM_s);
void Gr_to_IM(int n0, int prologue, int epilogue, bool* IM, vector<int> &Index_Var_IM, Equation_set *Equation,bool mixing, bool* IM_s);
void Free_Other(Variable_set* Variable);
void Free_All(int n, Equation_set* Equation, Variable_set* Variable);
int eq, eex;