Commit 25ca5d3e authored by ferhat's avatar ferhat
Browse files

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
parent 2eed8a5f
This diff is collapsed.
......@@ -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;
......
......@@ -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:
......
......@@ -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
......
......@@ -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];
......
......@@ -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:
......
......@@ -28,6 +28,7 @@ MAIN_OBJS = \
ExprNode.o \
ModelNormalization.o \
ModelBlocks.o \
MinimumFeedbackSet.o \
IncidenceMatrix.o \
BlockTriangular.o \
ModelGraph.o \
......
This diff is collapsed.
#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
......@@ -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);
}
......
......@@ -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;
......
......@@ -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);
......
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment