diff --git a/BlockTriangular.cc b/BlockTriangular.cc
index 63fc314c..d7119571 100644
--- a/BlockTriangular.cc
+++ b/BlockTriangular.cc
@@ -17,7 +17,6 @@
* along with Dynare. If not, see .
*/
-
#include
#include
#include
@@ -38,25 +37,27 @@ 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
void
-BlockTriangular::Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM0)
+BlockTriangular::Prologue_Epilogue(bool *IM, int &prologue, int &epilogue, int n, vector &Index_Var_IM, vector &Index_Equ_IM, bool *IM0)
{
bool modifie = 1;
int i, j, k, l = 0;
@@ -65,10 +66,10 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n
while (modifie)
{
modifie = 0;
- for (i = prologue;i < n;i++)
+ for (i = prologue; i < n; i++)
{
k = 0;
- for (j = prologue;j < n;j++)
+ for (j = prologue; j < n; j++)
{
if (IM[i*n + j])
{
@@ -89,10 +90,10 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n
while (modifie)
{
modifie = 0;
- for (i = prologue;i < n - epilogue;i++)
+ for (i = prologue; i < n - epilogue; i++)
{
k = 0;
- for (j = prologue;j < n - epilogue;j++)
+ for (j = prologue; j < n - epilogue; j++)
{
if (IM[j*n + i])
{
@@ -110,148 +111,199 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n
}
}
-
//------------------------------------------------------------------------------
// Find a matching between equations and endogenous variables
bool
BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector &Index_Equ_IM) const
- {
- int n = equation_number - prologue - epilogue;
+{
+ int n = equation_number - prologue - epilogue;
+ typedef adjacency_list BipartiteGraph;
- typedef adjacency_list BipartiteGraph;
+ /*
+ Vertices 0 to n-1 are for endogenous (using type specific ID)
+ Vertices n to 2*n-1 are for equations (using equation no.)
+ */
+ BipartiteGraph g(2 * n);
+ // Fill in the graph
+ for (int i = 0; i < n; i++)
+ for (int j = 0; j < n; j++)
+ if (IM0[(i+prologue) * equation_number+j+prologue])
+ add_edge(i + n, j, g);
- /*
- Vertices 0 to n-1 are for endogenous (using type specific ID)
- Vertices n to 2*n-1 are for equations (using equation no.)
- */
- BipartiteGraph g(2 * n);
- // Fill in the graph
- for (int i = 0; i < n; i++)
- for (int j = 0; j < n; j++)
- if (IM0[(i+prologue) * equation_number+j+prologue])
- add_edge(i + n, j, g);
+ // Compute maximum cardinality matching
+ typedef vector::vertex_descriptor> mate_map_t;
+ mate_map_t mate_map(2*n);
+ bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
+ //cout << "check = " << check << "\n";
+ if (check)
+ {
+ // Check if all variables are normalized
+ mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits::null_vertex());
+ if (it != mate_map.begin() + n)
+ {
+ if (verbose)
+ {
+ cerr << "ERROR: Could not normalize dynamic model. Variable "
+ << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
+ << " is not in the maximum cardinality matching." << endl;
+ exit(EXIT_FAILURE);
+ }
+ return false;
+ }
+ vector Index_Equ_IM_tmp(Index_Equ_IM);
+ bool *SIM;
+ SIM = (bool *) malloc(equation_number*equation_number*sizeof(bool));
+ memcpy(SIM, IM, equation_number*equation_number*sizeof(bool));
+ for (int i = 0; i < n; i++)
+ {
+ Index_Equ_IM[i + prologue] = Index_Equ_IM_tmp[mate_map[i] - n + prologue];
+ for (int k = 0; k < n; k++)
+ IM[(i+prologue)*equation_number+k +prologue] = SIM[(mate_map[i]-n+prologue)*equation_number+k +prologue];
+ }
+ free(SIM);
+ }
+ return check;
+}
- // Compute maximum cardinality matching
- typedef vector::vertex_descriptor> mate_map_t;
- mate_map_t mate_map(2*n);
-
- bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
- //cout << "check = " << check << "\n";
- if (check)
- {
- // Check if all variables are normalized
- mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits::null_vertex());
- if (it != mate_map.begin() + n)
- {
- if (verbose)
- {
- cerr << "ERROR: Could not normalize dynamic model. Variable "
- << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
- << " is not in the maximum cardinality matching." << endl;
- exit(EXIT_FAILURE);
- }
- return false;
- }
- vector Index_Equ_IM_tmp(Index_Equ_IM);
- bool *SIM;
- SIM = (bool*)malloc(equation_number*equation_number*sizeof(bool));
- memcpy(SIM, IM, equation_number*equation_number*sizeof(bool));
- for (int i = 0; i < n; i++)
- {
- Index_Equ_IM[i + prologue] = Index_Equ_IM_tmp[mate_map[i] - n + prologue];
- for (int k=0; k &components_set, int nb_blck_sim, int prologue, int epilogue) const
+{
+ int nb_endo = symbol_table.endo_nbr();
+ vector 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 &Index_Equ_IM, vector &Index_Var_IM, vector > &blocks, t_etype &Equation_Type, bool verbose_) const
- {
- int n = nb_var - prologue - epilogue;
- bool *AMp;
- AMp = (bool*) malloc(n*n*sizeof(bool));
- //transforms the incidence matrix of the complet model into an adjancent matrix of the non-recursive part of the model
- for (int i=prologue; i component(num_vertices(G2)), discover_time(num_vertices(G2));
+ //In a first step we compute the strong components of the graph representation of the static model.
+ // This insures that block are dynamically recursives.
+ GraphvizDigraph G2 = AM_2_GraphvizDigraph(AMp, n);
+ vector component(num_vertices(G2)), discover_time(num_vertices(G2));
- int num = strong_components(G2, &component[0]);
+ int num = strong_components(G2, &component[0]);
- blocks = vector >(num , make_pair(0,0));
+ blocks = vector >(num, make_pair(0, 0));
- //This vector contains for each block:
- // - first set = equations belonging to the block,
- // - second set = the feeback variables,
- // - third vector = the reordered non-feedback variables.
- vector, pair, vector > > > components_set(num);
+ //This vector contains for each block:
+ // - first set = equations belonging to the block,
+ // - second set = the feeback variables,
+ // - third vector = the reordered non-feedback variables.
+ vector, pair, vector > > > components_set(num);
- for (unsigned int i=0; i tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
- int order = prologue;
- bool *SIM;
- SIM = (bool*)malloc(nb_var*nb_var*sizeof(bool));
- memcpy(SIM, IM, nb_var*nb_var*sizeof(bool));
+ vector tmp_Index_Equ_IM(Index_Equ_IM), tmp_Index_Var_IM(Index_Var_IM);
+ int order = prologue;
+ bool *SIM;
+ SIM = (bool *) malloc(nb_var*nb_var*sizeof(bool));
+ memcpy(SIM, IM, nb_var*nb_var*sizeof(bool));
- //Add a loop on vertices which could not be normalized => force those vvertices to belong to the feedback set
- for(int i=0; i 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 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
- // and the non-feedback variables are reordered to get
- // a sub-recursive block without feedback variables
- for (int i=0; i feed_back_vertices;
- //Print(G);
- AdjacencyList_type G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
- property_map::type v_index = get(vertex_index, G);
- components_set[i].second.first = feed_back_vertices;
- blocks[i].second = feed_back_vertices.size();
- vector Reordered_Vertice;
- Reordered_Vertice = Reorder_the_recursive_variables(G, feed_back_vertices);
- //First we have the recursive equations conditional on feedback variables
- for (vector::iterator its=Reordered_Vertice.begin();its != Reordered_Vertice.end(); its++)
- {
- Index_Equ_IM[order] = tmp_Index_Equ_IM[*its+prologue];
- Index_Var_IM[order] = tmp_Index_Var_IM[*its+prologue];
- order++;
- }
- components_set[i].second.second = Reordered_Vertice;
- //Second we have the equations related to the feedback variables
- for (set::iterator its=feed_back_vertices.begin();its != feed_back_vertices.end(); its++)
- {
- Index_Equ_IM[order] = tmp_Index_Equ_IM[v_index[vertex(*its, G)]+prologue];
- Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue];
- order++;
- }
- }
- free(AMp);
- free(SIM);
- }
+ //For each block, the minimum set of feedback variable is computed
+ // and the non-feedback variables are reordered to get
+ // a sub-recursive block without feedback variables
+ for (int i = 0; i < num; i++)
+ {
+ AdjacencyList_type G = GraphvizDigraph_2_AdjacencyList(G2, components_set[i].first);
+ set feed_back_vertices;
+ //Print(G);
+ AdjacencyList_type G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
+ property_map::type v_index = get(vertex_index, G);
+ components_set[i].second.first = feed_back_vertices;
+ blocks[i].second = feed_back_vertices.size();
+ vector Reordered_Vertice;
+ Reordered_Vertice = Reorder_the_recursive_variables(G, feed_back_vertices);
+ //First we have the recursive equations conditional on feedback variables
+ for (vector::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
+ {
+ Index_Equ_IM[order] = tmp_Index_Equ_IM[*its+prologue];
+ Index_Var_IM[order] = tmp_Index_Var_IM[*its+prologue];
+ order++;
+ }
+ components_set[i].second.second = Reordered_Vertice;
+ //Second we have the equations related to the feedback variables
+ for (set::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
+ {
+ Index_Equ_IM[order] = tmp_Index_Equ_IM[v_index[vertex(*its, G)]+prologue];
+ Index_Var_IM[order] = tmp_Index_Var_IM[v_index[vertex(*its, G)]+prologue];
+ order++;
+ }
+ }
+ free(AMp);
+ free(SIM);
+}
void
-BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block * ModelBlock, t_etype &Equation_Type, int recurs_Size)
+BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size)
{
int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li;
int *tmp_size, *tmp_size_other_endo, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_other_endo, *tmp_exo, tmp_nb_other_endo, tmp_nb_exo, nb_lead_lag_endo;
@@ -260,36 +312,28 @@ 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].is_linear = true;
ModelBlock->Block_List[count_Block].Size = size;
ModelBlock->Block_List[count_Block].Type = type;
ModelBlock->Block_List[count_Block].Nb_Recursives = recurs_Size;
- ModelBlock->Block_List[count_Block].Temporary_InUse=new temporary_terms_inuse_type ();
+ ModelBlock->Block_List[count_Block].Temporary_InUse = new temporary_terms_inuse_type();
ModelBlock->Block_List[count_Block].Temporary_InUse->clear();
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].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));
+ 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_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));
Lead = Lag = 0;
first_count_equ = *count_Equ;
- tmp_var = (int*)malloc(size * sizeof(int));
- 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));
+ tmp_var = (int *) malloc(size * sizeof(int));
+ 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));
+ 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));
memset(tmp_size_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
@@ -299,18 +343,25 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
Lag_Endo = Lead_Endo = Lag_Other_Endo = Lead_Other_Endo = Lag_Exo = Lead_Exo = 0;
//Variable by variable looking for all leads and lags its occurence in each equation of the block
- tmp_variable_evaluated = (bool*)malloc(symbol_table.endo_nbr()*sizeof(bool));
+ tmp_variable_evaluated = (bool *) malloc(symbol_table.endo_nbr()*sizeof(bool));
memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr()*sizeof(bool));
- for (i = 0;i < size;i++)
+ for (i = 0; i < size; i++)
{
- ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]=new temporary_terms_type ();
+ ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i] = new temporary_terms_type();
ModelBlock->Block_List[count_Block].Temporary_Terms_in_Equation[i]->clear();
ModelBlock->Block_List[count_Block].Equation[i] = Index_Equ_IM[*count_Equ];
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++)
+ for (k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if (Cur_IM)
@@ -318,7 +369,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
OK = false;
if (k >= 0)
{
- for (j = 0;j < size;j++)
+ for (j = 0; j < size; j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j]*symbol_table.endo_nbr()])
{
@@ -337,7 +388,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
}
else
{
- for (j = 0;j < size;j++)
+ for (j = 0; j < size; j++)
{
if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j]*symbol_table.endo_nbr()])
{
@@ -363,15 +414,15 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
Lead_Endo = Lead;
tmp_nb_other_endo = 0;
- for (i = 0;i < size;i++)
+ for (i = 0; i < size; i++)
{
- for (k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
+ for (k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if (Cur_IM)
{
i_1 = Index_Equ_IM[first_count_equ+i] * symbol_table.endo_nbr();
- for (j = 0;j < symbol_table.endo_nbr();j++)
+ for (j = 0; j < symbol_table.endo_nbr(); j++)
if (Cur_IM[i_1 + j])
{
if (!tmp_variable_evaluated[j])
@@ -379,13 +430,13 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
tmp_other_endo[j] = 1;
tmp_nb_other_endo++;
}
- if (k>0 && k>Lead_Other_Endo)
+ if (k > 0 && k > Lead_Other_Endo)
Lead_Other_Endo = k;
- else if (k<0 && (-k)>Lag_Other_Endo)
+ else if (k < 0 && (-k) > Lag_Other_Endo)
Lag_Other_Endo = -k;
- if (k>0 && k>Lead)
+ if (k > 0 && k > Lead)
Lead = k;
- else if (k<0 && (-k)>Lag)
+ else if (k < 0 && (-k) > Lag)
Lag = -k;
tmp_size_other_endo[k+incidencematrix.Model_Max_Lag_Endo]++;
}
@@ -393,21 +444,20 @@ 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));
+ 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));
+ tmp_exo = (int *) malloc(symbol_table.exo_nbr() * sizeof(int));
memset(tmp_exo, 0, symbol_table.exo_nbr() * sizeof(int));
tmp_nb_exo = 0;
- for (i = 0;i < size;i++)
+ for (i = 0; i < size; i++)
{
- for (k = -incidencematrix.Model_Max_Lag_Exo; k<=incidencematrix.Model_Max_Lead_Exo; k++)
+ for (k = -incidencematrix.Model_Max_Lag_Exo; k <= incidencematrix.Model_Max_Lead_Exo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eExogenous);
if (Cur_IM)
{
i_1 = Index_Equ_IM[first_count_equ+i] * symbol_table.exo_nbr();
- for (j=0;j0 && k>Lead_Exo)
+ if (k > 0 && k > Lead_Exo)
Lead_Exo = k;
- else if (k<0 && (-k)>Lag_Exo)
+ else if (k < 0 && (-k) > Lag_Exo)
Lag_Exo = -k;
- if (k>0 && k>Lead)
+ if (k > 0 && k > Lead)
Lead = k;
- else if (k<0 && (-k)>Lag)
+ else if (k < 0 && (-k) > Lag)
Lag = -k;
tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++;
}
@@ -429,11 +479,10 @@ 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));
+ ModelBlock->Block_List[count_Block].Exogenous = (int *) malloc(tmp_nb_exo * sizeof(int));
k = 0;
- for (j=0;jBlock_List[count_Block].Exogenous[k] = j;
@@ -450,52 +499,52 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].Max_Lead_Other_Endo = Lead_Other_Endo;
ModelBlock->Block_List[count_Block].Max_Lag_Exo = Lag_Exo;
ModelBlock->Block_List[count_Block].Max_Lead_Exo = Lead_Exo;
- ModelBlock->Block_List[count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact));
+ ModelBlock->Block_List[count_Block].IM_lead_lag = (IM_compact *) malloc((Lead + Lag + 1) * sizeof(IM_compact));
ls = l = li = size;
i1 = 0;
ModelBlock->Block_List[count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo;
- for (i = 0;i < Lead + Lag + 1;i++)
+ for (i = 0; i < Lead + Lag + 1; i++)
{
- if (incidencematrix.Model_Max_Lag_Endo-Lag+i>=0)
+ if (incidencematrix.Model_Max_Lag_Endo-Lag+i >= 0)
{
ModelBlock->Block_List[count_Block].IM_lead_lag[i].size = tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i];
ModelBlock->Block_List[count_Block].IM_lead_lag[i].nb_endo = tmp_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].u = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].us = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].u = (int *) malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].us = (int *) malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var = (int *) malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ = (int *) malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index = (int *) malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index = (int *) malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_other_endo = tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
ModelBlock->Block_List[count_Block].IM_lead_lag[i].nb_other_endo = tmp_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i];
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_other_endo = (int *) malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_other_endo = (int *) malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_other_endo = (int *) malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index_other_endo = (int *) malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index_other_endo = (int *) malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int));
}
else
ModelBlock->Block_List[count_Block].IM_lead_lag[i].size = 0;
/*if (incidencematrix.Model_Max_Lag_Exo-Lag+i>=0)
- {
+ {
ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_exo = tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i];
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Exogenous = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Exogenous_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X_Index = (int*)malloc(tmp_size_exo[incidencematrix.Model_Max_Lag_Exo - Lag + i] * sizeof(int));
- }
- else
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_exo = 0;*/
+ }
+ else
+ ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_exo = 0;*/
ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_init = l;
memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr()*sizeof(bool));
IM = incidencematrix.Get_IM(i - Lag, eEndogenous);
if (IM)
{
- for (j = first_count_equ;j < size + first_count_equ;j++)
+ for (j = first_count_equ; j < size + first_count_equ; j++)
{
i_1 = Index_Var_IM[j];
m = 0;
- for (k = first_count_equ;k < size + first_count_equ;k++)
+ for (k = first_count_equ; k < size + first_count_equ; k++)
if (IM[i_1 + Index_Equ_IM[k] * symbol_table.endo_nbr()])
m++;
if (m > 0)
@@ -505,10 +554,10 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
}
}
m = 0;
- for (j = first_count_equ;j < size + first_count_equ;j++)
+ for (j = first_count_equ; j < size + first_count_equ; j++)
{
i_1 = Index_Equ_IM[j] * symbol_table.endo_nbr();
- for (k = first_count_equ;k < size + first_count_equ;k++)
+ for (k = first_count_equ; k < size + first_count_equ; k++)
if (IM[Index_Var_IM[k] + i_1])
{
if (i == Lag)
@@ -529,10 +578,10 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
}
ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_finish = li - 1;
m = 0;
- for (j = first_count_equ;j < size + first_count_equ;j++)
+ for (j = first_count_equ; j < size + first_count_equ; j++)
{
i_1 = Index_Equ_IM[j] * symbol_table.endo_nbr();
- for (k = 0;k < symbol_table.endo_nbr();k++)
+ for (k = 0; k < symbol_table.endo_nbr(); k++)
if ((!tmp_variable_evaluated[Index_Var_IM[k]]) && IM[Index_Var_IM[k] + i_1])
{
ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_other_endo[m] = l;
@@ -547,8 +596,8 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_other_endo = m;
}
/*IM = incidencematrix.Get_IM(i - Lag, eExogenous);
- if (IM)
- {
+ if (IM)
+ {
m = 0;
for (j = first_count_equ;j < size + first_count_equ;j++)
{
@@ -565,7 +614,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
}
}
}
- }*/
+ }*/
}
free(tmp_size);
free(tmp_size_other_endo);
@@ -577,55 +626,52 @@ 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);
- for (i = 0;i < ModelBlock->Block_List[blk].Max_Lag + ModelBlock->Block_List[blk].Max_Lead + 1;i++)
- {
- if (incidencematrix.Model_Max_Lag_Endo-ModelBlock->Block_List[blk].Max_Lag+i>=0 /*&& ModelBlock->Block_List[blk].IM_lead_lag[i].size*/)
- {
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].u_other_endo);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_other_endo);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_other_endo);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index_other_endo);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index_other_endo);
- }
- /*if (incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 )
- {
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index);
- free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X);
- }*/
- }
- free(ModelBlock->Block_List[blk].IM_lead_lag);
- for (i=0; iBlock_List[blk].Size; i++)
- delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
- free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
- delete(ModelBlock->Block_List[blk].Temporary_InUse);
- }
- free(ModelBlock->Block_List);
- free(ModelBlock);
- }
+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_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*/)
+ {
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].u);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].us);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].u_other_endo);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_other_endo);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_other_endo);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Var_Index_other_endo);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_Index_other_endo);
+ }
+ /*if (incidencematrix.Model_Max_Lag_Exo-ModelBlock->Block_List[blk].Max_Lag+i>=0 )
+ {
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Exogenous_Index);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X_Index);
+ free(ModelBlock->Block_List[blk].IM_lead_lag[i].Equ_X);
+ }*/
+ }
+ free(ModelBlock->Block_List[blk].IM_lead_lag);
+ for (i = 0; i < ModelBlock->Block_List[blk].Size; i++)
+ delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
+ free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
+ delete (ModelBlock->Block_List[blk].Temporary_InUse);
+ }
+ free(ModelBlock->Block_List);
+ free(ModelBlock);
+}
t_etype
BlockTriangular::Equation_Type_determination(vector &equations, map, NodeID> &first_cur_endo_derivatives, vector &Index_Var_IM, vector &Index_Equ_IM)
@@ -637,7 +683,7 @@ BlockTriangular::Equation_Type_determination(vector &equations,
temporary_terms_type temporary_terms;
EquationType Equation_Simulation_Type;
t_etype V_Equation_Simulation_Type(equations.size());
- for(unsigned int i = 0; i < equations.size(); i++)
+ for (unsigned int i = 0; i < equations.size(); i++)
{
temporary_terms.clear();
int eq = Index_Equ_IM[i];
@@ -646,59 +692,57 @@ BlockTriangular::Equation_Type_determination(vector &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, 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 > result;
- //result.clear();
+ pair res;
derivative->second->collectEndogenous(result);
+ /*for(set >::const_iterator iitt = result.begin(); iitt!=result.end(); iitt++)
+ cout << "result = " << iitt->first << ", " << iitt->second << "\n";*/
set >::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;
+ {
+ 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())
+ {
+ //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_R;
- //cout << "Equation " << eq << " is reversed\n";
- }
- else
- { //the equation could be normalized using the derivative independant of the endogenous variable
- if (d_endo_variable == result.end())
- {
- Equation_Simulation_Type = E_EVALUATE_S;
- Var_To_Derivate = var;
- }
+ Equation_Simulation_Type = E_EVALUATE_S;
+ //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(res.second));
}
- return(V_Equation_Simulation_Type);
+ return (V_Equation_Simulation_Type);
}
t_type
BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector > &blocks, vector &equations, t_etype &Equation_Type)
{
- int i=0;
- int count_equ = 0, blck_count_simult =0;
+ int i = 0;
+ int count_equ = 0, blck_count_simult = 0;
int Blck_Size, Recurs_Size;
int Lead, Lag;
t_type Type;
bool *Cur_IM;
- BlockSimulationType Simulation_Type , prev_Type=UNKNOWN;
+ BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
int eq = 0;
- for ( i=0; i 0)
{
if ((prev_Type == EVALUATE_FORWARD and Simulation_Type == EVALUATE_FORWARD)
- or (prev_Type == EVALUATE_BACKWARD and Simulation_Type == EVALUATE_BACKWARD))
+ or (prev_Type == EVALUATE_BACKWARD and Simulation_Type == EVALUATE_BACKWARD))
{
BlockSimulationType c_Type = (Type[Type.size()-1]).first;
int c_Size = (Type[Type.size()-1]).second.first;
- Type[Type.size()-1]=make_pair(c_Type, make_pair(++c_Size, Type[Type.size()-1].second.second));
+ Type[Type.size()-1] = make_pair(c_Type, make_pair(++c_Size, Type[Type.size()-1].second.second));
}
else
Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size)));
@@ -792,53 +836,53 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
prev_Type = Simulation_Type;
eq += Blck_Size;
}
- return(Type);
+ return (Type);
}
-
void
-BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0, jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map, NodeID> &first_cur_endo_derivatives)
+BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool *IM_0, jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map, NodeID> &first_cur_endo_derivatives)
{
int i, j, Nb_TotalBlocks, Nb_RecursBlocks, Nb_SimulBlocks;
BlockType Btype;
int count_Block, count_Equ;
- bool* SIM0, *SIM00;
+ bool *SIM0, *SIM00;
- SIM0 = (bool*)malloc(n * n * sizeof(bool));
- memcpy(SIM0,IM_0,n*n*sizeof(bool));
+ SIM0 = (bool *) malloc(n * n * sizeof(bool));
+ memcpy(SIM0, IM_0, n*n*sizeof(bool));
Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0);
free(SIM0);
- int counted=0;
- if (prologue+epilogue, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
+ double *max_val = (double *) malloc(n*sizeof(double));
+ memset(max_val, 0, n*sizeof(double));
+ for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
{
- if (fabs(iter->second)>max_val[iter->first.first])
- max_val[iter->first.first]=fabs(iter->second);
+ if (fabs(iter->second) > max_val[iter->first.first])
+ max_val[iter->first.first] = fabs(iter->second);
}
- for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
- iter->second/=max_val[iter->first.first];
+ for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
+ iter->second /= max_val[iter->first.first];
free(max_val);
- bool OK=false;
- double bi=0.99999999;
- int suppressed=0;
+ bool OK = false;
+ double bi = 0.99999999;
+ //double bi=1e-13;
+ int suppressed = 0;
vector Index_Equ_IM_save(Index_Equ_IM);
- while (!OK && bi>1e-14)
+ while (!OK && bi > 1e-14)
{
- int suppress=0;
+ int suppress = 0;
Index_Equ_IM = Index_Equ_IM_save;
- SIM0 = (bool*)malloc(n * n * sizeof(bool));
- memset(SIM0,0,n*n*sizeof(bool));
- SIM00 = (bool*)malloc(n * n * sizeof(bool));
- memset(SIM00,0,n*n*sizeof(bool));
- for ( map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++ )
+ SIM0 = (bool *) malloc(n * n * sizeof(bool));
+ memset(SIM0, 0, n*n*sizeof(bool));
+ SIM00 = (bool *) malloc(n * n * sizeof(bool));
+ memset(SIM00, 0, n*n*sizeof(bool));
+ for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
{
- if (fabs(iter->second)>bi)
+ if (fabs(iter->second) > bi)
{
- SIM0[iter->first.first*n+iter->first.second]=1;
+ SIM0[iter->first.first*n+iter->first.second] = 1;
if (!IM_0[iter->first.first*n+iter->first.second])
{
cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << " " << iter->second << "\n";
@@ -847,32 +891,31 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock,
else
suppress++;
}
- for (i = 0;i < n;i++)
- for (j = 0;j < n;j++)
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
{
SIM00[i*n + j] = SIM0[Index_Equ_IM[i] * n + Index_Var_IM[j]];
}
free(SIM0);
- if (suppress!=suppressed)
+ if (suppress != suppressed)
OK = Compute_Normalization(IM, n, prologue, epilogue, false, SIM00, Index_Equ_IM);
- suppressed=suppress;
+ suppressed = suppress;
if (!OK)
//bi/=1.07;
- bi/=3;
+ bi /= 2;
counted++;
- if (bi>1e-14)
+ 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);
cout << "Finding the optimal block decomposition of the model ...\n";
vector > blocks;
- if (prologue+epiloguefirst==SOLVE_FORWARD_COMPLETE || it->first==SOLVE_BACKWARD_COMPLETE || it->first==SOLVE_TWO_BOUNDARIES_COMPLETE)
+ if (it->first == SOLVE_FORWARD_COMPLETE || it->first == SOLVE_BACKWARD_COMPLETE || it->first == SOLVE_TWO_BOUNDARIES_COMPLETE)
{
Nb_SimulBlocks++;
- if (it->second.first>j)
+ if (it->second.first > j)
{
- j=it->second.first;
+ j = it->second.first;
Nb_fv = blocks[Nb_SimulBlocks-1].second;
}
}
@@ -898,20 +941,21 @@ 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);
+ 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++)
+
+ for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
{
- if (count_Equsecond.first==1)
+ else if (count_Equ < n-epilogue)
+ if (it->second.first == 1)
Btype = PROLOGUE;
else
Btype = SIMULTANS;
@@ -921,26 +965,25 @@ 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
void
BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector &equations, t_etype &equation_simulation_type, map, NodeID> &first_cur_endo_derivatives)
{
- bool* SIM, *SIM_0;
- bool* Cur_IM;
+ bool *SIM, *SIM_0;
+ bool *Cur_IM;
int i, k, size;
//First create a static model incidence matrix
size = symbol_table.endo_nbr() * symbol_table.endo_nbr() * sizeof(*SIM);
- SIM = (bool*)malloc(size);
- for (i = 0; i< symbol_table.endo_nbr() * symbol_table.endo_nbr(); i++) SIM[i] = 0;
- for (k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++)
+ SIM = (bool *) malloc(size);
+ for (i = 0; i < symbol_table.endo_nbr() * symbol_table.endo_nbr(); i++) SIM[i] = 0;
+ for (k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
{
Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
if (Cur_IM)
{
- for (i = 0;i < symbol_table.endo_nbr()*symbol_table.endo_nbr();i++)
+ for (i = 0; i < symbol_table.endo_nbr()*symbol_table.endo_nbr(); i++)
{
SIM[i] = (SIM[i]) || (Cur_IM[i]);
}
@@ -952,21 +995,21 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vec
incidencematrix.Print_SIM(SIM, eEndogenous);
}
Index_Equ_IM = vector(symbol_table.endo_nbr());
- for (i = 0;i < symbol_table.endo_nbr();i++)
+ for (i = 0; i < symbol_table.endo_nbr(); i++)
{
Index_Equ_IM[i] = i;
}
Index_Var_IM = vector(symbol_table.endo_nbr());
- for (i = 0;i < symbol_table.endo_nbr();i++)
+ for (i = 0; i < symbol_table.endo_nbr(); i++)
{
Index_Var_IM[i] = i;
}
if (ModelBlock != NULL)
Free_Block(ModelBlock);
- ModelBlock = (Model_Block*)malloc(sizeof(*ModelBlock));
+ ModelBlock = (Model_Block *) malloc(sizeof(*ModelBlock));
Cur_IM = incidencematrix.Get_IM(0, eEndogenous);
- SIM_0 = (bool*)malloc(symbol_table.endo_nbr() * symbol_table.endo_nbr() * sizeof(*SIM_0));
- for (i = 0;i < symbol_table.endo_nbr()*symbol_table.endo_nbr();i++)
+ SIM_0 = (bool *) malloc(symbol_table.endo_nbr() * symbol_table.endo_nbr() * sizeof(*SIM_0));
+ for (i = 0; i < symbol_table.endo_nbr()*symbol_table.endo_nbr(); i++)
SIM_0[i] = Cur_IM[i];
Normalize_and_BlockDecompose(SIM, ModelBlock, symbol_table.endo_nbr(), prologue, epilogue, Index_Var_IM, Index_Equ_IM, SIM_0, j_m, equations, equation_simulation_type, first_cur_endo_derivatives);
free(SIM_0);
diff --git a/BlockTriangular.hh b/BlockTriangular.hh
index 87761639..009741f5 100644
--- a/BlockTriangular.hh
+++ b/BlockTriangular.hh
@@ -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,double> jacob_map;
typedef vector > > t_type;
-typedef vector > t_etype;
+//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
+typedef vector > t_etype;
+
+//! Vector describing variables: max_lag in the block, max_lead in the block
+typedef vector > 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 &Index_Var_IM) const;
//! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables
void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector &Index_Equ_IM, vector &Index_Var_IM, vector > &blocks, 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 &equations, map, NodeID> &first_cur_endo_derivatives, vector &Index_Var_IM, vector &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 > &blocks, vector &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 &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 &equations, t_etype &V_Equation_Type, map, NodeID> &first_cur_endo_derivatives);
void Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM_0 , jacob_map &j_m, vector &equations, t_etype &equation_simulation_type, map, NodeID> &first_cur_endo_derivatives);
vector Index_Equ_IM;
vector 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 "
};
diff --git a/CodeInterpreter.hh b/CodeInterpreter.hh
index 1ecc837c..cb495131 100644
--- a/CodeInterpreter.hh
+++ b/CodeInterpreter.hh
@@ -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
diff --git a/DynamicModel.cc b/DynamicModel.cc
index 5e40bfc3..08ccf074 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -35,15 +35,15 @@
DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg) :
- ModelTree(symbol_table_arg, num_constants_arg),
- max_lag(0), max_lead(0),
- max_endo_lag(0), max_endo_lead(0),
- max_exo_lag(0), max_exo_lead(0),
- max_exo_det_lag(0), max_exo_det_lead(0),
- dynJacobianColsNbr(0),
- cutoff(1e-15),
- markowitz(0.7),
- block_triangular(symbol_table_arg)
+ ModelTree(symbol_table_arg, num_constants_arg),
+ max_lag(0), max_lead(0),
+ max_endo_lag(0), max_endo_lead(0),
+ max_exo_lag(0), max_exo_lead(0),
+ max_exo_det_lag(0), max_exo_det_lead(0),
+ dynJacobianColsNbr(0),
+ cutoff(1e-15),
+ markowitz(0.7),
+ block_triangular(symbol_table_arg, num_constants_arg)
{
}
@@ -55,13 +55,14 @@ 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)));
- if (it != first_derivatives.end())
- (it->second)->compile(code_file, false, temporary_terms, map_idx);
- else
- code_file.write(&FLDZ, sizeof(FLDZ));
-}
+ {
+ //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
+ code_file.write(&FLDZ, sizeof(FLDZ));
+ }
void
DynamicModel::BuildIncidenceMatrix()
@@ -77,7 +78,7 @@ DynamicModel::BuildIncidenceMatrix()
Id->collectEndogenous(endogenous);
for (set >::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 >::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++)
{
@@ -216,1838 +224,2109 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
void
DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &dynamic_basename) const
-{
- int i,j,k,m;
- string tmp_s, sps;
- ostringstream tmp_output, tmp1_output, global_output;
- NodeID lhs=NULL, rhs=NULL;
- BinaryOpNode *eq_node;
- ostringstream Uf[symbol_table.endo_nbr()];
- map reference_count;
- 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();
- int nze, nze_exo, nze_other_endo;
- //----------------------------------------------------------------------
- //For each block
- for (j = 0;j < ModelBlock->Size;j++)
- {
- //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++)
- 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);
- output << "%\n";
- output << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare\n";
- output << "%\n";
- output << "% Warning : this file is generated automatically by Dynare\n";
- output << "% from model file (.mod)\n\n";
- 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)
- {
- 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";
- 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";
- 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 << " % ////////////////////////////////////////////////////////////////////////" << endl
- << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type)
- << " //" << endl
- << " % // Simulation type "
- << BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //" << endl
- << " % ////////////////////////////////////////////////////////////////////////" << endl;
- //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)
- {
- 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 << " 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_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";
+ {
+ int i,j,k,m;
+ string tmp_s, sps;
+ ostringstream tmp_output, tmp1_output, global_output;
+ NodeID lhs=NULL, rhs=NULL;
+ BinaryOpNode *eq_node;
+ ostringstream Uf[symbol_table.endo_nbr()];
+ map reference_count;
+ 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();
+ int nze, nze_exo, nze_other_endo;
+ map recursive_variables;
+ vector 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+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);
+ output << "%\n";
+ output << "% " << tmp1_output.str() << " : Computes dynamic model for Dynare\n";
+ output << "%\n";
+ output << "% Warning : this file is generated automatically by Dynare\n";
+ output << "% from model file (.mod)\n\n";
+ 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*/)
+ {
+ 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, 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, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, it_, jacobian_eval)\n";
+ else
+ 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
+ << " % // Simulation type "
+ << BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //" << endl
+ << " % ////////////////////////////////////////////////////////////////////////" << endl;
+ //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*/)
+ {
+ output << " if(jacobian_eval)\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_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->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].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";
+ }
+
+ output << " g2=0;g3=0;\n";
+ if (ModelBlock->Block_List[j].Temporary_InUse->size())
+ {
+ tmp_output.str("");
+ for (temporary_terms_inuse_type::const_iterator it = ModelBlock->Block_List[j].Temporary_InUse->begin();
+ it != ModelBlock->Block_List[j].Temporary_InUse->end(); it++)
+ tmp_output << " T" << *it;
+ output << " global" << tmp_output.str() << ";\n";
+ }
+ 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 = 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";
+ output << " Per_K_=(it_-1)*y_size;\n";
+ sps=" ";
+ }
+ else
+ if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD )
+ sps = " ";
else
- output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size << ", " << ModelBlock->Block_List[j].Size << ", " << nze << ");\n";
- output << " end;\n";
- }
-
- output << " g2=0;g3=0;\n";
- if(ModelBlock->Block_List[j].Temporary_InUse->size())
- {
- tmp_output.str("");
- for (temporary_terms_inuse_type::const_iterator it = ModelBlock->Block_List[j].Temporary_InUse->begin();
- it != ModelBlock->Block_List[j].Temporary_InUse->end(); it++)
- 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)
- 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 << " 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";
- output << " Per_K_=(it_-1)*y_size;\n";
- 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)
- sps = " ";
+ sps="";
+ // The equations
+ for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ {
+ temporary_terms_type tt2;
+ tt2.clear();
+ if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size())
+ output << " " << sps << "% //Temporary variables" << endl;
+ 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++)
+ {
+ output << " " << sps;
+ (*it)->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ output << " = ";
+ (*it)->writeOutput(output, oMatlabDynamicModelSparse, tt2);
+ // Insert current node into tt2
+ tt2.insert(*it);
+ output << ";" << endl;
+ }
+ 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();
+ tmp_output.str("");
+ lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
+ switch (ModelBlock->Block_List[j].Simulation_Type)
+ {
+ case EVALUATE_BACKWARD:
+ case EVALUATE_FORWARD:
+evaluation:
+ 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 << " ";
+ if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
+ {
+ output << tmp_output.str();
+ output << " = ";
+ rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ }
+ /*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(EXIT_FAILURE);
+ }
+ output << ";\n";
+ break;
+ case SOLVE_BACKWARD_SIMPLE:
+ case SOLVE_FORWARD_SIMPLE:
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FORWARD_COMPLETE:
+ if (iBlock_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 << ") = (";
+ goto end;
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ if (iBlock_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-ModelBlock->Block_List[j].Nb_Recursives << ", it_)";
+ output << " residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ", it_) = (";
+ goto end;
+ default:
+end:
+ output << tmp_output.str();
+ output << ") - (";
+ rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ output << ");\n";
+#ifdef CONDITION
+ if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+ output << " condition(" << i+1 << ")=0;\n";
+#endif
+ }
+ }
+ // The Jacobian if we have to solve the block
+ if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE
+ || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
+ output << " " << sps << "% Jacobian " << endl;
else
- sps="";
- // The equations
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- temporary_terms_type tt2;
- tt2.clear();
- if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size())
- output << " " << sps << "% //Temporary variables" << endl;
- 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++)
- {
- output << " " << sps;
- (*it)->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << " = ";
- (*it)->writeOutput(output, oMatlabDynamicModelSparse, tt2);
- // Insert current node into tt2
- tt2.insert(*it);
- output << ";" << endl;
- }
- string sModel = symbol_table.getName(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();
- tmp_output.str("");
- lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
- switch (ModelBlock->Block_List[j].Simulation_Type)
- {
- case EVALUATE_BACKWARD:
- case EVALUATE_FORWARD:
- evaluation:
- 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 << " ";
- if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
- {
- output << tmp_output.str();
- output << " = ";
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- }
- 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
- {
- cerr << "Type missmatch for equation " << ModelBlock->Block_List[j].Equation[i]+1 << "\n";
- exit(-1);
- }
- 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;*/
- case SOLVE_BACKWARD_SIMPLE:
- case SOLVE_FORWARD_SIMPLE:
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FORWARD_COMPLETE:
- if(iBlock_List[j].Nb_Recursives)
- goto evaluation;
- 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 << ") = (";
- goto end;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- if(iBlock_List[j].Nb_Recursives)
- goto evaluation;
- 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_)";
- output << " residual(" << i+1-ModelBlock->Block_List[j].Nb_Recursives << ", it_) = (";
- goto end;
- default:
- end:
- output << tmp_output.str();
- output << ") - (";
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << ");\n";
-#ifdef CONDITION
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
- output << " condition(" << i+1 << ")=0;\n";
-#endif
- }
- }
- // The Jacobian if we have to solve the block
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE
- || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
- output << " " << sps << "% Jacobian " << endl;
- else
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE ||
- ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
- output << " % Jacobian " << endl << " if jacobian_eval" << endl;
- else
- output << " % Jacobian " << endl << " if jacobian_eval" << endl;
- switch (ModelBlock->Block_List[j].Simulation_Type)
- {
- case EVALUATE_BACKWARD:
- case EVALUATE_FORWARD:
- case EVALUATE_BACKWARD_R:
- case EVALUATE_FORWARD_R:
- count_derivates++;
- for (m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_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+(m+variable_table.max_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr*/
+ if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE ||
+ ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
+ output << " % Jacobian " << endl << " if jacobian_eval" << endl;
+ else
+ output << " % Jacobian " << endl << " if jacobian_eval" << endl;
+ switch (ModelBlock->Block_List[j].Simulation_Type)
+ {
+ case EVALUATE_BACKWARD:
+ case EVALUATE_FORWARD:
+ /*case EVALUATE_BACKWARD_R:
+ case EVALUATE_FORWARD_R:*/
+ count_derivates++;
+ for (m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
+ {
+ k=m-ModelBlock->Block_List[j].Max_Lag;
+ for (i=0;iBlock_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+(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)
- << "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
- << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }
- //jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr;
- /*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;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
- int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
- output << " g1_x(" << eqr+1 << ", "
- << varr+1+(m+max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr() << ") = ";
- writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
- output << "; % variable=" << symbol_table.getName(var)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }*/
- 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)
- {
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
- 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 << ", "
- << 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)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }
- }
- output << " varargout{1}=g1_x;\n";
- output << " varargout{2}=g1_o;\n";
- output << " end;" << endl;
- output << " end;" << endl;
- break;
- case SOLVE_BACKWARD_SIMPLE:
- case SOLVE_FORWARD_SIMPLE:
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FORWARD_COMPLETE:
- count_derivates++;
- for (m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_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+(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)
- << "(" << k
- << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }
- /*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;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
- int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
- output << " g1_x(" << eqr+1 << ", " << varr+1+(m+max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = ";
- writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
- output << "; % variable=" << symbol_table.getName(var)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }*/
- 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)
- {
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
- 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 << ", "
- << 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)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }
- }
- output << " varargout{1}=g1_x;\n";
- output << " varargout{2}=g1_o;\n";
- output << " else" << endl;
+ writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
+ 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;
+ }
+ }
+ //jacobian_max_endo_col=(variable_table.max_endo_lag+variable_table.max_endo_lead+1)*symbol_table.endo_nbr;
+ /*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;
+ for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
+ {
+ int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
+ int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
+ int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
+ output << " g1_x(" << eqr+1 << ", "
+ << varr+1+(m+max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr() << ") = ";
+ writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
+ output << "; % variable=" << symbol_table.getName(var)
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }*/
+ 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)
+ {
+ for (i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++)
+ {
+ int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
+ 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 << ", "
+ << 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(symbol_table.getID(eEndogenous, var))
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }
+ }
+ 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:
+ case SOLVE_FORWARD_SIMPLE:
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FORWARD_COMPLETE:
+ count_derivates++;
+ for (m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
+ {
+ k=m-ModelBlock->Block_List[j].Max_Lag;
+ for (i=0;iBlock_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 + m*(ModelBlock->Block_List[j].Size) << ") = ";
+ writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }
+ /*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;
+ for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
+ {
+ int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
+ int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
+ int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
+ output << " g1_x(" << eqr+1 << ", " << varr+1+(m+max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*ModelBlock->Block_List[j].nb_exo << ") = ";
+ writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
+ output << "; % variable=" << symbol_table.getName(var)
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }*/
+ 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)
+ {
+ for (i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++)
+ {
+ int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
+ 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-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(symbol_table.getID(eEndogenous, var))
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }
+ }
+ output << " varargout{1}=g1_x;\n";
+ output << " varargout{2}=g1_o;\n";
+ output << " else" << endl;
+
+ m=ModelBlock->Block_List[j].Max_Lag;
+ //cout << "\nDerivatives in Block " << j << "\n";
+ for (i=0;iBlock_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];
+ bool derivative_exist;
+ ostringstream tmp_output;
+ if (eqrBlock_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 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 (eqrBlock_List[j].Nb_Recursives or varrBlock_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;*/
+ }
+ output << " end;\n";
+ //output << " ya = y;\n";
+ break;
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ output << " if ~jacobian_eval" << endl;
+ 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;
+ for (i=0;iBlock_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];
+ bool derivative_exist;
+ ostringstream tmp_output;
+ //cout << "ModelBlock->Block_List[" << j << "].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n";
+ if (eqrBlock_List[j].Nb_Recursives)
+ {
+ /*if (varrBlock_List[j].Nb_Recursives)
+ {
+ if (k==0)
+ 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_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_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)
+ tmp_output << " g1_tmp_r(" << eqr+1 << ", "
+ << varr+1 << "+" << ModelBlock->Block_List[j].Size << "*" << ModelBlock->Block_List[j].Max_Lag << ") = ";
+ else if (k>0)
+ 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)
+ 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 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 << ";";
+
+ 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;
+ }*/
+ }
- m=ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_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 << ") = ";
- 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;
- }
- output << " end;\n";
- break;
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- output << " if ~jacobian_eval" << endl;
- 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;
- for (i=0;iBlock_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];
- 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 << ")";
- 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 << ")";
- 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 << ")";
- 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_) = ";
- else if (k>0)
- output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << ")) = ";
- else if (k<0)
- output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << ")) = ";
- writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
- output << "; % variable=" << symbol_table.getName(var)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
#ifdef CONDITION
- output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++)
- {
- output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
+ }
+ }
+ 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 << "))Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
- }
- }
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
+ 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;
+ for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+ output << " u(" << u+1 << "+Per_u_) = u(" << u+1 << "+Per_u_) / condition(" << eqr+1 << ");\n";
+ }
+ }
+ for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ output << " u(" << i+1 << "+Per_u_) = u(" << i+1 << "+Per_u_) / condition(" << i+1 << ");\n";
#endif
- output << " else" << endl;
- 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;
- for (i=0;iBlock_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+(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)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }
- jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size;
- /*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;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
- int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
- output << " g1_x(" << eqr+1 << ", "
- << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = ";
- writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
- output << "; % variable (exogenous)=" << symbol_table.getName(var)
- << "(" << k << ") " << var+1 << " " << varr+1
- << ", equation=" << eq+1 << endl;
- }
- }*/
- 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)
- {
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
- 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 << ", "
- << 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)
- << "(" << k << ") " << var+1
- << ", equation=" << eq+1 << endl;
- }
- }
- }
- output << " varargout{1}=g1_x;\n";
- output << " varargout{2}=g1_o;\n";
- output << " end;\n";
- output << " end;\n";
- break;
- default:
- break;
- }
- prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
- output.close();
- }
-}
+ output << " else" << endl;
+ 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;
+ for (i=0;iBlock_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+(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(symbol_table.getID(eEndogenous, var))
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }
+ jacobian_max_endo_col=(ModelBlock->Block_List[j].Max_Lead_Endo+ModelBlock->Block_List[j].Max_Lag_Endo+1)*ModelBlock->Block_List[j].Size;
+ /*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;
+ for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
+ {
+ int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X[i];
+ int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous[i];
+ int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
+ output << " g1_x(" << eqr+1 << ", "
+ << jacobian_max_endo_col+(m-(ModelBlock->Block_List[j].Max_Lag-ModelBlock->Block_List[j].Max_Lag_Exo))*ModelBlock->Block_List[j].nb_exo+varr+1 << ") = ";
+ writeDerivative(output, eq, symbol_table.getID(eExogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
+ output << "; % variable (exogenous)=" << symbol_table.getName(var)
+ << "(" << k << ") " << var+1 << " " << varr+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }*/
+ 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)
+ {
+ for (i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++)
+ {
+ int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
+ 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 << ", "
+ << 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(symbol_table.getID(eEndogenous, var))
+ << "(" << k << ") " << var+1
+ << ", equation=" << eq+1 << endl;
+ }
+ }
+ }
+ output << " varargout{1}=g1_x;\n";
+ output << " varargout{2}=g1_o;\n";
+ output << " end;\n";
+ //output << " ya = y;\n";
+ output << " end;\n";
+ break;
+ default:
+ break;
+ }
+ prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
+ output.close();
+ }
+ }
void
DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, map_idx_type map_idx) const
-{
- struct Uff_l
{
- int u, var, lag;
- Uff_l *pNext;
- };
+ struct Uff_l
+ {
+ int u, var, lag;
+ Uff_l *pNext;
+ };
- struct Uff
- {
- Uff_l *Ufl, *Ufl_First;
- int eqr;
- };
+ struct Uff
+ {
+ Uff_l *Ufl, *Ufl_First;
+ int eqr;
+ };
- int i,j,k,m, v, ModelBlock_Aggregated_Count, k0, k1;
- string tmp_s;
- ostringstream tmp_output;
- ofstream code_file;
- NodeID lhs=NULL, rhs=NULL;
- BinaryOpNode *eq_node;
- bool lhs_rhs_done;
- Uff Uf[symbol_table.endo_nbr()];
- map reference_count;
- map ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
- int prev_Simulation_Type=-1;
- //SymbolicGaussElimination SGE;
- bool file_open=false;
- temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
- //----------------------------------------------------------------------
- string main_name=file_name;
- main_name+=".cod";
- code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate );
- if (!code_file.is_open())
- {
- cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
- exit(EXIT_FAILURE);
- }
- //Temporary variables declaration
- code_file.write(&FDIMT, sizeof(FDIMT));
- k=temporary_terms.size();
- code_file.write(reinterpret_cast(&k),sizeof(k));
- //search for successive and identical blocks
- i=k=k0=0;
- ModelBlock_Aggregated_Count=-1;
- for (j = 0;j < ModelBlock->Size;j++)
- {
- 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 ))
- {
- }
- else
- {
- k=k0=0;
- ModelBlock_Aggregated_Count++;
- }
- k0+=ModelBlock->Block_List[j].Size;
- ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0;
- ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k;
- prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
- }
- ModelBlock_Aggregated_Count++;
- //For each block
- j=0;
- for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++)
- {
- k1=j;
- if (k0>0)
- code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
- code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK));
- v=ModelBlock_Aggregated_Number[k0];
- code_file.write(reinterpret_cast(&v),sizeof(v));
- v=ModelBlock->Block_List[j].Simulation_Type;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- for (k=0; kBlock_List[j].Size;i++)
- {
- code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
- code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
- code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i]));
- }
- j++;
- }
- j=k1;
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
- ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
- {
- code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].is_linear),sizeof(ModelBlock->Block_List[j].is_linear));
- v=block_triangular.ModelBlock->Block_List[j].IM_lead_lag[block_triangular.ModelBlock->Block_List[j].Max_Lag + block_triangular.ModelBlock->Block_List[j].Max_Lead].u_finish + 1;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- v=symbol_table.endo_nbr();
- code_file.write(reinterpret_cast(&v),sizeof(v));
- v=block_triangular.ModelBlock->Block_List[j].Max_Lag;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- v=block_triangular.ModelBlock->Block_List[j].Max_Lead;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
- //{
- int u_count_int=0;
- Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
- ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
- v=u_count_int;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- file_open=true;
- //}
- }
- for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++)
- {
- //For a block composed of a single equation determines whether we have to evaluate or to solve the equation
- if (ModelBlock->Block_List[j].Size==1)
- {
- lhs_rhs_done=true;
- eq_node = equations[ModelBlock->Block_List[j].Equation[0]];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- }
- else
- lhs_rhs_done=false;
- // The equations
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
- //The Temporary terms
- temporary_terms_type tt2;
+ int i,j,k,m, v, ModelBlock_Aggregated_Count, k0, k1;
+ string tmp_s;
+ ostringstream tmp_output;
+ ofstream code_file;
+ NodeID lhs=NULL, rhs=NULL;
+ BinaryOpNode *eq_node;
+ bool lhs_rhs_done;
+ Uff Uf[symbol_table.endo_nbr()];
+ map reference_count;
+ map ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
+ int prev_Simulation_Type=-1;
+ //SymbolicGaussElimination SGE;
+ bool file_open=false;
+ //temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
+ //----------------------------------------------------------------------
+ string main_name=file_name;
+ main_name+=".cod";
+ code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate );
+ if (!code_file.is_open())
+ {
+ cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
+ exit(EXIT_FAILURE);
+ }
+ //Temporary variables declaration
+ code_file.write(&FDIMT, sizeof(FDIMT));
+ k=temporary_terms.size();
+ code_file.write(reinterpret_cast(&k),sizeof(k));
+ //search for successive and identical blocks
+ i=k=k0=0;
+ ModelBlock_Aggregated_Count=-1;
+ for (j = 0;j < ModelBlock->Size;j++)
+ {
+ 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 */))
+ {
+ }
+ else
+ {
+ k=k0=0;
+ ModelBlock_Aggregated_Count++;
+ }
+ k0+=ModelBlock->Block_List[j].Size;
+ ModelBlock_Aggregated_Number[ModelBlock_Aggregated_Count]=k0;
+ ModelBlock_Aggregated_Size[ModelBlock_Aggregated_Count]=++k;
+ prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
+ }
+ ModelBlock_Aggregated_Count++;
+ //For each block
+ j=0;
+ for (k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++)
+ {
+ k1=j;
+ if (k0>0)
+ code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
+ code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK));
+ v=ModelBlock_Aggregated_Number[k0];
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ v=ModelBlock->Block_List[j].Simulation_Type;
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ for (k=0; kBlock_List[j].Size;i++)
+ {
+ code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
+ code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
+ code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].Own_Derivative[i]),sizeof(ModelBlock->Block_List[j].Own_Derivative[i]));
+ }
+ j++;
+ }
+ j=k1;
+ if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
+ ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
+ {
+ code_file.write(reinterpret_cast(&ModelBlock->Block_List[j].is_linear),sizeof(ModelBlock->Block_List[j].is_linear));
+ v=block_triangular.ModelBlock->Block_List[j].IM_lead_lag[block_triangular.ModelBlock->Block_List[j].Max_Lag + block_triangular.ModelBlock->Block_List[j].Max_Lead].u_finish + 1;
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ v=symbol_table.endo_nbr();
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ v=block_triangular.ModelBlock->Block_List[j].Max_Lag;
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ v=block_triangular.ModelBlock->Block_List[j].Max_Lead;
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ //if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+ //{
+ int u_count_int=0;
+ Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
+ ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
+ v=u_count_int;
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ file_open=true;
+ //}
+ }
+ for (k1 = 0; k1 < ModelBlock_Aggregated_Size[k0]; k1++)
+ {
+ //For a block composed of a single equation determines whether we have to evaluate or to solve the equation
+ if (ModelBlock->Block_List[j].Size==1)
+ {
+ lhs_rhs_done=true;
+ eq_node = equations[ModelBlock->Block_List[j].Equation[0]];
+ lhs = eq_node->get_arg1();
+ rhs = eq_node->get_arg2();
+ }
+ else
+ lhs_rhs_done=false;
+ // The equations
+ for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ {
+ //ModelBlock->Block_List[j].Variable_Sorted[i] = variable_table.getID(eEndogenous, ModelBlock->Block_List[j].Variable[i], 0);
+ //The Temporary terms
+ temporary_terms_type tt2;
#ifdef DEBUGC
- k=0;
+ k=0;
#endif
- 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)->compile(code_file, false, tt2, map_idx);
- code_file.write(&FSTPT, sizeof(FSTPT));
- map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
- v=(int)ii->second;
- code_file.write(reinterpret_cast(&v), sizeof(v));
- // Insert current node into tt2
- tt2.insert(*it);
+ 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)->compile(code_file, false, tt2, map_idx);
+ code_file.write(&FSTPT, sizeof(FSTPT));
+ map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
+ v=(int)ii->second;
+ code_file.write(reinterpret_cast(&v), sizeof(v));
+ // Insert current node into tt2
+ tt2.insert(*it);
#ifdef DEBUGC
- cout << "FSTPT " << v << "\n";
- code_file.write(&FOK, sizeof(FOK));
- code_file.write(reinterpret_cast(&k), sizeof(k));
- ki++;
+ cout << "FSTPT " << v << "\n";
+ code_file.write(&FOK, sizeof(FOK));
+ code_file.write(reinterpret_cast(&k), sizeof(k));
+ ki++;
#endif
- }
+ }
#ifdef DEBUGC
- for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin();
- it != ModelBlock->Block_List[j].Temporary_terms->end(); it++)
- {
- map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
- cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
- }
+ for (temporary_terms_type::const_iterator it = ModelBlock->Block_List[j].Temporary_terms->begin();
+ it != ModelBlock->Block_List[j].Temporary_terms->end(); it++)
+ {
+ map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
+ cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
+ }
#endif
- if (!lhs_rhs_done)
- {
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- }
- switch (ModelBlock->Block_List[j].Simulation_Type)
- {
- case EVALUATE_BACKWARD:
- case EVALUATE_FORWARD:
- 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_FORWARD_R:
- lhs->compile(code_file, false, temporary_terms, map_idx);
- rhs->compile(code_file, true, temporary_terms, map_idx);
- break;
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FORWARD_COMPLETE:
- v=ModelBlock->Block_List[j].Equation[i];
- Uf[v].eqr=i;
- Uf[v].Ufl=NULL;
- goto end;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- v=ModelBlock->Block_List[j].Equation[i];
- Uf[v].eqr=i;
- Uf[v].Ufl=NULL;
- goto end;
- default:
- end:
- lhs->compile(code_file, false, temporary_terms, map_idx);
- rhs->compile(code_file, false, temporary_terms, map_idx);
- code_file.write(&FBINARY, sizeof(FBINARY));
- int v=oMinus;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- code_file.write(&FSTPR, sizeof(FSTPR));
- code_file.write(reinterpret_cast(&i), sizeof(i));
+ if (!lhs_rhs_done)
+ {
+ eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
+ lhs = eq_node->get_arg1();
+ rhs = eq_node->get_arg2();
+ }
+ switch (ModelBlock->Block_List[j].Simulation_Type)
+ {
+ 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_FORWARD_R:
+ lhs->compile(code_file, false, temporary_terms, map_idx);
+ rhs->compile(code_file, true, temporary_terms, map_idx);
+ break;*/
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FORWARD_COMPLETE:
+ v=ModelBlock->Block_List[j].Equation[i];
+ Uf[v].eqr=i;
+ Uf[v].Ufl=NULL;
+ goto end;
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ v=ModelBlock->Block_List[j].Equation[i];
+ Uf[v].eqr=i;
+ Uf[v].Ufl=NULL;
+ goto end;
+ default:
+end:
+ lhs->compile(code_file, false, temporary_terms, map_idx);
+ rhs->compile(code_file, false, temporary_terms, map_idx);
+ code_file.write(&FBINARY, sizeof(FBINARY));
+ int v=oMinus;
+ code_file.write(reinterpret_cast(&v),sizeof(v));
+ code_file.write(&FSTPR, sizeof(FSTPR));
+ code_file.write(reinterpret_cast(&i), sizeof(i));
#ifdef CONDITION
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
- output << " condition[" << i << "]=0;\n";
+ if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
+ output << " condition[" << i << "]=0;\n";
#endif
- }
- }
- code_file.write(&FENDEQU, sizeof(FENDEQU));
- // 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)
- {
- switch (ModelBlock->Block_List[j].Simulation_Type)
- {
- case SOLVE_BACKWARD_SIMPLE:
- case SOLVE_FORWARD_SIMPLE:
- compileDerivative(code_file, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, map_idx);
- code_file.write(&FSTPG, sizeof(FSTPG));
- v=0;
- code_file.write(reinterpret_cast(&v), sizeof(v));
- break;
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FORWARD_COMPLETE:
- m=ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- int v=ModelBlock->Block_List[j].Equation[eqr];
- if (!Uf[v].Ufl)
- {
- Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
- Uf[v].Ufl_First=Uf[v].Ufl;
- }
- else
- {
- Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
- Uf[v].Ufl=Uf[v].Ufl->pNext;
- }
- Uf[v].Ufl->pNext=NULL;
- Uf[v].Ufl->u=u;
- Uf[v].Ufl->var=var;
- compileDerivative(code_file, eq, var, 0, map_idx);
- code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast(&u), sizeof(u));
- }
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- code_file.write(&FLDR, sizeof(FLDR));
- code_file.write(reinterpret_cast(&i), sizeof(i));
- code_file.write(&FLDZ, sizeof(FLDZ));
- int v=ModelBlock->Block_List[j].Equation[i];
- for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
- {
- code_file.write(&FLDU, sizeof(FLDU));
- code_file.write(reinterpret_cast(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
- code_file.write(&FLDV, sizeof(FLDV));
- char vc=eEndogenous;
- code_file.write(reinterpret_cast(&vc), sizeof(vc));
- code_file.write(reinterpret_cast(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var));
- int v1=0;
- code_file.write(reinterpret_cast(&v1), sizeof(v1));
- code_file.write(&FBINARY, sizeof(FBINARY));
- v1=oTimes;
- code_file.write(reinterpret_cast(&v1), sizeof(v1));
- code_file.write(&FCUML, sizeof(FCUML));
- }
- Uf[v].Ufl=Uf[v].Ufl_First;
- while (Uf[v].Ufl)
- {
- Uf[v].Ufl_First=Uf[v].Ufl->pNext;
- free(Uf[v].Ufl);
- Uf[v].Ufl=Uf[v].Ufl_First;
- }
- code_file.write(&FBINARY, sizeof(FBINARY));
- v=oMinus;
- code_file.write(reinterpret_cast(&v), sizeof(v));
- code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast(&i), sizeof(i));
- }
- break;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- 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;
- for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- int v=ModelBlock->Block_List[j].Equation[eqr];
- if (!Uf[v].Ufl)
- {
- Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
- Uf[v].Ufl_First=Uf[v].Ufl;
- }
- else
- {
- Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
- Uf[v].Ufl=Uf[v].Ufl->pNext;
- }
- Uf[v].Ufl->pNext=NULL;
- Uf[v].Ufl->u=u;
- Uf[v].Ufl->var=var;
- Uf[v].Ufl->lag=k;
- compileDerivative(code_file, eq, var, k, map_idx);
- code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast(&u), sizeof(u));
+ }
+ }
+ code_file.write(&FENDEQU, sizeof(FENDEQU));
+ // 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*/)
+ {
+ switch (ModelBlock->Block_List[j].Simulation_Type)
+ {
+ case SOLVE_BACKWARD_SIMPLE:
+ case SOLVE_FORWARD_SIMPLE:
+ compileDerivative(code_file, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, map_idx);
+ code_file.write(&FSTPG, sizeof(FSTPG));
+ v=0;
+ code_file.write(reinterpret_cast(&v), sizeof(v));
+ break;
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FORWARD_COMPLETE:
+ m=ModelBlock->Block_List[j].Max_Lag;
+ for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+ int v=ModelBlock->Block_List[j].Equation[eqr];
+ if (!Uf[v].Ufl)
+ {
+ Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
+ Uf[v].Ufl_First=Uf[v].Ufl;
+ }
+ else
+ {
+ Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
+ Uf[v].Ufl=Uf[v].Ufl->pNext;
+ }
+ Uf[v].Ufl->pNext=NULL;
+ Uf[v].Ufl->u=u;
+ Uf[v].Ufl->var=var;
+ compileDerivative(code_file, eq, var, 0, map_idx);
+ code_file.write(&FSTPU, sizeof(FSTPU));
+ code_file.write(reinterpret_cast(&u), sizeof(u));
+ }
+ for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ {
+ code_file.write(&FLDR, sizeof(FLDR));
+ code_file.write(reinterpret_cast(&i), sizeof(i));
+ code_file.write(&FLDZ, sizeof(FLDZ));
+ int v=ModelBlock->Block_List[j].Equation[i];
+ for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
+ {
+ code_file.write(&FLDU, sizeof(FLDU));
+ code_file.write(reinterpret_cast(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
+ code_file.write(&FLDV, sizeof(FLDV));
+ char vc=eEndogenous;
+ code_file.write(reinterpret_cast(&vc), sizeof(vc));
+ code_file.write(reinterpret_cast(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var));
+ int v1=0;
+ code_file.write(reinterpret_cast(&v1), sizeof(v1));
+ code_file.write(&FBINARY, sizeof(FBINARY));
+ v1=oTimes;
+ code_file.write(reinterpret_cast(&v1), sizeof(v1));
+ code_file.write(&FCUML, sizeof(FCUML));
+ }
+ Uf[v].Ufl=Uf[v].Ufl_First;
+ while (Uf[v].Ufl)
+ {
+ Uf[v].Ufl_First=Uf[v].Ufl->pNext;
+ free(Uf[v].Ufl);
+ Uf[v].Ufl=Uf[v].Ufl_First;
+ }
+ code_file.write(&FBINARY, sizeof(FBINARY));
+ v=oMinus;
+ code_file.write(reinterpret_cast(&v), sizeof(v));
+ code_file.write(&FSTPU, sizeof(FSTPU));
+ code_file.write(reinterpret_cast(&i), sizeof(i));
+ }
+ break;
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ 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;
+ for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+ int v=ModelBlock->Block_List[j].Equation[eqr];
+ if (!Uf[v].Ufl)
+ {
+ Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
+ Uf[v].Ufl_First=Uf[v].Ufl;
+ }
+ else
+ {
+ Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
+ Uf[v].Ufl=Uf[v].Ufl->pNext;
+ }
+ Uf[v].Ufl->pNext=NULL;
+ Uf[v].Ufl->u=u;
+ Uf[v].Ufl->var=var;
+ Uf[v].Ufl->lag=k;
+ compileDerivative(code_file, eq, var, k, map_idx);
+ code_file.write(&FSTPU, sizeof(FSTPU));
+ code_file.write(reinterpret_cast(&u), sizeof(u));
#ifdef CONDITION
- output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++)
- {
- code_file.write(&FLDR, sizeof(FLDR));
- code_file.write(reinterpret_cast(&i), sizeof(i));
- code_file.write(&FLDZ, sizeof(FLDZ));
- int v=ModelBlock->Block_List[j].Equation[i];
- for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
- {
- code_file.write(&FLDU, sizeof(FLDU));
- code_file.write(reinterpret_cast(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
- code_file.write(&FLDV, sizeof(FLDV));
- char vc=eEndogenous;
- code_file.write(reinterpret_cast(&vc), sizeof(vc));
- int v1=Uf[v].Ufl->var;
- code_file.write(reinterpret_cast(&v1), sizeof(v1));
- v1=Uf[v].Ufl->lag;
- code_file.write(reinterpret_cast(&v1), sizeof(v1));
- code_file.write(&FBINARY, sizeof(FBINARY));
- v1=oTimes;
- code_file.write(reinterpret_cast(&v1), sizeof(v1));
- code_file.write(&FCUML, sizeof(FCUML));
- }
- Uf[v].Ufl=Uf[v].Ufl_First;
- while (Uf[v].Ufl)
- {
- Uf[v].Ufl_First=Uf[v].Ufl->pNext;
- free(Uf[v].Ufl);
- Uf[v].Ufl=Uf[v].Ufl_First;
- }
- code_file.write(&FBINARY, sizeof(FBINARY));
- v=oMinus;
- code_file.write(reinterpret_cast(&v), sizeof(v));
- code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast(&i), sizeof(i));
+ }
+ }
+ for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ {
+ code_file.write(&FLDR, sizeof(FLDR));
+ code_file.write(reinterpret_cast(&i), sizeof(i));
+ code_file.write(&FLDZ, sizeof(FLDZ));
+ int v=ModelBlock->Block_List[j].Equation[i];
+ for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
+ {
+ code_file.write(&FLDU, sizeof(FLDU));
+ code_file.write(reinterpret_cast(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
+ code_file.write(&FLDV, sizeof(FLDV));
+ char vc=eEndogenous;
+ code_file.write(reinterpret_cast(&vc), sizeof(vc));
+ int v1=Uf[v].Ufl->var;
+ code_file.write(reinterpret_cast(&v1), sizeof(v1));
+ v1=Uf[v].Ufl->lag;
+ code_file.write(reinterpret_cast(&v1), sizeof(v1));
+ code_file.write(&FBINARY, sizeof(FBINARY));
+ v1=oTimes;
+ code_file.write(reinterpret_cast(&v1), sizeof(v1));
+ code_file.write(&FCUML, sizeof(FCUML));
+ }
+ Uf[v].Ufl=Uf[v].Ufl_First;
+ while (Uf[v].Ufl)
+ {
+ Uf[v].Ufl_First=Uf[v].Ufl->pNext;
+ free(Uf[v].Ufl);
+ Uf[v].Ufl=Uf[v].Ufl_First;
+ }
+ code_file.write(&FBINARY, sizeof(FBINARY));
+ v=oMinus;
+ code_file.write(reinterpret_cast(&v), sizeof(v));
+ code_file.write(&FSTPU, sizeof(FSTPU));
+ code_file.write(reinterpret_cast(&i), sizeof(i));
#ifdef CONDITION
- output << " if (fabs(condition[" << i << "])Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- output << " u[" << u << "+Per_u_] /= condition[" << eqr << "];\n";
- }
- }
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n";
+ 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;
+ for (i=0;iBlock_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 u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
+ int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
+ output << " u[" << u << "+Per_u_] /= condition[" << eqr << "];\n";
+ }
+ }
+ for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n";
#endif
- break;
- default:
- break;
- }
+ break;
+ default:
+ break;
+ }
- prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
- }
- j++;
- }
- }
- code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
- code_file.write(&FEND, sizeof(FEND));
- code_file.close();
-}
+ prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
+ }
+ j++;
+ }
+ }
+ code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
+ code_file.write(&FEND, sizeof(FEND));
+ code_file.close();
+ }
void
DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
-{
- string filename = dynamic_basename + ".m";
+ {
+ string filename = dynamic_basename + ".m";
- ofstream mDynamicModelFile;
- mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
- if (!mDynamicModelFile.is_open())
- {
- cerr << "Error: Can't open file " << filename << " for writing" << endl;
- exit(EXIT_FAILURE);
- }
- mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, it_)" << endl
- << "%" << endl
- << "% Status : Computes dynamic model for Dynare" << endl
- << "%" << endl
- << "% Warning : this file is generated automatically by Dynare" << endl
- << "% from model file (.mod)" << endl << endl;
+ ofstream mDynamicModelFile;
+ mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
+ if (!mDynamicModelFile.is_open())
+ {
+ cerr << "Error: Can't open file " << filename << " for writing" << endl;
+ exit(EXIT_FAILURE);
+ }
+ mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, it_)" << endl
+ << "%" << endl
+ << "% Status : Computes dynamic model for Dynare" << endl
+ << "%" << endl
+ << "% Warning : this file is generated automatically by Dynare" << endl
+ << "% from model file (.mod)" << endl << endl;
- writeDynamicModel(mDynamicModelFile);
+ writeDynamicModel(mDynamicModelFile);
- mDynamicModelFile.close();
-}
+ mDynamicModelFile.close();
+ }
void
DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
-{
- string filename = dynamic_basename + ".c";
- ofstream mDynamicModelFile;
+ {
+ string filename = dynamic_basename + ".c";
+ ofstream mDynamicModelFile;
- mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
- if (!mDynamicModelFile.is_open())
- {
- cerr << "Error: Can't open file " << filename << " for writing" << endl;
- exit(EXIT_FAILURE);
- }
- mDynamicModelFile << "/*" << endl
- << " * " << filename << " : Computes dynamic model for Dynare" << endl
- << " *" << endl
- << " * Warning : this file is generated automatically by Dynare" << endl
- << " * from model file (.mod)" << endl
- << endl
- << " */" << endl
- << "#include " << endl
- << "#include \"mex.h\"" << endl;
+ mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
+ if (!mDynamicModelFile.is_open())
+ {
+ cerr << "Error: Can't open file " << filename << " for writing" << endl;
+ exit(EXIT_FAILURE);
+ }
+ mDynamicModelFile << "/*" << endl
+ << " * " << filename << " : Computes dynamic model for Dynare" << endl
+ << " *" << endl
+ << " * Warning : this file is generated automatically by Dynare" << endl
+ << " * from model file (.mod)" << endl
+ << endl
+ << " */" << endl
+ << "#include " << endl
+ << "#include \"mex.h\"" << endl;
- // Writing the function body
- writeDynamicModel(mDynamicModelFile);
+ // Writing the function body
+ writeDynamicModel(mDynamicModelFile);
- // Writing the gateway routine
- mDynamicModelFile << "/* The gateway routine */" << endl
- << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
- << "{" << endl
- << " double *y, *x, *params;" << endl
- << " double *residual, *g1, *g2;" << endl
- << " int nb_row_x, it_;" << endl
- << endl
- << " /* Create a pointer to the input matrix y. */" << endl
- << " y = mxGetPr(prhs[0]);" << endl
- << endl
- << " /* Create a pointer to the input matrix x. */" << endl
- << " x = mxGetPr(prhs[1]);" << endl
- << endl
- << " /* Create a pointer to the input matrix params. */" << endl
- << " params = mxGetPr(prhs[2]);" << endl
- << endl
- << " /* Fetch time index */" << endl
- << " it_ = (int) mxGetScalar(prhs[3]) - 1;" << endl
- << endl
- << " /* Gets number of rows of matrix x. */" << endl
- << " nb_row_x = mxGetM(prhs[1]);" << endl
- << endl
- << " residual = NULL;" << endl
- << " if (nlhs >= 1)" << endl
- << " {" << endl
- << " /* Set the output pointer to the output matrix residual. */" << endl
- << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
- << " /* Create a C pointer to a copy of the output matrix residual. */" << endl
- << " residual = mxGetPr(plhs[0]);" << endl
- << " }" << endl
- << endl
- << " g1 = NULL;" << endl
- << " if (nlhs >= 2)" << endl
- << " {" << endl
- << " /* Set the output pointer to the output matrix g1. */" << endl
+ // Writing the gateway routine
+ mDynamicModelFile << "/* The gateway routine */" << endl
+ << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
+ << "{" << endl
+ << " double *y, *x, *params;" << endl
+ << " double *residual, *g1, *g2;" << endl
+ << " int nb_row_x, it_;" << endl
+ << endl
+ << " /* Create a pointer to the input matrix y. */" << endl
+ << " y = mxGetPr(prhs[0]);" << endl
+ << endl
+ << " /* Create a pointer to the input matrix x. */" << endl
+ << " x = mxGetPr(prhs[1]);" << endl
+ << endl
+ << " /* Create a pointer to the input matrix params. */" << endl
+ << " params = mxGetPr(prhs[2]);" << endl
+ << endl
+ << " /* Fetch time index */" << endl
+ << " it_ = (int) mxGetScalar(prhs[3]) - 1;" << endl
+ << endl
+ << " /* Gets number of rows of matrix x. */" << endl
+ << " nb_row_x = mxGetM(prhs[1]);" << endl
+ << endl
+ << " residual = NULL;" << endl
+ << " if (nlhs >= 1)" << endl
+ << " {" << endl
+ << " /* Set the output pointer to the output matrix residual. */" << endl
+ << " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
+ << " /* Create a C pointer to a copy of the output matrix residual. */" << endl
+ << " residual = mxGetPr(plhs[0]);" << endl
+ << " }" << endl
+ << endl
+ << " g1 = NULL;" << endl
+ << " if (nlhs >= 2)" << endl
+ << " {" << endl
+ << " /* Set the output pointer to the output matrix g1. */" << endl
- << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl
- << " /* Create a C pointer to a copy of the output matrix g1. */" << endl
- << " g1 = mxGetPr(plhs[1]);" << endl
- << " }" << endl
- << endl
- << " g2 = NULL;" << endl
- << " if (nlhs >= 3)" << endl
- << " {" << endl
- << " /* Set the output pointer to the output matrix g2. */" << endl
- << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr*dynJacobianColsNbr
- << ", mxREAL);" << endl
- << " /* Create a C pointer to a copy of the output matrix g1. */" << endl
- << " g2 = mxGetPr(plhs[2]);" << endl
- << " }" << endl
- << endl
- << " /* Call the C subroutines. */" << endl
- << " Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);" << endl
- << "}" << endl;
- mDynamicModelFile.close();
-}
+ << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl
+ << " /* Create a C pointer to a copy of the output matrix g1. */" << endl
+ << " g1 = mxGetPr(plhs[1]);" << endl
+ << " }" << endl
+ << endl
+ << " g2 = NULL;" << endl
+ << " if (nlhs >= 3)" << endl
+ << " {" << endl
+ << " /* Set the output pointer to the output matrix g2. */" << endl
+ << " plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr*dynJacobianColsNbr
+ << ", mxREAL);" << endl
+ << " /* Create a C pointer to a copy of the output matrix g1. */" << endl
+ << " g2 = mxGetPr(plhs[2]);" << endl
+ << " }" << endl
+ << endl
+ << " /* Call the C subroutines. */" << endl
+ << " Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);" << endl
+ << "}" << endl;
+ mDynamicModelFile.close();
+ }
string
DynamicModel::reform(const string name1) const
-{
- string name=name1;
- int pos = name.find("\\", 0);
- while (pos >= 0)
- {
- if (name.substr(pos + 1, 1) != "\\")
- {
- name = name.insert(pos, "\\");
- pos++;
- }
- pos++;
- pos = name.find("\\", pos);
- }
- return (name);
-}
+ {
+ string name=name1;
+ int pos = name.find("\\", 0);
+ while (pos >= 0)
+ {
+ if (name.substr(pos + 1, 1) != "\\")
+ {
+ name = name.insert(pos, "\\");
+ pos++;
+ }
+ pos++;
+ pos = name.find("\\", pos);
+ }
+ return (name);
+ }
void
DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num,
int &u_count_int, bool &file_open, bool is_two_boundaries) const
-{
- int j;
- std::ofstream SaveCode;
- if (file_open)
- SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate );
- else
- SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::binary);
- if (!SaveCode.is_open())
- {
- cout << "Error : Can't open file \"" << bin_basename << ".bin\" for writing\n";
- exit(EXIT_FAILURE);
- }
- u_count_int=0;
- for (int m=0;m<=block_triangular.ModelBlock->Block_List[num].Max_Lead+block_triangular.ModelBlock->Block_List[num].Max_Lag;m++)
- {
- int k1=m-block_triangular.ModelBlock->Block_List[num].Max_Lag;
- for (j=0;jBlock_List[num].IM_lead_lag[m].size;j++)
- {
- int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size;
- int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j];
- int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j];
- SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1));
- SaveCode.write(reinterpret_cast(&varr), sizeof(varr));
- SaveCode.write(reinterpret_cast(&k1), sizeof(k1));
- SaveCode.write(reinterpret_cast(&u), sizeof(u));
- u_count_int++;
- }
- }
- if(is_two_boundaries)
- {
- for (j=0;jBlock_List[num].Size;j++)
- {
- int eqr1=j;
- int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods
- +block_triangular.incidencematrix.Model_Max_Lead_Endo);
- int k1=0;
- SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1));
- SaveCode.write(reinterpret_cast(&varr), sizeof(varr));
- SaveCode.write(reinterpret_cast(&k1), sizeof(k1));
- SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1));
- u_count_int++;
- }
- }
- //cout << "u_count_int=" << u_count_int << "\n";
- for (j=0;jBlock_List[num].Size;j++)
- {
- int varr=block_triangular.ModelBlock->Block_List[num].Variable[j];
- SaveCode.write(reinterpret_cast(&varr), sizeof(varr));
- }
- for (j=0;jBlock_List[num].Size;j++)
- {
- int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j];
- SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1));
- }
- SaveCode.close();
-}
+ {
+ int j;
+ std::ofstream SaveCode;
+ if (file_open)
+ SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate );
+ else
+ SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::binary);
+ if (!SaveCode.is_open())
+ {
+ cout << "Error : Can't open file \"" << bin_basename << ".bin\" for writing\n";
+ exit(EXIT_FAILURE);
+ }
+ u_count_int=0;
+ for (int m=0;m<=block_triangular.ModelBlock->Block_List[num].Max_Lead+block_triangular.ModelBlock->Block_List[num].Max_Lag;m++)
+ {
+ int k1=m-block_triangular.ModelBlock->Block_List[num].Max_Lag;
+ for (j=0;jBlock_List[num].IM_lead_lag[m].size;j++)
+ {
+ int varr=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Var[j]+k1*block_triangular.ModelBlock->Block_List[num].Size;
+ int u=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].u[j];
+ int eqr1=block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].Equ[j];
+ SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1));
+ SaveCode.write(reinterpret_cast(&varr), sizeof(varr));
+ SaveCode.write(reinterpret_cast(&k1), sizeof(k1));
+ SaveCode.write(reinterpret_cast