From 916e950df966de7865ce35f28f6cff02eab2d242 Mon Sep 17 00:00:00 2001 From: ferhat Date: Tue, 12 May 2009 23:03:40 +0000 Subject: [PATCH] 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 --- BlockTriangular.cc | 528 +++++++++++++++++++--------------- BlockTriangular.hh | 28 +- DynamicModel.cc | 20 +- DynamicModel.hh | 1 + IncidenceMatrix.cc | 14 +- IncidenceMatrix.hh | 2 +- Makefile.in | 1 + MinimumFeedbackSet.cc | 637 ++++++++++++++++++++++++++++++++++++++++++ MinimumFeedbackSet.hh | 50 ++++ ModelBlocks.cc | 16 +- ModelBlocks.hh | 2 +- ModelNormalization.cc | 19 +- ModelNormalization.hh | 10 +- 13 files changed, 1052 insertions(+), 276 deletions(-) create mode 100644 MinimumFeedbackSet.cc create mode 100644 MinimumFeedbackSet.hh diff --git a/BlockTriangular.cc b/BlockTriangular.cc index 1dad468b..6cd75e56 100644 --- a/BlockTriangular.cc +++ b/BlockTriangular.cc @@ -25,15 +25,26 @@ #include #include #include -using namespace std; +#include +#include +#include "MinimumFeedbackSet.hh" +#include +#include +#include //------------------------------------------------------------------------------ #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 &Index_Var_IM, vector &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 &Index_Equ_IM) const + { + int n = equation_number - prologue - epilogue; + + + typedef adjacency_list 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::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::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 &Index_Equ_IM, vector &Index_Var_IM, vector > &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 component(num_vertices(G2)), discover_time(num_vertices(G2)); + + //write_graphviz(cout, G2); + + int num = strong_components(G2, &component[0]); + + blocks = vector >(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, vector > > > components_set(num); + + for (unsigned int i=0; i 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 feed_back_vertices; + AdjacencyList_type G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G); + property_map::type v_index = get(vertex_index, G); + components_set[i].second.first = feed_back_vertices; + blocks[i].second = feed_back_vertices.size(); + vector Reordered_Vertice; + Reordered_Vertice = Reorder_the_recursive_variables(G, feed_back_vertices); + for (vector::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::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 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; kBlock_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; iBlock_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; iBlock_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 equations ) +BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector 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; in_sets+epilogue; i++) + for ( i=0; in_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 equations ) +void +BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0, jacob_map j_m, vector 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, 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+*epiloguen_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 > blocks; + if (prologue+epiloguefirst==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_Equsecond==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(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(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); } diff --git a/BlockTriangular.hh b/BlockTriangular.hh index 2434a1a8..9b3e2696 100644 --- a/BlockTriangular.hh +++ b/BlockTriangular.hh @@ -34,31 +34,39 @@ -//! Matrix of doubles for representing jacobian +//! Sparse matrix of double to store the values of the Jacobian typedef map,double> jacob_map; typedef vector > 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 &Index_Var_IM, vector &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 &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 &Index_Equ_IM, vector &Index_Var_IM, vector > &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 > &blocks, vector 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 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 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 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 &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0 , jacob_map j_m, vector equations); + vector Index_Equ_IM; + vector Index_Var_IM; int prologue, epilogue; bool bt_verbose; //int endo_nbr, exo_nbr; diff --git a/DynamicModel.cc b/DynamicModel.cc index 03a3b020..0ff40bed 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -21,6 +21,7 @@ #include #include + #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;iBlock_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;iBlock_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: diff --git a/DynamicModel.hh b/DynamicModel.hh index fd0d16ff..107db04d 100644 --- a/DynamicModel.hh +++ b/DynamicModel.hh @@ -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 diff --git a/IncidenceMatrix.cc b/IncidenceMatrix.cc index 9d3205a0..5d8de6ef 100644 --- a/IncidenceMatrix.cc +++ b/IncidenceMatrix.cc @@ -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 &Index_Var_IM, vector &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]; diff --git a/IncidenceMatrix.hh b/IncidenceMatrix.hh index aefca204..9e8b35b1 100644 --- a/IncidenceMatrix.hh +++ b/IncidenceMatrix.hh @@ -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 &Index_Var_IM, vector &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: diff --git a/Makefile.in b/Makefile.in index 2fc4d8ec..5c57d887 100644 --- a/Makefile.in +++ b/Makefile.in @@ -28,6 +28,7 @@ MAIN_OBJS = \ ExprNode.o \ ModelNormalization.o \ ModelBlocks.o \ + MinimumFeedbackSet.o \ IncidenceMatrix.o \ BlockTriangular.o \ ModelGraph.o \ diff --git a/MinimumFeedbackSet.cc b/MinimumFeedbackSet.cc new file mode 100644 index 00000000..b4c15998 --- /dev/null +++ b/MinimumFeedbackSet.cc @@ -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 &circuit_stack) + { + property_map::type v_index = get(vertex_index, g); + color[u] = gray_color; + graph_traits::out_edge_iterator vi, vi_end; + for (tie(vi, vi_end) = out_edges(u, g); vi != vi_end; ++vi) + if (color[target(*vi, g)] == white_color) + { + 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 &circuit_stack) + { + color_type color; + graph_traits::vertex_iterator vi, vi_end; + for (tie(vi, vi_end) = vertices(g); vi != vi_end; vi++) + color[*vi] = white_color; + property_map::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 &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::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::type v_index = get(vertex_index, G); + property_map::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::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::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 select_index) + { + unsigned int n = select_index.size(); + //cout << "n=" << n << "\n"; + AdjacencyList_type G(n); + property_map::type v_index = get(vertex_index, G); + property_map::type v_index1 = get(vertex_index1, G); + property_map::type v1_index = get(vertex_index, G1); + set::iterator it = select_index.begin(); + map 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 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 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::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::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::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 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::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::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 &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::type v_index = get(vertex_index, G); + cout << "store v[*it] = " << v_index[*it]+1 << "\n"; +#endif + property_map::type v_index1 = get(vertex_index1, G); + feed_back_vertices.insert(v_index1[*it] ); + /*property_map::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 &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 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::type v_index1 = get(vertex_index1, G); + feed_back_vertices.insert(v_index1[*max_degree_index]); + /*property_map::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::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 + Reorder_the_recursive_variables(const AdjacencyList_type& G1, set &feedback_vertices) + { + AdjacencyList_type G(G1); + property_map::type v_index = get(vertex_index, G); + set::iterator its, ita; + set 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; + } + + } + + diff --git a/MinimumFeedbackSet.hh b/MinimumFeedbackSet.hh new file mode 100644 index 00000000..eeec3e99 --- /dev/null +++ b/MinimumFeedbackSet.hh @@ -0,0 +1,50 @@ +#ifndef MFE_BOOST +#define MFE_BOOST + +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace boost; + +//#define verbose + + typedef property > > > > VertexProperty; + typedef adjacency_list AdjacencyList_type; + typedef map::vertex_descriptor,default_color_type> color_type; + typedef vector 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 > &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 select_index); + bool has_cycle_dfs(AdjacencyList_type& g, AdjacencyList_type::vertex_descriptor u, color_type& color, vector &circuit_stack); + bool has_cylce(AdjacencyList_type& g, vector &circuit_stack, int size); + bool has_cycle(vector &circuit_stack, AdjacencyList_type& G); + AdjacencyList_type Minimal_set_of_feedback_vertex(set &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 Reorder_the_recursive_variables(const AdjacencyList_type& G1, set &feed_back_vertices); +}; + +#endif diff --git a/ModelBlocks.cc b/ModelBlocks.cc index a7a6ea12..8a09ceca 100644 --- a/ModelBlocks.cc +++ b/ModelBlocks.cc @@ -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 &Index_Equ_IM, vector &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 Index_Equ_IM_tmp(Index_Equ_IM), Index_Var_IM_tmp(Index_Var_IM); for(i=0;in_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;kvertices[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;kvertices[j]+prologue)]; l++; } } - free(Index_Equ_IM_tmp); - free(Index_Var_IM_tmp); free(SIM); } diff --git a/ModelBlocks.hh b/ModelBlocks.hh index 94daf995..1cc7e234 100644 --- a/ModelBlocks.hh +++ b/ModelBlocks.hh @@ -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 &Index_Equ_IM,vector &Index_Var_IM); Equation_vertex *vertices; int *block_vertices, *sets_s, *sets_f; int *block_stack, *sp, tos; diff --git a/ModelNormalization.cc b/ModelNormalization.cc index 0c0cea0c..7aed8c7c 100644 --- a/ModelNormalization.cc +++ b/ModelNormalization.cc @@ -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 &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 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 &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); diff --git a/ModelNormalization.hh b/ModelNormalization.hh index 8feaf904..e3e99f57 100644 --- a/ModelNormalization.hh +++ b/ModelNormalization.hh @@ -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 &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 &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;