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

- Correction of several bugs 

- normalize an equation linear in its endogenous variable
- Chained rule derivatives (necessary to reduce a block to the feedback equations and variables)

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2726 ac1d8469-bf42-47a9-8791-bf33cf982152
issue#70
ferhat 2009-06-05 14:45:23 +00:00
parent 518c5fba93
commit 3737c1aa2e
11 changed files with 4263 additions and 3973 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
BlockTriangular.cc
View File

@ -17,7 +17,6 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <sstream>
#include <fstream>
@ -38,20 +37,22 @@ 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)
BlockTriangular::BlockTriangular(SymbolTable &symbol_table_arg, NumericalConstants &num_const_arg)
: symbol_table(symbol_table_arg),
//normalization(symbol_table_arg),
incidencematrix(symbol_table_arg),
num_const(num_const_arg)
{
bt_verbose = 0;
ModelBlock = NULL;
periods = 0;
Normalized_Equation = new DataTree(symbol_table, num_const);
}
BlockTriangular::~BlockTriangular()
{
delete Normalized_Equation;
}
//------------------------------------------------------------------------------
// Find the prologue and the epilogue of the model
@ -110,7 +111,6 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n
}
}
//------------------------------------------------------------------------------
// Find a matching between equations and endogenous variables
bool
@ -118,7 +118,6 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
{
int n = equation_number - prologue - epilogue;
typedef adjacency_list<vecS, vecS, undirectedS> BipartiteGraph;
/*
@ -132,7 +131,6 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
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);
@ -169,9 +167,62 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
return check;
}
t_vtype
BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const
{
int nb_endo = symbol_table.endo_nbr();
vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
t_vtype Variable_Type(nb_endo);
for (int i = 0; i < nb_endo; i++)
{
if (i < prologue)
{
variable_blck[Index_Var_IM[i]] = i;
equation_blck[Index_Equ_IM[i]] = i;
}
else if (i < components_set.size() + prologue)
{
variable_blck[Index_Var_IM[i]] = components_set[i-prologue] + prologue;
equation_blck[Index_Equ_IM[i]] = components_set[i-prologue] + prologue;
}
else
{
variable_blck[Index_Var_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
equation_blck[Index_Equ_IM[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
}
//cout << "equation_blck[" << Index_Equ_IM[i] << "]=" << equation_blck[Index_Equ_IM[i]] << " variable_blck[" << Index_Var_IM[i] << "] = " << variable_blck[Index_Var_IM[i]] << "\n";
Variable_Type[i] = make_pair(0, 0);
}
for (int k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{
bool *Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if (Cur_IM)
{
for (int i = 0; i < nb_endo; i++)
{
int i_1 = Index_Var_IM[i];
for (int j = 0; j < nb_endo; j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[Index_Equ_IM[ j]])
{
if (k > Variable_Type[i_1].second)
Variable_Type[i_1] = make_pair(Variable_Type[i_1].first, k);
if (k < -Variable_Type[i_1].first)
Variable_Type[i_1] = make_pair(-k, Variable_Type[i_1].second);
}
}
}
}
}
/*for(int i = 0; i < nb_endo; i++)
cout << "Variable_Type[" << Index_Equ_IM[i] << "].first = " << Variable_Type[Index_Equ_IM[i]].first << " Variable_Type[" << Index_Equ_IM[i] << "].second=" << Variable_Type[Index_Equ_IM[i]].second << "\n";*/
return (Variable_Type);
}
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, t_etype &Equation_Type, bool verbose_) const
{
t_vtype V_Variable_Type;
int n = nb_var - prologue - epilogue;
bool *AMp;
AMp = (bool *) malloc(n*n*sizeof(bool));
@ -204,6 +255,7 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
components_set[component[i]].first.insert(i);
}
V_Variable_Type = Get_Variable_LeadLag_By_Block(component, num, prologue, epilogue);
vector<int> tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
int order = prologue;
@ -211,9 +263,9 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
SIM = (bool *) malloc(nb_var*nb_var*sizeof(bool));
memcpy(SIM, IM, nb_var*nb_var*sizeof(bool));
//Add a loop on vertices which could not be normalized => force those vvertices to belong to the feedback set
//Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
for (int i = 0; i < n; i++)
if(Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or Equation_Type[Index_Equ_IM[i+prologue]].first == E_EVALUATE_S)
if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or V_Variable_Type[Index_Var_IM[i+prologue]].second >= 0 or V_Variable_Type[Index_Var_IM[i+prologue]].first > 0)
add_edge(i, i, G2);
//For each block, the minimum set of feedback variable is computed
@ -260,8 +312,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
bool *IM, OK;
int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo;
cout << "Allocate Block=" << count_Block << " recurs_Size=" << recurs_Size << "\n";
ModelBlock->Periods = periods;
ModelBlock->Block_List[count_Block].is_linear = true;
ModelBlock->Block_List[count_Block].Size = size;
@ -272,12 +322,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].Simulation_Type = SimType;
ModelBlock->Block_List[count_Block].Equation = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
ModelBlock->Block_List[count_Block].Equation_Type = (EquationType *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(EquationType));
ModelBlock->Block_List[count_Block].Equation_Type_Var = (int*)malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
/*cout << "ModelBlock->Block_List[" << count_Block << "].Size = " << ModelBlock->Block_List[count_Block].Size << " E_UNKNOWN=" << E_UNKNOWN << "\n";
t_etype eType(size);
cout << "eType[0].first=" << eType[0].first << " eType[0].second=" << eType[0].second << " eType.size()=" << eType.size() << "\n";
ModelBlock->Block_List[count_Block].Equation_Type(eType);
cout << "after that\n";*/
ModelBlock->Block_List[count_Block].Equation_Normalized = (NodeID *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(NodeID));
ModelBlock->Block_List[count_Block].Variable = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation = (temporary_terms_type **) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(temporary_terms_type));
ModelBlock->Block_List[count_Block].Own_Derivative = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
@ -287,7 +332,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
tmp_endo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_other_endo = (int *) malloc(symbol_table.endo_nbr() * sizeof(int));
tmp_size = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
//cout << "tmp_size = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1= " << incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1 << ") * sizeof(int))\n";
tmp_size_other_endo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size_exo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
@ -308,7 +352,14 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].Equation[i] = Index_Equ_IM[*count_Equ];
ModelBlock->Block_List[count_Block].Variable[i] = Index_Var_IM[*count_Equ];
ModelBlock->Block_List[count_Block].Equation_Type[i] = Equation_Type[Index_Equ_IM[*count_Equ]].first;
ModelBlock->Block_List[count_Block].Equation_Type_Var[i] = Equation_Type[Index_Equ_IM[*count_Equ]].second;
ModelBlock->Block_List[count_Block].Equation_Normalized[i] = Equation_Type[Index_Equ_IM[*count_Equ]].second;
/*if(Equation_Type[Index_Equ_IM[*count_Equ]].second)
{
temporary_terms_type temporary_terms;
cout << "Equation_Type[Index_Equ_IM[*count_Equ]].second->get_op_code()=" << Equation_Type[Index_Equ_IM[*count_Equ]].second->get_op_code() << "\n";
cout << "ModelBlock->Block_List[" << count_Block << "].Equation_Normalized[" << i << "]->get_op_code()=" << ModelBlock->Block_List[count_Block].Equation_Normalized[i]->get_op_code() << "\n";
ModelBlock->Block_List[count_Block].Equation_Normalized[i]->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
}*/
i_1 = Index_Var_IM[*count_Equ];
for (k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{
@ -395,7 +446,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].nb_other_endo = tmp_nb_other_endo;
ModelBlock->Block_List[count_Block].Other_Endogenous = (int *) malloc(tmp_nb_other_endo * sizeof(int));
tmp_exo = (int *) malloc(symbol_table.exo_nbr() * sizeof(int));
memset(tmp_exo, 0, symbol_table.exo_nbr() * sizeof(int));
tmp_nb_exo = 0;
@ -429,7 +479,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
}
}
ModelBlock->Block_List[count_Block].nb_exo = tmp_nb_exo;
ModelBlock->Block_List[count_Block].Exogenous = (int *) malloc(tmp_nb_exo * sizeof(int));
k = 0;
@ -577,22 +626,19 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
free(tmp_variable_evaluated);
}
void
BlockTriangular::Free_Block(Model_Block *ModelBlock) const
{
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);
free(ModelBlock->Block_List[blk].Equation_Type);
free(ModelBlock->Block_List[blk].Equation_Type_Var);
free(ModelBlock->Block_List[blk].Equation_Normalized);
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*/)
@ -646,43 +692,41 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations,
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
Equation_Simulation_Type = E_SOLVE;
int Var_To_Derivate = -1;
tmp_s.str("");
tmp_output.str("");
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
tmp_s << "y(it_, " << Index_Var_IM[i]+1 << ")";
map<pair<int, int>, NodeID>::iterator derivative = first_cur_endo_derivatives.find(make_pair(eq, var));
/*cout << "eq=" << eq << " var=" << var << "=";
first_cur_endo_derivatives[make_pair(eq, var)]->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
cout << "\n";
cout << "equation : ";
eq_node->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
cout << "\n";*/
set<pair<int, int> > result;
//result.clear();
pair<bool, NodeID> res;
derivative->second->collectEndogenous(result);
/*for(set<pair<int, int> >::const_iterator iitt = result.begin(); iitt!=result.end(); iitt++)
cout << "result = " << iitt->first << ", " << iitt->second << "\n";*/
set<pair<int, int> >::const_iterator d_endo_variable = result.find(make_pair(var, 0));
//Determine whether the equation could be evaluated rather than to be solved
if (tmp_output.str() == tmp_s.str() and d_endo_variable == result.end())
Equation_Simulation_Type = E_EVALUATE;
else
{ //the equation could be normalized by a permutation of the rhs and the lhs
tmp_output.str("");
rhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
if (tmp_output.str() == tmp_s.str() and d_endo_variable == result.end())
{
Equation_Simulation_Type = E_EVALUATE_R;
//cout << "Equation " << eq << " is reversed\n";
Equation_Simulation_Type = E_EVALUATE;
}
else
{ //the equation could be normalized using the derivative independant of the endogenous variable
if (d_endo_variable == result.end())
{
//the equation could be normalized by a permutation of the rhs and the lhs
if (d_endo_variable == result.end()) //the equation is linear in the endogenous and could be normalized using the derivative
{
Equation_Simulation_Type = E_EVALUATE_S;
Var_To_Derivate = var;
//cout << " gone normalized : ";
res = equations[eq]->normalizeLinearInEndoEquation(var, derivative->second);
/*res.second->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
cout << " done\n";*/
}
}
}
/*cout << "-----------------------------------------------------------\n";
lhs->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
cout << " = ";
rhs->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
cout << "% " << c_Equation_Type(Equation_Simulation_Type) << " " << var+1 << "\n";*/
V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, Var_To_Derivate);
V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second));
}
return (V_Equation_Simulation_Type);
}
@ -763,7 +807,7 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
}
if (Blck_Size == 1)
{
if(Equation_Type[Index_Equ_IM[eq]].first==E_EVALUATE or Equation_Type[Index_Equ_IM[eq]].first==E_EVALUATE_R)
if (Equation_Type[Index_Equ_IM[eq]].first == E_EVALUATE /*or Equation_Type[Index_Equ_IM[eq]].first==E_EVALUATE_R*/ or Equation_Type[Index_Equ_IM[eq]].first == E_EVALUATE_S)
{
if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
Simulation_Type = EVALUATE_BACKWARD;
@ -795,7 +839,6 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
return (Type);
}
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, t_etype &V_Equation_Type, map<pair<int, int >, NodeID> &first_cur_endo_derivatives)
{
@ -824,6 +867,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
free(max_val);
bool OK = false;
double bi = 0.99999999;
//double bi=1e-13;
int suppressed = 0;
vector<int> Index_Equ_IM_save(Index_Equ_IM);
while (!OK && bi > 1e-14)
@ -858,14 +902,13 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
suppressed = suppress;
if (!OK)
//bi/=1.07;
bi/=3;
bi /= 2;
counted++;
if (bi > 1e-14)
free(SIM00);
}
if (!OK)
Compute_Normalization(IM, n, prologue, epilogue, true, SIM00, Index_Equ_IM);
}
V_Equation_Type = Equation_Type_determination(equations, first_cur_endo_derivatives, Index_Var_IM, Index_Equ_IM);
@ -898,14 +941,15 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
cout << Nb_TotalBlocks << " block(s) found:\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";
cout << " the largest simultaneous block has " << j << " equation(s)\n"
<<" and " << Nb_fv << " feedback variable(s).\n";
ModelBlock->Size = Nb_TotalBlocks;
ModelBlock->Periods = periods;
ModelBlock->Block_List = (Block *) malloc(sizeof(ModelBlock->Block_List[0]) * Nb_TotalBlocks);
count_Equ = count_Block = 0;
for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
{
if (count_Equ < prologue)
@ -921,7 +965,6 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
}
}
//------------------------------------------------------------------------------
// normalize each equation of the dynamic model
// and find the optimal block triangular decomposition of the static model

