diff --git a/BlockTriangular.cc b/BlockTriangular.cc
deleted file mode 100644
index ad0dd471..00000000
--- a/BlockTriangular.cc
+++ /dev/null
@@ -1,1166 +0,0 @@
-/*
- * Copyright (C) 2007-2009 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-#include "MinimumFeedbackSet.hh"
-#include
-#include
-#include
-#include
-
-//------------------------------------------------------------------------------
-#include "BlockTriangular.hh"
-//------------------------------------------------------------------------------
-using namespace std;
-using namespace boost;
-using namespace MFS;
-
-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;
- prologue = epilogue = 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)
-{
- bool modifie = 1;
- int i, j, k, l = 0;
- /*Looking for a prologue */
- prologue = 0;
- while (modifie)
- {
- modifie = 0;
- for (i = prologue; i < n; i++)
- {
- k = 0;
- for (j = prologue; j < n; j++)
- {
- if (IM[i*n + j])
- {
- k++;
- l = j;
- }
- }
- if ((k == 1) && IM0[Index_Equ_IM[i]*n + Index_Var_IM[l]])
- {
- modifie = 1;
- incidencematrix.swap_IM_c(IM, prologue, i, l, Index_Var_IM, Index_Equ_IM, n);
- prologue++;
- }
- }
- }
- epilogue = 0;
- modifie = 1;
- while (modifie)
- {
- modifie = 0;
- for (i = prologue; i < n - epilogue; i++)
- {
- k = 0;
- for (j = prologue; j < n - epilogue; j++)
- {
- if (IM[j*n + i])
- {
- k++;
- l = j;
- }
- }
- if ((k == 1) && IM0[Index_Equ_IM[l]*n + Index_Var_IM[i]])
- {
- modifie = 1;
- incidencematrix.swap_IM_c(IM, n - (1 + epilogue), l, i, Index_Var_IM, Index_Equ_IM, n);
- epilogue++;
- }
- }
- }
-}
-
-//------------------------------------------------------------------------------
-// Find a matching between equations and endogenous variables
-bool
-BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, int verbose, bool *IM0, vector &Index_Equ_IM) const
-{
- int n = equation_number - prologue - epilogue;
-
- typedef adjacency_list BipartiteGraph;
-
-
- if(verbose == 2)
- cout << "trying to normlized even in singular case\n";
- /*
- 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])
- {
- //printf("equation=%3d variable=%3d\n",i,j);
- 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 == 1)
- {
- cout << "WARNING: Could not normalize dynamic model. Variable "
- << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
- << " is not in the maximum cardinality matching. Trying to find a singular normalization." << endl;
- //exit(EXIT_FAILURE);
- return false;
- }
- else if (verbose == 2)
- {
- cerr << "ERROR: Could not normalize dynamic model (even with a singularity). 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++)
- {
- //printf("match equation %4d with variable %4d \n", mate_map[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;
-}
-
-
-
-
-t_vtype
-BlockTriangular::Get_Variable_LeadLag_By_Block(vector &components_set, int nb_blck_sim, int prologue, int epilogue, t_vtype &equation_lead_lag) 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 < (int)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);
- }
- Variable_Type[i] = make_pair(0, 0);
- }
- equation_lead_lag = Variable_Type;
- 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++)
- {
- int j_l = Index_Equ_IM[ j];
- if (Cur_IM[i_1 + Index_Equ_IM[ j] * nb_endo] and variable_blck[i_1] == equation_blck[j_l])
- {
- 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);
- if (k > equation_lead_lag[j_l].second)
- equation_lead_lag[j_l] = make_pair(equation_lead_lag[j_l].first, k);
- if (k < -equation_lead_lag[j_l].first)
- equation_lead_lag[j_l] = make_pair(-k, equation_lead_lag[j_l].second);
- }
- }
- }
- }
- }
- 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_, bool select_feedback_variable, int mfs) const
-{
- t_vtype V_Variable_Type;
- int n = nb_var - prologue - epilogue;
- bool *AMp;
- AMp = (bool *) malloc(n*n*sizeof(bool));
- //transforms the incidence matrix of the complet model into an adjancent matrix of the non-recursive part of the model
- for (int i = prologue; i < nb_var - epilogue; i++)
- for (int j = prologue; j < nb_var - epilogue; j++)
- if (j != i)
- AMp[(i-prologue)*n+j-prologue] = IM[i*nb_var + j];
- else
- AMp[(i-prologue)*n+j-prologue] = 0;
-
- //In a first step we compute the strong components of the graph representation of the static model.
- // This insures that block are dynamically recursives.
- GraphvizDigraph G2 = AM_2_GraphvizDigraph(AMp, n);
- vector endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
-
- int num = strong_components(G2, &endo2block[0]);
-
- blocks = vector >(num, make_pair(0, 0));
-
-
-
-/*New*/
- // Compute strongly connected components
- // Create directed acyclic graph associated to the strongly connected components
- typedef adjacency_list DirectedGraph;
- DirectedGraph dag(num);
- /*graph_traits::edge_iterator ei, ei_end;
- for(tie(ei, ei_end) = edges(G2); ei != ei_end; ++ei)
- {
- int s = endo2block[source(*ei, G2)];
- int t = endo2block[target(*ei, G2)];
- if (s != t)
- add_edge(s, t, dag);
- }*/
- for (unsigned int i = 0;i < num_vertices(G2);i++)
- {
- GraphvizDigraph::out_edge_iterator it_out, out_end;
- GraphvizDigraph::vertex_descriptor vi = vertex(i, G2);
- for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
- {
- int t_b = endo2block[target(*it_out, G2)];
- int s_b = endo2block[source(*it_out, G2)];
- if (s_b != t_b)
- add_edge(s_b, t_b, dag);
- }
- }
- // Compute topological sort of DAG (ordered list of unordered SCC)
- deque ordered2unordered;
- topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
- // Construct mapping from unordered SCC to ordered SCC
- vector unordered2ordered(num);
- for(int i = 0; i < num; i++)
- unordered2ordered[ordered2unordered[i]] = i;
-/*EndNew*/
-
-
- //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 < endo2block.size(); i++)
- {
- endo2block[i] = unordered2ordered[endo2block[i]];
- blocks[endo2block[i]].first++;
- components_set[endo2block[i]].first.insert(i);
- }
-
-
- t_vtype equation_lead_lag;
- V_Variable_Type = Get_Variable_LeadLag_By_Block(endo2block, num, prologue, epilogue, equation_lead_lag);
-
- 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 or vertices related to lead variables => force those vertices to belong to the feedback set
- if(select_feedback_variable)
- 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
- or equation_lead_lag[Index_Equ_IM[i+prologue]].second > 0 or equation_lead_lag[Index_Equ_IM[i+prologue]].first > 0
- or mfs == 0)
- add_edge(i, i, G2);
- else
- for (int i = 0; i < n; i++)
- if (Equation_Type[Index_Equ_IM[i+prologue]].first == E_SOLVE or mfs == 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 < 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;
- Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
- //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, vector &Index_Var_IM, vector &Index_Equ_IM)
-{
- 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;
- bool *tmp_variable_evaluated;
- bool *Cur_IM;
- bool *IM, OK;
- int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo;
- //cout << "block " << count_Block << " size " << size << " SimType=" << BlockSim(SimType) << "\n";
-
- ModelBlock->Periods = periods;
- 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].Chain_Rule_Derivatives = new chain_rule_derivatives_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_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((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * 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));
- memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
- memset(tmp_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
- nb_lead_lag_endo = 0;
- 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));
- memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr()*sizeof(bool));
- 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]->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_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++)
- {
- Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
- if (Cur_IM)
- {
- OK = false;
- if (k >= 0)
- {
- for (j = 0; j < size; j++)
- {
- if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j]*symbol_table.endo_nbr()])
- {
- tmp_variable_evaluated[i_1] = true;
- tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
- if (!OK)
- {
- tmp_endo[incidencematrix.Model_Max_Lag + k]++;
- nb_lead_lag_endo++;
- OK = true;
- }
- if (k > Lead)
- Lead = k;
- }
- }
- }
- else
- {
- for (j = 0; j < size; j++)
- {
- if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j]*symbol_table.endo_nbr()])
- {
- tmp_variable_evaluated[i_1] = true;
- tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++;
- if (!OK)
- {
- tmp_variable_evaluated[i_1] = true;
- tmp_endo[incidencematrix.Model_Max_Lag + k]++;
- nb_lead_lag_endo++;
- OK = true;
- }
- if (-k > Lag)
- Lag = -k;
- }
- }
- }
- }
- }
- (*count_Equ)++;
- }
- Lag_Endo = Lag;
- Lead_Endo = Lead;
-
- tmp_nb_other_endo = 0;
- for (i = 0; i < size; i++)
- {
- 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++)
- if (Cur_IM[i_1 + j])
- {
- if (!tmp_variable_evaluated[j])
- {
- tmp_other_endo[incidencematrix.Model_Max_Lag + k]++;
- tmp_nb_other_endo++;
- }
- if (k > 0 && k > Lead_Other_Endo)
- Lead_Other_Endo = k;
- else if (k < 0 && (-k) > Lag_Other_Endo)
- Lag_Other_Endo = -k;
- if (k > 0 && k > Lead)
- Lead = k;
- else if (k < 0 && (-k) > Lag)
- Lag = -k;
- tmp_size_other_endo[k+incidencematrix.Model_Max_Lag_Endo]++;
- }
- }
- }
- }
- ModelBlock->Block_List[count_Block].nb_other_endo = tmp_nb_other_endo;
- ModelBlock->Block_List[count_Block].Other_Endogenous = (int *) malloc(tmp_nb_other_endo * sizeof(int));
-
- tmp_exo = (int *) malloc(symbol_table.exo_nbr() * sizeof(int));
- memset(tmp_exo, 0, symbol_table.exo_nbr() * sizeof(int));
- tmp_nb_exo = 0;
- for (i = 0; i < size; i++)
- {
- 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; j < symbol_table.exo_nbr(); j++)
- if (Cur_IM[i_1 + j])
- {
- if (!tmp_exo[j])
- {
- tmp_exo[j] = 1;
- tmp_nb_exo++;
- }
- if (k > 0 && k > Lead_Exo)
- Lead_Exo = k;
- else if (k < 0 && (-k) > Lag_Exo)
- Lag_Exo = -k;
- if (k > 0 && k > Lead)
- Lead = k;
- else if (k < 0 && (-k) > Lag)
- Lag = -k;
- tmp_size_exo[k+incidencematrix.Model_Max_Lag_Exo]++;
- }
- }
- }
- }
-
- ModelBlock->Block_List[count_Block].nb_exo = tmp_nb_exo;
- ModelBlock->Block_List[count_Block].Exogenous = (int *) malloc(tmp_nb_exo * sizeof(int));
- k = 0;
- for (j = 0; j < symbol_table.exo_nbr(); j++)
- if (tmp_exo[j])
- {
- ModelBlock->Block_List[count_Block].Exogenous[k] = j;
- k++;
- }
-
- ModelBlock->Block_List[count_Block].nb_exo_det = 0;
-
- ModelBlock->Block_List[count_Block].Max_Lag = Lag;
- ModelBlock->Block_List[count_Block].Max_Lead = Lead;
- ModelBlock->Block_List[count_Block].Max_Lag_Endo = Lag_Endo;
- ModelBlock->Block_List[count_Block].Max_Lead_Endo = Lead_Endo;
- ModelBlock->Block_List[count_Block].Max_Lag_Other_Endo = Lag_Other_Endo;
- 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));
- 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++)
- {
- 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].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));
- }
- 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;*/
- 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++)
- {
- i_1 = Index_Var_IM[j];
- m = 0;
- for (k = first_count_equ; k < size + first_count_equ; k++)
- if (IM[i_1 + Index_Equ_IM[k] * symbol_table.endo_nbr()])
- m++;
- if (m > 0)
- {
- tmp_var[j - first_count_equ] = i1;
- i1++;
- }
- }
- m = 0;
- 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++)
- if (IM[Index_Var_IM[k] + i_1])
- {
- if (i == Lag)
- {
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].us[m] = ls;
- ls++;
- }
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].u[m] = li;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j];
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k];
- tmp_variable_evaluated[Index_Var_IM[k]] = true;
- l++;
- m++;
- li++;
- }
- }
- 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++)
- {
- i_1 = Index_Equ_IM[j] * symbol_table.endo_nbr();
- for (k = 0; k < symbol_table.endo_nbr(); k++)
- if ((!tmp_variable_evaluated[Index_Var_IM[k]]) && IM[Index_Var_IM[k] + i_1])
- {
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].u_other_endo[m] = l;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_other_endo[m] = j - first_count_equ;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_other_endo[m] = k;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_Index_other_endo[m] = Index_Equ_IM[j];
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Var_Index_other_endo[m] = Index_Var_IM[k];
- l++;
- m++;
- }
- }
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].size_other_endo = m;
- }
- /*IM = incidencematrix.Get_IM(i - Lag, eExogenous);
- if (IM)
- {
- m = 0;
- for (j = first_count_equ;j < size + first_count_equ;j++)
- {
- i_1 = Index_Equ_IM[j] * symbol_table.exo_nbr();
- for (k = 0; kBlock_List[count_Block].Exogenous[k]+i_1])
- {
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Exogenous[m] = k;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Exogenous_Index[m] = ModelBlock->Block_List[count_Block].Exogenous[k];
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X[m] = j - first_count_equ;
- ModelBlock->Block_List[count_Block].IM_lead_lag[i].Equ_X_Index[m] = Index_Equ_IM[j];
- m++;
- }
- }
- }
- }*/
- }
- free(tmp_size);
- free(tmp_size_other_endo);
- free(tmp_size_exo);
- free(tmp_endo);
- free(tmp_other_endo);
- free(tmp_exo);
- free(tmp_var);
- 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_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);
- delete ModelBlock->Block_List[blk].Chain_Rule_Derivatives;
- }
- free(ModelBlock->Block_List);
- free(ModelBlock);
-}
-
-t_etype
-BlockTriangular::Equation_Type_determination(vector &equations, map >, NodeID> &first_order_endo_derivatives, vector &Index_Var_IM, vector &Index_Equ_IM, int mfs)
-{
- NodeID lhs, rhs;
- ostringstream tmp_output;
- BinaryOpNode *eq_node;
- ostringstream tmp_s;
- 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++)
- {
- temporary_terms.clear();
- int eq = Index_Equ_IM[i];
- int var = Index_Var_IM[i];
- eq_node = equations[eq];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- Equation_Simulation_Type = E_SOLVE;
- 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_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
- pair res;
- if(derivative != first_order_endo_derivatives.end())
- {
- set > result;
- derivative->second->collectEndogenous(result);
- set >::const_iterator d_endo_variable = result.find(make_pair(var, 0));
- //Determine whether the equation could be evaluated rather than to be solved
- ostringstream tt("");
- derivative->second->writeOutput(tt, oMatlabDynamicModelSparse, temporary_terms);
- if (tmp_output.str() == tmp_s.str() and tt.str()=="1")
- {
- Equation_Simulation_Type = E_EVALUATE;
- }
- else
- {
- vector > > List_of_Op_RHS;
- res = equations[eq]->normalizeEquation(var, List_of_Op_RHS);
- if(mfs==2)
- {
- if(d_endo_variable == result.end() && res.second)
- Equation_Simulation_Type = E_EVALUATE_S;
- }
- else if(mfs==3)
- {
- if(res.second) // The equation could be solved analytically
- Equation_Simulation_Type = E_EVALUATE_S;
- }
- }
- }
- V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast(res.second));
- }
- 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, vector &Index_Var_IM, vector &Index_Equ_IM)
-{
- 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;
- int eq = 0;
- for (i = 0; i < prologue+(int) blocks.size()+epilogue; i++)
- {
- int first_count_equ = count_equ;
- if (i < prologue)
- {
- Blck_Size = 1;
- Recurs_Size = 0;
- }
- else if (i < prologue+(int) blocks.size())
- {
- Blck_Size = blocks[blck_count_simult].first;
- Recurs_Size = Blck_Size - blocks[blck_count_simult].second;
- blck_count_simult++;
- }
- else if (i < prologue+(int) blocks.size()+epilogue)
- {
- Blck_Size = 1;
- Recurs_Size = 0;
- }
-
- Lag = Lead = 0;
- for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
- {
- int i_1 = Index_Var_IM[count_equ];
- for (int k = -incidencematrix.Model_Max_Lag_Endo; k <= incidencematrix.Model_Max_Lead_Endo; k++)
- {
- Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
- if (Cur_IM)
- {
- for (int j = 0; j < Blck_Size; j++)
- {
- if (Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j] * symbol_table.endo_nbr()])
- {
- if (k > Lead)
- Lead = k;
- else if (-k > Lag)
- Lag = -k;
- }
- }
- }
- }
- }
- if ((Lag > 0) && (Lead > 0))
- {
- if (Blck_Size == 1)
- Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
- else
- Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
- }
- else if (Blck_Size > 1)
- {
- if (Lead > 0)
- Simulation_Type = SOLVE_BACKWARD_COMPLETE;
- else
- Simulation_Type = SOLVE_FORWARD_COMPLETE;
- }
- else
- {
- if (Lead > 0)
- Simulation_Type = SOLVE_BACKWARD_SIMPLE;
- else
- Simulation_Type = SOLVE_FORWARD_SIMPLE;
- }
- if (Blck_Size == 1)
- {
- if (Equation_Type[Index_Equ_IM[eq]].first == E_EVALUATE /*or Equation_Type[Index_Equ_IM[eq]].first==E_EVALUATE_R*/ or Equation_Type[Index_Equ_IM[eq]].first == E_EVALUATE_S)
- {
- if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
- Simulation_Type = EVALUATE_BACKWARD;
- else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
- Simulation_Type = EVALUATE_FORWARD;
- }
- if (i > 0)
- {
- if ((prev_Type == EVALUATE_FORWARD and Simulation_Type == EVALUATE_FORWARD)
- 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));
- }
- else
- Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size)));
- }
- else
- Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size)));
- }
- else
- {
- Type.push_back(make_pair(Simulation_Type, make_pair(Blck_Size, Recurs_Size)));
- }
- prev_Type = Simulation_Type;
- eq += Blck_Size;
- }
- return (Type);
-}
-
-
-map >, pair >, int>
-BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
-{
- map >, pair >, int> Derivatives;
- Derivatives.clear();
- int nb_endo = symbol_table.endo_nbr();
- /*ModelBlock.Block_List[Blck].first_order_determinstic_simulation_derivatives = new*/
- for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++)
- {
- bool *IM=incidencematrix.Get_IM(lag, eEndogenous);
- if(IM)
- {
- for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++)
- {
- int eqr = ModelBlock->Block_List[blck].Equation[eq];
- for(int var = 0; var < ModelBlock->Block_List[blck].Size; var++)
- {
- int varr = ModelBlock->Block_List[blck].Variable[var];
- /*cout << "IM=" << IM << "\n";
- cout << "varr=" << varr << " eqr=" << eqr << " lag=" << lag << "\n";*/
- if(IM[varr+eqr*nb_endo])
- {
- bool OK = true;
- map >, pair >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr)));
- if(its!=Derivatives.end())
- {
- if(its->second == 2)
- OK=false;
- }
-
- if(OK)
- {
- if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eqBlock_List[blck].Nb_Recursives)
- //It's a normalized equation, we have to recompute the derivative using chain rule derivative function*/
- Derivatives[make_pair(make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr))] = 1;
- else
- //It's a feedback equation we can use the derivatives
- Derivatives[make_pair(make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr))] = 0;
- }
- if(varBlock_List[blck].Nb_Recursives)
- {
- int eqs = ModelBlock->Block_List[blck].Equation[var];
- for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; varsBlock_List[blck].Size; vars++)
- {
- int varrs = ModelBlock->Block_List[blck].Variable[vars];
- //A new derivative need to be computed using the chain rule derivative function (a feedback variable appear in a recursive equation)
- if(Derivatives.find(make_pair(make_pair(lag, make_pair(var, vars)), make_pair(eqs, varrs)))!=Derivatives.end())
- Derivatives[make_pair(make_pair(lag, make_pair(eq, vars)), make_pair(eqr, varrs))] = 2;
- }
- }
- }
- }
- }
- }
- }
- return(Derivatives);
-}
-
-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_order_endo_derivatives, bool dynamic, int mfs, double cutoff)
-{
- int i, j, Nb_TotalBlocks, Nb_RecursBlocks, Nb_SimulBlocks;
- BlockType Btype;
- int count_Block, count_Equ;
- bool *SIM0, *SIM00;
-
-
-
- int counted = 0;
- if (prologue+epilogue < n)
- {
- cout << "Normalizing the model ...\n";
- double *max_val = (double *) malloc(n*sizeof(double));
- memset(max_val, 0, n*sizeof(double));
- for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
- {
- if (fabs(iter->second) > max_val[iter->first.first])
- max_val[iter->first.first] = fabs(iter->second);
- }
- for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
- iter->second /= max_val[iter->first.first];
- free(max_val);
- bool OK = false;
- double bi = 0.99999999;
- //double bi=1e-13;
- int suppressed = 0;
- vector Index_Equ_IM_save(Index_Equ_IM);
- while (!OK && bi > 1e-19)
- {
- 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));
- //cout << "---------------------------------\n";
- for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
- {
- //printf("iter->second=% 1.10f iter->first.first=%3d iter->first.second=%3d bi=%f\n", iter->second, iter->first.first, iter->first.second, bi);
- if (fabs(iter->second) > max(bi, cutoff))
- {
- SIM0[iter->first.first*n+iter->first.second] = 1;
- if (!IM_0[iter->first.first*n+iter->first.second])
- {
- cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << " " << iter->second << "\n";
- }
- }
- else
- suppress++;
- }
- for (i = 0; i < n; i++)
- for (j = 0; j < n; j++)
- {
- SIM00[i*n + j] = SIM0[Index_Equ_IM[i] * n + Index_Var_IM[j]];
- }
- free(SIM0);
- if (suppress != suppressed)
- OK = Compute_Normalization(IM, n, prologue, epilogue, 0, SIM00, Index_Equ_IM);
- suppressed = suppress;
- if (!OK)
- //bi/=1.07;
- bi /= 2;
- counted++;
- if (bi > 1e-19)
- free(SIM00);
- }
- if (!OK)
- {
- Compute_Normalization(IM, n, prologue, epilogue, 1, SIM00, Index_Equ_IM);
- Compute_Normalization(IM, n, prologue, epilogue, 2, IM_0, Index_Equ_IM);
- }
- }
- 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);
-
- V_Equation_Type = Equation_Type_determination(equations, first_order_endo_derivatives, Index_Var_IM, Index_Equ_IM, mfs);
-
- cout << "Finding the optimal block decomposition of the model ...\n";
- vector > blocks;
- if (prologue+epilogue < n)
- {
- if(dynamic)
- Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, V_Equation_Type, false, true, mfs);
- else
- Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, V_Equation_Type, false, false, mfs);
- }
-
- t_type Type = Reduce_Blocks_and_type_determination(prologue, epilogue, blocks, equations, V_Equation_Type, Index_Var_IM, Index_Equ_IM);
-
- i = 0;
- j = 0;
- Nb_SimulBlocks = 0;
- int Nb_feedback_variable = 0;
- for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
- {
- if (it->first == SOLVE_FORWARD_COMPLETE || it->first == SOLVE_BACKWARD_COMPLETE || it->first == SOLVE_TWO_BOUNDARIES_COMPLETE)
- {
- Nb_SimulBlocks++;
- if (it->second.first > j)
- {
- j = it->second.first;
- Nb_feedback_variable = blocks[Nb_SimulBlocks-1].second;
- }
- }
- }
-
- Nb_TotalBlocks = Type.size();
- Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
- cout << Nb_TotalBlocks << " block(s) found:\n";
- cout << " " << Nb_RecursBlocks << " recursive block(s) and " << blocks.size() << " simultaneous block(s). \n";
- cout << " the largest simultaneous block has " << j << " equation(s)\n"
- <<" and " << Nb_feedback_variable << " feedback variable(s).\n";
-
- ModelBlock->Size = Nb_TotalBlocks;
- ModelBlock->Periods = periods;
- ModelBlock->Block_List = (Block *) malloc(sizeof(ModelBlock->Block_List[0]) * Nb_TotalBlocks);
-
- count_Equ = count_Block = 0;
-
-
-
- for (t_type::const_iterator it = Type.begin(); it != Type.end(); it++)
- {
- if (count_Equ < prologue)
- Btype = PROLOGUE;
- else if (count_Equ < n-epilogue)
- if (it->second.first == 1)
- Btype = PROLOGUE;
- else
- {
- Btype = SIMULTANS;
- }
- else
- Btype = EPILOGUE;
- Allocate_Block(it->second.first, &count_Equ, count_Block++, Btype, it->first, ModelBlock, V_Equation_Type, it->second.second, Index_Var_IM, Index_Equ_IM);
- }
-}
-
-//------------------------------------------------------------------------------
-// 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_order_endo_derivatives, int mfs, double cutoff)
-{
- 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++)
- {
- Cur_IM = incidencematrix.Get_IM(k, eEndogenous);
- if (Cur_IM)
- {
- for (i = 0; i < symbol_table.endo_nbr()*symbol_table.endo_nbr(); i++)
- {
- SIM[i] = (SIM[i]) || (Cur_IM[i]);
- }
- }
- }
- if (bt_verbose)
- {
- cout << "incidence matrix for the static model (unsorted) \n";
- incidencematrix.Print_SIM(SIM, eEndogenous);
- }
- Index_Equ_IM = vector(symbol_table.endo_nbr());
- 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++)
- {
- Index_Var_IM[i] = i;
- }
- if (ModelBlock != NULL)
- Free_Block(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[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_order_endo_derivatives, true, mfs, cutoff);
- free(SIM_0);
- free(SIM);
-}
diff --git a/BlockTriangular.hh b/BlockTriangular.hh
deleted file mode 100644
index c0ebfb4b..00000000
--- a/BlockTriangular.hh
+++ /dev/null
@@ -1,200 +0,0 @@
-/*
- * Copyright (C) 2007-2008 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see .
- */
-
-#ifndef _BLOCKTRIANGULAR_HH
-#define _BLOCKTRIANGULAR_HH
-
-#include
-#include "CodeInterpreter.hh"
-#include "ExprNode.hh"
-#include "SymbolTable.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;
-
-//! 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;
-
-typedef set temporary_terms_inuse_type;
-
-typedef vector >, pair > > chain_rule_derivatives_type;
-
-//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
-struct IM_compact
-{
- int size, u_init, u_finish, nb_endo, nb_other_endo, size_exo, size_other_endo;
- int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index;
- int *u_other_endo, *Var_other_endo, *Equ_other_endo, *Var_Index_other_endo, *Equ_Index_other_endo;
-};
-
-//! One block of the model
-struct Block
-{
- int Size, Sized, nb_exo, nb_exo_det, nb_other_endo, Nb_Recursives;
- BlockType Type;
- BlockSimulationType Simulation_Type;
- int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo;
- int Max_Lag_Endo, Max_Lead_Endo;
- int Max_Lag_Other_Endo, Max_Lead_Other_Endo;
- int Max_Lag_Exo, Max_Lead_Exo;
- bool is_linear;
- int *Equation, *Own_Derivative;
- EquationType *Equation_Type;
- NodeID *Equation_Normalized;
- int *Variable, *Other_Endogenous, *Exogenous;
- temporary_terms_type **Temporary_Terms_in_Equation;
- //temporary_terms_type *Temporary_terms;
- temporary_terms_inuse_type *Temporary_InUse;
- IM_compact *IM_lead_lag;
- chain_rule_derivatives_type *Chain_Rule_Derivatives;
- int Code_Start, Code_Length;
-};
-
-
-
-//! The set of all blocks of the model
-struct Model_Block
-{
- int Size, Periods;
- Block* Block_List;
- //int *in_Block_Equ, *in_Block_Var, *in_Equ_of_Block, *in_Var_of_Block;
-};
-
-
-//! Creates the incidence matrix, computes prologue & epilogue, normalizes the model and computes the block decomposition
-class BlockTriangular
-{
-private:
- //! Find equations and endogenous variables belonging to the prologue and epilogue of the model
- void Prologue_Epilogue(bool* IM, int &prologue, int &epilogue, int n, vector &Index_Var_IM, vector &Index_Equ_IM, bool* IM0);
- //! Allocates and fills the Model structure describing the content of each block
- void Allocate_Block(int size, int *count_Equ, int count_Block, BlockType type, BlockSimulationType SimType, Model_Block *ModelBlock, t_etype &Equation_Type, int recurs_Size, vector &Index_Var_IM, vector &Index_Equ_IM);
- //! Finds a matching between equations and endogenous variables
- bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, int 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_, bool select_feedback_variable, int mfs) const;
- //! 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_order_endo_derivatives, vector &Index_Var_IM, vector &Index_Equ_IM, int mfs);
- //! 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, vector &Index_Var_IM, vector &Index_Equ_IM);
- //! 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, t_vtype &equation_lead_lag) const;
-public:
- 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;
-
- map >, pair >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
-
-
- void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map >, NodeID> &first_order_endo_derivatives, int mfs, double cutoff);
- 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_order_endo_derivatives, bool dynamic, int mfs, double cutoff);
- vector Index_Equ_IM;
- vector Index_Var_IM;
- int prologue, epilogue;
- bool bt_verbose;
- Model_Block* ModelBlock;
- int periods;
- inline static std::string BlockType0(int type)
- {
- switch (type)
- {
- case 0:
- return ("SIMULTANEOUS TIME SEPARABLE ");
- break;
- case 1:
- return ("PROLOGUE ");
- break;
- case 2:
- return ("EPILOGUE ");
- break;
- case 3:
- return ("SIMULTANEOUS TIME UNSEPARABLE");
- break;
- default:
- return ("UNKNOWN ");
- break;
- }
- };
- inline static std::string BlockSim(int type)
- {
- switch (type)
- {
- case EVALUATE_FORWARD:
- //case EVALUATE_FORWARD_R:
- return ("EVALUATE FORWARD ");
- break;
- case EVALUATE_BACKWARD:
- //case EVALUATE_BACKWARD_R:
- return ("EVALUATE BACKWARD ");
- break;
- case SOLVE_FORWARD_SIMPLE:
- return ("SOLVE FORWARD SIMPLE ");
- break;
- case SOLVE_BACKWARD_SIMPLE:
- return ("SOLVE BACKWARD SIMPLE ");
- break;
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- return ("SOLVE TWO BOUNDARIES SIMPLE ");
- break;
- case SOLVE_FORWARD_COMPLETE:
- return ("SOLVE FORWARD COMPLETE ");
- break;
- case SOLVE_BACKWARD_COMPLETE:
- return ("SOLVE BACKWARD COMPLETE ");
- break;
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- return ("SOLVE TWO BOUNDARIES COMPLETE");
- break;
- default:
- return ("UNKNOWN ");
- break;
- }
- };
- inline static std::string c_Equation_Type(int type)
- {
- char c_Equation_Type[4][13]=
- {
- "E_UNKNOWN ",
- "E_EVALUATE ",
- "E_EVALUATE_S",
- "E_SOLVE "
- };
- return(c_Equation_Type[type]);
- };
-};
-#endif
diff --git a/CodeInterpreter.hh b/CodeInterpreter.hh
index abfa44d1..ae925158 100644
--- a/CodeInterpreter.hh
+++ b/CodeInterpreter.hh
@@ -100,10 +100,10 @@ enum Tags
enum BlockType
{
- SIMULTANS = 0, //!< Simultaneous time separable block
- PROLOGUE = 1, //!< Prologue block (one equation at the beginning, later merged)
- EPILOGUE = 2, //!< Epilogue block (one equation at the beginning, later merged)
- SIMULTAN = 3 //!< Simultaneous time unseparable block
+ SIMULTANS, //!< Simultaneous time separable block
+ PROLOGUE, //!< Prologue block (one equation at the beginning, later merged)
+ EPILOGUE, //!< Epilogue block (one equation at the beginning, later merged)
+ SIMULTAN //!< Simultaneous time unseparable block
};
enum EquationType
@@ -126,7 +126,7 @@ enum BlockSimulationType
SOLVE_TWO_BOUNDARIES_SIMPLE, //!< Block of one equation, newton solver needed, forward & ackward
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
+ SOLVE_TWO_BOUNDARIES_COMPLETE //!< Block of several equations, newton solver needed, forward and backwar
};
//! Enumeration of possible symbol types
@@ -475,9 +475,8 @@ private:
uint8_t op_code;
int size;
uint8_t type;
- int *variable;
- int *equation;
- int *own_derivatives;
+ vector variable;
+ vector equation;
bool is_linear;
vector Block_Contain_;
int endo_nbr;
@@ -485,11 +484,14 @@ private:
int Max_Lead;
int u_count_int;
public:
- inline FBEGINBLOCK_(){ op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; variable = NULL; equation = NULL; own_derivatives = NULL;
+ inline FBEGINBLOCK_(){ op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; /*variable = NULL; equation = NULL;*/
is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0;};
- inline FBEGINBLOCK_(const int size_arg, const BlockSimulationType type_arg, int *variable_arg, int *equation_arg, int *own_derivatives_arg,
- bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg)
- { op_code = FBEGINBLOCK; size = size_arg; type = type_arg; variable = variable_arg; equation = equation_arg; own_derivatives = own_derivatives_arg;
+ inline FBEGINBLOCK_(unsigned int &size_arg, BlockSimulationType &type_arg, int unsigned first_element, int unsigned &block_size,
+ const vector &variable_arg, const vector &equation_arg,
+ bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg)
+ { op_code = FBEGINBLOCK; size = size_arg; type = type_arg;
+ variable = vector(variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size));
+ equation = vector(equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size));
is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg;/*Block_Contain.clear();*/};
inline unsigned int get_size() { return size;};
inline uint8_t get_type() { return type;};
@@ -508,7 +510,6 @@ public:
{
CompileCode.write(reinterpret_cast(&variable[i]), sizeof(variable[0]));
CompileCode.write(reinterpret_cast(&equation[i]), sizeof(equation[0]));
- CompileCode.write(reinterpret_cast(&own_derivatives[i]), sizeof(own_derivatives[0]));
}
if (type==SOLVE_TWO_BOUNDARIES_SIMPLE || type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
type==SOLVE_BACKWARD_COMPLETE || type==SOLVE_FORWARD_COMPLETE)
@@ -532,7 +533,6 @@ public:
Block_contain_type bc;
memcpy(&bc.Variable, code, sizeof(bc.Variable)); code += sizeof(bc.Variable);
memcpy(&bc.Equation, code, sizeof(bc.Equation)); code += sizeof(bc.Equation);
- memcpy(&bc.Own_Derivative, code, sizeof(bc.Own_Derivative)); code += sizeof(bc.Own_Derivative);
Block_Contain_.push_back(bc);
}
if (type==SOLVE_TWO_BOUNDARIES_SIMPLE || type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index dc109d5c..949c21f0 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -876,13 +876,13 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct)
void
PlannerObjectiveStatement::computingPass()
{
- model_tree->computingPass(false, true, false);
+ model_tree->computingPass(eval_context_type(), false, true, false);
}
void
PlannerObjectiveStatement::writeOutput(ostream &output, const string &basename) const
{
- model_tree->writeStaticFile(basename + "_objective", false);
+ model_tree->writeStaticFile(basename + "_objective", false, false);
}
BVARDensityStatement::BVARDensityStatement(int maxnlags_arg, const OptionsList &options_list_arg) :
diff --git a/DynamicModel.cc b/DynamicModel.cc
index ce95029d..90f3d229 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -23,6 +23,7 @@
#include
#include
#include
+#include
#include "DynamicModel.hh"
// For mkdir() and chdir()
@@ -42,9 +43,9 @@ DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
max_exo_lag(0), max_exo_lead(0),
max_exo_det_lag(0), max_exo_det_lead(0),
dynJacobianColsNbr(0),
+ global_temporary_terms(true),
cutoff(1e-15),
- mfs(0),
- block_triangular(symbol_table_arg, num_constants_arg)
+ mfs(0)
{
}
@@ -57,12 +58,10 @@ DynamicModel::AddVariable(int symb_id, int lag)
void
DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const
{
- //first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag)));
if (it != first_derivatives.end())
(it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
else
- /*code_file.write(&FLDZ, sizeof(FLDZ));*/
{
FLDZ_ fldz;
fldz.write(code_file);
@@ -81,240 +80,178 @@ DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr,
FLDZ_ fldz;
fldz.write(code_file);
}
- //code_file.write(&FLDZ, sizeof(FLDZ));
}
void
-DynamicModel::BuildIncidenceMatrix()
-{
- set > endogenous, exogenous;
- for (int eq = 0; eq < (int) equations.size(); eq++)
- {
- BinaryOpNode *eq_node = equations[eq];
- endogenous.clear();
- NodeID Id = eq_node->get_arg1();
- Id->collectEndogenous(endogenous);
- Id = eq_node->get_arg2();
- Id->collectEndogenous(endogenous);
- for (set >::iterator it_endogenous=endogenous.begin();it_endogenous!=endogenous.end();it_endogenous++)
- {
- block_triangular.incidencematrix.fill_IM(eq, it_endogenous->first, it_endogenous->second, eEndogenous);
- }
- exogenous.clear();
- Id = eq_node->get_arg1();
- Id->collectExogenous(exogenous);
- Id = eq_node->get_arg2();
- Id->collectExogenous(exogenous);
- for (set >::iterator it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++)
- {
- block_triangular.incidencematrix.fill_IM(eq, it_exogenous->first, it_exogenous->second, eExogenous);
- }
- }
-}
-
-void
-DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
+DynamicModel::computeTemporaryTermsOrdered()
{
map > first_occurence;
map reference_count;
- int i, j, m, eq, var, eqr, varr, lag;
- temporary_terms_type vect;
- ostringstream tmp_output;
BinaryOpNode *eq_node;
first_derivatives_type::const_iterator it;
first_chain_rule_derivatives_type::const_iterator it_chr;
ostringstream tmp_s;
-
- temporary_terms.clear();
+ v_temporary_terms.clear();
map_idx.clear();
- for (j = 0;j < ModelBlock->Size;j++)
- {
- // Compute the temporary terms reordered
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S && iBlock_List[j].Nb_Recursives && 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);
- else
- {
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx);
- }
- }
- for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
- lag=it.first.first;
- int eqr=it.second.first;
- int varr=it.second.second;
- it_chr=first_chain_rule_derivatives.find(make_pair(eqr, make_pair( varr, lag)));
- it_chr->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
- }
- for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+ unsigned int nb_blocks = getNbBlocks();
+ v_temporary_terms = vector >(nb_blocks);
+ v_temporary_terms_inuse = vector (nb_blocks);
+ temporary_terms.clear();
+
+ if(!global_temporary_terms)
+ {
+ for (unsigned int block = 0; block < nb_blocks; block++)
{
- lag=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++)
+ reference_count.clear();
+ temporary_terms.clear();
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_nb_mfs = getBlockMfs(block);
+ unsigned int block_nb_recursives = block_size - block_nb_mfs;
+ v_temporary_terms[block] = vector(block_size);
+ for (unsigned int i = 0; i < block_size; i++)
{
- eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
- it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
- }
- }
- /*for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
- lag=it.first.first;
- eqr=it.second.first;
- varr=it.second.second;
- it_chr=first_chain_rule_derivatives.find(make_pair(eqr, make_pair( varr, lag)));
- it_chr->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
- }*/
- /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- lag=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
- {
- eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
- var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
- it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eExogenous, var), lag)));
- it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
- }
- }*/
- //jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
- for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- lag=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++)
+ if (icomputeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
+ else
{
- eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
- var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
- it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
- it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
+ eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
+ eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
}
}
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+ {
+ NodeID id=it->second.second;
+ id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+ }
+ for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+ for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+
+ set temporary_terms_in_use;
+ temporary_terms_in_use.clear();
+ v_temporary_terms_inuse[block] = temporary_terms_in_use;
}
}
- for (j = 0;j < ModelBlock->Size;j++)
+ else
{
- // Collecte the temporary terms reordered
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ for (unsigned int block = 0; block < nb_blocks; block++)
{
- if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S && iBlock_List[j].Nb_Recursives && ModelBlock->Block_List[j].Equation_Normalized[i])
- ModelBlock->Block_List[j].Equation_Normalized[i]->collectTemporary_terms(temporary_terms, ModelBlock, j);
- else
+ // Compute the temporary terms reordered
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_nb_mfs = getBlockMfs(block);
+ unsigned int block_nb_recursives = block_size - block_nb_mfs;
+ v_temporary_terms[block] = vector(block_size);
+ for (unsigned int i = 0; i < block_size; i++)
{
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- eq_node->collectTemporary_terms(temporary_terms, ModelBlock, j);
- }
-
- /*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++)
- {
- lag=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++)
- {
- eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
- //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
- //if(it!=first_derivatives.end())
- it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
- }
- }
- for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
- lag=it.first.first;
- eqr=it.second.first;
- varr=it.second.second;
- it_chr=first_chain_rule_derivatives.find(make_pair(eqr, make_pair( varr, lag)));
- it_chr->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
- }
- /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- lag=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size_exo;i++)
- {
- eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i];
- var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i];
- it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eExogenous, var), lag)));
- //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
- it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
- }
- }*/
- //jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr;
- for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- lag=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++)
+ if (icomputeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
+ else
{
- eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i];
- var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i];
- it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var), lag)));
- //it=first_derivatives.find(make_pair(eq,variable_table.getID(var, lag)));
- //if(it!=first_derivatives.end())
- it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
+ eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
+ eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
}
}
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+ {
+ NodeID id=it->second.second;
+ id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+ }
+ for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
+ for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
+ it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
}
+ for (unsigned int block = 0; block < nb_blocks; block++)
+ {
+ // Collect the temporary terms reordered
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_nb_mfs = getBlockMfs(block);
+ unsigned int block_nb_recursives = block_size - block_nb_mfs;
+ set temporary_terms_in_use;
+ for (unsigned int i = 0; i < block_size; i++)
+ {
+ if (icollectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+ else
+ {
+ eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
+ eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+ }
+ }
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+ {
+ NodeID id=it->second.second;
+ id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+ }
+ for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
+ it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+ for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
+ it->second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+ v_temporary_terms_inuse[block] = temporary_terms_in_use;
+ }
+ // Add a mapping form node ID to temporary terms order
+ int j=0;
+ for (temporary_terms_type::const_iterator it = temporary_terms.begin();
+ it != temporary_terms.end(); it++)
+ map_idx[(*it)->idx]=j++;
}
- // Add a mapping form node ID to temporary terms order
- j=0;
- for (temporary_terms_type::const_iterator it = temporary_terms.begin();
- it != temporary_terms.end(); it++)
- map_idx[(*it)->idx]=j++;
}
void
-DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string &dynamic_basename) const
+DynamicModel::writeModelEquationsOrdered_M(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;
+ temporary_terms_type local_temporary_terms;
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;
+ ExprNodeOutputType local_output_type;
+
+ if(global_temporary_terms)
+ {
+ local_output_type = oMatlabDynamicModelSparse;
+ local_temporary_terms = temporary_terms;
+ }
+ else
+ local_output_type = oMatlabDynamicModelSparseLocalTemporaryTerms;
+
//----------------------------------------------------------------------
//For each block
- for (j = 0;j < ModelBlock->Size;j++)
+ for (unsigned int block = 0; block < getNbBlocks(); block++)
{
+
//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++)
+ nze = derivative_endo[block].size();
+ nze_other_endo = derivative_other_endo[block].size();
+ nze_exo = derivative_exo[block].size();
+ BlockSimulationType simulation_type = getBlockSimulationType(block);
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_mfs = getBlockMfs(block);
+ unsigned int block_recursive = block_size - block_mfs;
+ unsigned int block_exo_size = exo_block[block].size();
+ unsigned int block_exo_det_size = exo_det_block[block].size();
+ unsigned int block_other_endo_size = other_endo_block[block].size();
+ int block_max_lag=max_leadlag_block[block].first;
+ if(global_temporary_terms)
{
- 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;
+ local_output_type = oMatlabDynamicModelSparse;
+ local_temporary_terms = temporary_terms;
}
+ else
+ local_output_type = oMatlabDynamicModelSparseLocalTemporaryTerms;
+
tmp1_output.str("");
- tmp1_output << dynamic_basename << "_" << j+1 << ".m";
+ tmp1_output << dynamic_basename << "_" << block+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";
@@ -322,183 +259,196 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
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*/)
+ if (simulation_type ==EVALUATE_BACKWARD || simulation_type ==EVALUATE_FORWARD)
{
- output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << j+1 << "(y, x, params, jacobian_eval, y_kmin, periods)\n";
+ output << "function [y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+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 if (simulation_type ==SOLVE_FORWARD_COMPLETE || simulation_type ==SOLVE_BACKWARD_COMPLETE)
+ output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, it_, jacobian_eval)\n";
+ else if (simulation_type ==SOLVE_BACKWARD_SIMPLE || simulation_type ==SOLVE_FORWARD_SIMPLE)
+ output << "function [residual, y, g1, g2, g3, varargout] = " << dynamic_basename << "_" << block+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 << "function [residual, y, g1, g2, g3, b, varargout] = " << dynamic_basename << "_" << block+1 << "(y, x, params, periods, jacobian_eval, y_kmin, y_size)\n";
+ BlockType block_type;
+ if(simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE ||simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
+ block_type = SIMULTAN;
+ else if(simulation_type == SOLVE_FORWARD_COMPLETE ||simulation_type == SOLVE_BACKWARD_COMPLETE)
+ block_type = SIMULTANS;
+ else if((simulation_type == SOLVE_FORWARD_SIMPLE ||simulation_type == SOLVE_BACKWARD_SIMPLE ||
+ simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
+ && getBlockFirstEquation(block) < prologue)
+ block_type = PROLOGUE;
+ else if((simulation_type == SOLVE_FORWARD_SIMPLE ||simulation_type == SOLVE_BACKWARD_SIMPLE ||
+ simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
+ && getBlockFirstEquation(block) >= equations.size() - epilogue)
+ block_type = EPILOGUE;
+ else
+ block_type = SIMULTANS;
output << " % ////////////////////////////////////////////////////////////////////////" << endl
- << " % //" << string(" Block ").substr(int(log10(j + 1))) << j + 1 << " " << BlockTriangular::BlockType0(ModelBlock->Block_List[j].Type)
+ << " % //" << string(" Block ").substr(int(log10(block + 1))) << block + 1 << " " << BlockType0(block_type)
<< " //" << endl
<< " % // Simulation type "
- << BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //" << endl
+ << BlockSim(simulation_type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << endl;
output << " global options_;" << endl;
//The Temporary terms
- //output << " relax = 1;\n";
- if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD
- ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
+ if (simulation_type ==EVALUATE_BACKWARD || simulation_type ==EVALUATE_FORWARD)
{
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 << " g1 = spalloc(" << block_mfs << ", " << block_mfs*(1+getBlockMaxLag(block)+getBlockMaxLead(block)) << ", " << nze << ");\n";
+ output << " g1_x=spalloc(" << block_size << ", " << (block_exo_size + block_exo_det_size)*
+ (1+max(exo_det_max_leadlag_block[block].first, exo_max_leadlag_block[block].first)+max(exo_det_max_leadlag_block[block].second, exo_max_leadlag_block[block].second))
+ << ", " << nze_exo << ");\n";
+ output << " g1_o=spalloc(" << block_size << ", " << block_other_endo_size*
+ (1+other_endo_max_leadlag_block[block].first+other_endo_max_leadlag_block[block].second)
+ << ", " << 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 << " g1 = spalloc(" << block_size << ", " << block_size*(1+getBlockMaxLag(block)+getBlockMaxLead(block)) << ", " << nze << ");\n";
+ output << " g1_x=spalloc(" << block_size << ", " << (block_exo_size + block_exo_det_size)*
+ (1+max(exo_det_max_leadlag_block[block].first, exo_max_leadlag_block[block].first)+max(exo_det_max_leadlag_block[block].second, exo_max_leadlag_block[block].second))
+ << ", " << nze_exo << ");\n";
+ output << " g1_o=spalloc(" << block_size << ", " << block_other_endo_size*
+ (1+other_endo_max_leadlag_block[block].first+other_endo_max_leadlag_block[block].second)
+ << ", " << 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)
+ if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
- output << " g1 = spalloc(" << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives) << "*options_.periods, "
- << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives) << "*(options_.periods+" << ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1 << ")"
+ output << " g1 = spalloc(" << block_mfs << "*options_.periods, "
+ << block_mfs << "*(options_.periods+" << max_leadlag_block[block].first+max_leadlag_block[block].second+1 << ")"
<< ", " << nze << "*options_.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 << " g1 = spalloc(" << block_mfs
+ << ", " << block_mfs << ", " << nze << ");\n";
+ output << " g1_tmp_r = spalloc(" << block_recursive
+ << ", " << block_size << ", " << nze << ");\n";
+ output << " g1_tmp_b = spalloc(" << block_mfs
+ << ", " << block_size << ", " << nze << ");\n";
}
output << " end;\n";
}
output << " g2=0;g3=0;\n";
- if (ModelBlock->Block_List[j].Temporary_InUse->size())
+ if (v_temporary_terms_inuse[block].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++)
+ for (temporary_terms_inuse_type::const_iterator it = v_temporary_terms_inuse[block].begin();
+ it != v_temporary_terms_inuse[block].end(); it++)
tmp_output << " T" << *it;
output << " global" << tmp_output.str() << ";\n";
}
- if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD && 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)
+ if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+ {
+ temporary_terms_type tt2;
+ tt2.clear();
+ for(int i=0; i< (int) block_size; i++)
+ {
+ if (v_temporary_terms[block][i].size() && global_temporary_terms)
+ {
+ output << " " << "% //Temporary variables initialization" << endl
+ << " " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
+ for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
+ it != v_temporary_terms[block][i].end(); it++)
+ {
+ output << " ";
+ (*it)->writeOutput(output, oMatlabDynamicModel, local_temporary_terms);
+ output << " = T_zeros;" << endl;
+ }
+ }
+ }
+ }
+ if (simulation_type==SOLVE_BACKWARD_SIMPLE || simulation_type==SOLVE_FORWARD_SIMPLE || simulation_type==SOLVE_BACKWARD_COMPLETE || simulation_type==SOLVE_FORWARD_COMPLETE)
+ output << " residual=zeros(" << block_mfs << ",1);\n";
+ else if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+ output << " residual=zeros(" << block_mfs << ",y_kmin+periods);\n";
+ if (simulation_type==EVALUATE_BACKWARD)
output << " for it_ = (y_kmin+periods):y_kmin+1\n";
- if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
+ if (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)
+ if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || 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";
+ output << " b = zeros(periods*y_size,1);" << endl
+ << " for it_ = y_kmin+1:(periods+y_kmin)" << endl
+ << " Per_y_=it_*y_size;" << endl
+ << " Per_J_=(it_-y_kmin-1)*y_size;" << endl
+ << " Per_K_=(it_-1)*y_size;" << endl;
sps=" ";
}
else
- if (ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD || ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD )
+ if (simulation_type==EVALUATE_BACKWARD || simulation_type==EVALUATE_FORWARD )
sps = " ";
else
sps="";
// The equations
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ for (unsigned int i = 0; i < block_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++)
+ if (v_temporary_terms[block].size())
{
- output << " " << sps;
- (*it)->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << " = ";
- (*it)->writeOutput(output, oMatlabDynamicModelSparse, tt2);
- // Insert current node into tt2
- tt2.insert(*it);
- output << ";" << endl;
+ output << " " << "% //Temporary variables" << endl;
+ for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
+ it != v_temporary_terms[block][i].end(); it++)
+ {
+ output << " " << sps;
+ (*it)->writeOutput(output, local_output_type, local_temporary_terms);
+ output << " = ";
+ (*it)->writeOutput(output, local_output_type, 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]];
+
+ int variable_ID = getBlockVariableID(block, i);
+ int equation_ID = getBlockEquationID(block, i);
+ EquationType equ_type = getBlockEquationType(block, i);
+ string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, variable_ID)) ;
+ eq_node = (BinaryOpNode*)getBlockEquationNodeID(block,i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
tmp_output.str("");
- /*if((ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) and (iBlock_List[j].Nb_Recursives))
- lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
- else*/
- lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
- switch (ModelBlock->Block_List[j].Simulation_Type)
+ lhs->writeOutput(tmp_output, local_output_type, local_temporary_terms);
+ switch (simulation_type)
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
-evaluation: if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
- 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;
+evaluation: if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE)
+ output << " % equation " << getBlockEquationID(block, i)+1 << " variable : " << sModel
+ << " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
output << " ";
- if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
+ if (equ_type == E_EVALUATE)
{
output << tmp_output.str();
output << " = ";
- /*if(!(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD))
- {
- lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << "-relax*(";
- lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << "-(";
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << "))";
- }
- else*/
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ rhs->writeOutput(output, local_output_type, local_temporary_terms);
}
- else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
+ else if (equ_type == E_EVALUATE_S)
{
output << "%" << tmp_output.str();
output << " = ";
- if (ModelBlock->Block_List[j].Equation_Normalized[i])
+ if (isBlockEquationRenormalized(block, i))
{
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ rhs->writeOutput(output, local_output_type, local_temporary_terms);
output << "\n ";
tmp_output.str("");
- eq_node = (BinaryOpNode *)ModelBlock->Block_List[j].Equation_Normalized[i];
+ eq_node = (BinaryOpNode *)getBlockEquationRenormalizedNodeID(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ lhs->writeOutput(output, local_output_type, local_temporary_terms);
output << " = ";
- /*if(!(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD))
- {
- lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << "-relax*(";
- lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << "-(";
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << "))";
- }
- else*/
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ rhs->writeOutput(output, local_output_type, local_temporary_terms);
}
}
else
{
- cerr << "Type missmatch for equation " << ModelBlock->Block_List[j].Equation[i]+1 << "\n";
+ cerr << "Type missmatch for equation " << equation_ID+1 << "\n";
exit(EXIT_FAILURE);
}
output << ";\n";
@@ -507,342 +457,251 @@ evaluation: if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDAR
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 << ") = (";
+ 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_) = (";
+ if (iwriteOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ rhs->writeOutput(output, local_output_type, local_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)
+ if (simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || 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)
+ if (simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE || 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)
+ if (simulation_type==SOLVE_BACKWARD_SIMPLE || simulation_type==SOLVE_FORWARD_SIMPLE ||
+ simulation_type==SOLVE_BACKWARD_COMPLETE || 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)
+ switch (simulation_type)
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
- for (m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
+ for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
{
- 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(symbol_table.getID(eEndogenous, var))
- << "(" << k//variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0]))
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ int varr = getBlockInitialVariableID(block, var);
+
+ NodeID id = it->second;
+
+ output << " g1(" << eqr+1 << ", " << varr+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << lag
<< ") " << 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++)
+ for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
- 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;
- }
- }
+ int lag = it->first.first;
+ int eq = it->first.second.first;
+ int var = it->first.second.second;
+ int eqr = getBlockInitialEquationID(block, eq);
+ NodeID id = it->second;
+
+ output << " g1_o(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << lag
+ << ") " << 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:
- for (m=0;mBlock_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
+ for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
{
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++)
- {
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
- output << " g1(" << eq+1 << ", "
- << var+1 + m*(ModelBlock->Block_List[j].Size) << ") = ";
- writeDerivative(output, eqr, symbol_table.getID(eEndogenous, varr), k, oMatlabDynamicModelSparse, temporary_terms);
- output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
- << "(" << k << ") " << varr+1
- << ", equation=" << eqr+1 << endl;
- }
+ int lag = it->first.first;
+ unsigned int eq = it->first.second.first;
+ unsigned int var = it->first.second.second;
+ NodeID id = it->second;
+
+ output << " g1(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
}
- /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+
+ for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
- 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;
- }
- }
+ int lag = it->first.first;
+ unsigned int eq = it->first.second.first;
+ unsigned int var = it->first.second.second;
+ NodeID id = it->second;
+
+ output << " g1_o(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << lag
+ << ") " << var+1
+ << ", equation=" << eq+1 << endl;
}
output << " varargout{1}=g1_x;\n";
output << " varargout{2}=g1_o;\n";
output << " else" << endl;
-
- for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
{
- pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
- k=it.first.first;
- int eq=it.first.second.first;
- int var=it.first.second.second;
- int eqr=it.second.first;
- int varr=it.second.second;
- output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
- << var+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
- writeChainRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
- output << "; %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
- << "(" << k << ") " << varr+1 << ", equation=" << eqr+1 << endl;
+ unsigned int eq = it->first.first;
+ unsigned int var = it->first.second;
+ unsigned int eqr = getBlockEquationID(block, eq);
+ unsigned int varr = getBlockVariableID(block, var);
+ NodeID id = it->second.second;
+ int lag = it->second.first;
+ output << " g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
+ << "(" << lag
+ << ") " << varr+1
+ << ", equation=" << eqr+1 << endl;
}
output << " end;\n";
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
output << " if ~jacobian_eval" << endl;
- for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
{
- pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
- k=it.first.first;
- int eq=it.first.second.first;
- int var=it.first.second.second;
- int eqr=it.second.first;
- int varr=it.second.second;
+ unsigned int eq = it->first.first;
+ unsigned int var = it->first.second;
+ unsigned int eqr = getBlockEquationID(block, eq);
+ unsigned int varr = getBlockVariableID(block, var);
ostringstream tmp_output;
- if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
+ NodeID id = it->second.second;
+ int lag = it->second.first;
+ if(eq>=block_recursive and var>=block_recursive)
{
- if (k==0)
- Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
- << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
+ if (lag==0)
+ Uf[eqr] << "+g1(" << eq+1-block_recursive
+ << "+Per_J_, " << var+1-block_recursive
<< "+Per_K_)*y(it_, " << varr+1 << ")";
- else if (k==1)
- Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
- << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
+ else if (lag==1)
+ Uf[eqr] << "+g1(" << eq+1-block_recursive
+ << "+Per_J_, " << var+1-block_recursive
<< "+Per_y_)*y(it_+1, " << varr+1 << ")";
- else if (k>0)
- Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
- << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
- << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << varr+1 << ")";
- else if (k<0)
- Uf[ModelBlock->Block_List[j].Equation[eq]] << "+g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives
- << "+Per_J_, " << var+1-ModelBlock->Block_List[j].Nb_Recursives
- << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << varr+1 << ")";
- if (k==0)
- tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
- << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_K_) = ";
- else if (k==1)
- tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
- << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_y_) = ";
- else if (k>0)
- tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
- << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_+" << k-1 << ")) = ";
- else if (k<0)
- tmp_output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << "+Per_J_, "
- << var+1-ModelBlock->Block_List[j].Nb_Recursives << "+y_size*(it_" << k-1 << ")) = ";
+ else if (lag>0)
+ Uf[eqr] << "+g1(" << eq+1-block_recursive
+ << "+Per_J_, " << var+1-block_recursive
+ << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")";
+ else if (lag<0)
+ Uf[eqr] << "+g1(" << eq+1-block_recursive
+ << "+Per_J_, " << var+1-block_recursive
+ << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")";
+ if (lag==0)
+ tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, "
+ << var+1-block_recursive << "+Per_K_) = ";
+ else if (lag==1)
+ tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, "
+ << var+1-block_recursive << "+Per_y_) = ";
+ else if (lag>0)
+ tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, "
+ << var+1-block_recursive << "+y_size*(it_+" << lag-1 << ")) = ";
+ else if (lag<0)
+ tmp_output << " g1(" << eq+1-block_recursive << "+Per_J_, "
+ << var+1-block_recursive << "+y_size*(it_" << lag-1 << ")) = ";
output << " " << tmp_output.str();
-
- writeChainRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
-
+ id->writeOutput(output, local_output_type, local_temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
- << "(" << k << ") " << varr+1
+ << "(" << lag << ") " << varr+1
<< ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
}
+
#ifdef CONDITION
output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++)
+ for (unsigned int i = 0; i < block_size; i++)
{
- if (i>=ModelBlock->Block_List[j].Nb_Recursives)
- output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n";
+ if (i>=block_recursive)
+ output << " " << Uf[getBlockEquationID(block, i)].str() << ";\n";
#ifdef CONDITION
output << " if (fabs(condition(" << i+1 << "))Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+ for (m=0;m<=ModelBlock->Block_List[block].Max_Lead+ModelBlock->Block_List[block].Max_Lag;m++)
{
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[m].size;i++)
+ k=m-ModelBlock->Block_List[block].Max_Lag;
+ for (i=0;iBlock_List[block].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];
+ unsigned int eq=ModelBlock->Block_List[block].IM_lead_lag[m].Equ_Index[i];
+ unsigned int var=ModelBlock->Block_List[block].IM_lead_lag[m].Var_Index[i];
+ unsigned int u=ModelBlock->Block_List[block].IM_lead_lag[m].u[i];
+ unsigned int eqr=ModelBlock->Block_List[block].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++)
+ for (i = 0;i < ModelBlock->Block_List[block].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++)
+
+ for (t_derivative::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
{
- 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
+ int lag = it->first.first;
+ unsigned int eq = it->first.second.first;
+ unsigned int var = it->first.second.second;
+ NodeID id = it->second;
+ output << " g1(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << lag
+ << ") " << 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++)
+ for (t_derivative::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
- 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;
- }
- }
+ int lag = it->first.first;
+ unsigned int eq = it->first.second.first;
+ unsigned int var = it->first.second.second;
+ NodeID id = it->second;
+
+ output << " g1_o(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
+ id->writeOutput(output, local_output_type, local_temporary_terms);
+ output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
+ << "(" << lag
+ << ") " << 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:
@@ -853,7 +712,7 @@ end:
}
void
-DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model_Block *ModelBlock, const string bin_basename, map_idx_type map_idx) const
+DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const
{
struct Uff_l
{
@@ -866,7 +725,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
Uff_l *Ufl, *Ufl_First;
};
- int i,j,k,v;
+ int i,v;
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
@@ -876,6 +735,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
map reference_count;
vector feedback_variables;
bool file_open=false;
+
string main_name=file_name;
main_name+=".cod";
code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate );
@@ -885,93 +745,61 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
exit(EXIT_FAILURE);
}
//Temporary variables declaration
- /*code_file.write(&FDIMT, sizeof(FDIMT));
- k=temporary_terms.size();
- code_file.write(reinterpret_cast(&k),sizeof(k));*/
+
FDIMT_ fdimt(temporary_terms.size());
fdimt.write(code_file);
- for (j = 0; j < ModelBlock->Size ;j++)
+ for (unsigned int block = 0; block < getNbBlocks(); block++)
{
feedback_variables.clear();
- if (j>0)
+ if (block>0)
{
FENDBLOCK_ fendblock;
fendblock.write(code_file);
- //code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
}
int count_u;
int u_count_int=0;
- 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)
+ BlockSimulationType simulation_type = getBlockSimulationType(block);
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_mfs = getBlockMfs(block);
+ unsigned int block_recursive = block_size - block_mfs;
+ int block_max_lag=max_leadlag_block[block].first;
+
+ if (simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
+ simulation_type==SOLVE_BACKWARD_COMPLETE || simulation_type==SOLVE_FORWARD_COMPLETE)
{
- //cout << "ModelBlock->Block_List[j].Nb_Recursives = " << ModelBlock->Block_List[j].Nb_Recursives << "\n";
- 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);
- //cout << "u_count_int=" << u_count_int << "\n";
-
-
- /*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;
- 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));
-
- v=u_count_int;
- code_file.write(reinterpret_cast(&v),sizeof(v));*/
+ Write_Inf_To_Bin_File(file_name, bin_basename, block, u_count_int,file_open,
+ simulation_type==SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type==SOLVE_TWO_BOUNDARIES_SIMPLE);
file_open=true;
}
- FBEGINBLOCK_ fbeginblock(ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives,
- ModelBlock->Block_List[j].Simulation_Type,
- ModelBlock->Block_List[j].Variable,
- ModelBlock->Block_List[j].Equation,
- ModelBlock->Block_List[j].Own_Derivative,
- ModelBlock->Block_List[j].is_linear,
+ FBEGINBLOCK_ fbeginblock(block_mfs,
+ simulation_type,
+ getBlockFirstEquation(block),
+ block_size,
+ variable_reordered,
+ equation_reordered,
+ blocks_linear[block],
symbol_table.endo_nbr(),
- ModelBlock->Block_List[j].Max_Lag,
- ModelBlock->Block_List[j].Max_Lead,
+ block_max_lag,
+ block_max_lag,
u_count_int
);
fbeginblock.write(code_file);
- /*code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK));
- v=ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives;
- //cout << "v (Size) = " << v << "\n";
- code_file.write(reinterpret_cast(&v),sizeof(v));
- v=ModelBlock->Block_List[j].Simulation_Type;
- code_file.write(reinterpret_cast(&v),sizeof(v));
- for (i=ModelBlock->Block_List[j].Nb_Recursives; i < ModelBlock->Block_List[j].Size;i++)
+ // The equations
+ for (i = 0;i < (int) block_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]));
- }*/
-
- // The equations
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
+ //The Temporary terms
+ temporary_terms_type tt2;
+ tt2.clear();
+ if (v_temporary_terms[block][i].size())
{
- //The Temporary terms
- temporary_terms_type tt2;
- tt2.clear();
-#ifdef DEBUGC
- 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++)
+ for (temporary_terms_type::const_iterator it = v_temporary_terms[block][i].begin();
+ it != v_temporary_terms[block][i].end(); it++)
{
(*it)->compile(code_file, false, tt2, map_idx, true, false);
-
FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
fstpt.write(code_file);
-
- /*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
@@ -982,204 +810,167 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
#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 = v_temporary_terms[block][i].begin();
+ it != v_temporary_terms[block][i].end(); it++)
+ {
+ map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
+ cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
+ }
#endif
- switch (ModelBlock->Block_List[j].Simulation_Type)
- {
+
+ int variable_ID, equation_ID;
+ EquationType equ_type;
+
+
+ switch (simulation_type)
+ {
evaluation:
- case EVALUATE_BACKWARD:
- case EVALUATE_FORWARD:
- if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
- {
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
- lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
- }
- else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
- {
- eq_node = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[i];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
- lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
- }
- break;
- case SOLVE_BACKWARD_COMPLETE:
- case SOLVE_FORWARD_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- if (iBlock_List[j].Nb_Recursives)
- goto evaluation;
- feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
- v=ModelBlock->Block_List[j].Equation[i];
- Uf[v].Ufl=NULL;
- goto end;
- default:
-end:
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
+ case EVALUATE_BACKWARD:
+ case EVALUATE_FORWARD:
+ equ_type = getBlockEquationType(block, i);
+ if (equ_type == E_EVALUATE)
+ {
+ eq_node = (BinaryOpNode*)getBlockEquationNodeID(block,i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- lhs->compile(code_file, false, temporary_terms, map_idx, true, false);
rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
-
- FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
- /*code_file.write(&FBINARY, sizeof(FBINARY));
- int v=oMinus;
- code_file.write(reinterpret_cast(&v),sizeof(v));*/
- FSTPR_ fstpr(i - ModelBlock->Block_List[j].Nb_Recursives);
- fstpr.write(code_file);
- /*code_file.write(&FSTPR, sizeof(FSTPR));
- v = i - ModelBlock->Block_List[j].Nb_Recursives;
- code_file.write(reinterpret_cast(&v), sizeof(v));*/
+ lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
}
- }
- FENDEQU_ fendequ;
- fendequ.write(code_file);
- //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)
- {
- switch (ModelBlock->Block_List[j].Simulation_Type)
+ else if (equ_type == E_EVALUATE_S)
{
- 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);
- {
- FSTPG_ fstpg(0);
- fstpg.write(code_file);
- }
- /*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:
- case SOLVE_TWO_BOUNDARIES_COMPLETE:
- case SOLVE_TWO_BOUNDARIES_SIMPLE:
- //count_u=ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives;
- count_u = feedback_variables.size();
- for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- pair< pair >, pair > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
- k=it.first.first;
- int eq=it.first.second.first;
- int var=it.first.second.second;
- int eqr=it.second.first;
- int varr=it.second.second;
- //cout << "k=" << k << " eq=" << eq << " (" << eq-ModelBlock->Block_List[j].Nb_Recursives << ") var=" << var << " (" << var-ModelBlock->Block_List[j].Nb_Recursives << ") eqr=" << eqr << " varr=" << varr << " count_u=" << count_u << "\n";
- int v=ModelBlock->Block_List[j].Equation[eq];
- /*m = ModelBlock->Block_List[j].Max_Lag + k;
- int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];*/
- if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
- {
- 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=count_u;
- Uf[v].Ufl->var=varr;
- Uf[v].Ufl->lag=k;
- compileChainRuleDerivative(code_file, eqr, varr, k, map_idx);
-
- FSTPU_ fstpu(count_u);
- fstpu.write(code_file);
-
- /*code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast(&count_u), sizeof(count_u));*/
- count_u++;
- }
- }
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- if(i>=ModelBlock->Block_List[j].Nb_Recursives)
- {
- FLDR_ fldr(i-ModelBlock->Block_List[j].Nb_Recursives);
- fldr.write(code_file);
- /*code_file.write(&FLDR, sizeof(FLDR));
- v = i-ModelBlock->Block_List[j].Nb_Recursives;
- code_file.write(reinterpret_cast(&v), sizeof(v));*/
-
- FLDZ_ fldz;
- fldz.write(code_file);
- //code_file.write(&FLDZ, sizeof(FLDZ));
-
- 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)
- {
- FLDU_ fldu(Uf[v].Ufl->u);
- fldu.write(code_file);
- /*code_file.write(&FLDU, sizeof(FLDU));
- code_file.write(reinterpret_cast(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));*/
- FLDV_ fldv(eEndogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag);
- fldv.write(code_file);
-
- /*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));*/
- FBINARY_ fbinary(oTimes);
- fbinary.write(code_file);
- /*code_file.write(&FBINARY, sizeof(FBINARY));
- v1=oTimes;
- code_file.write(reinterpret_cast(&v1), sizeof(v1));*/
-
- FCUML_ fcuml;
- fcuml.write(code_file);
- //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;
- }
- FBINARY_ fbinary(oMinus);
- fbinary.write(code_file);
- /*code_file.write(&FBINARY, sizeof(FBINARY));
- v=oMinus;
- code_file.write(reinterpret_cast(&v), sizeof(v));*/
-
- FSTPU_ fstpu(i - ModelBlock->Block_List[j].Nb_Recursives);
- fstpu.write(code_file);
- /*code_file.write(&FSTPU, sizeof(FSTPU));
- v = i - ModelBlock->Block_List[j].Nb_Recursives;
- code_file.write(reinterpret_cast(&v), sizeof(v));*/
- }
- }
- break;
- default:
- break;
+ eq_node = (BinaryOpNode*)getBlockEquationRenormalizedNodeID(block,i);
+ lhs = eq_node->get_arg1();
+ rhs = eq_node->get_arg2();
+ rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
+ lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
}
+ break;
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FORWARD_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ if (i< (int) block_recursive)
+ goto evaluation;
+ variable_ID = getBlockVariableID(block, i);
+ equation_ID = getBlockEquationID(block, i);
+ feedback_variables.push_back(variable_ID);
+ Uf[equation_ID].Ufl=NULL;
+ goto end;
+ default:
+end:
+ eq_node = (BinaryOpNode*)getBlockEquationNodeID(block, i);
+ lhs = eq_node->get_arg1();
+ rhs = eq_node->get_arg2();
+ lhs->compile(code_file, false, temporary_terms, map_idx, true, false);
+ rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
+
+ FBINARY_ fbinary(oMinus);
+ fbinary.write(code_file);
+ FSTPR_ fstpr(i - block_recursive);
+ fstpr.write(code_file);
}
+ }
+ FENDEQU_ fendequ;
+ fendequ.write(code_file);
+ // The Jacobian if we have to solve the block
+ if (simulation_type != EVALUATE_BACKWARD
+ && simulation_type != EVALUATE_FORWARD)
+ {
+ switch (simulation_type)
+ {
+ case SOLVE_BACKWARD_SIMPLE:
+ case SOLVE_FORWARD_SIMPLE:
+ compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
+ {
+ FSTPG_ fstpg(0);
+ fstpg.write(code_file);
+ }
+ break;
+
+ case SOLVE_BACKWARD_COMPLETE:
+ case SOLVE_FORWARD_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_COMPLETE:
+ case SOLVE_TWO_BOUNDARIES_SIMPLE:
+ count_u = feedback_variables.size();
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
+ {
+ unsigned int eq = it->first.first;
+ unsigned int var = it->first.second;
+ unsigned int eqr = getBlockEquationID(block, eq);
+ unsigned int varr = getBlockVariableID(block, var);
+ int lag = it->second.first;
+ if(eq>=block_recursive and var>=block_recursive)
+ {
+ if (!Uf[eqr].Ufl)
+ {
+ Uf[eqr].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
+ Uf[eqr].Ufl_First=Uf[eqr].Ufl;
+ }
+ else
+ {
+ Uf[eqr].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
+ Uf[eqr].Ufl=Uf[eqr].Ufl->pNext;
+ }
+ Uf[eqr].Ufl->pNext=NULL;
+ Uf[eqr].Ufl->u=count_u;
+ Uf[eqr].Ufl->var=varr;
+ Uf[eqr].Ufl->lag=lag;
+ compileChainRuleDerivative(code_file, eqr, varr, lag, map_idx);
+ FSTPU_ fstpu(count_u);
+ fstpu.write(code_file);
+ count_u++;
+ }
+ }
+ for (i = 0;i < (int) block_size;i++)
+ {
+ if(i>= (int) block_recursive)
+ {
+ FLDR_ fldr(i-block_recursive);
+ fldr.write(code_file);
+
+ FLDZ_ fldz;
+ fldz.write(code_file);
+
+ v=getBlockEquationID(block, i);
+ for (Uf[v].Ufl=Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl=Uf[v].Ufl->pNext)
+ {
+ FLDU_ fldu(Uf[v].Ufl->u);
+ fldu.write(code_file);
+ FLDV_ fldv(eEndogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag);
+ fldv.write(code_file);
+
+ FBINARY_ fbinary(oTimes);
+ fbinary.write(code_file);
+
+ FCUML_ fcuml;
+ fcuml.write(code_file);
+ }
+ 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;
+ }
+ FBINARY_ fbinary(oMinus);
+ fbinary.write(code_file);
+
+ FSTPU_ fstpu(i - block_recursive);
+ fstpu.write(code_file);
+ }
+ }
+ break;
+ default:
+ break;
+ }
+ }
}
FENDBLOCK_ fendblock;
fendblock.write(code_file);
- //code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
FEND_ fend;
fend.write(code_file);
- //code_file.write(&FEND, sizeof(FEND));
code_file.close();
}
@@ -1201,7 +992,7 @@ DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
<< "%" << endl
<< "% Warning : this file is generated automatically by Dynare" << endl
<< "% from model file (.mod)" << endl << endl;
-
+
if (containsSteadyStateOperator())
mDynamicModelFile << "global oo_;" << endl << endl;
@@ -1215,7 +1006,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
{
string filename = dynamic_basename + ".c";
ofstream mDynamicModelFile;
-
+
mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
if (!mDynamicModelFile.is_open())
{
@@ -1336,73 +1127,38 @@ DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string
exit(EXIT_FAILURE);
}
u_count_int=0;
- int Size = block_triangular.ModelBlock->Block_List[num].Size - block_triangular.ModelBlock->Block_List[num].Nb_Recursives;
- for(int i=0; i<(int)block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->size();i++)
- {
- //Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
- pair< pair >, pair > it = block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->at(i);
- int k=it.first.first;
- int eq=it.first.second.first;
-
- int var_init=it.first.second.second;
- /*int eqr=it.second.first;
- int varr=it.second.second;*/
- if(eq>=block_triangular.ModelBlock->Block_List[num].Nb_Recursives and var_init>=block_triangular.ModelBlock->Block_List[num].Nb_Recursives)
+ unsigned int block_size = getBlockSize(num);
+ unsigned int block_mfs = getBlockMfs(num);
+ unsigned int block_recursive = block_size - block_mfs;
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[num].begin(); it != (blocks_derivatives[num]).end(); it++)
+ {
+ unsigned int eq = it->first.first;
+ unsigned int var = it->first.second;
+ int lag = it->second.first;
+ if(eq>=block_recursive and var>=block_recursive)
{
- int v=eq-block_triangular.ModelBlock->Block_List[num].Nb_Recursives;
+ int v = eq - block_recursive;
SaveCode.write(reinterpret_cast(&v), sizeof(v));
- int var=it.first.second.second-block_triangular.ModelBlock->Block_List[num].Nb_Recursives + k * Size;
- SaveCode.write(reinterpret_cast(&var), sizeof(var));
- SaveCode.write(reinterpret_cast(&k), sizeof(k));
- int u = u_count_int + Size;
+ int varr = var - block_recursive + lag * block_mfs;
+ SaveCode.write(reinterpret_cast(&varr), sizeof(varr));
+ SaveCode.write(reinterpret_cast(&lag), sizeof(lag));
+ int u = u_count_int + block_mfs;
SaveCode.write(reinterpret_cast(&u), sizeof(u));
- //cout << "eq=" << eq << " var=" << var << " k=" << k << " u=" << u << "\n";
u_count_int++;
}
- }
-
-
- /*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;j(&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++;
- }*/
- u_count_int+=Size;
}
- //cout << "u_count_int=" << u_count_int << "\n";
- for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;jBlock_List[num].Size;j++)
+
+ if (is_two_boundaries)
+ u_count_int+=block_mfs;
+ for (j = block_recursive; j < (int) block_size; j++)
{
- int varr=block_triangular.ModelBlock->Block_List[num].Variable[j];
+ unsigned int varr=getBlockVariableID(num, j);
SaveCode.write(reinterpret_cast(&varr), sizeof(varr));
}
- for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;jBlock_List[num].Size;j++)
+ for (j = block_recursive; j < (int) block_size; j++)
{
- int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j];
- SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1));
+ unsigned int eqr=getBlockEquationID(num, j);
+ SaveCode.write(reinterpret_cast(&eqr), sizeof(eqr));
}
SaveCode.close();
}
@@ -1413,8 +1169,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
string sp;
ofstream mDynamicModelFile;
ostringstream tmp, tmp1, tmp_eq;
- int prev_Simulation_Type, tmp_i;
- //SymbolicGaussElimination SGE;
+ int prev_Simulation_Type;
bool OK;
chdir(basename.c_str());
string filename = dynamic_basename + ".m";
@@ -1431,7 +1186,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << "% from model file (.mod)\n\n";
mDynamicModelFile << "%/\n";
- int i, k, Nb_SGE=0;
+ int Nb_SGE=0;
bool skip_head, open_par=false;
mDynamicModelFile << "function [varargout] = " << dynamic_basename << "(varargin)\n";
@@ -1464,95 +1219,73 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
if (tmp_output.str().length()>0)
mDynamicModelFile << tmp_output.str();
- mDynamicModelFile << " y_kmin=M_.maximum_lag;\n";
- mDynamicModelFile << " y_kmax=M_.maximum_lead;\n";
- mDynamicModelFile << " y_size=M_.endo_nbr;\n";
- mDynamicModelFile << " if(length(varargin)>0)\n";
- mDynamicModelFile << " %it is a simple evaluation of the dynamic model for time _it\n";
- mDynamicModelFile << " params=varargin{3};\n";
- mDynamicModelFile << " it_=varargin{4};\n";
- /*i = symbol_table.endo_nbr*(variable_table.max_endo_lag+variable_table.max_endo_lead+1)+
- symbol_table.exo_nbr*(variable_table.max_exo_lag+variable_table.max_exo_lead+1);
- mDynamicModelFile << " g1=spalloc(" << symbol_table.endo_nbr << ", " << i << ", " << i*symbol_table.endo_nbr << ");\n";*/
- mDynamicModelFile << " Per_u_=0;\n";
- mDynamicModelFile << " Per_y_=it_*y_size;\n";
- mDynamicModelFile << " y=varargin{1};\n";
- mDynamicModelFile << " ys=y(it_,:);\n";
- mDynamicModelFile << " x=varargin{2};\n";
+ mDynamicModelFile << " y_kmin=M_.maximum_lag;" << endl
+ << " y_kmax=M_.maximum_lead;" << endl
+ << " y_size=M_.endo_nbr;" << endl
+ << " if(length(varargin)>0)" << endl
+ << " %it is a simple evaluation of the dynamic model for time _it" << endl
+ << " params=varargin{3};" << endl
+ << " it_=varargin{4};" << endl
+ << " Per_u_=0;" << endl
+ << " Per_y_=it_*y_size;" << endl
+ << " y=varargin{1};" << endl
+ << " ys=y(it_,:);" << endl
+ << " x=varargin{2};" << endl;
prev_Simulation_Type=-1;
tmp.str("");
tmp_eq.str("");
- for (int count_call=1, i = 0;i < block_triangular.ModelBlock->Size;i++, count_call++)
+ unsigned int nb_blocks = getNbBlocks();
+ unsigned int block = 0;
+ for (int count_call=1; block < nb_blocks; block++, count_call++)
{
- mDynamicModelFile << " %block_triangular.ModelBlock->Block_List[i].Nb_Recursives=" << block_triangular.ModelBlock->Block_List[i].Nb_Recursives << " block_triangular.ModelBlock->Block_List[i].Size=" << block_triangular.ModelBlock->Block_List[i].Size << "\n";
- k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD)
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_mfs = getBlockMfs(block);
+ unsigned int block_recursive = block_size - block_mfs;
+ BlockSimulationType simulation_type = getBlockSimulationType(block);
+
+ if(simulation_type==EVALUATE_FORWARD || simulation_type==EVALUATE_BACKWARD)
{
- for (int ik=0 ;ikBlock_List[i].Size;ik++)
+ for (unsigned int ik=0 ;ikBlock_List[i].Variable[ik]+1;
- tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
+ tmp << " " << getBlockVariableID(block, ik)+1;
+ tmp_eq << " " << getBlockEquationID(block, ik)+1;
}
}
else
{
- for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_List[i].Size;ik++)
+ for (unsigned int ik = block_recursive; ik < block_size; ik++)
{
- tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
- tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
+ tmp << " " << getBlockVariableID(block, ik)+1;
+ tmp_eq << " " << getBlockEquationID(block, ik)+1;
}
}
mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
- switch (k)
+ switch (simulation_type)
{
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD:
- mDynamicModelFile << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n";
+ mDynamicModelFile << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << block + 1 << "(y, x, params, 1, it_-1, 1);\n";
mDynamicModelFile << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
break;
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_SIMPLE:
- //mDynamicModelFile << " y_index_eq = " << block_triangular.ModelBlock->Block_List[i].Equation[0]+1 << ";\n";
- mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n";
+ mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << block + 1 << "(y, x, params, it_, 1);\n";
mDynamicModelFile << " residual(y_index_eq)=r;\n";
break;
case SOLVE_FORWARD_COMPLETE:
case SOLVE_BACKWARD_COMPLETE:
- //mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
- mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_, 1);\n";
+ mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << block + 1 << "(y, x, params, it_, 1);\n";
mDynamicModelFile << " residual(y_index_eq)=r;\n";
break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
- int j;
- /*mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
- tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
- mDynamicModelFile << " y_index = [";
- for (j=0;jBlock_List[i].Size;ik++)
- {
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1+j*symbol_table.endo_nbr();
- }
- int tmp_ix=block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1;
- for (j=0;jBlock_List[i].nb_exo;ik++)
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Exogenous[ik]+1+j*symbol_table.exo_nbr()+symbol_table.endo_nbr()*tmp_i;
- mDynamicModelFile << " ];\n";*/
- //mDynamicModelFile << " ga = [];\n";
- j = block_triangular.ModelBlock->Block_List[i].Size*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1)
- + block_triangular.ModelBlock->Block_List[i].nb_exo*(block_triangular.ModelBlock->Block_List[i].Max_Lag_Exo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Exo+1);
- /*mDynamicModelFile << " ga=spalloc(" << block_triangular.ModelBlock->Block_List[i].Size << ", " << j << ", " <<
- block_triangular.ModelBlock->Block_List[i].Size*j << ");\n";*/
- tmp_i=block_triangular.ModelBlock->Block_List[i].Max_Lag_Endo+block_triangular.ModelBlock->Block_List[i].Max_Lead_Endo+1;
- mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << max_lag << ", " << block_triangular.ModelBlock->Block_List[i].Size-block_triangular.ModelBlock->Block_List[i].Nb_Recursives << ");\n";
- /*if(block_triangular.ModelBlock->Block_List[i].Max_Lag==variable_table.max_lag && block_triangular.ModelBlock->Block_List[i].Max_Lead==variable_table.max_lead)
- mDynamicModelFile << " g1(y_index_eq,y_index) = ga;\n";
- else
- mDynamicModelFile << " g1(y_index_eq,y_index) = ga(:," << 1+(variable_table.max_lag-block_triangular.ModelBlock->Block_List[i].Max_Lag)*block_triangular.ModelBlock->Block_List[i].Size << ":" << (variable_table.max_lag+1+block_triangular.ModelBlock->Block_List[i].Max_Lead)*block_triangular.ModelBlock->Block_List[i].Size << ");\n";*/
+ mDynamicModelFile << " [r, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, b, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << block + 1 << "(y, x, params, it_-" << max_lag << ", 1, " << max_lag << ", " << block_recursive << ");\n";
mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
break;
+ default:
+ break;
}
tmp_eq.str("");
tmp.str("");
@@ -1562,43 +1295,47 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << tmp1.str();
tmp1.str("");
}
- mDynamicModelFile << " varargout{1}=residual;\n";
- mDynamicModelFile << " varargout{2}=dr;\n";
- mDynamicModelFile << " return;\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " %it is the deterministic simulation of the block decomposed dynamic model\n";
- mDynamicModelFile << " if(options_.stack_solve_algo==1)\n";
- mDynamicModelFile << " mthd='Sparse LU';\n";
- mDynamicModelFile << " elseif(options_.stack_solve_algo==2)\n";
- mDynamicModelFile << " mthd='GMRES';\n";
- mDynamicModelFile << " elseif(options_.stack_solve_algo==3)\n";
- mDynamicModelFile << " mthd='BICGSTAB';\n";
- mDynamicModelFile << " elseif(options_.stack_solve_algo==4)\n";
- mDynamicModelFile << " mthd='OPTIMPATH';\n";
- mDynamicModelFile << " else\n";
- mDynamicModelFile << " mthd='UNKNOWN';\n";
- mDynamicModelFile << " end;\n";
- mDynamicModelFile << " disp (['-----------------------------------------------------']) ;\n";
- mDynamicModelFile << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;\n";
- mDynamicModelFile << " fprintf('\\n') ;\n";
- mDynamicModelFile << " periods=options_.periods;\n";
- mDynamicModelFile << " maxit_=options_.maxit_;\n";
- mDynamicModelFile << " solve_tolf=options_.solve_tolf;\n";
- mDynamicModelFile << " y=oo_.endo_simul';\n";
- mDynamicModelFile << " x=oo_.exo_simul;\n";
+ mDynamicModelFile << " varargout{1}=residual;" << endl
+ << " varargout{2}=dr;" << endl
+ << " return;" << endl
+ << " end;" << endl
+ << " %it is the deterministic simulation of the block decomposed dynamic model" << endl
+ << " if(options_.stack_solve_algo==1)" << endl
+ << " mthd='Sparse LU';" << endl
+ << " elseif(options_.stack_solve_algo==2)" << endl
+ << " mthd='GMRES';" << endl
+ << " elseif(options_.stack_solve_algo==3)" << endl
+ << " mthd='BICGSTAB';" << endl
+ << " elseif(options_.stack_solve_algo==4)" << endl
+ << " mthd='OPTIMPATH';" << endl
+ << " else" << endl
+ << " mthd='UNKNOWN';" << endl
+ << " end;" << endl
+ << " disp (['-----------------------------------------------------']) ;" << endl
+ << " disp (['MODEL SIMULATION: (method=' mthd ')']) ;" << endl
+ << " fprintf('\\n') ;" << endl
+ << " periods=options_.periods;" << endl
+ << " maxit_=options_.maxit_;" << endl
+ << " solve_tolf=options_.solve_tolf;" << endl
+ << " y=oo_.endo_simul';" << endl
+ << " x=oo_.exo_simul;" << endl;
prev_Simulation_Type=-1;
mDynamicModelFile << " params=M_.params;\n";
mDynamicModelFile << " oo_.deterministic_simulation.status = 0;\n";
- for (i = 0;i < block_triangular.ModelBlock->Size;i++)
+ for (block = 0;block < nb_blocks; block++)
{
- k = block_triangular.ModelBlock->Block_List[i].Simulation_Type;
- if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
- (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/))
+ unsigned int block_size = getBlockSize(block);
+ unsigned int block_mfs = getBlockMfs(block);
+ unsigned int block_recursive = block_size - block_mfs;
+ BlockSimulationType simulation_type = getBlockSimulationType(block);
+
+ if (BlockSim(prev_Simulation_Type)==BlockSim(simulation_type) &&
+ (simulation_type==EVALUATE_FORWARD || simulation_type==EVALUATE_BACKWARD ))
skip_head=true;
else
skip_head=false;
- if ((k == EVALUATE_FORWARD /*|| k == EVALUATE_FORWARD_R*/) && (block_triangular.ModelBlock->Block_List[i].Size))
+ if ((simulation_type == EVALUATE_FORWARD ) && (block_size))
{
if (!skip_head)
{
@@ -1618,17 +1355,15 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n";
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n";
mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
- //mDynamicModelFile << " for it_ = y_kmin+1:(periods+y_kmin)\n";
- mDynamicModelFile << " y=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n";
- mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
+ mDynamicModelFile << " y=" << dynamic_basename << "_" << block + 1 << "(y, x, params, 0, y_kmin, periods);\n";
+ mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
- mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << i <<"']);\n";
+ mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << block <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
- //open_par=true;
}
- else if ((k == EVALUATE_BACKWARD /*|| k == EVALUATE_BACKWARD_R*/) && (block_triangular.ModelBlock->Block_List[i].Size))
+ else if ((simulation_type == EVALUATE_BACKWARD ) && (block_size))
{
if (!skip_head)
{
@@ -1648,15 +1383,15 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).error = 0;\n";
mDynamicModelFile << " oo_.deterministic_simulation.block(blck_num).iterations = 0;\n";
mDynamicModelFile << " g1=[];g2=[];g3=[];\n";
- mDynamicModelFile << " " << dynamic_basename << "_" << i + 1 << "(y, x, params, 0, y_kmin, periods);\n";
- mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
+ mDynamicModelFile << " " << dynamic_basename << "_" << block + 1 << "(y, x, params, 0, y_kmin, periods);\n";
+ mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
- mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << i <<"']);\n";
+ mDynamicModelFile << " disp(['Inf or Nan value during the evaluation of block " << block <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
}
- else if ((k == SOLVE_FORWARD_COMPLETE || k == SOLVE_FORWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ else if ((simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_FORWARD_SIMPLE) && (block_size))
{
if (open_par)
mDynamicModelFile << " end\n";
@@ -1664,30 +1399,28 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " g1=0;\n";
mDynamicModelFile << " r=0;\n";
tmp.str("");
- for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_List[i].Size;ik++)
+ for (unsigned int ik = block_recursive; ik < block_size; ik++)
{
- tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp << " " << getBlockVariableID(block, ik)+1;
}
mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
- int nze, m;
- for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
- nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
+ int nze = blocks_derivatives[block].size();
mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n";
mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n";
mDynamicModelFile << " else\n";
mDynamicModelFile << " blck_num = 1;\n";
mDynamicModelFile << " end;\n";
- mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" <<
+ mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << block + 1 << "'" <<
", y, x, params, y_index, " << nze <<
- ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
+ ", options_.periods, " << blocks_linear[block] <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);\n";
- mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
+ mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
- mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n";
+ mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << block <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
- else if ((k == SOLVE_BACKWARD_COMPLETE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ else if ((simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_SIMPLE) && (block_size))
{
if (open_par)
mDynamicModelFile << " end\n";
@@ -1695,42 +1428,39 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " g1=0;\n";
mDynamicModelFile << " r=0;\n";
tmp.str("");
- for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_List[i].Size;ik++)
+ for (unsigned int ik = block_recursive; ik < block_size; ik++)
{
- tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp << " " << getBlockVariableID(block, ik)+1;
}
mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
- int nze, m;
- for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
- nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
+ int nze = blocks_derivatives[block].size();
+
mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n";
mDynamicModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n";
mDynamicModelFile << " else\n";
mDynamicModelFile << " blck_num = 1;\n";
mDynamicModelFile << " end;\n";
- mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" <<
+ mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << block + 1 << "'" <<
", y, x, params, y_index, " << nze <<
- ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear <<
+ ", options_.periods, " << blocks_linear[block] <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);\n";
- mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
+ mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
- mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n";
+ mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << block <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
- else if ((k == SOLVE_TWO_BOUNDARIES_COMPLETE || k == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size))
+ else if ((simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE) && (block_size))
{
if (open_par)
mDynamicModelFile << " end\n";
open_par=false;
Nb_SGE++;
- int nze, m;
- for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
- nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
+ int nze = blocks_derivatives[block].size();
mDynamicModelFile << " y_index=[";
- for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_List[i].Size;ik++)
+ for (unsigned int ik = block_recursive; ik < block_size; ik++)
{
- mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ mDynamicModelFile << " " << getBlockVariableID(block, ik)+1;
}
mDynamicModelFile << " ];\n";
mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n";
@@ -1738,19 +1468,19 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " else\n";
mDynamicModelFile << " blck_num = 1;\n";
mDynamicModelFile << " end;\n";
- mDynamicModelFile << " y = solve_two_boundaries('" << dynamic_basename << "_" << i + 1 << "'" <<
+ mDynamicModelFile << " y = solve_two_boundaries('" << dynamic_basename << "_" << block + 1 << "'" <<
", y, x, params, y_index, " << nze <<
- ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].Max_Lag <<
- ", " << block_triangular.ModelBlock->Block_List[i].Max_Lead <<
- ", " << block_triangular.ModelBlock->Block_List[i].is_linear <<
+ ", options_.periods, " << max_leadlag_block[block].first <<
+ ", " << max_leadlag_block[block].second <<
+ ", " << blocks_linear[block] <<
", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo);\n";
- mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << i + 1 << ").variable);\n";
+ mDynamicModelFile << " tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);\n";
mDynamicModelFile << " if(isnan(tmp) | isinf(tmp))\n";
- mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << i <<"']);\n";
+ mDynamicModelFile << " disp(['Inf or Nan value during the resolution of block " << block <<"']);\n";
mDynamicModelFile << " return;\n";
mDynamicModelFile << " end;\n";
}
- prev_Simulation_Type=k;
+ prev_Simulation_Type=simulation_type;
}
if (open_par)
mDynamicModelFile << " end;\n";
@@ -1760,7 +1490,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile.close();
- writeModelEquationsOrdered_M( block_triangular.ModelBlock, dynamic_basename);
+ writeModelEquationsOrdered_M(dynamic_basename);
chdir("..");
}
@@ -1993,7 +1723,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
}
void
-DynamicModel::writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll) const
+DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll) const
{
/* Writing initialisation for M_.lead_lag_incidence matrix
M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -2035,58 +1765,50 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block, b
<< equation_tags[i].second.second << "' ;" << endl;
output << "};" << endl;
- //In case of sparse model, writes the block structure of the model
- if (block)
+ //In case of sparse model, writes the block_decomposition structure of the model
+ if (block_decomposition)
{
- //int prev_Simulation_Type=-1;
- //bool skip_the_head;
- int k=0;
int count_lead_lag_incidence = 0;
int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo;
- for (int j = 0;j < block_triangular.ModelBlock->Size;j++)
+ unsigned int nb_blocks = getNbBlocks();
+ for (unsigned int block = 0; block < nb_blocks; block++)
{
//For a block composed of a single equation determines wether we have to evaluate or to solve the equation
- //skip_the_head=false;
- k++;
count_lead_lag_incidence = 0;
- int Block_size=block_triangular.ModelBlock->Block_List[j].Size;
- max_lag =block_triangular.ModelBlock->Block_List[j].Max_Lag ;
- max_lead=block_triangular.ModelBlock->Block_List[j].Max_Lead;
- max_lag_endo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Endo ;
- max_lead_endo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Endo;
- max_lag_exo =block_triangular.ModelBlock->Block_List[j].Max_Lag_Exo ;
- max_lead_exo=block_triangular.ModelBlock->Block_List[j].Max_Lead_Exo;
- bool evaluate=false;
- vector exogenous;
+ BlockSimulationType simulation_type = getBlockSimulationType(block);
+ int block_size = getBlockSize(block);
+ max_lag = max_leadlag_block[block].first;
+ max_lead = max_leadlag_block[block].second;
+ max_lag_endo = endo_max_leadlag_block[block].first;
+ max_lead_endo= endo_max_leadlag_block[block].second;
+ max_lag_exo = max(exo_max_leadlag_block[block].first, exo_det_max_leadlag_block[block].first);
+ max_lead_exo= max(exo_max_leadlag_block[block].second, exo_det_max_leadlag_block[block].second);
+ vector exogenous(symbol_table.exo_nbr(), -1);
vector::iterator it_exogenous;
exogenous.clear();
ostringstream tmp_s, tmp_s_eq;
tmp_s.str("");
tmp_s_eq.str("");
- for (int i=0;iBlock_List[j].Size;i++)
+ for (int i=0;iBlock_List[j].Variable[i]+1;
- tmp_s_eq << " " << block_triangular.ModelBlock->Block_List[j].Equation[i]+1;
+ tmp_s << " " << getBlockVariableID(block, i)+1;
+ tmp_s_eq << " " << getBlockEquationID(block, i)+1;
}
- for (int i=0;iBlock_List[j].nb_exo;i++)
- {
- int ii=block_triangular.ModelBlock->Block_List[j].Exogenous[i];
- for (it_exogenous=exogenous.begin();it_exogenous!=exogenous.end() && *it_exogenous!=ii;it_exogenous++) /*cout << "*it_exogenous=" << *it_exogenous << "\n"*/;
- if (it_exogenous==exogenous.end() || exogenous.begin()==exogenous.end())
- exogenous.push_back(ii);
- }
- output << "M_.block_structure.block(" << k << ").num = " << j+1 << ";\n";
- output << "M_.block_structure.block(" << k << ").Simulation_Type = " << block_triangular.ModelBlock->Block_List[j].Simulation_Type << ";\n";
- output << "M_.block_structure.block(" << k << ").maximum_lag = " << max_lag << ";\n";
- output << "M_.block_structure.block(" << k << ").maximum_lead = " << max_lead << ";\n";
- output << "M_.block_structure.block(" << k << ").maximum_endo_lag = " << max_lag_endo << ";\n";
- output << "M_.block_structure.block(" << k << ").maximum_endo_lead = " << max_lead_endo << ";\n";
- output << "M_.block_structure.block(" << k << ").maximum_exo_lag = " << max_lag_exo << ";\n";
- output << "M_.block_structure.block(" << k << ").maximum_exo_lead = " << max_lead_exo << ";\n";
- output << "M_.block_structure.block(" << k << ").endo_nbr = " << Block_size << ";\n";
- output << "M_.block_structure.block(" << k << ").equation = [" << tmp_s_eq.str() << "];\n";
- output << "M_.block_structure.block(" << k << ").variable = [" << tmp_s.str() << "];\n";
- output << "M_.block_structure.block(" << k << ").exogenous = [";
+ it_exogenous = exogenous.begin();
+ for(t_lag_var::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
+ it_exogenous = set_union(it->second.begin(), it->second.end(), exogenous.begin(), exogenous.end(), it_exogenous);
+ output << "M_.block_structure.block(" << block+1 << ").num = " << block+1 << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_endo_lag = " << max_lag_endo << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
+ output << "M_.block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];\n";
+ output << "M_.block_structure.block(" << block+1 << ").exogenous = [";
int i=0;
for (it_exogenous=exogenous.begin();it_exogenous!=exogenous.end();it_exogenous++)
if (*it_exogenous>=0)
@@ -2095,112 +1817,79 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block, b
i++;
}
output << "];\n";
- output << "M_.block_structure.block(" << k << ").exo_nbr = " << i << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").exo_nbr = " << i << ";\n";
- output << "M_.block_structure.block(" << k << ").exo_det_nbr = " << block_triangular.ModelBlock->Block_List[j].nb_exo_det << ";\n";
+ output << "M_.block_structure.block(" << block+1 << ").exo_det_nbr = " << exo_det_block.size() << ";\n";
tmp_s.str("");
-
- bool done_IM=false;
- if (!evaluate)
+ count_lead_lag_incidence = 0;
+ dynamic_jacob_map reordered_dynamic_jacobian;
+ for (t_block_derivatives_equation_variable_laglead_nodeid::const_iterator it = blocks_derivatives[block].begin(); it != blocks_derivatives[block].end(); it++)
+ reordered_dynamic_jacobian[make_pair(it->second.first, make_pair(it->first.second, it->first.first))] = it->second.second;
+ output << "M_.block_structure.block(" << block+1 << ").lead_lag_incidence = [];\n";
+ int last_var = -1;
+ for (int lag=-max_lag_endo;lagfirst.first && last_var != it->first.second.first)
{
- for (int l_var=0;l_varBlock_List[j].Size;l_var++)
- {
- for (int l_equ=0;l_equBlock_List[j].Size;l_equ++)
- if (tmp_IM[block_triangular.ModelBlock->Block_List[j].Equation[l_equ]*symbol_table.endo_nbr()+block_triangular.ModelBlock->Block_List[j].Variable[l_var]])
- {
- count_lead_lag_incidence++;
- if (tmp_s.str().length())
- tmp_s << " ";
- tmp_s << count_lead_lag_incidence;
- done_IM=true;
- break;
- }
- if (!done_IM)
- tmp_s << " 0";
- done_IM=false;
- }
- output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [ M_.block_structure.block(" << k << ").lead_lag_incidence; " << tmp_s.str() << "];\n";
- tmp_s.str("");
+ count_lead_lag_incidence++;
+ for( int i = last_var; i < it->first.second.first-1; i++)
+ tmp_s << " 0";
+ if (tmp_s.str().length())
+ tmp_s << " ";
+ tmp_s << count_lead_lag_incidence;
+ last_var = it->first.second.first;
}
}
- }
- else
- {
- bool done_some_where;
- output << "M_.block_structure.block(" << k << ").lead_lag_incidence = [\n";
- for (int l=-max_lag_endo;lBlock_List[ii].Size;l_var++)
- {
- for (int l_equ=0;l_equBlock_List[ii].Size;l_equ++)
- if (tmp_IM[block_triangular.ModelBlock->Block_List[ii].Equation[l_equ]*symbol_table.endo_nbr()+block_triangular.ModelBlock->Block_List[ii].Variable[l_var]])
- {
- //if(not_increm && l==-max_lag)
- count_lead_lag_incidence++;
- not_increm=false;
- if (tmp_s.str().length())
- tmp_s << " ";
- //tmp_s << count_lead_lag_incidence+(l+max_lag)*Block_size;
- tmp_s << count_lead_lag_incidence;
- done_IM=true;
- break;
- }
- if (!done_IM)
- tmp_s << " 0";
- else
- done_some_where = true;
- done_IM=false;
- }
- ii++;
- }
- output << tmp_s.str() << "\n";
- tmp_s.str("");
- }
- }
- output << "];\n";
- }
-
- }
- for (int j=-block_triangular.incidencematrix.Model_Max_Lag_Endo;j<=block_triangular.incidencematrix.Model_Max_Lead_Endo;j++)
- {
- bool* IM = block_triangular.incidencematrix.Get_IM(j, eEndogenous);
- if (IM)
- {
- bool new_entry=true;
- output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").lead_lag = " << j << ";\n";
- output << "M_.block_structure.incidence(" << block_triangular.incidencematrix.Model_Max_Lag_Endo+j+1 << ").sparse_IM = [";
- for (int i=0;i >, int> lag_row_incidence;
+ for (first_derivatives_type::const_iterator it = first_derivatives.begin();
+ it != first_derivatives.end(); it++)
+ {
+ int deriv_id = it->first.second;
+ if (getTypeByDerivID(deriv_id) == eEndogenous)
+ {
+ int eq = it->first.first;
+ int symb = getSymbIDByDerivID(deriv_id);
+ int var = symbol_table.getTypeSpecificID(symb);
+ int lag = getLagByDerivID(deriv_id);
+ int eqr = inv_equation_reordered[eq];
+ int varr = inv_variable_reordered[var];
+ lag_row_incidence[make_pair(lag, make_pair(eqr, varr))] = 1;
+ }
+ }
+ int prev_lag=-1000000;
+ for(map >, int> ::const_iterator it=lag_row_incidence.begin(); it != lag_row_incidence.end(); it++)
+ {
+ if(prev_lag != it->first.first)
+ {
+ if(prev_lag != -1000000)
+ output << "];\n";
+ prev_lag = it->first.first;
+ output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").lead_lag = " << prev_lag << ";\n";
+ output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").sparse_IM = [";
+ }
+ output << it->first.second.first << " " << it->first.second.second << ";\n";
+ }
+ output << "];\n";
}
// Writing initialization for some other variables
output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];" << endl
@@ -2234,139 +1923,6 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block, b
<< "M_.NNZDerivatives(3) = " << NNZDerivatives[2] << ";" << endl;
}
-void
-DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic)
-{
- int i=0;
- int j=0;
- bool *IM=NULL;
- int a_variable_lag=-9999;
- for (first_derivatives_type::iterator it = first_derivatives.begin();
- it != first_derivatives.end(); it++)
- {
- //cout << "it->first.second=" << it->first.second << " variable_table.getSymbolID(it->first.second)=" << variable_table.getSymbolID(it->first.second) << " Type=" << variable_table.getType(it->first.second) << " eEndogenous=" << eEndogenous << " eExogenous=" << eExogenous << " variable_table.getLag(it->first.second)=" << variable_table.getLag(it->first.second) << "\n";
- if (getTypeByDerivID(it->first.second) == eEndogenous)
- {
- NodeID Id = it->second;
- double val = 0;
- try
- {
- val = Id->eval(eval_context);
- }
- catch (ExprNode::EvalException &e)
- {
- cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ") [" << getSymbIDByDerivID(it->first.second) << "] !" << endl;
- Id->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
- cout << "\n";
- cerr << "DynamicModel::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getName(getSymbIDByDerivID(it->first.second)) << "(" << getLagByDerivID(it->first.second) << ")!" << endl;
- }
- int eq=it->first.first;
- int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(it->first.second));///symbol_table.getID(eEndogenous,it->first.second);//variable_table.getSymbolID(it->first.second);
- int k1 = getLagByDerivID(it->first.second);
- if (a_variable_lag!=k1)
- {
- IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous);
- a_variable_lag=k1;
- }
- if (k1==0 or !dynamic)
- {
- j++;
- (*j_m)[make_pair(eq,var)]+=val;
- }
- if (IM[eq*symbol_table.endo_nbr()+var] && (fabs(val) < cutoff))
- {
- if (block_triangular.bt_verbose)
- cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")\n";
- block_triangular.incidencematrix.unfill_IM(eq, var, k1, eEndogenous);
- i++;
- }
- }
- }
- //Get ride of the elements of the incidence matrix equal to Zero
- IM=block_triangular.incidencematrix.Get_IM(0, eEndogenous);
- for (int i=0;i0)
- {
- cout << i << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded\n";
- cout << "the contemporaneous incidence matrix has " << j << " elements\n";
- }
-}
-
-void
-DynamicModel::BlockLinear(Model_Block *ModelBlock)
-{
- int i,j,l,m,ll;
- for (j = 0;j < ModelBlock->Size;j++)
- {
- if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE ||
- ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
- {
- ll=ModelBlock->Block_List[j].Max_Lag;
- for (i=0;iBlock_List[j].IM_lead_lag[ll].size;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[ll].Equ_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[ll].Var_Index[i];
- //first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,variable_table.getID(var,0)));
- first_derivatives_type::const_iterator it=first_derivatives.find(make_pair(eq,getDerivID(symbol_table.getID(eEndogenous, var),0)));
- if (it!= first_derivatives.end())
- {
- NodeID Id = it->second;
- set > endogenous;
- Id->collectEndogenous(endogenous);
- if (endogenous.size() > 0)
- {
- for (l=0;lBlock_List[j].Size;l++)
- {
- if (endogenous.find(make_pair(ModelBlock->Block_List[j].Variable[l], 0)) != endogenous.end())
- {
- ModelBlock->Block_List[j].is_linear=false;
- goto follow;
- }
- }
- }
- }
- }
- }
- else if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
- {
- for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- int k1=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;i