41
BlockTriangular.hh
View File

@ -24,21 +24,23 @@
#include "CodeInterpreter.hh"
#include "ExprNode.hh"
#include "SymbolTable.hh"
#include "ModelNormalization.hh"
#include "ModelBlocks.hh"
//#include "ModelNormalization.hh"
//#include "ModelBlocks.hh"
#include "IncidenceMatrix.hh"
#include "ModelTree.hh"
//! Sparse matrix of double to store the values of the Jacobian
typedef map<pair<int ,int >,double> jacob_map;
typedef vector<pair<BlockSimulationType, pair<int, int> > > t_type;
typedef vector<pair<EquationType, int> > t_etype;
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
typedef vector<pair<EquationType, NodeID > > t_etype;
//! Vector describing variables: max_lag in the block, max_lead in the block
typedef vector<pair< int, int> > t_vtype;
//! Creates the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
class BlockTriangular
@ -52,27 +54,32 @@ private:
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, t_etype &Equation_Type, bool verbose_) const;
//! determine the type of each equation of the model (couble evaluated or need to be solved)
//! determines the type of each equation of the model (could be evaluated or need to be solved)
t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, int >, NodeID> &first_cur_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! 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, t_etype &Equation_Type);
//! Compute for each variable its maximum lead and lag in its block
t_vtype Get_Variable_LeadLag_By_Block(vector<int > &components_set, int nb_blck_sim, int prologue, int epilogue) const;
public:
const SymbolTable &symbol_table;
BlockTriangular(const SymbolTable &symbol_table_arg);
SymbolTable &symbol_table;
/*Blocks blocks;
Normalization normalization;*/
IncidenceMatrix incidencematrix;
NumericalConstants &num_const;
DataTree *Normalized_Equation;
BlockTriangular(SymbolTable &symbol_table_arg, NumericalConstants &num_const_arg);
~BlockTriangular();
//! 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(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, int >, NodeID> &first_cur_endo_derivatives);
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, t_etype &equation_simulation_type, map<pair<int, int >, NodeID> &first_cur_endo_derivatives);
vector<int> Index_Equ_IM;
vector<int> Index_Var_IM;
int prologue, epilogue;
bool bt_verbose;
//int endo_nbr, exo_nbr;
Model_Block* ModelBlock;
int periods;
inline static std::string BlockType0(int type)
@ -101,11 +108,11 @@ public:
switch (type)
{
case EVALUATE_FORWARD:
case EVALUATE_FORWARD_R:
//case EVALUATE_FORWARD_R:
return ("EVALUATE FORWARD ");
break;
case EVALUATE_BACKWARD:
case EVALUATE_BACKWARD_R:
//case EVALUATE_BACKWARD_R:
return ("EVALUATE BACKWARD ");
break;
case SOLVE_FORWARD_SIMPLE:
@ -137,7 +144,7 @@ public:
{
"E_UNKNOWN ",
"E_EVALUATE ",
"E_EVALUATE_R",
//"E_EVALUATE_R",
"E_EVALUATE_S",
"E_SOLVE "
};

8
CodeInterpreter.hh
View File

@ -53,8 +53,8 @@ enum EquationType
{
E_UNKNOWN, //!< Unknown equation type
E_EVALUATE, //!< Simple evaluation, normalized variable on left-hand side
E_EVALUATE_R, //!< Simple evaluation, normalized variable on right-hand side
E_EVALUATE_S, //!< Simple evaluation, normalize using the first order derivative which does not involve the normalized variable
//E_EVALUATE_R, //!< Simple evaluation, normalized variable on right-hand side
E_EVALUATE_S, //!< Simple evaluation, normalize using the first order derivative
E_SOLVE //!< No simple evaluation of the equation, it has to be solved
};
@ -71,8 +71,8 @@ enum BlockSimulationType
SOLVE_FORWARD_COMPLETE, //!< Block of several equations, newton solver needed, forward
SOLVE_BACKWARD_COMPLETE, //!< Block of several equations, newton solver needed, backward
SOLVE_TWO_BOUNDARIES_COMPLETE, //!< Block of several equations, newton solver needed, forward and backwar
EVALUATE_FORWARD_R, //!< Simple evaluation, normalized variable on right-hand side, forward
EVALUATE_BACKWARD_R //!< Simple evaluation, normalized variable on right-hand side, backward
//EVALUATE_FORWARD_R, //!< Simple evaluation, normalized variable on right-hand side, forward
//EVALUATE_BACKWARD_R //!< Simple evaluation, normalized variable on right-hand side, backward
};
//! Enumeration of possible symbol types

461
DynamicModel.cc
View File

@ -43,7 +43,7 @@ DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
dynJacobianColsNbr(0),
cutoff(1e-15),
markowitz(0.7),
block_triangular(symbol_table_arg)
block_triangular(symbol_table_arg, num_constants_arg)
{
}
@ -56,7 +56,8 @@ DynamicModel::AddVariable(const string &name, int lag)
void
DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const
{
first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
//first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag)));
if (it != first_derivatives.end())
(it->second)->compile(code_file, false, temporary_terms, map_idx);
else
@ -77,7 +78,7 @@ DynamicModel::BuildIncidenceMatrix()
Id->collectEndogenous(endogenous);
for (set<pair<int, int> >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++)
{
block_triangular.incidencematrix.fill_IM(eq, symbol_table.getTypeSpecificID(it_endogenous->first), it_endogenous->second, eEndogenous);
block_triangular.incidencematrix.fill_IM(eq, it_endogenous->first, it_endogenous->second, eEndogenous);
}
exogenous.clear();
Id = eq_node->get_arg1();
@ -86,7 +87,7 @@ DynamicModel::BuildIncidenceMatrix()
Id->collectExogenous(exogenous);
for (set<pair<int, int> >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++)
{
block_triangular.incidencematrix.fill_IM(eq, symbol_table.getTypeSpecificID(it_exogenous->first), it_exogenous->second, eExogenous);
block_triangular.incidencematrix.fill_IM(eq, it_exogenous->first, it_exogenous->second, eExogenous);
}
}
}
@ -107,12 +108,14 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
map_idx.clear();
for (j = 0;j < ModelBlock->Size;j++)
{
//cerr << "block=" << j+1 << "\n";
// Compute the temporary terms reordered
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx);
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
if(ModelBlock->Block_List[j].Equation_Normalized[i])
ModelBlock->Block_List[j].Equation_Normalized[i]->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx);
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -158,11 +161,16 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
}
for (j = 0;j < ModelBlock->Size;j++)
{
// Compute the temporary terms reordered
// Collecte the temporary terms reordered
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
eq_node->collectTemporary_terms(temporary_terms, ModelBlock, j);
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
if(ModelBlock->Block_List[j].Equation_Normalized[i])
ModelBlock->Block_List[j].Equation_Normalized[i]->collectTemporary_terms(temporary_terms, ModelBlock, j);
for(temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->begin(); it!= ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->end(); it++)
(*it)->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@ -227,20 +235,28 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
int prev_Simulation_Type=-1, count_derivates=0;
int jacobian_max_endo_col;
ofstream output;
temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
int nze, nze_exo, nze_other_endo;
map<int, NodeID> recursive_variables;
vector<int> feedback_variables;
//----------------------------------------------------------------------
//For each block
for (j = 0;j < ModelBlock->Size;j++)
{
recursive_variables.clear();
feedback_variables.clear();
//For a block composed of a single equation determines wether we have to evaluate or to solve the equation
nze = nze_exo = nze_other_endo = 0;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
nze+=ModelBlock->Block_List[j].IM_lead_lag[m].size;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead_Exo+ModelBlock->Block_List[j].Max_Lag_Exo;m++)
nze_exo+=ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead_Other_Endo+ModelBlock->Block_List[j].Max_Lag_Other_Endo;m++)
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead_Exo+ModelBlock->Block_List[j].Max_Lag_Exo;m++)
nze_exo+=ModelBlock->Block_List[j].IM_lead_lag[m].size_exo;*/
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
if (block_triangular.incidencematrix.Model_Max_Lag_Endo - ModelBlock->Block_List[j].Max_Lag +m >=0)
nze_other_endo+=ModelBlock->Block_List[j].IM_lead_lag[m].size_other_endo;
}
tmp1_output.str("");
tmp1_output << dynamic_basename << "_" << j+1 << ".m";
output.open(tmp1_output.str().c_str(), ios::out | ios::binary);
@ -252,19 +268,19 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << "%/\n";
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
/*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R*/)
{
output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, jacobian_eval, y_kmin, periods)\n";
}
else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE)
output << "function [residual, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE
|| ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE)
output << "function [residual, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
else
output << "function [residual, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, periods, jacobian_eval, y_kmin, y_size)\n";
output << "function [residual, y, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, periods, jacobian_eval, y_kmin, y_size)\n";
output << " % ////////////////////////////////////////////////////////////////////////" << endl
<< " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type)
<< " //" << endl
@ -274,26 +290,45 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
//The Temporary terms
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
/*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R*/)
{
output << " if(jacobian_eval)\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size*(1+ModelBlock->Block_List[j].Max_Lag_Endo+ModelBlock->Block_List[j].Max_Lead_Endo) << ", " << nze << ");\n";
//output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
//output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)*(1+ModelBlock->Block_List[j].Max_Lag_Endo+ModelBlock->Block_List[j].Max_Lead_Endo)
<< ", " << nze << ");\n";
output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
output << " end;\n";
}
else
{
output << " if(jacobian_eval)\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size*(1+ModelBlock->Block_List[j].Max_Lag_Endo+ModelBlock->Block_List[j].Max_Lead_Endo) << ", " << nze << ");\n";
//output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
output << " g1_x=spalloc(" << ModelBlock->Block_List[j].Size << ", " << (ModelBlock->Block_List[j].nb_exo + ModelBlock->Block_List[j].nb_exo_det)*(1+ModelBlock->Block_List[j].Max_Lag_Exo+ModelBlock->Block_List[j].Max_Lead_Exo) << ", " << nze_exo << ");\n";
output << " g1_o=spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].nb_other_endo*(1+ModelBlock->Block_List[j].Max_Lag_Other_Endo+ModelBlock->Block_List[j].Max_Lead_Other_Endo) << ", " << nze_other_endo << ");\n";
output << " else\n";
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size*ModelBlock->Periods << ", " << ModelBlock->Block_List[j].Size*(ModelBlock->Periods+ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead) << ", " << nze*ModelBlock->Periods << ");\n";
{
output << " g1 = spalloc(" << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)*ModelBlock->Periods
<< ", " << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)*(ModelBlock->Periods+ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1)
<< ", " << nze*ModelBlock->Periods << ");\n";
output << " g1_tmp_r = spalloc(" << (ModelBlock->Block_List[j].Nb_Recursives)
<< ", " << (ModelBlock->Block_List[j].Size)*(ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1)
<< ", " << nze << ");\n";
output << " g1_tmp_b = spalloc(" << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)
<< ", " << (ModelBlock->Block_List[j].Size)*(ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1)
<< ", " << nze << ");\n";
}
else
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
{
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives << ", " << nze << ");\n";
output << " g1_tmp_r = spalloc(" << ModelBlock->Block_List[j].Nb_Recursives
<< ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
output << " g1_tmp_b = spalloc(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
}
output << " end;\n";
}
@ -306,17 +341,16 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
tmp_output << " T" << *it;
output << " global" << tmp_output.str() << ";\n";
}
output << " residual=zeros(" << ModelBlock->Block_List[j].Size << ",1);\n";
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD and ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD)
output << " residual=zeros(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives << ",1);\n";
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD)
output << " for it_ = (y_kmin+periods):y_kmin+1\n";
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
output << " for it_ = y_kmin+1:(y_kmin+periods)\n";
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
output << " b = [];\n";
output << " b = zeros(periods*y_size,1);\n";
output << " for it_ = y_kmin+1:(periods+y_kmin)\n";
output << " Per_y_=it_*y_size;\n";
output << " Per_J_=(it_-y_kmin-1)*y_size;\n";
@ -324,8 +358,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
sps=" ";
}
else
if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD ||
ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R)
if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD )
sps = " ";
else
sps="";
@ -347,7 +380,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
tt2.insert(*it);
output << ";" << endl;
}
string sModel = symbol_table.getName(ModelBlock->Block_List[j].Variable[i]) ;
string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i])) ;
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
@ -367,35 +400,46 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " = ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
}
else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_R)
/*else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_R)
{
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " = ";
output << tmp_output.str();
output << "; %reversed " << ModelBlock->Block_List[j].Equation_Type[i] << " \n";
}*/
else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
{
output << "%" << tmp_output.str();
output << " = ";
if (ModelBlock->Block_List[j].Equation_Normalized[i])
{
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << "\n ";
temporary_terms_type tt2;
tt2.clear();
ModelBlock->Block_List[j].Equation_Normalized[i]->writeOutput(output , oMatlabDynamicModelSparse, temporary_terms/*tt2*/);
}
}
else
{
cerr << "Type missmatch for equation " << ModelBlock->Block_List[j].Equation[i]+1 << "\n";
exit(-1);
exit(EXIT_FAILURE);
}
output << ";\n";
break; /*case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ")" << endl;
output << " ";
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " = ";
lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << ";\n";
break;*/
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
if (i<ModelBlock->Block_List[j].Nb_Recursives)
{
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i];
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];
goto evaluation;
}
feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
output << " " << "residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ") = (";
@ -403,10 +447,17 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
if (i<ModelBlock->Block_List[j].Nb_Recursives)
{
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = ModelBlock->Block_List[j].Equation_Normalized[i];
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i]), 0)] = equations[ModelBlock->Block_List[j].Equation[i]];
goto evaluation;
}
feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
output << " % equation " << ModelBlock->Block_List[j].Equation[i]+1 << " variable : " << sModel
<< " (" << ModelBlock->Block_List[j].Variable[i]+1 << ") " << block_triangular.c_Equation_Type(ModelBlock->Block_List[j].Equation_Type[i]) << endl;
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_) = -residual(" << i+1 << ", it_)";
Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_) = -residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ", it_)";
output << " residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ", it_) = (";
goto end;
default:
@ -435,8 +486,8 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
/*case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:*/
count_derivates++;
for (m=0;m<ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
{
@ -450,7 +501,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " g1(" << eqr+1 << ", " << /*varr+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr*/
varr+1+m*ModelBlock->Block_List[j].Size << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
@ -488,7 +539,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " g1_o(" << eqr+1 << ", "
<< varr+1+(m+max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
@ -497,6 +548,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " varargout{1}=g1_x;\n";
output << " varargout{2}=g1_o;\n";
output << " end;" << endl;
//output << " ya = y;\n";
output << " end;" << endl;
break;
case SOLVE_BACKWARD_SIMPLE:
@ -513,12 +565,11 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
output << " g1(" << eqr+1 << ", " << /*varr+1+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr*/
varr+1+m*ModelBlock->Block_List[j].Size << ") = ";
output << " g1(" << eqr+1 << ", "
<< varr+1 + m*(ModelBlock->Block_List[j].Size) << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
<< "(" << k
<< ") " << var+1
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
}
@ -549,10 +600,10 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_other_endo[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_other_endo[i];
output << " g1_o(" << eqr+1 << ", "
output << " g1_o(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1+(m+max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
@ -563,19 +614,77 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " else" << endl;
m=ModelBlock->Block_List[j].Max_Lag;
//cout << "\nDerivatives in Block " << j << "\n";
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
output << " g1(" << eqr+1 << ", " << varr+1 << ") = ";
bool derivative_exist;
ostringstream tmp_output;
if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{
/*tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
NodeID tmp_n;
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE)
tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
else
tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];
int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0);
NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables);
ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;*/
}
}
else
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{
output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
/*writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);*/
/*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables);
else
derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, 0, 0, recursive_variables, feedback_variables);
//if (derivative_exist)
output << tmp_output.str() << ";";*/
NodeID tmp_n;
//if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
/*else
tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/
//cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
//cout << "derivaive eq=" << eq << " var=" << var << " k0=" << k << "\n";
int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),0);
map<int, NodeID> recursive_variables_save(recursive_variables);
NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables, var, 0);
recursive_variables = recursive_variables_save;
ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
}
/*if (eqr<ModelBlock->Block_List[j].Nb_Recursives or varr<ModelBlock->Block_List[j].Nb_Recursives)
output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), 0, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
<< "(0) " << var+1
<< ", equation=" << eq+1 << endl;
<< ", equation=" << eq+1 << endl;*/
}
output << " end;\n";
//output << " ya = y;\n";
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
@ -589,26 +698,157 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
bool derivative_exist;
ostringstream tmp_output;
//cout << "ModelBlock->Block_List[" << j << "].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n";
if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
{
/*if (varr<ModelBlock->Block_List[j].Nb_Recursives)
{
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_)*y(it_, " << var+1 << ")";
else if (k==1)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_)*y(it_+1, " << var+1 << ")";
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
<< ", " << varr+1
<< "+" << ModelBlock->Block_List[j].Size*ModelBlock->Block_List[j].Max_Lag << ")*y(it_, " << var+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")";
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
<< ", " << varr+1
<< "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_+" << k << ", " << var+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")";
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_r(" << eqr+1
<< ", " << varr+1
<< "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ")*y(it_" << k << ", " << var+1 << ")";
if (k==0)
output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = ";
else if (k==1)
output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_y_) = ";
tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k>0)
output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = ";
tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k<0)
output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = ";
tmp_output << " g1_tmp_r(" << eqr+1 << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
else
{
BinaryOpNode* tt = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[eqr];
derivative_exist = tt->get_arg2()->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
}
if (derivative_exist)
output << tmp_output.str() << ";";
else
{
//output << "1" << ";";
if (ModelBlock->Block_List[j].Equation_Type[eqr] != E_EVALUATE)
{
}
}
//writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << " %1 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << " derivative_exist=" << derivative_exist << " varr+1=" << varr+1 << endl;
}*/
}
else
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_K_)*y(it_, " << var+1 << ")";
else if (k==1)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_y_)*y(it_+1, " << var+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+Per_J_, " << varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives
<< "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")";
if (k==0)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = ";
else if (k==1)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = ";
else if (k>0)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = ";
else if (k<0)
tmp_output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = ";
/*NodeID tmp_n;
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
tmp_n = equations[ModelBlock->Block_List[j].Equation[eqr]];
else
tmp_n = ModelBlock->Block_List[j].Equation_Normalized[eqr];*/
/*int deriv_id = getDerivID(symbol_table.getID(eEndogenous, var),k);
//cout << "-------------------------------------------------------------------------------------\n";
//cout << "derivaive eq=" << eq << " var=" << var << " k=" << k << "\n";
//output << " " << tmp_output.str();
map<int, NodeID> recursive_variables_save(recursive_variables);
NodeID ChaineRule_Derivative = tmp_n->getChaineRuleDerivative(deriv_id ,recursive_variables, var, k);
recursive_variables = recursive_variables_save;
//ChaineRule_Derivative->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
*/
/*if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
else
derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
if (derivative_exist)
output << tmp_output.str() << ";";*/
/*output << tmp_output.str();*/
//output << ";";
//output << "\n%";
output << tmp_output.str();
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
/*else
{
if (k==0)
tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k>0)
tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
else if (k<0)
tmp_output << " g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag << ") = ";
if (ModelBlock->Block_List[j].Equation_Type[eqr] == E_EVALUATE or ModelBlock->Block_List[j].Equation_Type[eqr] == E_SOLVE)
derivative_exist = equations[ModelBlock->Block_List[j].Equation[eqr]]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
else
derivative_exist = ModelBlock->Block_List[j].Equation_Normalized[eqr]->writeOutputDerivativesRespectToFeedBackVariables(tmp_output, oMatlabDynamicModelSparse, temporary_terms, eq, var, varr, k, ModelBlock->Block_List[j].Max_Lag, recursive_variables, feedback_variables);
if (derivative_exist)
output << tmp_output.str() << ";";
if (k==0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag
<< ")*y(it_, " << var+1 << ")";
else if (k>0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag
<< ")*y(it_+" << k << ", " << var+1 << ")";
else if (k<0)
Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1_tmp_b(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives
<< ", " << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << k+ModelBlock->Block_List[j].Max_Lag
<< ")*y(it_" << k << ", " << var+1 << ")";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}*/
}
#ifdef CONDITION
output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
@ -617,6 +857,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
}
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
if (i>=ModelBlock->Block_List[j].Nb_Recursives)
output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
#ifdef CONDITION
output << " if (fabs(condition(" << i+1 << "))<fabs(u(" << i << "+Per_u_)))\n";
@ -652,7 +893,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
output << " g1(" << eqr+1 << ", " << varr+1+(m-ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lag_Endo)*ModelBlock->Block_List[j].Size << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
@ -689,7 +930,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " g1_o(" << eqr+1 << ", "
<< varr+1+(m+max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr() << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(var)
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
}
@ -698,6 +939,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " varargout{1}=g1_x;\n";
output << " varargout{2}=g1_o;\n";
output << " end;\n";
//output << " ya = y;\n";
output << " end;\n";
break;
default:
@ -736,7 +978,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
int prev_Simulation_Type=-1;
//SymbolicGaussElimination SGE;
bool file_open=false;
temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
//----------------------------------------------------------------------
string main_name=file_name;
main_name+=".cod";
@ -758,8 +1000,8 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type)
&& (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R ))
/*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R */))
{
}
else
@ -875,14 +1117,26 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
{
rhs->compile(code_file, false, temporary_terms, map_idx);
lhs->compile(code_file, true, temporary_terms, map_idx);
}
else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
{
eq_node = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[i];
//cout << "EVALUATE_S var " << ModelBlock->Block_List[j].Variable[i] << "\n";
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
rhs->compile(code_file, false, temporary_terms, map_idx);
lhs->compile(code_file, true, temporary_terms, map_idx);
}
break;
case EVALUATE_BACKWARD_R:
/*case EVALUATE_BACKWARD_R:
case EVALUATE_FORWARD_R:
lhs->compile(code_file, false, temporary_terms, map_idx);
rhs->compile(code_file, true, temporary_terms, map_idx);
break;
break;*/
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
v=ModelBlock->Block_List[j].Equation[i];
@ -914,8 +1168,8 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
// The Jacobian if we have to solve the block
if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD_R)
/*&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD_R*/)
{
switch (ModelBlock->Block_List[j].Simulation_Type)
{
@ -1357,8 +1611,8 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
{
k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
if ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD || prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R)
|| (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)))
((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD /*|| prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R*/)
|| (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/)))
{
mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
tmp_eq.str("");
@ -1372,7 +1626,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
}
if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)
if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/)
{
if (i==block_triangular.ModelBlock->Size-1)
{
@ -1385,7 +1639,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
}
}
if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))
(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/))
skip_head=true;
else
skip_head=false;
@ -1393,8 +1647,8 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
{
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD_R:
case EVALUATE_BACKWARD_R:
/*case EVALUATE_FORWARD_R:
case EVALUATE_BACKWARD_R:*/
if (!skip_head)
{
tmp1 << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n";
@ -1439,7 +1693,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
/*mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " <<
block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";*/
tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size << ");\n";
mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size-block_triangular.ModelBlock->Block_List[i].Nb_Recursives << ");\n";
/*if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n";
else
@ -1484,11 +1738,11 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
{
k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R))
(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/))
skip_head=true;
else
skip_head=false;
if ((k == EVALUATE_FORWARD || k == EVALUATE_FORWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
if ((k == EVALUATE_FORWARD /*|| k == EVALUATE_FORWARD_R*/) && (block_triangular.ModelBlock->Block_List[i].Size))
{
if (!skip_head)
{
@ -1510,10 +1764,15 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
//mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n";
mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << i <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
//open_par=true;
}
else if ((k == EVALUATE_BACKWARD || k == EVALUATE_BACKWARD_R) && (block_triangular.ModelBlock->Block_List[i].Size))
else if ((k == EVALUATE_BACKWARD /*|| k == EVALUATE_BACKWARD_R*/) && (block_triangular.ModelBlock->Block_List[i].Size))
{
if (!skip_head)
{
@ -1534,6 +1793,11 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n";
mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n";
mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << i <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
}
else if ((k == SOLVE_FORWARD_COMPLETE || k == SOLVE_FORWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
@ -1561,7 +1825,11 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
", y, x, params, y_index, " << nze <<
", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n";
mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
else if ((k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
{
@ -1573,6 +1841,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
tmp.str("");
for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
if (ik>=block_triangular.ModelBlock->Block_List[i].Nb_Recursives)
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
@ -1588,6 +1857,11 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
", y, x, params, y_index, " << nze <<
", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n";
mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE || k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
{
@ -1601,6 +1875,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " y_index=[";
for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
if (ik>=block_triangular.ModelBlock->Block_List[i].Nb_Recursives)
mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
mDynamicModelFile << " ];\n";
@ -1615,7 +1890,11 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
", " << block_triangular.ModelBlock->Block_List[i].Max_Lead <<
", " << block_triangular.ModelBlock->Block_List[i].is_linear <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method);\n";
mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
prev_Simulation_Type=k;
}

585
ExprNode.cc
View File

@ -60,6 +60,32 @@ ExprNode::getDerivative(int deriv_id)
}
}
NodeID
ExprNode::getChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_)
{
// Return zero if derivative is necessarily null (using symbolic a priori)
/*set<int>::const_iterator it = non_null_derivatives.find(deriv_id);
if (it == non_null_derivatives.end())
{
cout << "0\n";
return datatree.Zero;
}
*/
// If derivative is stored in cache, use the cached value, otherwise compute it (and cache it)
/*map<int, NodeID>::const_iterator it2 = derivatives.find(deriv_id);
if (it2 != derivatives.end())
return it2->second;
else*/
{
NodeID d = computeChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
//derivatives[deriv_id] = d;
return d;
}
}
int
ExprNode::precedence(ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const
{
@ -94,12 +120,25 @@ ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
// Nothing to do for a terminal node
}
void
ExprNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
//writeOutput(output, oMatlabOutsideModel, temporary_terms_type());
}
void
ExprNode::writeOutput(ostream &output)
{
writeOutput(output, oMatlabOutsideModel, temporary_terms_type());
}
pair<bool, NodeID>
ExprNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
return(make_pair(true, (NodeID)NULL));
}
NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
ExprNode(datatree_arg),
@ -129,6 +168,7 @@ void
NumConstNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
//cout << "writeOutput constante\n";
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<NumConstNode *>(this));
if (it != temporary_terms.end())
if (output_type == oMatlabDynamicModelSparse)
@ -166,6 +206,27 @@ NumConstNode::collectExogenous(set<pair<int, int> > &result) const
{
}
pair<bool, NodeID>
NumConstNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
return(make_pair(false, datatree.AddNumConstant(datatree.num_constants.get(id))));
}
NodeID
NumConstNode::computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_)
{
return datatree.Zero;
}
/*
pair<bool, NodeID>
NumConstNode::computeDerivativeRespectToFeedbackVariable(int equ, int var, int varr, int lag_, int max_lag, vector<int> &recursive_variables, vector<int> &feeback_variables) const
{
return(make_pair(false, datatree.Zero));
}
*/
NodeID
NumConstNode::toStatic(DataTree &static_datatree) const
{
@ -251,7 +312,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
{
// If node is a temporary term
#ifdef DEBUGC
cout << "write_ouput output_type=" << output_type << "\n";
cout << "write_ouput Variable output_type=" << output_type << "\n";
#endif
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this));
if (it != temporary_terms.end())
@ -428,15 +489,20 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms
cout << "output_type=" << output_type << "\n";
#endif
if (!lhs_rhs)
{
CompileCode.write(&FLDV, sizeof(FLDV));
}
else
{
CompileCode.write(&FSTPV, sizeof(FSTPV));
}
char typel=(char)type;
CompileCode.write(&typel, sizeof(typel));
int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
switch (type)
{
case eParameter:
//cout << "Parameter=" << tsid << "\n";
i = tsid;
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
#ifdef DEBUGC
@ -444,12 +510,14 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms
#endif
break;
case eEndogenous :
i = symb_id;
//cout << "Endogenous=" << symb_id << "\n";
i = tsid;//symb_id;
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
lagl=lag;
CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
break;
case eExogenous :
//cout << "Exogenous=" << tsid << "\n";
i = tsid;
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
lagl=lag;
@ -457,12 +525,14 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms
break;
case eExogenousDet:
i = tsid + datatree.symbol_table.exo_nbr();
//cout << "ExogenousDet=" << i << "\n";
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
lagl=lag;
CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
break;
case eModelLocalVariable:
case eModFileLocalVariable:
//cout << "eModelLocalVariable=" << symb_id << "\n";
datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
break;
case eUnknownFunction:
@ -488,7 +558,7 @@ void
VariableNode::collectEndogenous(set<pair<int, int> > &result) const
{
if (type == eEndogenous)
result.insert(make_pair(symb_id, lag));
result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(symb_id), lag));
else if (type == eModelLocalVariable)
datatree.local_variables_table[symb_id]->collectEndogenous(result);
}
@ -497,11 +567,87 @@ void
VariableNode::collectExogenous(set<pair<int, int> > &result) const
{
if (type == eExogenous)
result.insert(make_pair(symb_id, lag));
result.insert(make_pair(datatree.symbol_table.getTypeSpecificID(symb_id), lag));
else if (type == eModelLocalVariable)
datatree.local_variables_table[symb_id]->collectExogenous(result);
}
pair<bool, NodeID>
VariableNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
if (type ==eEndogenous)
{
if (datatree.symbol_table.getTypeSpecificID(symb_id)==var_endo and lag==0)
return(make_pair(true, (NodeID)NULL));
else
return(make_pair(false, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), lag)));
}
else
{
if (type == eParameter)
return(make_pair(false, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), 0)));
else
return(make_pair(false, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), lag)));
}
}
template<class InputIterator, class T>
pair<InputIterator, int> find_r ( InputIterator first, InputIterator last, const T& value )
{
int i=0;
for (; first!=last; first++)
{
if ( *first==value ) break;
i++;
}
return make_pair(first, i);
}
NodeID
VariableNode::computeChaineRuleDerivative(int deriv_id_arg, map<int, NodeID> &recursive_variables, int var, int lag_)
{
switch (type)
{
case eEndogenous:
case eExogenous:
case eExogenousDet:
case eParameter:
//cout << "deriv_id=" << deriv_id << " deriv_id_arg=" << deriv_id_arg << " symb_id=" << symb_id << " type=" << type << " lag=" << lag << " var=" << var << " lag_ = " << lag_ << "\n";
if (deriv_id == deriv_id_arg)
return datatree.One;
else
{
//if there is in the equation a recursive variable we could use a chaine rule derivation
if(lag == lag_)
{
map<int, NodeID>::const_iterator it = recursive_variables.find(deriv_id);
if (it != recursive_variables.end())
{
recursive_variables.erase(it->first);
return datatree.AddUMinus(it->second->getChaineRuleDerivative(deriv_id_arg, recursive_variables, var, lag_));
}
else
return datatree.Zero;
}
else
return datatree.Zero;
}
case eModelLocalVariable:
return datatree.local_variables_table[symb_id]->getChaineRuleDerivative(deriv_id_arg, recursive_variables, var, lag_);
case eModFileLocalVariable:
cerr << "ModFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
case eUnknownFunction:
cerr << "Impossible case!" << endl;
exit(EXIT_FAILURE);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
}
NodeID
VariableNode::toStatic(DataTree &static_datatree) const
{
@ -739,6 +885,7 @@ void
UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
//cout << "writeOutput unary\n";
// If node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
if (it != temporary_terms.end())
@ -918,6 +1065,127 @@ UnaryOpNode::collectExogenous(set<pair<int, int> > &result) const
arg->collectExogenous(result);
}
pair<bool, NodeID>
UnaryOpNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
pair<bool, NodeID> res = arg->normalizeLinearInEndoEquation(var_endo, Derivative);
bool is_endogenous_present = res.first;
NodeID New_NodeID = res.second;
if (not is_endogenous_present)
{
switch (op_code)
{
case oUminus:
return(make_pair(false, /*tmp_*/datatree.AddUMinus(New_NodeID)));
case oExp:
return(make_pair(false, /*tmp_*/datatree.AddExp(New_NodeID)));
case oLog:
return(make_pair(false, /*tmp_*/datatree.AddLog(New_NodeID)));
case oLog10:
return(make_pair(false, /*tmp_*/datatree.AddLog10(New_NodeID)));
case oCos:
return(make_pair(false, /*tmp_*/datatree.AddCos(New_NodeID)));
case oSin:
return(make_pair(false, /*tmp_*/datatree.AddSin(New_NodeID)));
case oTan:
return(make_pair(false, /*tmp_*/datatree.AddTan(New_NodeID)));
case oAcos:
return(make_pair(false, /*tmp_*/datatree.AddAcos(New_NodeID)));
case oAsin:
return(make_pair(false, /*tmp_*/datatree.AddAsin(New_NodeID)));
case oAtan:
return(make_pair(false, /*tmp_*/datatree.AddAtan(New_NodeID)));
case oCosh:
return(make_pair(false, /*tmp_*/datatree.AddCosh(New_NodeID)));
case oSinh:
return(make_pair(false, /*tmp_*/datatree.AddSinh(New_NodeID)));
case oTanh:
return(make_pair(false, /*tmp_*/datatree.AddTanh(New_NodeID)));
case oAcosh:
return(make_pair(false, /*tmp_*/datatree.AddAcosh(New_NodeID)));
case oAsinh:
return(make_pair(false, /*tmp_*/datatree.AddAsinh(New_NodeID)));
case oAtanh:
return(make_pair(false, /*tmp_*/datatree.AddAtanh(New_NodeID)));
case oSqrt:
return(make_pair(false, /*tmp_*/datatree.AddSqrt(New_NodeID)));
}
}
return(make_pair(true, (NodeID)NULL));
}
NodeID
UnaryOpNode::computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_)
{
NodeID darg = arg->getChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
NodeID t11, t12, t13;
switch (op_code)
{
case oUminus:
return datatree.AddUMinus(darg);
case oExp:
return datatree.AddTimes(darg, this);
case oLog:
return datatree.AddDivide(darg, arg);
case oLog10:
t11 = datatree.AddExp(datatree.One);
t12 = datatree.AddLog10(t11);
t13 = datatree.AddDivide(darg, arg);
return datatree.AddTimes(t12, t13);
case oCos:
t11 = datatree.AddSin(arg);
t12 = datatree.AddUMinus(t11);
return datatree.AddTimes(darg, t12);
case oSin:
t11 = datatree.AddCos(arg);
return datatree.AddTimes(darg, t11);
case oTan:
t11 = datatree.AddTimes(this, this);
t12 = datatree.AddPlus(t11, datatree.One);
return datatree.AddTimes(darg, t12);
case oAcos:
t11 = datatree.AddSin(this);
t12 = datatree.AddDivide(darg, t11);
return datatree.AddUMinus(t12);
case oAsin:
t11 = datatree.AddCos(this);
return datatree.AddDivide(darg, t11);
case oAtan:
t11 = datatree.AddTimes(arg, arg);
t12 = datatree.AddPlus(datatree.One, t11);
return datatree.AddDivide(darg, t12);
case oCosh:
t11 = datatree.AddSinh(arg);
return datatree.AddTimes(darg, t11);
case oSinh:
t11 = datatree.AddCosh(arg);
return datatree.AddTimes(darg, t11);
case oTanh:
t11 = datatree.AddTimes(this, this);
t12 = datatree.AddMinus(datatree.One, t11);
return datatree.AddTimes(darg, t12);
case oAcosh:
t11 = datatree.AddSinh(this);
return datatree.AddDivide(darg, t11);
case oAsinh:
t11 = datatree.AddCosh(this);
return datatree.AddDivide(darg, t11);
case oAtanh:
t11 = datatree.AddTimes(arg, arg);
t12 = datatree.AddMinus(datatree.One, t11);
return datatree.AddTimes(darg, t12);
case oSqrt:
t11 = datatree.AddPlus(this, this);
return datatree.AddDivide(darg, t11);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
}
NodeID
UnaryOpNode::toStatic(DataTree &static_datatree) const
{
@ -1001,11 +1269,16 @@ BinaryOpNode::computeDerivative(int deriv_id)
t12 = datatree.AddTimes(darg2, arg1);
return datatree.AddPlus(t11, t12);
case oDivide:
if (darg2!=datatree.Zero)
{
t11 = datatree.AddTimes(darg1, arg2);
t12 = datatree.AddTimes(darg2, arg1);
t13 = datatree.AddMinus(t11, t12);
t14 = datatree.AddTimes(arg2, arg2);
return datatree.AddDivide(t13, t14);
}
else
return datatree.AddDivide(darg1, arg2);
case oLess:
case oGreater:
case oLessEqual:
@ -1303,6 +1576,7 @@ void
BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
//cout << "writeOutput binary\n";
// If current node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
if (it != temporary_terms.end())
@ -1474,6 +1748,213 @@ BinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
arg2->collectExogenous(result);
}
pair<bool, NodeID>
BinaryOpNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
pair<bool, NodeID> res = arg1->normalizeLinearInEndoEquation(var_endo, Derivative);
bool is_endogenous_present_1 = res.first;
NodeID NodeID_1 = res.second;
res = arg2->normalizeLinearInEndoEquation(var_endo, Derivative);
bool is_endogenous_present_2 = res.first;
NodeID NodeID_2 = res.second;
switch (op_code)
{
case oPlus:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddPlus(NodeID_1, NodeID_2)));
else if (is_endogenous_present_1 and is_endogenous_present_2)
return(make_pair(true, (NodeID)NULL));
else if (not is_endogenous_present_1 and is_endogenous_present_2)
return(make_pair(false, NodeID_1));
else if (is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, NodeID_2));
break;
case oMinus:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddMinus(NodeID_1, NodeID_2)));
else if (is_endogenous_present_1 and is_endogenous_present_2)
return(make_pair(true, (NodeID)NULL));
else if (not is_endogenous_present_1 and is_endogenous_present_2)
return(make_pair(false, NodeID_1));
else if (is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddUMinus(NodeID_2)));
break;
case oTimes:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddTimes(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oDivide:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, datatree.AddDivide(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oPower:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, datatree.AddPower(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oEqual:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
{
if (Derivative!=datatree.One)
return( make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddDivide(datatree.AddMinus(NodeID_2, NodeID_1), Derivative))) );
else
return( make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddMinus(NodeID_2, NodeID_1))) );
}
else if (is_endogenous_present_1 and is_endogenous_present_2)
{
return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.Zero)));
}
else if (not is_endogenous_present_1 and is_endogenous_present_2)
{
if (Derivative!=datatree.One)
return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddDivide(datatree.AddUMinus(NodeID_1), Derivative))));
else
return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddUMinus(NodeID_1))));
}
else if (is_endogenous_present_1 and not is_endogenous_present_2)
{
if (Derivative!=datatree.One)
return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddDivide(NodeID_2, Derivative))));
else
return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), NodeID_2)));
}
break;
case oMax:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, datatree.AddMax(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oMin:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, datatree.AddMin(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oLess:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddLess(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oGreater:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddGreater(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oLessEqual:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddLessEqual(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oGreaterEqual:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddGreaterEqual(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oEqualEqual:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddEqualEqual(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
case oDifferent:
if (not is_endogenous_present_1 and not is_endogenous_present_2)
return(make_pair(false, /*tmp_*/datatree.AddDifferent(NodeID_1, NodeID_2)));
else
return(make_pair(true, (NodeID)NULL));
break;
}
// Suppress GCC warning
exit(EXIT_FAILURE);
}
NodeID
BinaryOpNode::computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_)
{
NodeID darg1 = arg1->getChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
NodeID darg2 = arg2->getChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
NodeID t11, t12, t13, t14, t15;
switch (op_code)
{
case oPlus:
return datatree.AddPlus(darg1, darg2);
case oMinus:
return datatree.AddMinus(darg1, darg2);
case oTimes:
t11 = datatree.AddTimes(darg1, arg2);
t12 = datatree.AddTimes(darg2, arg1);
return datatree.AddPlus(t11, t12);
case oDivide:
if (darg2!=datatree.Zero)
{
t11 = datatree.AddTimes(darg1, arg2);
t12 = datatree.AddTimes(darg2, arg1);
t13 = datatree.AddMinus(t11, t12);
t14 = datatree.AddTimes(arg2, arg2);
return datatree.AddDivide(t13, t14);
}
else
return datatree.AddDivide(darg1, arg2);
case oLess:
case oGreater:
case oLessEqual:
case oGreaterEqual:
case oEqualEqual:
case oDifferent:
return datatree.Zero;
case oPower:
if (darg2 == datatree.Zero)
{
if (darg1 == datatree.Zero)
return datatree.Zero;
else
{
t11 = datatree.AddMinus(arg2, datatree.One);
t12 = datatree.AddPower(arg1, t11);
t13 = datatree.AddTimes(arg2, t12);
return datatree.AddTimes(darg1, t13);
}
}
else
{
t11 = datatree.AddLog(arg1);
t12 = datatree.AddTimes(darg2, t11);
t13 = datatree.AddTimes(darg1, arg2);
t14 = datatree.AddDivide(t13, arg1);
t15 = datatree.AddPlus(t12, t14);
return datatree.AddTimes(t15, this);
}
case oMax:
t11 = datatree.AddGreater(arg1,arg2);
t12 = datatree.AddTimes(t11,darg1);
t13 = datatree.AddMinus(datatree.One,t11);
t14 = datatree.AddTimes(t13,darg2);
return datatree.AddPlus(t14,t12);
case oMin:
t11 = datatree.AddGreater(arg2,arg1);
t12 = datatree.AddTimes(t11,darg1);
t13 = datatree.AddMinus(datatree.One,t11);
t14 = datatree.AddTimes(t13,darg2);
return datatree.AddPlus(t14,t12);
case oEqual:
return datatree.AddMinus(darg1, darg2);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
}
NodeID
BinaryOpNode::toStatic(DataTree &static_datatree) const
{
@ -1796,6 +2277,75 @@ TrinaryOpNode::collectExogenous(set<pair<int, int> > &result) const
arg3->collectExogenous(result);
}
pair<bool, NodeID>
TrinaryOpNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
pair<bool, NodeID> res = arg1->normalizeLinearInEndoEquation(var_endo, Derivative);
bool is_endogenous_present_1 = res.first;
NodeID NodeID_1 = res.second;
res = arg2->normalizeLinearInEndoEquation(var_endo, Derivative);
bool is_endogenous_present_2 = res.first;
NodeID NodeID_2 = res.second;
res = arg3->normalizeLinearInEndoEquation(var_endo, Derivative);
bool is_endogenous_present_3 = res.first;
NodeID NodeID_3 = res.second;
if (not is_endogenous_present_1 and not is_endogenous_present_2 and not is_endogenous_present_3)
return(make_pair(false, /*tmp_*/datatree.AddNormcdf(NodeID_1, NodeID_2, NodeID_3)));
else
return(make_pair(true, (NodeID)NULL));
}
NodeID
TrinaryOpNode::computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_)
{
NodeID darg1 = arg1->getChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
NodeID darg2 = arg2->getChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
NodeID darg3 = arg3->getChaineRuleDerivative(deriv_id, recursive_variables, var, lag_);
NodeID t11, t12, t13, t14, t15;
switch (op_code)
{
case oNormcdf:
// normal pdf is inlined in the tree
NodeID y;
// sqrt(2*pi)
t14 = datatree.AddSqrt(datatree.AddTimes(datatree.Two, datatree.Pi));
// x - mu
t12 = datatree.AddMinus(arg1,arg2);
// y = (x-mu)/sigma
y = datatree.AddDivide(t12,arg3);
// (x-mu)^2/sigma^2
t12 = datatree.AddTimes(y,y);
// -(x-mu)^2/sigma^2
t13 = datatree.AddUMinus(t12);
// -((x-mu)^2/sigma^2)/2
t12 = datatree.AddDivide(t13, datatree.Two);
// exp(-((x-mu)^2/sigma^2)/2)
t13 = datatree.AddExp(t12);
// derivative of a standardized normal
// t15 = (1/sqrt(2*pi))*exp(-y^2/2)
t15 = datatree.AddDivide(t13,t14);
// derivatives thru x
t11 = datatree.AddDivide(darg1,arg3);
// derivatives thru mu
t12 = datatree.AddDivide(darg2,arg3);
// intermediary sum
t14 = datatree.AddMinus(t11,t12);
// derivatives thru sigma
t11 = datatree.AddDivide(y,arg3);
t12 = datatree.AddTimes(t11,darg3);
//intermediary sum
t11 = datatree.AddMinus(t14,t12);
// total derivative:
// (darg1/sigma - darg2/sigma - darg3*(x-mu)/sigma^2) * t15
// where t15 is the derivative of a standardized normal
return datatree.AddTimes(t11, t15);
}
// Suppress GCC warning
exit(EXIT_FAILURE);
}
NodeID
TrinaryOpNode::toStatic(DataTree &static_datatree) const
{
@ -1828,6 +2378,14 @@ UnknownFunctionNode::computeDerivative(int deriv_id)
exit(EXIT_FAILURE);
}
NodeID
UnknownFunctionNode::computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_)
{
cerr << "UnknownFunctionNode::computeDerivative: operation impossible!" << endl;
exit(EXIT_FAILURE);
}
void
UnknownFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
temporary_terms_type &temporary_terms,
@ -1907,6 +2465,25 @@ UnknownFunctionNode::compile(ofstream &CompileCode, bool lhs_rhs, const temporar
exit(EXIT_FAILURE);
}
pair<bool, NodeID>
UnknownFunctionNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
{
vector<pair<bool, NodeID> > V_arguments;
vector<NodeID> V_NodeID;
bool present = false;
for (vector<NodeID>::const_iterator it = arguments.begin();
it != arguments.end(); it++)
{
V_arguments.push_back((*it)->normalizeLinearInEndoEquation(var_endo, Derivative));
present = present or V_arguments[V_arguments.size()-1].first;
V_NodeID.push_back(V_arguments[V_arguments.size()-1].second);
}
if (not present)
return(make_pair(false, datatree.AddUnknownFunction(datatree.symbol_table.getName(symb_id), V_NodeID)));
else
return(make_pair(true, (NodeID)NULL));
}
NodeID
UnknownFunctionNode::toStatic(DataTree &static_datatree) const
{

30
ExprNode.hh
View File

@ -109,6 +109,10 @@ private:
//! Computes derivative w.r. to a derivation ID (but doesn't store it in derivatives map)
/*! You shoud use getDerivative() to get the benefit of symbolic a priori and of caching */
virtual NodeID computeDerivative(int deriv_id) = 0;
//! Computes derivative w.r. to a derivation ID and use chaine rule derivatives (but doesn't store it in derivatives map)
/*! You shoud use getDerivative() to get the benefit of symbolic a priori and of caching */
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_) = 0;
protected:
//! Reference to the enclosing DataTree
@ -136,6 +140,10 @@ public:
For an equal node, returns the derivative of lhs minus rhs */
NodeID getDerivative(int deriv_id);
//! Returns derivative w.r. to derivation ID and use if it possible chaine rule derivatives
NodeID getChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
//! Returns precedence of node
/*! Equals 100 for constants, variables, unary ops, and temporary terms */
virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
@ -145,7 +153,7 @@ public:
virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_type &temporary_terms, bool is_matlab) const;
//! Writes output of node, using a Txxx notation for nodes in temporary_terms
virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const = 0;
virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const /*= 0*/;
//! Writes output of node (with no temporary terms and with "outside model" output type)
void writeOutput(ostream &output);
@ -174,6 +182,7 @@ public:
map_idx_type &map_idx) const;
class EvalException
{
};
@ -185,6 +194,8 @@ public:
adds the result in the static_datatree argument (and not in the original datatree), and returns it.
*/
virtual NodeID toStatic(DataTree &static_datatree) const = 0;
//! Try to normalize an equation linear in its endogenous variable
virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! Object used to compare two nodes (using their indexes)
@ -204,6 +215,7 @@ private:
//! Id from numerical constants table
const int id;
virtual NodeID computeDerivative(int deriv_id);
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
public:
NumConstNode(DataTree &datatree_arg, int id_arg);
virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
@ -213,6 +225,7 @@ public:
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! Symbol or variable node
@ -226,6 +239,7 @@ private:
//! Derivation ID
const int deriv_id;
virtual NodeID computeDerivative(int deriv_id_arg);
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
public:
VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg, int deriv_id_arg);
virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms = temporary_terms_type()) const;
@ -243,6 +257,7 @@ public:
virtual void compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
int get_symb_id() const { return symb_id; };
virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! Unary operator node
@ -252,6 +267,7 @@ private:
const NodeID arg;
const UnaryOpcode op_code;
virtual NodeID computeDerivative(int deriv_id);
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const;
public:
@ -276,6 +292,7 @@ public:
//! Returns op code
UnaryOpcode get_op_code() const { return(op_code); };
virtual NodeID toStatic(DataTree &static_datatree) const;
virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! Binary operator node
@ -285,6 +302,8 @@ private:
const NodeID arg1, arg2;
const BinaryOpcode op_code;
virtual NodeID computeDerivative(int deriv_id);
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const;
public:
BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
@ -312,6 +331,7 @@ public:
//! Returns op code
BinaryOpcode get_op_code() const { return(op_code); };
virtual NodeID toStatic(DataTree &static_datatree) const;
pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! Trinary operator node
@ -322,6 +342,8 @@ private:
const NodeID arg1, arg2, arg3;
const TrinaryOpcode op_code;
virtual NodeID computeDerivative(int deriv_id);
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
virtual int cost(const temporary_terms_type &temporary_terms, bool is_matlab) const;
public:
TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
@ -343,6 +365,7 @@ public:
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! Unknown function node
@ -352,6 +375,7 @@ private:
const int symb_id;
const vector<NodeID> arguments;
virtual NodeID computeDerivative(int deriv_id);
virtual NodeID computeChaineRuleDerivative(int deriv_id, map<int, NodeID> &recursive_variables, int var, int lag_);
public:
UnknownFunctionNode(DataTree &datatree_arg, int symb_id_arg,
const vector<NodeID> &arguments_arg);
@ -370,6 +394,7 @@ public:
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
virtual void compile(ofstream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
};
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
@ -393,7 +418,8 @@ struct Block
bool is_linear;
int *Equation, *Own_Derivative;
EquationType *Equation_Type;
int *Variable, *Other_Endogenous, *Exogenous, *Equation_Type_Var;
NodeID *Equation_Normalized;
int *Variable, *Other_Endogenous, *Exogenous;
temporary_terms_type **Temporary_Terms_in_Equation;
//temporary_terms_type *Temporary_terms;
temporary_terms_inuse_type *Temporary_InUse;

123
MinimumFeedbackSet.cc
View File

@ -186,7 +186,6 @@ namespace MFS
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);
@ -208,13 +207,7 @@ namespace MFS
{
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;
@ -367,110 +360,8 @@ namespace MFS
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)
@ -537,21 +428,9 @@ namespace MFS
#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";
cout << "3 something_has_been_done=" << something_has_been_done << "\n";
#endif
}
vector<int> circuit;

31
MinimumFeedbackSet.hh
View File

@ -25,8 +25,6 @@
#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;
@ -44,25 +42,42 @@ using namespace boost;
namespace MFS
{
//! Eliminate a vertex i:
//! For a vertex i replace all edges e_k_i and e_i_j by a shorcut e_k_j and then Suppress the vertex i
void Eliminate(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, 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
vector_vertex_descriptor Collect_Doublet(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
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);
//! Graph reduction: eliminating purely intermediate variables or variables outside of any circuit
bool Elimination_of_Vertex_With_One_or_Less_Indegree_or_Outdegree_Step(AdjacencyList_type& G);
//! Graphe reduction: elimination of a vertex inside a clique
bool Elimination_of_Vertex_belonging_to_a_clique_Step(AdjacencyList_type& G);
//! A vertex belong to the feedback vertex set if the vertex loop on itself.
//! We have to suppress this vertex and store it into the feedback set.
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);
//! Print the Graph
void Print(GraphvizDigraph& G);
void Print(AdjacencyList_type& G);
//! Create a GraphvizDigraph from a Adjacency Matrix (an incidence Matrix without the diagonal terms)
GraphvizDigraph AM_2_GraphvizDigraph(bool* AM, unsigned int n);
//! Create an adjacency graph from a Adjacency Matrix (an incidence Matrix without the diagonal terms)
AdjacencyList_type AM_2_AdjacencyList(bool* AMp,unsigned int n);
//! Create an adjacency graph from a GraphvizDigraph
AdjacencyList_type GraphvizDigraph_2_AdjacencyList(GraphvizDigraph& G1, set<int> select_index);
//! Check if the graph contains any cycle (true if the model contains at least one cycle, false otherwise)
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);
//! Return the feedback set
AdjacencyList_type Minimal_set_of_feedback_vertex(set<int> &feed_back_vertices, const AdjacencyList_type& G);
//! clear all in and out edges of vertex_to_eliminate
//! and remove vertex_to_eliminate from the graph
void Suppress(AdjacencyList_type::vertex_descriptor vertex_to_eliminate, AdjacencyList_type& G);
void Suppress(int vertex_num, AdjacencyList_type& G);
//! reorder the recursive variable:
//! They appear first in a quasi triangular form and they are followed by the feedback variables
vector<int> Reorder_the_recursive_variables(const AdjacencyList_type& G1, set<int> &feed_back_vertices);
};

479
ModelGraph.cc
View File

@ -1,479 +0,0 @@
/*
* Copyright (C) 2007-2008 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <string>
#include <ctime>
#include <stack>
#include <cmath>
#include "ModelTree.hh"
#include "ModelGraph.hh"
#include "BlockTriangular.hh"
using namespace std;
void
free_model_graph(t_model_graph* model_graph)
{
int i;
for(i = 0;i < model_graph->nb_vertices;i++)
{
free(model_graph->vertex[i].in_degree_edge);
free(model_graph->vertex[i].out_degree_edge);
}
free(model_graph->vertex);
free(model_graph);
}
void
print_Graph(t_model_graph* model_graph)
{
int i, j;
for(i = 0;i < model_graph->nb_vertices;i++)
{
cout << "vertex " << model_graph->vertex[i].index << "(" << i << " ," << model_graph->vertex[i].nb_out_degree_edges << ")\n";
cout << " -> ";
for(j = 0;j < model_graph->vertex[i].nb_out_degree_edges;j++)
cout << model_graph->vertex[model_graph->vertex[i].out_degree_edge[j].index].index << /*" -" << model_graph->vertex[i].out_degree_edge[j].index << "-*/" (" << model_graph->vertex[i].out_degree_edge[j].u_count << "), ";
cout << "\n";
cout << " <- ";
for(j = 0;j < model_graph->vertex[i].nb_in_degree_edges;j++)
cout << model_graph->vertex[model_graph->vertex[i].in_degree_edge[j].index].index << /*" -" << model_graph->vertex[i].in_degree_edge[j].index << "-*/" (" << model_graph->vertex[i].in_degree_edge[j].u_count << "), ";
cout << "\n";
}
}
void Check_Graph(t_model_graph* model_graph)
{
int i, j, k, i1, i2;
bool OK, OK_u_count;
for(i = 0;i < model_graph->nb_vertices;i++)
{
for(j = 0;j < model_graph->vertex[i].nb_in_degree_edges;j++)
{
i1 = model_graph->vertex[i].in_degree_edge[j].index;
i2 = model_graph->vertex[i].in_degree_edge[j].u_count;
OK = 0;
OK_u_count = 0;
for(k = 0;(k < model_graph->vertex[i1].nb_out_degree_edges) && (!OK);k++)
{
if(model_graph->vertex[i1].out_degree_edge[k].index == i)
{
OK = 1;
if(model_graph->vertex[i1].out_degree_edge[k].u_count == i2)
OK_u_count = 1;
}
}
if(!OK)
{
cout << "not symetric for edge between vertices " << model_graph->vertex[i1].index << " and " << model_graph->vertex[i].index << " (in_degree)\n";
print_Graph(model_graph);
system("pause");
exit(EXIT_FAILURE);
}
if(!OK_u_count)
{
cout << "valeur de u_count non symétrique sur l'arc entre " << model_graph->vertex[i1].index << " et " << model_graph->vertex[i].index << " (in_degree)\n";
print_Graph(model_graph);
system("pause");
exit(EXIT_FAILURE);
}
}
for(j = 0;j < model_graph->vertex[i].nb_out_degree_edges;j++)
{
i1 = model_graph->vertex[i].out_degree_edge[j].index;
i2 = model_graph->vertex[i].out_degree_edge[j].u_count;
OK = 0;
OK_u_count = 0;
for(k = 0;(k < model_graph->vertex[i1].nb_in_degree_edges) && (!OK);k++)
{
if(model_graph->vertex[i1].in_degree_edge[k].index == i)
{
OK = 1;
if(model_graph->vertex[i1].in_degree_edge[k].u_count == i2)
OK_u_count = 1;
}
}
if(!OK)
{
cout << "pas symétrique sur l'arc entre " << model_graph->vertex[i1].index << " et " << model_graph->vertex[i].index << " (out_degree)\n";
print_Graph(model_graph);
system("pause");
exit(EXIT_FAILURE);
}
if(!OK_u_count)
{
cout << "valeur de u_count non symétrique sur l'arc entre " << model_graph->vertex[i1].index << " et " << model_graph->vertex[i].index << " (out_degree)\n";
print_Graph(model_graph);
system("pause");
exit(EXIT_FAILURE);
}
}
}
}
int
ModelBlock_Graph(Model_Block *ModelBlock, int Blck_num, bool dynamic, t_model_graph* model_graph, int nb_endo, int* block_u_count, int *starting_vertex, int *periods, int *nb_table_y, int *mean_var_in_equ)
{
int i, j, k, l, m, lag, per, lag1, k2, complete_size = 0, u_count;
int max_lead, max_lag, size, Lead, Lag;
int *Used, *todo_lag, *todo_lead, *vertex_ref, *vertex_index, *todo_lag1, *todo_lead1 ;
max_lag = ModelBlock->Block_List[Blck_num].Max_Lag;
max_lead = ModelBlock->Block_List[Blck_num].Max_Lead;
if(!dynamic)
{
/*It's a static model that have to be solved at each period*/
/*size=ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].size;*/
size = ModelBlock->Block_List[Blck_num].Size;
/*We add an extra vertex to take into account of the f(x0) constant term in f(x)=0 approximated by f(x0) + (x-x0) f'(x0) = 0*/
//cout << "Static, Blck_num= " << Blck_num << "size= " << size << "\n";
model_graph->nb_vertices = size + 1;
*starting_vertex = 0;
model_graph->vertex = (t_vertex*)malloc(model_graph->nb_vertices * sizeof(*model_graph->vertex));
for(i = 0;i < size;i++)
{
/*It's not f(x0) vertex*/
model_graph->vertex[i].in_degree_edge = (t_edge*)malloc((size + 1) * sizeof(t_edge));
model_graph->vertex[i].out_degree_edge = (t_edge*)malloc((size + 1) * sizeof(t_edge));
model_graph->vertex[i].nb_in_degree_edges = 0;
model_graph->vertex[i].nb_out_degree_edges = 0;
model_graph->vertex[i].index = ModelBlock->Block_List[Blck_num].Variable[i];
model_graph->vertex[i].lag_lead = 0;
}
/*It's f(x0) vertex*/
model_graph->vertex[size].in_degree_edge = (t_edge*)malloc(0 * sizeof(t_edge));
model_graph->vertex[size].out_degree_edge = (t_edge*)malloc((size) * sizeof(t_edge));
model_graph->vertex[size].nb_in_degree_edges = 0;
model_graph->vertex[size].index = -1;
model_graph->vertex[size].lag_lead = 0;
for(i = 0;i < ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].size;i++)
{
k = ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].Equ[i];
m = ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].Var[i];
j = model_graph->vertex[k].nb_in_degree_edges++;
l = model_graph->vertex[m].nb_out_degree_edges++;
model_graph->vertex[k].in_degree_edge[j].index = m;
model_graph->vertex[m].out_degree_edge[l].index = k;
model_graph->vertex[k].in_degree_edge[j].u_count = ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].us[i];
model_graph->vertex[m].out_degree_edge[l].u_count = ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].us[i];
}
model_graph->vertex[size].nb_out_degree_edges = size;
for(i = 0;i < size;i++)
{
j = model_graph->vertex[i].nb_in_degree_edges++;
model_graph->vertex[i].in_degree_edge[j].index = size;
model_graph->vertex[i].in_degree_edge[j].u_count = i;
model_graph->vertex[size].out_degree_edge[i].index = i;
model_graph->vertex[size].out_degree_edge[i].u_count = i;
}
u_count = ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].u_finish - ModelBlock->Block_List[Blck_num].IM_lead_lag[max_lag].u_init + 1
+ ModelBlock->Block_List[Blck_num].Size;
*block_u_count = u_count;
*nb_table_y = size;
return (u_count);
}
else
{
int sup;
Lead = ModelBlock->Block_List[Blck_num].Max_Lead;
Lag = ModelBlock->Block_List[Blck_num].Max_Lag;
cout << "---> *periods=" << *periods << "\n";
if(*periods>3)
{
sup = Lead + Lag +3;
*periods = Lead + Lag + sup;
}
#ifdef PRINT_OUT
cout << "Lag=" << Lag << " Lead=" << Lead << "\n";
cout << "periods=Lead+2*Lag+2= " << *periods << "\n";
#endif
size = ModelBlock->Block_List[Blck_num].Size;
/*It's a dynamic model that have to be solved for all periods.
So we consider the incidence matrice for all lead and lags plus the current value*/
model_graph->nb_vertices = 0;
vertex_ref = (int*)malloc(size * (Lag + Lead + *periods) * sizeof(int));
memset(vertex_ref, -1, size*(Lag + Lead + *periods)*sizeof(int));
vertex_index = (int*)malloc(size * (Lag + Lead + *periods) * sizeof(int));
complete_size = ModelBlock->Block_List[Blck_num].IM_lead_lag[Lag].size * (*periods);
if(Lag > 0)
{
todo_lag = (int*)malloc(size * Lag * sizeof(int));
todo_lag1 = (int*)malloc(size * Lag * sizeof(int));
memset(todo_lag, -1, size*Lag*sizeof(int));
memset(todo_lag1, -1, size*Lag*sizeof(int));
Used = (int*)malloc(size * Lag * sizeof(int));
for(lag = 0;lag < Lag;lag++)
{
memset(Used, -1, size*Lag*sizeof(int));
complete_size += ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size;
for(i = 0;i < ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size;i++)
{
if(Used[ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i]] < 0)
{
k = ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i];
todo_lag[lag*size + k] = k;
vertex_ref[lag*size + k] = model_graph->nb_vertices;
vertex_index[model_graph->nb_vertices] = lag * nb_endo + ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var_Index[i];
todo_lag1[lag*size + k] = i;
model_graph->nb_vertices++;
Used[k] = i;
}
}
if(lag > 0)
{
for(lag1 = 0;lag1 < lag;lag1++)
for(i = 0;i < size;i++)
if(todo_lag[(lag1)*size + i] >= 0)
{
if(Used[i] < 0)
{
todo_lag[lag*size + i] = i;
k = todo_lag[(lag1) * size + i];
vertex_ref[lag*size + k] = model_graph->nb_vertices;
j = todo_lag1[(lag1) * size + i];
vertex_index[model_graph->nb_vertices] = lag * nb_endo + ModelBlock->Block_List[Blck_num].IM_lead_lag[lag1].Var_Index[k];
model_graph->nb_vertices++;
}
}
}
}
*starting_vertex = model_graph->nb_vertices;
free(Used);
free(todo_lag);
free(todo_lag1);
}
int nb_vertices_1=model_graph->nb_vertices;
#ifdef PRINT_OUT
cout << "nb_vertices in the first part: " << nb_vertices_1 << "\n";
#endif
for(per = Lag;per < Lag + *periods;per++)
for(i = 0;i < size;i++)
{
vertex_ref[per*size + i] = model_graph->nb_vertices;
vertex_index[model_graph->nb_vertices] = (per) * nb_endo + ModelBlock->Block_List[Blck_num].Variable[i];
model_graph->nb_vertices++;
}
int nb_vertices_2=model_graph->nb_vertices-nb_vertices_1;
#ifdef PRINT_OUT
cout << "nb_vertices in the second part: " << nb_vertices_2 << "\n";
#endif
if(Lead > 0)
{
todo_lead = (int*)malloc(size * Lead * sizeof(int));
todo_lead1 = (int*)malloc(size * Lead * sizeof(int));
memset(todo_lead, -1, size*Lead*sizeof(int));
memset(todo_lead1, -1, size*Lead*sizeof(int));
Used = (int*)malloc(size * Lead * sizeof(int));
k2 = model_graph->nb_vertices;
for(lag = Lag + Lead;lag > Lag;lag--)
{
complete_size += ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size;
memset(Used, -1, size*Lead*sizeof(int));
for(i = 0;i < ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size;i++)
{
if(Used[ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i]] < 0)
{
k = ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i];
todo_lead[(lag - Lag - 1)*size + k] = k;
todo_lead1[(lag - Lag - 1)*size + k] = i;
Used[k] = i;
model_graph->nb_vertices++;
}
}
if(lag < Lag + Lead)
{
for(lag1 = Lag + Lead;lag1 > lag;lag1--)
for(i = 0;i < size;i++)
{
if(todo_lead[(lag1 - Lag - 1)*size + i] >= 0)
{
if(Used[i] < 0)
{
k = todo_lead[(lag1 - Lag - 1) * size + i];
model_graph->nb_vertices++;
}
}
}
}
}
k2 = model_graph->nb_vertices;
memset(todo_lead, -1, size*Lead*sizeof(int));
for(lag = Lag + Lead;lag > Lag;lag--)
{
complete_size += ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size;
memset(Used, -1, size*Lead*sizeof(int));
for(i = ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size - 1;i >= 0;i--)
{
if(Used[ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i]] < 0)
{
k2--;
k = ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i];
todo_lead[(lag - Lag - 1)*size + k] = k;
todo_lead1[(lag - Lag - 1)*size + k] = i;
vertex_ref[(lag + *periods - 1)*size + k] = k2;
vertex_index[k2] = (lag + *periods - 1) * nb_endo + ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var_Index[i];
Used[k] = i;
}
}
if(lag < Lag + Lead)
{
for(lag1 = Lag + Lead;lag1 > lag;lag1--)
{
for(i = size - 1;i >= 0;i--)
{
if(todo_lead[(lag1 - Lag - 1)*size + i] >= 0)
{
if(Used[i] < 0)
{
k2--;
todo_lead[(lag - Lag - 1)*size + i] = i;
todo_lead1[(lag - Lag - 1)*size + i] = todo_lead1[(lag1 - Lag - 1)*size + i];
k = todo_lead[(lag1 - Lag - 1) * size + i];
vertex_ref[(lag + *periods - 1)*size + k] = k2;
//#ifdef PRINT_OUT
//#endif
j = todo_lead1[(lag1 - Lag - 1) * size + i];
//#ifdef PRINT_OUT
if(j>ModelBlock->Block_List[Blck_num].IM_lead_lag[lag1].size||j==-1)
{
cout << "Error in model graph construction (lead part): j (" << j << ")>size (" << ModelBlock->Block_List[Blck_num].IM_lead_lag[lag1].size << ")\n";
system("pause");
exit(EXIT_FAILURE);
}
//#endif
vertex_index[k2] = (lag + *periods - 1) * nb_endo + ModelBlock->Block_List[Blck_num].IM_lead_lag[lag1].Var_Index[j];
}
}
}
}
}
}
free(Used);
free(todo_lead);
free(todo_lead1);
}
int nb_vertices_3=model_graph->nb_vertices-nb_vertices_2-nb_vertices_1;
#ifdef PRINT_OUT
cout << "nb_vertices in the last part: " << nb_vertices_3 << "\n";
#endif
/*We add an extra vertex to take into account of the f(x0) constant term in f(x)=0 approx f(x0) + (x-x0) f'(x0) = 0*/
model_graph->nb_vertices++;
model_graph->vertex = (t_vertex*)malloc(model_graph->nb_vertices * sizeof(*model_graph->vertex));
vertex_index[model_graph->nb_vertices - 1] = -1;
#ifdef PRINT_OUT
cout << "ok0\n";
cout << "model_graph->nb_vertices=" << model_graph->nb_vertices << " Lag=" << Lag << " Lead=" << Lead << "\n";
cout << "*periods=" << *periods << " size=" << size << "\n";
cout << "allocated / vertex = " << (size + nb_vertices_1 + nb_vertices_3+ 1) << "\n";
#endif
int nb_table_u= size+nb_vertices_1+nb_vertices_3+2;
for(k = 0;k < model_graph->nb_vertices-1;k++)
{
model_graph->vertex[k].index = vertex_index[k];
model_graph->vertex[k].in_degree_edge = (t_edge*)malloc(nb_table_u * sizeof(t_edge));
model_graph->vertex[k].out_degree_edge = (t_edge*)malloc(nb_table_u * sizeof(t_edge));
model_graph->vertex[k].nb_in_degree_edges = 0;
model_graph->vertex[k].nb_out_degree_edges = 0;
model_graph->vertex[k].max_nb_in_degree_edges = nb_table_u;
model_graph->vertex[k].max_nb_out_degree_edges = nb_table_u;
#ifdef PRINT_OUT
//if(k==8)
{
cout << " model_graph->vertex[" << k << "].in_degree_edge=" << model_graph->vertex[k].in_degree_edge << "\n";
}
#endif
}
model_graph->vertex[model_graph->nb_vertices-1].index = vertex_index[model_graph->nb_vertices-1];
model_graph->vertex[model_graph->nb_vertices-1].in_degree_edge = (t_edge*)malloc(/*model_graph->nb_vertices **/ sizeof(t_edge));
model_graph->vertex[model_graph->nb_vertices-1].out_degree_edge = (t_edge*)malloc(model_graph->nb_vertices * sizeof(t_edge));
model_graph->vertex[model_graph->nb_vertices-1].nb_in_degree_edges = 0;
model_graph->vertex[model_graph->nb_vertices-1].nb_out_degree_edges = 0;
model_graph->vertex[model_graph->nb_vertices-1].max_nb_in_degree_edges = 0;
model_graph->vertex[model_graph->nb_vertices-1].max_nb_out_degree_edges = model_graph->nb_vertices;
#ifdef PRINT_OUT
cout << "ok1\n";
system("pause");
#endif
u_count = 0;
*mean_var_in_equ = 0;
for(per = 0;per < *periods;per++)
{
j = model_graph->nb_vertices - 1;
for(i = 0;i < size;i++)
{
k = vertex_ref[(Lag + per) * size + i];
model_graph->vertex[k].in_degree_edge[model_graph->vertex[k].nb_in_degree_edges].index = j;
model_graph->vertex[j].out_degree_edge[model_graph->vertex[j].nb_out_degree_edges].index = k;
model_graph->vertex[k].in_degree_edge[model_graph->vertex[k].nb_in_degree_edges].u_count = u_count;
model_graph->vertex[j].out_degree_edge[model_graph->vertex[j].nb_out_degree_edges].u_count = u_count;
model_graph->vertex[k].nb_in_degree_edges++;
model_graph->vertex[j].nb_out_degree_edges++;
u_count++;
}
for(lag = 0;lag < Lag + Lead + 1;lag++)
{
#ifdef PRINT_OUT
cout << "ModelBlock->Block_List[" << Blck_num << "].IM_lead_lag[" << lag << "].size = " << ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size << "\n";
#endif
for(i = 0;i < ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].size;i++)
{
j = vertex_ref[(lag + per) * size + ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Var[i]];
k = vertex_ref[(Lag + per) * size + ModelBlock->Block_List[Blck_num].IM_lead_lag[lag].Equ[i]];
#ifdef PRINT_OUT
cout << "per=" << per << " lag=" << lag << " i=" << i << " j=" << j << " k=" << k << "\n";
#endif
model_graph->vertex[k].in_degree_edge[model_graph->vertex[k].nb_in_degree_edges].index = j;
model_graph->vertex[j].out_degree_edge[model_graph->vertex[j].nb_out_degree_edges].index = k;
model_graph->vertex[k].in_degree_edge[model_graph->vertex[k].nb_in_degree_edges].u_count = u_count;
model_graph->vertex[j].out_degree_edge[model_graph->vertex[j].nb_out_degree_edges].u_count = u_count;
if(per==(Lag+2))/*&&(lag==(Lag+1))*/
(*mean_var_in_equ)++;
model_graph->vertex[k].nb_in_degree_edges++;
model_graph->vertex[j].nb_out_degree_edges++;
u_count++;
}
}
}
(*mean_var_in_equ) += size;
//cout << "Total variables=" << *mean_var_in_equ << " nb_endo=" << size << "\n";
i=*mean_var_in_equ ;
i =int(ceil(double(i)/size));
*mean_var_in_equ = i;
//cout << "Mean var in equation=" << *mean_var_in_equ << "\n";
*block_u_count = u_count / (*periods);
free(vertex_index);
free(vertex_ref);
if(nb_vertices_1+nb_vertices_3+1>size)
*nb_table_y = nb_vertices_1+nb_vertices_3+1;
else
*nb_table_y = size;
return (u_count);
}
}

57
ModelGraph.hh
View File

@ -1,57 +0,0 @@
/*
* Copyright (C) 2007-2008 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MODEL_GRAPH
#define MODEL_GRAPH
#define DIRECT_COMPUTE
#define SORTED
#define SIMPLIFY
#define SIMPLIFYS
#define SAVE
#define COMPUTE
//#define PRINT_OUT_OUT
//#define PRINT_OUT_1
#define DIRECT_SAVE
#include "ModelTree.hh"
#include "BlockTriangular.hh"
struct t_edge
{
int index, u_count;
};
struct t_vertex
{
t_edge *out_degree_edge, *in_degree_edge;
int nb_out_degree_edges, nb_in_degree_edges;
int max_nb_out_degree_edges, max_nb_in_degree_edges;
int index, lag_lead;
};
struct t_model_graph
{
int nb_vertices;
t_vertex* vertex;
};
void free_model_graph(t_model_graph* model_graph);
void print_Graph(t_model_graph* model_graph);
void Check_Graph(t_model_graph* model_graph);
int ModelBlock_Graph(Model_Block *ModelBlock, int Blck_num,bool dynamic, t_model_graph* model_graph, int nb_endo, int *block_u_count, int *starting_vertex, int* periods, int *nb_table_y, int *mean_var_in_equ);
#endif