diff --git a/BlockTriangular.cc b/BlockTriangular.cc
index 3260e546959ed08a60f0f440d1efbe68d1c18d0a..14c8063e275682d2a7cac176ade4a6aae20421b0 100644
--- a/BlockTriangular.cc
+++ b/BlockTriangular.cc
@@ -20,7 +20,6 @@
#include <iostream>
#include <sstream>
#include <fstream>
-#include <ctime>
#include <cstdlib>
#include <cstring>
#include <cmath>
@@ -46,6 +45,7 @@ BlockTriangular::BlockTriangular(SymbolTable &symbol_table_arg, NumericalConstan
bt_verbose = 0;
ModelBlock = NULL;
periods = 0;
+ prologue = epilogue = 0;
Normalized_Equation = new DataTree(symbol_table, num_const);
}
@@ -114,7 +114,7 @@ BlockTriangular::Prologue_Epilogue(bool *IM, int &prologue, int &epilogue, int n
//------------------------------------------------------------------------------
// Find a matching between equations and endogenous variables
bool
-BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector<int> &Index_Equ_IM) const
+BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, int verbose, bool *IM0, vector<int> &Index_Equ_IM) const
{
int n = equation_number - prologue - epilogue;
@@ -143,9 +143,17 @@ BlockTriangular::Compute_Normalization(bool *IM, int equation_number, int prolog
mate_map_t::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
if (it != mate_map.begin() + n)
{
- if (verbose)
+ if (verbose == 1)
{
- cerr << "ERROR: Could not normalize dynamic model. Variable "
+ cerr << "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);
@@ -183,7 +191,7 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
variable_blck[Index_Var_IM[i]] = i;
equation_blck[Index_Equ_IM[i]] = i;
}
- else if (i < components_set.size() + prologue)
+ 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;
@@ -226,7 +234,7 @@ BlockTriangular::Get_Variable_LeadLag_By_Block(vector<int > &components_set, int
}
void
-BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_etype &Equation_Type, bool verbose_) const
+BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_etype &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs) const
{
t_vtype V_Variable_Type;
int n = nb_var - prologue - epilogue;
@@ -272,11 +280,16 @@ BlockTriangular::Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Blo
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
- 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)
- add_edge(i, i, G2);
-
+ 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
@@ -686,7 +699,7 @@ BlockTriangular::Free_Block(Model_Block *ModelBlock) const
}
t_etype
-BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM)
+BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, int mfs)
{
NodeID lhs, rhs;
ostringstream tmp_output;
@@ -709,44 +722,40 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations,
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
tmp_s << "y(it_, " << Index_Var_IM[i]+1 << ")";
map<pair<int, pair<int, int> >, NodeID>::iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
- /*cout << "eq=" << eq << " var=" << var << "=";
- first_cur_endo_derivatives[make_pair(eq, var)]->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
- cout << "\n";
- cout << "equation : ";
- eq_node->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
- cout << "\n";*/
set<pair<int, int> > result;
pair<bool, NodeID> res;
derivative->second->collectEndogenous(result);
- /*for(set<pair<int, int> >::const_iterator iitt = result.begin(); iitt!=result.end(); iitt++)
- cout << "result = " << iitt->first << ", " << iitt->second << "\n";*/
set<pair<int, int> >::const_iterator d_endo_variable = result.find(make_pair(var, 0));
//Determine whether the equation could be evaluated rather than to be solved
- if (tmp_output.str() == tmp_s.str() and d_endo_variable == result.end())
+ ostringstream tt("");
+ derivative->second->writeOutput(tt, oMatlabDynamicModelSparse, temporary_terms);
+ //cout << tt.str().c_str() << " tmp_output.str()=" << tmp_output.str() << " tmp_s.str()=" << tmp_s.str();
+ if (tmp_output.str() == tmp_s.str() and tt.str()=="1")
{
Equation_Simulation_Type = E_EVALUATE;
+ //cout << " E_EVALUATE ";
}
else
{
- //vector<pair<int, NodeID> > List_of_Op_RHS;
- //the equation could be normalized by a permutation of the rhs and the lhs
- if (d_endo_variable == result.end()) //the equation is linear in the endogenous and could be normalized using the derivative
+ vector<pair<int, pair<NodeID, NodeID> > > List_of_Op_RHS;
+ res = equations[eq]->normalizeEquation(var, List_of_Op_RHS);
+ if(mfs==2)
{
- Equation_Simulation_Type = E_EVALUATE_S;
- //cout << " gone normalized : ";
- res = equations[eq]->normalizeLinearInEndoEquation(var, derivative->second/*, List_of_Op_RHS*/);
- /*res.second->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
- cout << " done\n";*/
+ 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;
}
}
+ //cout << " " << c_Equation_Type(Equation_Simulation_Type) << endl;
V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second));
}
return (V_Equation_Simulation_Type);
}
-/*void
-BlockTriangular::Recompute_Deriavtives_Respect_to_Feedback_Variables(
-*/
t_type
BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue, vector<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM)
{
@@ -859,15 +868,15 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>
BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
{
- map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> Derivatives;
- Derivatives.clear();
- int nb_endo = symbol_table.endo_nbr();
- /*ModelBlock.Block_List[Blck].first_order_determinstic_simulation_derivatives = new*/
+ map<pair<pair<int, pair<int, int> >, pair<int, int> >, 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)
- {
+ 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];
@@ -886,26 +895,26 @@ BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
OK=false;
}
- if(OK)
- {
- if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
- //It's a normalized equation, we have to recompute the derivative using chain rule derivative function*/
+ if(OK)
+ {
+ if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_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(var<ModelBlock->Block_List[blck].Nb_Recursives)
- {
- int eqs = ModelBlock->Block_List[blck].Equation[var];
- for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; vars<ModelBlock->Block_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;
- }
- }
+ 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(var<ModelBlock->Block_List[blck].Nb_Recursives)
+ {
+ int eqs = ModelBlock->Block_List[blck].Equation[var];
+ for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; vars<ModelBlock->Block_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;
+ }
+ }
}
}
}
@@ -915,18 +924,18 @@ BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
}
void
-BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool *IM_0, jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives)
+BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool *IM_0, jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, 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;
-
- SIM0 = (bool *) malloc(n * n * sizeof(bool));
+ SIM0 = (bool *) malloc(n * n * sizeof(bool));
memcpy(SIM0, IM_0, n*n*sizeof(bool));
Prologue_Epilogue(IM, prologue, epilogue, n, Index_Var_IM, Index_Equ_IM, SIM0);
free(SIM0);
+
int counted = 0;
if (prologue+epilogue < n)
{
@@ -956,7 +965,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
memset(SIM00, 0, n*n*sizeof(bool));
for (map< pair< int, int >, double >::iterator iter = j_m.begin(); iter != j_m.end(); iter++)
{
- if (fabs(iter->second) > bi)
+ if (fabs(iter->second) > max(bi, cutoff))
{
SIM0[iter->first.first*n+iter->first.second] = 1;
if (!IM_0[iter->first.first*n+iter->first.second])
@@ -974,7 +983,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
}
free(SIM0);
if (suppress != suppressed)
- OK = Compute_Normalization(IM, n, prologue, epilogue, false, SIM00, Index_Equ_IM);
+ OK = Compute_Normalization(IM, n, prologue, epilogue, 0, SIM00, Index_Equ_IM);
suppressed = suppress;
if (!OK)
//bi/=1.07;
@@ -984,15 +993,23 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
free(SIM00);
}
if (!OK)
- Compute_Normalization(IM, n, prologue, epilogue, true, SIM00, Index_Equ_IM);
+ {
+ Compute_Normalization(IM, n, prologue, epilogue, 1, SIM00, Index_Equ_IM);
+ Compute_Normalization(IM, n, prologue, epilogue, 2, IM_0, Index_Equ_IM);
+ }
}
- V_Equation_Type = Equation_Type_determination(equations, first_order_endo_derivatives, Index_Var_IM, Index_Equ_IM);
+ 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<pair<int, int> > blocks;
if (prologue+epilogue < n)
- Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(IM, n, prologue, epilogue, Index_Equ_IM, Index_Var_IM, blocks, V_Equation_Type, false);
+ {
+ 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);
@@ -1049,7 +1066,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool *IM, Model_Block *ModelBlock,
// normalize each equation of the dynamic model
// and find the optimal block triangular decomposition of the static model
void
-BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &equation_simulation_type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives)
+BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &equation_simulation_type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, int mfs, double cutoff)
{
bool *SIM, *SIM_0;
bool *Cur_IM;
@@ -1091,7 +1108,7 @@ BlockTriangular::Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vec
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);
+ 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
index 3d6e6dabf1b2165de81c577d46b20d1e7384d8af..c0ebfb4bb2cb96126c81ba6e7b9513a5ea633e31 100644
--- a/BlockTriangular.hh
+++ b/BlockTriangular.hh
@@ -97,11 +97,11 @@ private:
//! 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<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Finds a matching between equations and endogenous variables
- bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, bool verbose, bool *IM0, vector<int> &Index_Var_IM) const;
+ bool Compute_Normalization(bool *IM, int equation_number, int prologue, int epilogue, int verbose, bool *IM0, vector<int> &Index_Var_IM) const;
//! Decomposes into recurive blocks the non purely recursive equations and determines for each block the minimum feedback variables
- void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_etype &Equation_Type, bool verbose_) const;
+ void Compute_Block_Decomposition_and_Feedback_Variables_For_Each_Block(bool *IM, int nb_var, int prologue, int epilogue, vector<int> &Index_Equ_IM, vector<int> &Index_Var_IM, vector<pair<int, int> > &blocks, t_etype &Equation_Type, bool verbose_, 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<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
+ t_etype Equation_Type_determination(vector<BinaryOpNode *> &equations, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, vector<int> &Index_Var_IM, vector<int> &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<pair<int, int> > &blocks, vector<BinaryOpNode *> &equations, t_etype &Equation_Type, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM);
//! Compute for each variable its maximum lead and lag in its block
@@ -121,8 +121,8 @@ public:
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
- void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives);
- void Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, int n, int &prologue, int &epilogue, vector<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM_0 , jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &equation_simulation_type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives);
+ void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, 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<int> &Index_Var_IM, vector<int> &Index_Equ_IM, bool* IM_0 , jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &equation_simulation_type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, bool dynamic, int mfs, double cutoff);
vector<int> Index_Equ_IM;
vector<int> Index_Var_IM;
int prologue, epilogue;
@@ -187,11 +187,10 @@ public:
};
inline static std::string c_Equation_Type(int type)
{
- char c_Equation_Type[5][13]=
+ char c_Equation_Type[4][13]=
{
"E_UNKNOWN ",
"E_EVALUATE ",
- //"E_EVALUATE_R",
"E_EVALUATE_S",
"E_SOLVE "
};
diff --git a/CodeInterpreter.hh b/CodeInterpreter.hh
index cb495131eff5d239aac0103b20fae99a57160a42..a33301f5486287d8dfa2ec864eb72ea8e710db28 100644
--- a/CodeInterpreter.hh
+++ b/CodeInterpreter.hh
@@ -40,6 +40,15 @@ const char FDIMT=16;
const char FEND=17;
const char FOK=18;
const char FENDEQU=19;
+const char FLDSV=20;
+const char FSTPSV=21;
+const char FLDSU=22;
+const char FSTPSU=23;
+const char FLDST=24;
+const char FSTPST=25;
+const char FDIMST=26;
+
+
enum BlockType
{
@@ -53,7 +62,6 @@ enum EquationType
{
E_UNKNOWN, //!< Unknown equation type
E_EVALUATE, //!< Simple evaluation, normalized variable on left-hand side
- //E_EVALUATE_R, //!< Simple evaluation, normalized variable on right-hand side
E_EVALUATE_S, //!< Simple evaluation, normalize using the first order derivative
E_SOLVE //!< No simple evaluation of the equation, it has to be solved
};
diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index 533f280676028c354a3e8c0e199275841b0e14fe..22ba5950ad6be4c2fac22614154cfce0d3639139 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -27,8 +27,8 @@ using namespace std;
#include "ComputingTasks.hh"
#include "Statement.hh"
-SteadyStatement::SteadyStatement(const OptionsList &options_list_arg) :
- options_list(options_list_arg)
+SteadyStatement::SteadyStatement(const OptionsList &options_list_arg, StaticDllModel::mode_t mode_arg) :
+ options_list(options_list_arg), mode(mode_arg)
{
}
@@ -37,12 +37,17 @@ SteadyStatement::checkPass(ModFileStructure &mod_file_struct)
{
if (options_list.num_options.find("block_mfs") != options_list.num_options.end())
mod_file_struct.steady_block_mfs_option = true;
+ else if (options_list.num_options.find("block_mfs_dll") != options_list.num_options.end())
+ mod_file_struct.steady_block_mfs_dll_option = true;
}
void
SteadyStatement::writeOutput(ostream &output, const string &basename) const
{
options_list.writeOutput(output);
+ /*if (mode == StaticDllModel::eSparseDLLMode)
+ output << "oo_.steady_state=simulate('steady');" << endl;
+ else*/
output << "steady;\n";
}
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index fad2123a949c63be3816d08596a4d63a9da4270d..37f7ebc329ddd71e2ffa33b28f506e6412120ce9 100644
--- a/ComputingTasks.hh
+++ b/ComputingTasks.hh
@@ -32,8 +32,9 @@ class SteadyStatement : public Statement
{
private:
const OptionsList options_list;
+ const StaticDllModel::mode_t mode;
public:
- SteadyStatement(const OptionsList &options_list_arg);
+ SteadyStatement(const OptionsList &options_list_arg, StaticDllModel::mode_t mode_arg);
virtual void checkPass(ModFileStructure &mod_file_struct);
virtual void writeOutput(ostream &output, const string &basename) const;
};
diff --git a/DynamicModel.cc b/DynamicModel.cc
index 27471170ff861368ee237aac7f38dd9da4d49bca..89afb45582f2c5f5c38a5972766ad2309a95bf75 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -23,7 +23,6 @@
#include <cassert>
#include <cstdio>
#include <cerrno>
-
#include "DynamicModel.hh"
// For mkdir() and chdir()
@@ -46,6 +45,7 @@ DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
mode(eStandardMode),
cutoff(1e-15),
markowitz(0.7),
+ mfs(0),
block_triangular(symbol_table_arg, num_constants_arg)
{
}
@@ -62,7 +62,7 @@ DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int la
//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);
+ (it->second)->compile(code_file, false, temporary_terms, map_idx, true);
else
code_file.write(&FLDZ, sizeof(FLDZ));
}
@@ -73,7 +73,7 @@ DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr,
{
map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
if (it != first_chain_rule_derivatives.end())
- (it->second)->compile(code_file, false, temporary_terms, map_idx);
+ (it->second)->compile(code_file, false, temporary_terms, map_idx, true);
else
code_file.write(&FLDZ, sizeof(FLDZ));
}
@@ -112,7 +112,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
{
map<NodeID, pair<int, int> > first_occurence;
map<NodeID, int> reference_count;
- int i, j, m, eq, var, lag;
+ int i, j, m, eq, var, eqr, varr, lag;
temporary_terms_type vect;
ostringstream tmp_output;
BinaryOpNode *eq_node;
@@ -127,12 +127,24 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
// Compute the temporary terms reordered
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, i, map_idx);
- if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
- if(ModelBlock->Block_List[j].Equation_Normalized[i])
+ if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S && i<ModelBlock->Block_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<int, pair<int, int> >, pair<int, int> > 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++)
{
lag=m-ModelBlock->Block_List[j].Max_Lag;
@@ -141,21 +153,18 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
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)));
- //printf("it=%d eq=%d var=%s (%d)\n",it!=first_derivatives.end(), eq, symbol_table.getName(symbol_table.getID(eEndogenous, var)).c_str(), var);
- //if(it!=first_derivatives.end())
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
- for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
+ /*for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
lag=it.first.first;
- eq=it.first.second.first;
- var=it.first.second.second;
- it_chr=first_chain_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
- //it_chr->second->writeChainRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+ 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;
@@ -178,8 +187,6 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
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->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
@@ -190,13 +197,21 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
// Collecte the temporary terms reordered
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
+ if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S && i<ModelBlock->Block_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
+ {
+ 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);
+ (*it)->collectTemporary_terms(temporary_terms, ModelBlock, j);*/
}
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
@@ -211,14 +226,13 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
- for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
+ for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
lag=it.first.first;
- eq=it.first.second.first;
- var=it.first.second.second;
- it_chr=first_chain_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
- //it_chr->second->writeChainRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
+ 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++)
@@ -268,7 +282,7 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
BinaryOpNode *eq_node;
ostringstream Uf[symbol_table.endo_nbr()];
map<NodeID, int> reference_count;
- int prev_Simulation_Type=-1, count_derivates=0;
+ //int prev_Simulation_Type=-1, count_derivates=0;
int jacobian_max_endo_col;
ofstream output;
//temporary_terms_type::const_iterator it_temp=temporary_terms.begin();
@@ -324,10 +338,9 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
<< BlockTriangular::BlockSim(ModelBlock->Block_List[j].Simulation_Type) << " //" << endl
<< " % ////////////////////////////////////////////////////////////////////////" << 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
- /*||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD_R
- ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD_R*/)
+ ||ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
{
output << " if(jacobian_eval)\n";
output << " g1 = spalloc(" << ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives
@@ -437,15 +450,18 @@ evaluation: if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDAR
{
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);
}
- /*else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_R)
- {
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
- output << " = ";
- output << tmp_output.str();
- output << "; %reversed " << ModelBlock->Block_List[j].Equation_Type[i] << " \n";
- }*/
else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
{
output << "%" << tmp_output.str();
@@ -454,19 +470,23 @@ evaluation: if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDAR
{
rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << "\n ";
- /*temporary_terms_type tt2;
- tt2.clear();*/
tmp_output.str("");
eq_node = (BinaryOpNode *)ModelBlock->Block_List[j].Equation_Normalized[i];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- if(ModelBlock->Block_List[j].Simulation_Type==EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type==EVALUATE_FORWARD)
- lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
- else
- lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
- output << tmp_output.str();
+ lhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
output << " = ";
- rhs->writeOutput(output, oMatlabDynamicModelSparse, temporary_terms);
+ /*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);
}
}
else
@@ -535,9 +555,6 @@ end:
{
case EVALUATE_BACKWARD:
case EVALUATE_FORWARD:
- /*case EVALUATE_BACKWARD_R:
- case EVALUATE_FORWARD_R:*/
- count_derivates++;
for (m=0;m<ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -604,7 +621,6 @@ end:
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
- count_derivates++;
for (m=0;m<ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag+1;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
@@ -649,7 +665,7 @@ end:
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 << ", "
+ 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))
@@ -662,130 +678,76 @@ end:
output << " varargout{2}=g1_o;\n";
output << " else" << endl;
- m=ModelBlock->Block_List[j].Max_Lag;
- //cout << "\nDerivatives in Block " << j << "\n";
- for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- //Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
- pair< pair<int, pair<int, int> >, pair<int, int> > 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;
-
- /*for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
+ for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->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];
- ostringstream tmp_output;
- if (eqr<ModelBlock->Block_List[j].Nb_Recursives)
- {
- if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
- {*/
- 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 << ";";
- output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
- << "(" << k
- << ") " << varr+1
- << ", equation=" << eqr+1 << endl;
- }
- /*}
- }
- }*/
+ pair< pair<int, pair<int, int> >, pair<int, int> > 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;
+ }
output << " end;\n";
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
output << " if ~jacobian_eval" << endl;
- /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
+ for(i=0; i<(int)ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
+ pair< pair<int, pair<int, int> >, pair<int, int> > 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;
+ ostringstream tmp_output;
+ if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
{
- 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];*/
- for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- //Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
- pair< pair<int, pair<int, int> >, pair<int, int> > 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;
-
- //bool derivative_exist;
- ostringstream tmp_output;
- /*if (eqr>=ModelBlock->Block_List[j].Nb_Recursives)
- {
- if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
- {*/
- /*for(int equation = ModelBlock->Block_List[j].Nb_Recursives; equation<ModelBlock->Block_List[j].Size; equation++)
- {
- int eq = ModelBlock->Block_List[j].Equation[equation];
- int eqr = equation - ModelBlock->Block_List[j].Nb_Recursives;
- for(int variable = ModelBlock->Block_List[j].Nb_Recursives; variable<ModelBlock->Block_List[j].Size; variable++)
- {
- int var = ModelBlock->Block_List[j].Variable[variable];
- int varr = variable - ModelBlock->Block_List[j].Nb_Recursives;*/
- //cout << "eqr=" << eqr << " varr=" << varr;
- //cout << "k=" << k << " eq=" << eq << " var=" << var << " eqr=" << eqr << " varr=" << varr << " ModelBlock->Block_List[j].Equation[eq]=" << ModelBlock->Block_List[j].Equation[eq] << "\n";
- if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
- {
-
- 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
- << "+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
- << "+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 << ")) = ";
-
-
- output << " " << tmp_output.str();
-
- writeChainRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
-
- output << ";";
- output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
- << "(" << k << ") " << varr+1
- << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
- }
- //cout << " done\n";
- /* }
- }*/
-
+ 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
+ << "+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
+ << "+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 << ")) = ";
+ output << " " << tmp_output.str();
+
+ writeChainRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
+
+ output << ";";
+ output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
+ << "(" << k << ") " << varr+1
+ << ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
+ }
#ifdef CONDITION
- output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
- output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
+ output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
+ output << " condition(" << eqr << ")=u(" << u << "+Per_u_);\n";
#endif
//}
}
@@ -879,7 +841,6 @@ end:
default:
break;
}
- prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
output.close();
}
}
@@ -896,16 +857,14 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
struct Uff
{
Uff_l *Ufl, *Ufl_First;
- int eqr;
};
- int i,j,k,m, v;
+ int i,j,k,v;
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
NodeID lhs=NULL, rhs=NULL;
BinaryOpNode *eq_node;
- bool lhs_rhs_done;
Uff Uf[symbol_table.endo_nbr()];
map<NodeID, int> reference_count;
vector<int> feedback_variables;
@@ -948,7 +907,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
//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";
+ //cout << "u_count_int=" << u_count_int << "\n";
code_file.write(reinterpret_cast<char *>(&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 = u_count_int ;
@@ -964,30 +923,19 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
file_open=true;
}
- /*if (ModelBlock->Block_List[j].Size==1)
- {
- lhs_rhs_done=true;
- eq_node = equations[ModelBlock->Block_List[j].Equation[0]];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- }
- else
- lhs_rhs_done=false;*/
// The equations
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
//The Temporary terms
temporary_terms_type tt2;
tt2.clear();
- /*if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size())
- output << " " << sps << "% //Temporary variables" << endl;*/
#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++)
{
- (*it)->compile(code_file, false, tt2, map_idx);
+ (*it)->compile(code_file, false, tt2, map_idx, true);
code_file.write(&FSTPT, sizeof(FSTPT));
map_idx_type::const_iterator ii=map_idx.find((*it)->idx);
v=(int)ii->second;
@@ -1010,12 +958,6 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
}
#endif
- /*if (!lhs_rhs_done)
- {*/
- eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
- lhs = eq_node->get_arg1();
- rhs = eq_node->get_arg2();
- /*}*/
switch (ModelBlock->Block_List[j].Simulation_Type)
{
evaluation:
@@ -1023,17 +965,19 @@ evaluation:
case EVALUATE_FORWARD:
if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE)
{
- rhs->compile(code_file, false, temporary_terms, map_idx);
- lhs->compile(code_file, true, temporary_terms, map_idx);
+ 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);
+ lhs->compile(code_file, true, temporary_terms, map_idx, true);
}
else if (ModelBlock->Block_List[j].Equation_Type[i] == E_EVALUATE_S)
{
eq_node = (BinaryOpNode*)ModelBlock->Block_List[j].Equation_Normalized[i];
- //cout << "EVALUATE_S var " << ModelBlock->Block_List[j].Variable[i] << "\n";
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
- rhs->compile(code_file, false, temporary_terms, map_idx);
- lhs->compile(code_file, true, temporary_terms, map_idx);
+ rhs->compile(code_file, false, temporary_terms, map_idx, true);
+ lhs->compile(code_file, true, temporary_terms, map_idx, true);
}
break;
case SOLVE_BACKWARD_COMPLETE:
@@ -1044,29 +988,27 @@ evaluation:
goto evaluation;
feedback_variables.push_back(ModelBlock->Block_List[j].Variable[i]);
v=ModelBlock->Block_List[j].Equation[i];
- Uf[v].eqr=i;
Uf[v].Ufl=NULL;
goto end;
default:
end:
- lhs->compile(code_file, false, temporary_terms, map_idx);
- rhs->compile(code_file, false, temporary_terms, map_idx);
+ eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
+ lhs = eq_node->get_arg1();
+ rhs = eq_node->get_arg2();
+ lhs->compile(code_file, false, temporary_terms, map_idx, true);
+ rhs->compile(code_file, false, temporary_terms, map_idx, true);
code_file.write(&FBINARY, sizeof(FBINARY));
int v=oMinus;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
code_file.write(&FSTPR, sizeof(FSTPR));
v = i - ModelBlock->Block_List[j].Nb_Recursives;
- //cout << "residual[" << v << "]\n";
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
}
}
code_file.write(&FENDEQU, sizeof(FENDEQU));
// The Jacobian if we have to solve the block
- bool feedback_variable = (feedback_variables.size()>0);
if (ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD
- && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD
- /*&& ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD_R
- && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD_R*/)
+ && ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD)
{
switch (ModelBlock->Block_List[j].Simulation_Type)
{
@@ -1077,87 +1019,15 @@ end:
v=0;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
break;
+
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
- count_u = feedback_variables.size();
- //cout << "todo SOLVE_COMPLETE\n";
- for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- //Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
- pair< pair<int, pair<int, int> >, pair<int, int> > 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;
- int v=ModelBlock->Block_List[j].Equation[eq];
- 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;
- compileChainRuleDerivative(code_file, eqr, varr, 0, map_idx);
- code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast<char *>(&count_u), sizeof(count_u));
- count_u++;
- }
- //cout << "done SOLVE_COMPLETE\n";
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- {
- if(i>=ModelBlock->Block_List[j].Nb_Recursives)
- {
- code_file.write(&FLDR, sizeof(FLDR));
- v = i-ModelBlock->Block_List[j].Nb_Recursives;
- code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
- code_file.write(&FLDZ, sizeof(FLDZ));
- int v=ModelBlock->Block_List[j].Equation[i];
- for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
- {
- code_file.write(&FLDU, sizeof(FLDU));
- code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
- code_file.write(&FLDV, sizeof(FLDV));
- char vc=eEndogenous;
- code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc));
- code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var));
- int v1=0;
- code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
- code_file.write(&FBINARY, sizeof(FBINARY));
- v1=oTimes;
- code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
- code_file.write(&FCUML, sizeof(FCUML));
- }
- Uf[v].Ufl=Uf[v].Ufl_First;
- while (Uf[v].Ufl)
- {
- Uf[v].Ufl_First=Uf[v].Ufl->pNext;
- free(Uf[v].Ufl);
- Uf[v].Ufl=Uf[v].Ufl_First;
- }
- code_file.write(&FBINARY, sizeof(FBINARY));
- v=oMinus;
- code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
- code_file.write(&FSTPU, sizeof(FSTPU));
- v = i - ModelBlock->Block_List[j].Nb_Recursives;
- code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
- }
- }
- break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
- count_u=ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives;
- //cout << "todo SOLVE_TWO_BOUNDARIES\n";
- //cout << "ModelBlock->Block_List[j].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n";
- for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
- {
- //Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
+ //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<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
k=it.first.first;
int eq=it.first.second.first;
@@ -1166,11 +1036,11 @@ end:
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;
+ /*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)
+ {
+ if (!Uf[v].Ufl)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
@@ -1184,53 +1054,21 @@ end:
Uf[v].Ufl->u=count_u;
Uf[v].Ufl->var=varr;
Uf[v].Ufl->lag=k;
- //writeChainRuleDerivative(cout, eqr, varr, k, oMatlabDynamicModelSparse, /*map_idx*/temporary_terms);
- //cout <<"\n";
compileChainRuleDerivative(code_file, eqr, varr, k, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&count_u), sizeof(count_u));
count_u++;
}
}
- //cout << "done it SOLVE_TWO_BOUNDARIES\n";
- /*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- int v=ModelBlock->Block_List[j].Equation[eqr];
- if (!Uf[v].Ufl)
- {
- Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
- Uf[v].Ufl_First=Uf[v].Ufl;
- }
- else
- {
- Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
- Uf[v].Ufl=Uf[v].Ufl->pNext;
- }
- Uf[v].Ufl->pNext=NULL;
- Uf[v].Ufl->u=u;
- Uf[v].Ufl->var=var;
- Uf[v].Ufl->lag=k;
- compileDerivative(code_file, eq, var, k, map_idx);
- code_file.write(&FSTPU, sizeof(FSTPU));
- code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
- }
- }*/
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
- if(i>=ModelBlock->Block_List[j].Nb_Recursives)
- {
+ if(i>=ModelBlock->Block_List[j].Nb_Recursives)
+ {
code_file.write(&FLDR, sizeof(FLDR));
v = i-ModelBlock->Block_List[j].Nb_Recursives;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FLDZ, sizeof(FLDZ));
- int v=ModelBlock->Block_List[j].Equation[i];
+ v=ModelBlock->Block_List[j].Equation[i];
for (Uf[v].Ufl=Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl=Uf[v].Ufl->pNext)
{
code_file.write(&FLDU, sizeof(FLDU));
@@ -1260,24 +1098,8 @@ end:
code_file.write(&FSTPU, sizeof(FSTPU));
v = i - ModelBlock->Block_List[j].Nb_Recursives;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
- }
- }
-#ifdef CONDITION
- for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
- {
- k=m-ModelBlock->Block_List[j].Max_Lag;
- for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
- {
- int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
- int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
- int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];
- int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
- output << " u[" << u << "+Per_u_] /= condition[" << eqr << "];\n";
}
}
- for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
- output << " u[" << i << "+Per_u_] /= condition[" << i << "];\n";
-#endif
break;
default:
break;
@@ -1421,17 +1243,17 @@ DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string
int j;
std::ofstream SaveCode;
if (file_open)
- SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate );
+ SaveCode.open((bin_basename + "_dynamic.bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate );
else
- SaveCode.open((bin_basename + ".bin").c_str(), ios::out | ios::binary);
+ SaveCode.open((bin_basename + "_dynamic.bin").c_str(), ios::out | ios::binary);
if (!SaveCode.is_open())
{
- cout << "Error : Can't open file \"" << bin_basename << ".bin\" for writing\n";
+ cout << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing\n";
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<block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->size();i++)
+ 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<int, pair<int, int> >, pair<int, int> > it = block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->at(i);
@@ -1578,70 +1400,50 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
tmp_eq.str("");
for (int count_call=1, i = 0;i < block_triangular.ModelBlock->Size;i++, 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 ((BlockTriangular::BlockSim(prev_Simulation_Type)!=BlockTriangular::BlockSim(k)) &&
- ((prev_Simulation_Type==EVALUATE_FORWARD || prev_Simulation_Type==EVALUATE_BACKWARD /*|| prev_Simulation_Type==EVALUATE_FORWARD_R || prev_Simulation_Type==EVALUATE_BACKWARD_R*/)
- || (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/)))
- {
- mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
- tmp_eq.str("");
- mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
- tmp.str("");
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
- }
- for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+ if(k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD)
{
- tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
- tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
+ for (int ik=0 ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
+ {
+ tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
+ }
}
- if (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/)
+ else
{
- if (i==block_triangular.ModelBlock->Size-1)
+ for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
- mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
- tmp_eq.str("");
- mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
- tmp.str("");
- mDynamicModelFile << tmp1.str();
- tmp1.str("");
+ tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
+ tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
}
}
- if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) &&
- (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD /*|| k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R*/))
- skip_head=true;
- else
- skip_head=false;
+ mDynamicModelFile << " y_index_eq=[" << tmp_eq.str() << "];\n";
+ mDynamicModelFile << " y_index=[" << tmp.str() << "];\n";
+
switch (k)
{
case EVALUATE_FORWARD:
case EVALUATE_BACKWARD:
- /*case EVALUATE_FORWARD_R:
- case EVALUATE_BACKWARD_R:*/
- if (!skip_head)
- {
- tmp1 << " [y, dr(" << count_call << ").g1, dr(" << count_call << ").g2, dr(" << count_call << ").g3, dr(" << count_call << ").g1_x, dr(" << count_call << ").g1_o]=" << dynamic_basename << "_" << i + 1 << "(y, x, params, 1, it_-1, 1);\n";
- tmp1 << " residual(y_index_eq)=ys(y_index)-y(it_, y_index);\n";
- }
- break;
+ 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 << " 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 << " 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 << " residual(y_index_eq)=r;\n";
- tmp_eq.str("");
- tmp.str("");
break;
case SOLVE_FORWARD_COMPLETE:
case SOLVE_BACKWARD_COMPLETE:
- mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n";
+ //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 << " 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";
+ /*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;j<tmp_i;j++)
@@ -1653,9 +1455,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
for (j=0;j<tmp_ix;j++)
for (int ik=0;ik<block_triangular.ModelBlock->Block_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";
- tmp.str("");
- tmp_eq.str("");
+ 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);
@@ -1670,7 +1470,8 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " residual(y_index_eq)=r(:,M_.maximum_lag+1);\n";
break;
}
- prev_Simulation_Type=k;
+ tmp_eq.str("");
+ tmp.str("");
}
if (tmp1.str().length())
{
@@ -2343,7 +2144,7 @@ DynamicModel::writeOutputPostComputing(ostream &output, const string &basename)
}
void
-DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m)
+DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic)
{
int i=0;
int j=0;
@@ -2376,10 +2177,10 @@ DynamicModel::evaluateJacobian(const eval_context_type &eval_context, jacob_map
IM=block_triangular.incidencematrix.Get_IM(k1, eEndogenous);
a_variable_lag=k1;
}
- if (k1==0)
+ if (k1==0 or !dynamic)
{
j++;
- (*j_m)[make_pair(eq,var)]=val;
+ (*j_m)[make_pair(eq,var)]+=val;
}
if (IM[eq*symbol_table.endo_nbr()+var] && (fabs(val) < cutoff))
{
@@ -2547,8 +2348,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
BuildIncidenceMatrix();
jacob_map j_m;
-
- evaluateJacobian(eval_context, &j_m);
+ evaluateJacobian(eval_context, &j_m, true);
if (block_triangular.bt_verbose)
@@ -2559,12 +2359,15 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
t_etype equation_simulation_type;
map<pair<int, pair<int, int> >, NodeID> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
- block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations, equation_simulation_type, first_order_endo_derivatives);
+ block_triangular.Normalize_and_BlockDecompose_Static_0_Model(j_m, equations, equation_simulation_type, first_order_endo_derivatives, mfs, cutoff);
+
BlockLinear(block_triangular.ModelBlock);
+
+ computeChainRuleJacobian(block_triangular.ModelBlock);
+
if (!no_tmp_terms)
computeTemporaryTermsOrdered(block_triangular.ModelBlock);
- computeChainRuleJacobian(block_triangular.ModelBlock);
}
else
if (!no_tmp_terms)
@@ -2633,6 +2436,20 @@ DynamicModel::toStatic(StaticModel &static_model) const
static_model.addEquation((*it)->toStatic(static_model));
}
+void
+DynamicModel::toStaticDll(StaticDllModel &static_model) const
+ {
+ // Convert model local variables (need to be done first)
+ for (map<int, NodeID>::const_iterator it = local_variables_table.begin();
+ it != local_variables_table.end(); it++)
+ static_model.AddLocalVariable(symbol_table.getName(it->first), it->second->toStatic(static_model));
+
+ // Convert equations
+ for (vector<BinaryOpNode *>::const_iterator it = equations.begin();
+ it != equations.end(); it++)
+ static_model.addEquation((*it)->toStatic(static_model));
+ }
+
int
DynamicModel::computeDerivID(int symb_id, int lag)
{
@@ -2789,17 +2606,13 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti
void
DynamicModel::computeChainRuleJacobian(Model_Block *ModelBlock)
{
- //cout << "computeChainRuleJacobian\n";
- //clock_t t1 = clock();
map<int, NodeID> recursive_variables;
first_chain_rule_derivatives.clear();
for(int blck = 0; blck<ModelBlock->Size; blck++)
{
- //cout << "blck=" << blck << "\n";
recursive_variables.clear();
if (ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
- //cout << "SOLVE_TWO_BOUNDARIES_COMPLETE \n";
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->clear();
for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
{
@@ -2808,56 +2621,41 @@ DynamicModel::computeChainRuleJacobian(Model_Block *ModelBlock)
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]];
}
- //cout << "After recursive_alloc\n";
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> Derivatives = block_triangular.get_Derivatives(ModelBlock, blck);
- for(map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++)
- {
- int lag = it->first.first.first;
- int eq = it->first.first.second.first;
- int var = it->first.first.second.second;
- int eqr = it->first.second.first;
- int varr = it->first.second.second;
- int Deriv_type = it->second;
- if(Deriv_type == 0)
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
- else if (Deriv_type == 1)
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
- else if (Deriv_type == 2)
- {
- if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
- else
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
- }
- ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr)));
- }
-
-
- /*for(int lag = -ModelBlock->Block_List[blck].Max_Lag; lag <= ModelBlock->Block_List[blck].Max_Lead; lag++)
+ map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>::const_iterator it = Derivatives.begin();
+ //#pragma omp parallel for shared(it, blck)
+ for(int i=0; i<(int)Derivatives.size(); i++)
{
-
- for(int eq = 0; eq < ModelBlock->Block_List[blck].Size; eq++)
+ int Deriv_type = it->second;
+ pair<pair<int, pair<int, int> >, pair<int, int> > it_l(it->first);
+ it++;
+ int lag = it_l.first.first;
+ int eq = it_l.first.second.first;
+ int var = it_l.first.second.second;
+ int eqr = it_l.second.first;
+ int varr = it_l.second.second;
+ if(Deriv_type == 0)
{
- int eqr = ModelBlock->Block_List[blck].Equation[eq];
- for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++)
- {
- int varr = ModelBlock->Block_List[blck].Variable[var];
- NodeID d1 = equations[eqr]->getChainRuleDerivative(recursive_variables, varr, lag);
- if (d1 == Zero)
- continue;
- first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = d1;
- ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
- }
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
}
- }*/
-
-
+ else if (Deriv_type == 1)
+ {
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+ }
+ else if (Deriv_type == 2)
+ {
+ if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S && eq<ModelBlock->Block_List[blck].Nb_Recursives)
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+ else
+ first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
+ }
+ ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr)));
+ }
}
else if( ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_SIMPLE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_SIMPLE
or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_COMPLETE)
{
- //cout << "SOLVE_FORWARD_SIMPLE \n";
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->clear();
for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
{
@@ -2881,7 +2679,6 @@ DynamicModel::computeChainRuleJacobian(Model_Block *ModelBlock)
}
}
}
- //cout << "elapsed time in milliseconds = " << 1000.0*(double(clock()) - double(t1))/double(CLOCKS_PER_SEC) << "\n";
}
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 54263fbecc22841922dcecb1305ea7b45e705429..ee05ea9546316b56694b12dce96eaaf59aea6bd6 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -25,6 +25,7 @@ using namespace std;
#include <fstream>
#include "StaticModel.hh"
+#include "StaticDllModel.hh"
#include "BlockTriangular.hh"
//! Stores a dynamic model
@@ -100,7 +101,7 @@ private:
- computes the jacobian for the model w.r. to contemporaneous variables
- removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
*/
- void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m);
+ void evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_m, bool dynamic);
void BlockLinear(Model_Block *ModelBlock);
string reform(string name) const;
map_idx_type map_idx;
@@ -152,8 +153,15 @@ public:
virtual NodeID AddVariable(const string &name, int lag = 0);
//! Absolute value under which a number is considered to be zero
double cutoff;
- //! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 from simulate.cc)
+ //! The weight of the Markowitz criteria to determine the pivot in the linear solver (simul_NG1 and simul_NG from simulate.cc)
double markowitz;
+ //! Compute the minimum feedback set in the dynamic model:
+ /*! 0 : all endogenous variables are considered as feedback variables
+ 1 : the variables belonging to non normalized equation are considered as feedback variables
+ 2 : the variables belonging to a non linear equation are considered as feedback variables
+ 3 : the variables belonging to a non normalizable non linear equation are considered as feedback variables
+ default value = 0 */
+ int mfs;
//! the file containing the model and the derivatives code
ofstream code_file;
//! Execute computations (variable sorting + derivation)
@@ -184,6 +192,8 @@ public:
/*! It assumes that the static model given in argument has just been allocated */
void toStatic(StaticModel &static_model) const;
+ void toStaticDll(StaticDllModel &static_model) const;
+
//! Writes LaTeX file with the equations of the dynamic model
void writeLatexFile(const string &basename) const;
diff --git a/DynareBison.yy b/DynareBison.yy
index d7423e8288c4702c5395dc4b49cda9f81dcfe255..439c7532e706b37f9bbd450469803a6bf1b21c2e 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -87,7 +87,7 @@ class ParsingDriver;
%}
%token AR AUTOCORR
-%token BAYESIAN_IRF BETA_PDF BICGSTAB BLOCK_MFS
+%token BAYESIAN_IRF BETA_PDF BICGSTAB BLOCK_MFS BLOCK_MFS_DLL
%token BVAR_DENSITY BVAR_FORECAST
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
@@ -106,7 +106,7 @@ class ParsingDriver;
%token KALMAN_ALGO KALMAN_TOL
%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_MH_FILE LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LU
%token MARKOWITZ MARGINAL_DENSITY MAX
-%token METHOD MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN
+%token METHOD MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN
%token MODE_CHECK MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS
%token MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER
%token <string_val> NAME
@@ -449,6 +449,7 @@ model_sparse_options_list : model_sparse_options_list COMMA model_sparse_options
model_sparse_options : o_cutoff
| o_markowitz
+ | o_mfs
;
model : MODEL ';' { driver.begin_model(); }
@@ -641,6 +642,10 @@ steady_options : o_solve_algo
| o_homotopy_mode
| o_homotopy_steps
| o_block_mfs
+ | o_block_mfs_dll
+ | o_cutoff
+ | o_markowitz
+ | o_mfs
;
check : CHECK ';'
@@ -1496,6 +1501,7 @@ o_method : METHOD EQUAL INT_NUMBER { driver.option_num("simulation_method",$3);}
| METHOD EQUAL GMRES { driver.option_num("simulation_method", "2"); }
| METHOD EQUAL BICGSTAB { driver.option_num("simulation_method", "3"); };
o_markowitz : MARKOWITZ EQUAL number { driver.option_num("markowitz", $3); };
+o_mfs : MFS EQUAL number { driver.option_num("mfs", $3); };
o_simul : SIMUL { driver.option_num("simul", "1"); };
o_simul_seed : SIMUL_SEED EQUAL INT_NUMBER { driver.option_num("simul_seed", $3); } ;
o_qz_criterium : QZ_CRITERIUM EQUAL number { driver.option_num("qz_criterium", $3); };
@@ -1602,6 +1608,8 @@ o_gsa_trans_ident : TRANS_IDENT EQUAL INT_NUMBER { driver.option_num("trans_iden
o_homotopy_mode : HOMOTOPY_MODE EQUAL INT_NUMBER {driver.option_num("homotopy_mode",$3); };
o_homotopy_steps : HOMOTOPY_STEPS EQUAL INT_NUMBER {driver.option_num("homotopy_steps",$3); };
o_block_mfs : BLOCK_MFS { driver.option_num("block_mfs", "1"); }
+o_block_mfs_dll : BLOCK_MFS_DLL { driver.option_num("block_mfs_dll", "1"); }
+
o_parameters : PARAMETERS EQUAL symbol {driver.option_str("parameters",$3);};
o_shocks : SHOCKS EQUAL '(' list_of_symbol_lists ')' { driver.option_symbol_list("shocks"); };
o_labels : LABELS EQUAL '(' symbol_list ')' { driver.option_symbol_list("labels"); };
diff --git a/DynareFlex.ll b/DynareFlex.ll
index 0180d384f1f22e518794d729197f42d054782c53..395050c03e5871efbde0ccb3103ae381c34d5026 100644
--- a/DynareFlex.ll
+++ b/DynareFlex.ll
@@ -210,6 +210,7 @@ int sigma_e = 0;
<DYNARE_STATEMENT>periods {return token::PERIODS;}
<DYNARE_STATEMENT>cutoff {return token::CUTOFF;}
<DYNARE_STATEMENT>markowitz {return token::MARKOWITZ;}
+<DYNARE_STATEMENT>mfs {return token::MFS;}
<DYNARE_STATEMENT>marginal_density {return token::MARGINAL_DENSITY;}
<DYNARE_STATEMENT>laplace {return token::LAPLACE;}
<DYNARE_STATEMENT>modifiedharmonicmean {return token::MODIFIEDHARMONICMEAN;}
@@ -220,6 +221,8 @@ int sigma_e = 0;
<DYNARE_STATEMENT>diffuse_filter {return token::DIFFUSE_FILTER;}
<DYNARE_STATEMENT>plot_priors {return token::PLOT_PRIORS;}
<DYNARE_STATEMENT>block_mfs {return token::BLOCK_MFS;}
+<DYNARE_STATEMENT>block_mfs_dll {return token::BLOCK_MFS_DLL;}
+
<DYNARE_STATEMENT>freq {return token::FREQ;}
<DYNARE_STATEMENT>initial_year {return token::INITIAL_YEAR;}
<DYNARE_STATEMENT>initial_subperiod {return token::INITIAL_SUBPERIOD;}
@@ -302,6 +305,7 @@ int sigma_e = 0;
<DYNARE_BLOCK>periods {return token::PERIODS;}
<DYNARE_BLOCK>cutoff {return token::CUTOFF;}
<DYNARE_BLOCK>markowitz {return token::MARKOWITZ;}
+<DYNARE_BLOCK>mfs {return token::MFS;}
<DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;}
<DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;}
<DYNARE_BLOCK>normal_pdf {return token::NORMAL_PDF;}
diff --git a/ExprNode.cc b/ExprNode.cc
index 5c799fcc0ac50475a2f68988e2f85e3a7fe065dd..6acccd01ba96c261ff6cc3f9ed4ef9e3158d87ca 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -131,6 +131,13 @@ ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
}
+pair<int, NodeID >
+ExprNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+ {
+ return(make_pair(0, (NodeID)NULL));
+ }
+
+
void
ExprNode::writeOutput(ostream &output)
{
@@ -184,13 +191,10 @@ NumConstNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
}
void
-NumConstNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
+NumConstNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const
{
CompileCode.write(&FLDC, sizeof(FLDC));
double vard = datatree.num_constants.getDouble(id);
-#ifdef DEBUGC
- cout << "FLDC " << vard << "\n";
-#endif
CompileCode.write(reinterpret_cast<char *>(&vard),sizeof(vard));
}
@@ -199,10 +203,10 @@ NumConstNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result
{
}
-pair<bool, NodeID>
-NumConstNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
+pair<int, NodeID >
+NumConstNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
{
- return(make_pair(false, datatree.AddNumConstant(datatree.num_constants.get(id))));
+ return(make_pair(0, datatree.AddNumConstant(datatree.num_constants.get(id))));
}
NodeID
@@ -295,9 +299,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
// If node is a temporary term
-#ifdef DEBUGC
- cout << "write_ouput Variable output_type=" << output_type << "\n";
-#endif
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<VariableNode *>(this));
if (it != temporary_terms.end())
{
@@ -464,19 +465,22 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
}
void
-VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
+VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const
{
int i, lagl;
-#ifdef DEBUGC
- cout << "output_type=" << output_type << "\n";
-#endif
if (!lhs_rhs)
{
- CompileCode.write(&FLDV, sizeof(FLDV));
+ if(dynamic)
+ CompileCode.write(&FLDV, sizeof(FLDV));
+ else
+ CompileCode.write(&FLDSV, sizeof(FLDSV));
}
else
{
- CompileCode.write(&FSTPV, sizeof(FSTPV));
+ if(dynamic)
+ CompileCode.write(&FSTPV, sizeof(FSTPV));
+ else
+ CompileCode.write(&FSTPSV, sizeof(FSTPSV));
}
char typel=(char)type;
CompileCode.write(&typel, sizeof(typel));
@@ -487,35 +491,41 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
//cout << "Parameter=" << tsid << "\n";
i = tsid;
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
-#ifdef DEBUGC
- cout << "FLD Param[ " << i << ", symb_id=" << symb_id << "]\n";
-#endif
break;
case eEndogenous :
//cout << "Endogenous=" << symb_id << "\n";
i = tsid;//symb_id;
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
- lagl=lag;
- CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
+ if(dynamic)
+ {
+ lagl=lag;
+ CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
+ }
break;
case eExogenous :
//cout << "Exogenous=" << tsid << "\n";
i = tsid;
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
- lagl=lag;
- CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
+ if(dynamic)
+ {
+ lagl=lag;
+ CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
+ }
break;
case eExogenousDet:
i = tsid + datatree.symbol_table.exo_nbr();
//cout << "ExogenousDet=" << i << "\n";
CompileCode.write(reinterpret_cast<char *>(&i), sizeof(i));
- lagl=lag;
- CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
+ if(dynamic)
+ {
+ lagl=lag;
+ CompileCode.write(reinterpret_cast<char *>(&lagl), sizeof(lagl));
+ }
break;
case eModelLocalVariable:
case eModFileLocalVariable:
//cout << "eModelLocalVariable=" << symb_id << "\n";
- datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
+ datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
break;
case eUnknownFunction:
cerr << "Impossible case: eUnknownFuncion" << endl;
@@ -545,22 +555,22 @@ VariableNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result
datatree.local_variables_table[symb_id]->collectVariables(type_arg, result);
}
-pair<bool, NodeID>
-VariableNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
+pair<int, NodeID>
+VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
{
if (type ==eEndogenous)
{
if (datatree.symbol_table.getTypeSpecificID(symb_id)==var_endo && lag==0)
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
else
- return(make_pair(false, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), lag)));
+ return(make_pair(0, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), lag) ));
}
else
{
if (type == eParameter)
- return(make_pair(false, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), 0)));
+ return(make_pair(0, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), 0) ));
else
- return(make_pair(false, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), lag)));
+ return(make_pair(0, datatree.AddVariableInternal(datatree.symbol_table.getName(symb_id), lag) ));
}
}
@@ -581,9 +591,19 @@ VariableNode::getChainRuleDerivative(int deriv_id_arg, const map<int, NodeID> &r
map<int, NodeID>::const_iterator it = recursive_variables.find(deriv_id);
if (it != recursive_variables.end())
{
- map<int, NodeID> recursive_vars2(recursive_variables);
- recursive_vars2.erase(it->first);
- return datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id_arg, recursive_vars2));
+ map<int, NodeID>::const_iterator it2 = derivatives.find(deriv_id_arg);
+ if (it2 != derivatives.end())
+ return it2->second;
+ else
+ {
+ map<int, NodeID> recursive_vars2(recursive_variables);
+ recursive_vars2.erase(it->first);
+ //NodeID c = datatree.AddNumConstant("1");
+ NodeID d = datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id_arg, recursive_vars2));
+ //d = datatree.AddTimes(c, d);
+ derivatives[deriv_id_arg] = d;
+ return d;
+ }
}
else
return datatree.Zero;
@@ -997,17 +1017,20 @@ UnaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExcept
}
void
-UnaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
+UnaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const
{
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
if (it != temporary_terms.end())
{
- CompileCode.write(&FLDT, sizeof(FLDT));
+ if(dynamic)
+ CompileCode.write(&FLDT, sizeof(FLDT));
+ else
+ CompileCode.write(&FLDST, sizeof(FLDST));
int var=map_idx[idx];
CompileCode.write(reinterpret_cast<char *>(&var), sizeof(var));
return;
}
- arg->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
+ arg->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
CompileCode.write(&FUNARY, sizeof(FUNARY));
UnaryOpcode op_codel=op_code;
CompileCode.write(reinterpret_cast<char *>(&op_codel), sizeof(op_codel));
@@ -1019,53 +1042,101 @@ UnaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result)
arg->collectVariables(type_arg, result);
}
-pair<bool, NodeID>
-UnaryOpNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
+pair<int, NodeID>
+UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
{
- pair<bool, NodeID> res = arg->normalizeLinearInEndoEquation(var_endo, Derivative);
- bool is_endogenous_present = res.first;
+ pair<bool, NodeID > res = arg->normalizeEquation(var_endo, List_of_Op_RHS);
+ int is_endogenous_present = res.first;
NodeID New_NodeID = res.second;
- if (!is_endogenous_present)
+ /*if(res.second.second)*/
+ if(is_endogenous_present==2)
+ return(make_pair(2, (NodeID)NULL));
+ else if (is_endogenous_present)
+ {
+ switch (op_code)
+ {
+ case oUminus:
+ List_of_Op_RHS.push_back(make_pair(oUminus, make_pair((NodeID)NULL, (NodeID)NULL)));
+ return(make_pair(1, (NodeID)NULL));
+ case oExp:
+ List_of_Op_RHS.push_back(make_pair(oLog, make_pair((NodeID)NULL, (NodeID)NULL)));
+ return(make_pair(1, (NodeID)NULL));
+ case oLog:
+ List_of_Op_RHS.push_back(make_pair(oExp, make_pair((NodeID)NULL, (NodeID)NULL)));
+ return(make_pair(1, (NodeID)NULL));
+ case oLog10:
+ List_of_Op_RHS.push_back(make_pair(oPower, make_pair((NodeID)NULL, datatree.AddNumConstant("10"))));
+ return(make_pair(1, (NodeID)NULL));
+ case oCos:
+ return(make_pair(1, (NodeID)NULL));
+ case oSin:
+ return(make_pair(1, (NodeID)NULL));
+ case oTan:
+ return(make_pair(1, (NodeID)NULL));
+ case oAcos:
+ return(make_pair(1, (NodeID)NULL));
+ case oAsin:
+ return(make_pair(1, (NodeID)NULL));
+ case oAtan:
+ return(make_pair(1, (NodeID)NULL));
+ case oCosh:
+ return(make_pair(1, (NodeID)NULL));
+ case oSinh:
+ return(make_pair(1, (NodeID)NULL));
+ case oTanh:
+ return(make_pair(1, (NodeID)NULL));
+ case oAcosh:
+ return(make_pair(1, (NodeID)NULL));
+ case oAsinh:
+ return(make_pair(1, (NodeID)NULL));
+ case oAtanh:
+ return(make_pair(1, (NodeID)NULL));
+ case oSqrt:
+ List_of_Op_RHS.push_back(make_pair(oPower, make_pair((NodeID)NULL, datatree.AddNumConstant("2"))));
+ return(make_pair(1, (NodeID)NULL));
+ }
+ }
+ else
{
switch (op_code)
{
case oUminus:
- return(make_pair(false, /*tmp_*/datatree.AddUMinus(New_NodeID)));
+ return(make_pair(0, datatree.AddUMinus(New_NodeID)));
case oExp:
- return(make_pair(false, /*tmp_*/datatree.AddExp(New_NodeID)));
+ return(make_pair(0, datatree.AddExp(New_NodeID)));
case oLog:
- return(make_pair(false, /*tmp_*/datatree.AddLog(New_NodeID)));
+ return(make_pair(0, datatree.AddLog(New_NodeID)));
case oLog10:
- return(make_pair(false, /*tmp_*/datatree.AddLog10(New_NodeID)));
+ return(make_pair(0, datatree.AddLog10(New_NodeID)));
case oCos:
- return(make_pair(false, /*tmp_*/datatree.AddCos(New_NodeID)));
+ return(make_pair(0, datatree.AddCos(New_NodeID)));
case oSin:
- return(make_pair(false, /*tmp_*/datatree.AddSin(New_NodeID)));
+ return(make_pair(0, datatree.AddSin(New_NodeID)));
case oTan:
- return(make_pair(false, /*tmp_*/datatree.AddTan(New_NodeID)));
+ return(make_pair(0, datatree.AddTan(New_NodeID)));
case oAcos:
- return(make_pair(false, /*tmp_*/datatree.AddAcos(New_NodeID)));
+ return(make_pair(0, datatree.AddAcos(New_NodeID)));
case oAsin:
- return(make_pair(false, /*tmp_*/datatree.AddAsin(New_NodeID)));
+ return(make_pair(0, datatree.AddAsin(New_NodeID)));
case oAtan:
- return(make_pair(false, /*tmp_*/datatree.AddAtan(New_NodeID)));
+ return(make_pair(0, datatree.AddAtan(New_NodeID)));
case oCosh:
- return(make_pair(false, /*tmp_*/datatree.AddCosh(New_NodeID)));
+ return(make_pair(0, datatree.AddCosh(New_NodeID)));
case oSinh:
- return(make_pair(false, /*tmp_*/datatree.AddSinh(New_NodeID)));
+ return(make_pair(0, datatree.AddSinh(New_NodeID)));
case oTanh:
- return(make_pair(false, /*tmp_*/datatree.AddTanh(New_NodeID)));
+ return(make_pair(0, datatree.AddTanh(New_NodeID)));
case oAcosh:
- return(make_pair(false, /*tmp_*/datatree.AddAcosh(New_NodeID)));
+ return(make_pair(0, datatree.AddAcosh(New_NodeID)));
case oAsinh:
- return(make_pair(false, /*tmp_*/datatree.AddAsinh(New_NodeID)));
+ return(make_pair(0, datatree.AddAsinh(New_NodeID)));
case oAtanh:
- return(make_pair(false, /*tmp_*/datatree.AddAtanh(New_NodeID)));
+ return(make_pair(0, datatree.AddAtanh(New_NodeID)));
case oSqrt:
- return(make_pair(false, /*tmp_*/datatree.AddSqrt(New_NodeID)));
+ return(make_pair(0, datatree.AddSqrt(New_NodeID)));
}
}
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL));
}
@@ -1435,19 +1506,22 @@ BinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExcep
}
void
-BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
+BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const
{
// If current node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
if (it != temporary_terms.end())
{
- CompileCode.write(&FLDT, sizeof(FLDT));
+ if(dynamic)
+ CompileCode.write(&FLDT, sizeof(FLDT));
+ else
+ CompileCode.write(&FLDST, sizeof(FLDST));
int var=map_idx[idx];
CompileCode.write(reinterpret_cast<char *>(&var), sizeof(var));
return;
}
- arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
- arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
+ arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
+ arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
CompileCode.write(&FBINARY, sizeof(FBINARY));
BinaryOpcode op_codel=op_code;
CompileCode.write(reinterpret_cast<char *>(&op_codel),sizeof(op_codel));
@@ -1636,129 +1710,272 @@ BinaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result
arg2->collectVariables(type_arg, result);
}
-pair<bool, NodeID>
-BinaryOpNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
+NodeID
+BinaryOpNode::Compute_RHS(NodeID arg1, NodeID arg2, int op, int op_type) const
+{
+ temporary_terms_type temp;
+ switch(op_type)
+ {
+ case 0: /*Unary Operator*/
+ switch(op)
+ {
+ case oUminus:
+ return(datatree.AddUMinus(arg1));
+ break;
+ case oExp:
+ return(datatree.AddExp(arg1));
+ break;
+ case oLog:
+ return(datatree.AddLog(arg1));
+ break;
+ case oLog10:
+ return(datatree.AddLog10(arg1));
+ break;
+ }
+ break;
+ case 1: /*Binary Operator*/
+ switch(op)
+ {
+ case oPlus:
+ return(datatree.AddPlus(arg1, arg2));
+ break;
+ case oMinus:
+ return(datatree.AddMinus(arg1, arg2));
+ break;
+ case oTimes:
+ return(datatree.AddTimes(arg1, arg2));
+ break;
+ case oDivide:
+ return(datatree.AddDivide(arg1, arg2));
+ break;
+ case oPower:
+ return(datatree.AddPower(arg1, arg2));
+ break;
+ }
+ break;
+ }
+ return((NodeID)NULL);
+}
+
+pair<int, NodeID>
+BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
{
- pair<bool, NodeID> res = arg1->normalizeLinearInEndoEquation(var_endo, Derivative);
- bool is_endogenous_present_1 = res.first;
- NodeID NodeID_1 = res.second;
- res = arg2->normalizeLinearInEndoEquation(var_endo, Derivative);
- bool is_endogenous_present_2 = res.first;
- NodeID NodeID_2 = res.second;
+ vector<pair<int, pair<NodeID, NodeID> > > List_of_Op_RHS1, List_of_Op_RHS2;
+ int is_endogenous_present_1, is_endogenous_present_2;
+ pair<int, NodeID> res;
+ NodeID NodeID_1, NodeID_2;
+ res = arg1->normalizeEquation(var_endo, List_of_Op_RHS1);
+ is_endogenous_present_1 = res.first;
+ NodeID_1 = res.second;
+
+ res = arg2->normalizeEquation(var_endo, List_of_Op_RHS2);
+ is_endogenous_present_2 = res.first;
+ NodeID_2 = res.second;
+ if(is_endogenous_present_1==2 || is_endogenous_present_2==2)
+ return(make_pair(2,(NodeID)NULL));
+ else if(is_endogenous_present_1 && is_endogenous_present_2)
+ return(make_pair(2,(NodeID)NULL));
+ else if(is_endogenous_present_1)
+ {
+ if(op_code==oEqual)
+ {
+ pair<int, pair<NodeID, NodeID> > it;
+ int oo=List_of_Op_RHS1.size();
+ for(int i=0;i<oo;i++)
+ {
+ it = List_of_Op_RHS1.back();
+ List_of_Op_RHS1.pop_back();
+ if(it.second.first && !it.second.second) /*Binary operator*/
+ NodeID_2 = Compute_RHS(NodeID_2, (BinaryOpNode*)it.second.first, it.first, 1);
+ else if(it.second.second && !it.second.first) /*Binary operator*/
+ NodeID_2 = Compute_RHS(it.second.second, NodeID_2, it.first, 1);
+ else if(it.second.second && it.second.first) /*Binary operator*/
+ NodeID_2 = Compute_RHS(it.second.first, it.second.second, it.first, 1);
+ else /*Unary operator*/
+ NodeID_2 = Compute_RHS((UnaryOpNode*)NodeID_2, (UnaryOpNode*)it.second.first, it.first, 0);
+ }
+ }
+ else
+ List_of_Op_RHS = List_of_Op_RHS1;
+ }
+ else if(is_endogenous_present_2)
+ {
+ if(op_code==oEqual)
+ {
+ int oo=List_of_Op_RHS2.size();
+ for(int i=0;i<oo;i++)
+ {
+ pair<int, pair<NodeID, NodeID> > it;
+ it = List_of_Op_RHS2.back();
+ List_of_Op_RHS2.pop_back();
+ if(it.second.first && !it.second.second) /*Binary operator*/
+ NodeID_1 = Compute_RHS((BinaryOpNode*)NodeID_1, (BinaryOpNode*)it.second.first, it.first, 1);
+ else if(it.second.second && !it.second.first) /*Binary operator*/
+ NodeID_1 = Compute_RHS((BinaryOpNode*)it.second.second, (BinaryOpNode*)NodeID_1, it.first, 1);
+ else if(it.second.second && it.second.first) /*Binary operator*/
+ NodeID_1 = Compute_RHS(it.second.first, it.second.second, it.first, 1);
+ else
+ NodeID_1 = Compute_RHS((UnaryOpNode*)NodeID_1, (UnaryOpNode*)it.second.first, it.first, 0);
+ }
+ }
+ else
+ List_of_Op_RHS =List_of_Op_RHS2;
+ }
switch (op_code)
{
case oPlus:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddPlus(NodeID_1, NodeID_2)));
+ {
+ List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(datatree.AddPlus(NodeID_1, NodeID_2), (NodeID)NULL)));
+ return(make_pair(0, datatree.AddPlus(NodeID_1, NodeID_2)));
+ }
else if (is_endogenous_present_1 && is_endogenous_present_2)
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL));
else if (!is_endogenous_present_1 && is_endogenous_present_2)
- return(make_pair(false, NodeID_1));
+ {
+ List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(NodeID_1, (NodeID)NULL)));
+ return(make_pair(1, NodeID_1));
+ }
else if (is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, NodeID_2));
+ {
+ List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(NodeID_2, (NodeID)NULL) ));
+ return(make_pair(1, NodeID_2));
+ }
break;
case oMinus:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddMinus(NodeID_1, NodeID_2)));
+ {
+ List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(datatree.AddMinus(NodeID_1, NodeID_2), (NodeID)NULL) ));
+ return(make_pair(0, datatree.AddMinus(NodeID_1, NodeID_2)));
+ }
else if (is_endogenous_present_1 && is_endogenous_present_2)
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL));
else if (!is_endogenous_present_1 && is_endogenous_present_2)
- return(make_pair(false, NodeID_1));
+ {
+ List_of_Op_RHS.push_back(make_pair(oUminus, make_pair((NodeID)NULL, (NodeID)NULL)));
+ List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(NodeID_1, (NodeID)NULL) ));
+ return(make_pair(1, NodeID_1));
+ }
else if (is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddUMinus(NodeID_2)));
+ {
+ List_of_Op_RHS.push_back(make_pair(oPlus, make_pair(NodeID_2, (NodeID) NULL) ));
+ return(make_pair(1, datatree.AddUMinus(NodeID_2)));
+ }
break;
case oTimes:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddTimes(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddTimes(NodeID_1, NodeID_2)));
+ else if(!is_endogenous_present_1 && is_endogenous_present_2)
+ {
+ List_of_Op_RHS.push_back(make_pair(oDivide, make_pair(NodeID_1, (NodeID)NULL) ));
+ return(make_pair(1, NodeID_1));
+ }
+ else if(is_endogenous_present_1 && !is_endogenous_present_2)
+ {
+ List_of_Op_RHS.push_back(make_pair(oDivide, make_pair(NodeID_2, (NodeID)NULL) ));
+ return(make_pair(1, NodeID_2));
+ }
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL));
break;
case oDivide:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, datatree.AddDivide(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddDivide(NodeID_1, NodeID_2)));
+ else if(!is_endogenous_present_1 && is_endogenous_present_2)
+ {
+ List_of_Op_RHS.push_back(make_pair(oDivide, make_pair((NodeID)NULL, NodeID_1) ));
+ return(make_pair(1, NodeID_1));
+ }
+ else if(is_endogenous_present_1 && !is_endogenous_present_2)
+ {
+ List_of_Op_RHS.push_back(make_pair(oTimes, make_pair(NodeID_2, (NodeID)NULL) ));
+ return(make_pair(1, NodeID_2));
+ }
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL));
break;
case oPower:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, datatree.AddPower(NodeID_1, NodeID_2)));
- else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(0, datatree.AddPower(NodeID_1, NodeID_2)));
+ else if(is_endogenous_present_1 && !is_endogenous_present_2)
+ {
+ List_of_Op_RHS.push_back(make_pair(oPower, make_pair(datatree.AddDivide( datatree.AddNumConstant("1"), NodeID_2), (NodeID)NULL) ));
+ return(make_pair(1, (NodeID)NULL));
+ }
break;
case oEqual:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
{
- if (Derivative!=datatree.One)
- return( make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddDivide(datatree.AddMinus(NodeID_2, NodeID_1), Derivative))) );
- else
- return( make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddMinus(NodeID_2, NodeID_1))) );
+ return( make_pair(0,
+ datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddMinus(NodeID_2, NodeID_1))
+ ));
}
else if (is_endogenous_present_1 && is_endogenous_present_2)
{
- return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.Zero)));
+ return(make_pair(0,
+ datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.Zero)
+ ));
}
else if (!is_endogenous_present_1 && is_endogenous_present_2)
{
- if (Derivative!=datatree.One)
- return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddDivide(datatree.AddUMinus(NodeID_1), Derivative))));
- else
- return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddUMinus(NodeID_1))));
+ return(make_pair(0,
+ datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), /*datatree.AddUMinus(NodeID_1)*/NodeID_1)
+ ));
}
else if (is_endogenous_present_1 && !is_endogenous_present_2)
{
- if (Derivative!=datatree.One)
- return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), datatree.AddDivide(NodeID_2, Derivative))));
- else
- return(make_pair(false, datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), NodeID_2)));
+ return(make_pair(0,
+ datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getName(datatree.symbol_table.getID(eEndogenous, var_endo)), 0), NodeID_2)
+ ));
}
break;
case oMax:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, datatree.AddMax(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddMax(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oMin:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, datatree.AddMin(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddMin(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oLess:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddLess(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddLess(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oGreater:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddGreater(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddGreater(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oLessEqual:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddLessEqual(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddLessEqual(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oGreaterEqual:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddGreaterEqual(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddGreaterEqual(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oEqualEqual:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddEqualEqual(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddEqualEqual(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
case oDifferent:
if (!is_endogenous_present_1 && !is_endogenous_present_2)
- return(make_pair(false, /*tmp_*/datatree.AddDifferent(NodeID_1, NodeID_2)));
+ return(make_pair(0, datatree.AddDifferent(NodeID_1, NodeID_2) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
break;
}
// Suppress GCC warning
@@ -2023,20 +2240,23 @@ TrinaryOpNode::eval(const eval_context_type &eval_context) const throw (EvalExce
}
void
-TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
+TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const
{
// If current node is a temporary term
temporary_terms_type::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
if (it != temporary_terms.end())
{
- CompileCode.write(&FLDT, sizeof(FLDT));
+ if(dynamic)
+ CompileCode.write(&FLDT, sizeof(FLDT));
+ else
+ CompileCode.write(&FLDST, sizeof(FLDST));
int var=map_idx[idx];
CompileCode.write(reinterpret_cast<char *>(&var), sizeof(var));
return;
}
- arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
- arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
- arg3->compile(CompileCode, lhs_rhs, temporary_terms, map_idx);
+ arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
+ arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
+ arg3->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic);
CompileCode.write(&FBINARY, sizeof(FBINARY));
TrinaryOpcode op_codel=op_code;
CompileCode.write(reinterpret_cast<char *>(&op_codel),sizeof(op_codel));
@@ -2094,22 +2314,22 @@ TrinaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &resul
arg3->collectVariables(type_arg, result);
}
-pair<bool, NodeID>
-TrinaryOpNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
+pair<int, NodeID>
+TrinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
{
- pair<bool, NodeID> res = arg1->normalizeLinearInEndoEquation(var_endo, Derivative);
+ pair<int, NodeID> res = arg1->normalizeEquation(var_endo, List_of_Op_RHS);
bool is_endogenous_present_1 = res.first;
NodeID NodeID_1 = res.second;
- res = arg2->normalizeLinearInEndoEquation(var_endo, Derivative);
+ res = arg2->normalizeEquation(var_endo, List_of_Op_RHS);
bool is_endogenous_present_2 = res.first;
NodeID NodeID_2 = res.second;
- res = arg3->normalizeLinearInEndoEquation(var_endo, Derivative);
+ res = arg3->normalizeEquation(var_endo, List_of_Op_RHS);
bool is_endogenous_present_3 = res.first;
NodeID NodeID_3 = res.second;
if (!is_endogenous_present_1 && !is_endogenous_present_2 && !is_endogenous_present_3)
- return(make_pair(false, /*tmp_*/datatree.AddNormcdf(NodeID_1, NodeID_2, NodeID_3)));
+ return(make_pair(0, datatree.AddNormcdf(NodeID_1, NodeID_2, NodeID_3) ));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
}
NodeID
@@ -2226,14 +2446,14 @@ UnknownFunctionNode::eval(const eval_context_type &eval_context) const throw (Ev
}
void
-UnknownFunctionNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const
+UnknownFunctionNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const
{
cerr << "UnknownFunctionNode::compile: operation impossible!" << endl;
exit(EXIT_FAILURE);
}
-pair<bool, NodeID>
-UnknownFunctionNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivative) const
+pair<int, NodeID>
+UnknownFunctionNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
{
vector<pair<bool, NodeID> > V_arguments;
vector<NodeID> V_NodeID;
@@ -2241,14 +2461,14 @@ UnknownFunctionNode::normalizeLinearInEndoEquation(int var_endo, NodeID Derivati
for (vector<NodeID>::const_iterator it = arguments.begin();
it != arguments.end(); it++)
{
- V_arguments.push_back((*it)->normalizeLinearInEndoEquation(var_endo, Derivative));
+ V_arguments.push_back((*it)->normalizeEquation(var_endo, List_of_Op_RHS));
present = present || V_arguments[V_arguments.size()-1].first;
V_NodeID.push_back(V_arguments[V_arguments.size()-1].second);
}
if (!present)
- return(make_pair(false, datatree.AddUnknownFunction(datatree.symbol_table.getName(symb_id), V_NodeID)));
+ return(make_pair(0, datatree.AddUnknownFunction(datatree.symbol_table.getName(symb_id), V_NodeID)));
else
- return(make_pair(true, (NodeID)NULL));
+ return(make_pair(1, (NodeID)NULL ));
}
NodeID
diff --git a/ExprNode.hh b/ExprNode.hh
index cb2c751a386cef3de76c14ca8f21cbfb6dae39a3..b4d0b0505785c2c384b435146705e1fd06c7194a 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -94,6 +94,7 @@ class ExprNode
{
friend class DataTree;
friend class DynamicModel;
+ friend class StaticDllModel;
friend class ExprNodeLess;
friend class NumConstNode;
friend class VariableNode;
@@ -200,7 +201,7 @@ public:
};
virtual double eval(const eval_context_type &eval_context) const throw (EvalException) = 0;
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const = 0;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const = 0;
//! Creates a static version of this node
/*!
This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
@@ -208,7 +209,7 @@ public:
*/
virtual NodeID toStatic(DataTree &static_datatree) const = 0;
//! Try to normalize an equation linear in its endogenous variable
- virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const = 0;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const = 0;
};
//! Object used to compare two nodes (using their indexes)
@@ -234,9 +235,9 @@ public:
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
- virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const;
virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
};
@@ -264,10 +265,10 @@ public:
map_idx_type &map_idx) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
int get_symb_id() const { return symb_id; };
- virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const;
virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
};
@@ -296,13 +297,13 @@ public:
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const;
//! Returns operand
NodeID get_arg() const { return(arg); };
//! Returns op code
UnaryOpcode get_op_code() const { return(op_code); };
virtual NodeID toStatic(DataTree &static_datatree) const;
- virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const;
virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
};
@@ -333,7 +334,8 @@ public:
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const;
+ virtual NodeID Compute_RHS(NodeID arg1, NodeID arg2, int op, int op_type) const;
//! Returns first operand
NodeID get_arg1() const { return(arg1); };
//! Returns second operand
@@ -341,7 +343,7 @@ public:
//! Returns op code
BinaryOpcode get_op_code() const { return(op_code); };
virtual NodeID toStatic(DataTree &static_datatree) const;
- pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const;
virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
};
@@ -373,9 +375,9 @@ public:
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
- virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const;
virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
};
@@ -401,9 +403,9 @@ public:
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_type &temporary_terms, Model_Block *ModelBlock, int Curr_Block) const;
virtual double eval(const eval_context_type &eval_context) const throw (EvalException);
- virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx) const;
+ virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_type &temporary_terms, map_idx_type &map_idx, bool dynamic) const;
virtual NodeID toStatic(DataTree &static_datatree) const;
- virtual pair<bool, NodeID> normalizeLinearInEndoEquation(int symb_id_endo, NodeID Derivative) const;
+ virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const;
virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
};
diff --git a/Makefile.in b/Makefile.in
index a8d2af7c9cbef075de40ad50f6dea22c21d39810..6dc04bddf198c19e034955fd8002ccade599ca0c 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -14,6 +14,7 @@ MAIN_OBJS = \
ComputingTasks.o \
ModelTree.o \
StaticModel.o \
+ StaticDllModel.o \
DynamicModel.o \
NumericalConstants.o \
NumericalInitialization.o \
diff --git a/MinimumFeedbackSet.cc b/MinimumFeedbackSet.cc
index 0be2e17a9bcf0fab4994203aef3e1cc92d41a1d8..d9d50ee20542afd05306198f7588c0b929667d0e 100644
--- a/MinimumFeedbackSet.cc
+++ b/MinimumFeedbackSet.cc
@@ -382,6 +382,7 @@ namespace MFS
{
bool something_has_been_done = true;
int cut_ = 0;
+ feed_back_vertices.clear();
AdjacencyList_type G(G1);
while (num_vertices(G) > 0)
{
diff --git a/ModFile.cc b/ModFile.cc
index 138f3b1b665946b37c19046fd0eb5d4879a01e54..c5bb7914828ae47430f6ce92bfa9c5e385ed460f 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -25,6 +25,7 @@
ModFile::ModFile() : expressions_tree(symbol_table, num_constants),
static_model(symbol_table, num_constants),
+ static_dll_model(symbol_table, num_constants),
dynamic_model(symbol_table, num_constants),
linear(false)
{
@@ -143,9 +144,16 @@ ModFile::computingPass(bool no_tmp_terms)
if (dynamic_model.equation_number() > 0)
{
// Compute static model and its derivatives
- dynamic_model.toStatic(static_model);
- static_model.computingPass(mod_file_struct.steady_block_mfs_option, false, no_tmp_terms);
-
+ if(mod_file_struct.steady_block_mfs_dll_option)
+ {
+ dynamic_model.toStaticDll(static_dll_model);
+ static_dll_model.computingPass(global_eval_context, no_tmp_terms);
+ }
+ else
+ {
+ dynamic_model.toStatic(static_model);
+ static_model.computingPass(mod_file_struct.steady_block_mfs_option, false, no_tmp_terms);
+ }
// Set things to compute for dynamic model
if (mod_file_struct.simul_present)
@@ -228,7 +236,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
if (dynamic_model.equation_number() > 0)
{
dynamic_model.writeOutput(mOutputFile, basename);
- static_model.writeOutput(mOutputFile);
+ if(mod_file_struct.steady_block_mfs_dll_option)
+ static_dll_model.writeOutput(mOutputFile, basename);
+ else
+ static_model.writeOutput(mOutputFile);
}
// Print statements
@@ -248,7 +259,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all) const
// Create static and dynamic files
if (dynamic_model.equation_number() > 0)
{
- static_model.writeStaticFile(basename);
+ if(mod_file_struct.steady_block_mfs_dll_option)
+ static_dll_model.writeStaticFile(basename);
+ else
+ static_model.writeStaticFile(basename);
dynamic_model.writeDynamicFile(basename);
dynamic_model.writeParamsDerivativesFile(basename);
}
diff --git a/ModFile.hh b/ModFile.hh
index b43f408fbb36fe5064c30fc7eca30300d6282ed1..c7f4647eca871e2bcd0da161710b82176e4db2f4 100644
--- a/ModFile.hh
+++ b/ModFile.hh
@@ -29,6 +29,7 @@ using namespace std;
#include "NumericalConstants.hh"
#include "NumericalInitialization.hh"
#include "StaticModel.hh"
+#include "StaticDllModel.hh"
#include "DynamicModel.hh"
#include "Statement.hh"
@@ -46,6 +47,8 @@ public:
DataTree expressions_tree;
//! Static model
StaticModel static_model;
+ //! Static Dll model
+ StaticDllModel static_dll_model;
//! Dynamic model
DynamicModel dynamic_model;
//! Option linear
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index c3641a3308b166c5ef2135d70573b2af900c735b..fd8a2258351b24afd2fca0a51dd7e992777e3842 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -586,7 +586,7 @@ ParsingDriver::add_to_row(NodeID v)
void
ParsingDriver::steady()
{
- mod_file->addStatement(new SteadyStatement(options_list));
+ mod_file->addStatement(new SteadyStatement(options_list, mod_file->static_dll_model.mode));
options_list.clear();
}
@@ -620,7 +620,17 @@ ParsingDriver::option_num(const string &name_option, const string &opt)
if ((name_option == "periods") && (mod_file->dynamic_model.mode == DynamicModel::eSparseDLLMode || mod_file->dynamic_model.mode == DynamicModel::eSparseMode))
mod_file->dynamic_model.block_triangular.periods = atoi(opt.c_str());
else if (name_option == "cutoff")
- mod_file->dynamic_model.cutoff = atof(opt.c_str());
+ {
+ mod_file->dynamic_model.cutoff = atof(opt.c_str());
+ mod_file->static_dll_model.cutoff = atof(opt.c_str());
+ }
+ else if (name_option == "mfs")
+ {
+ mod_file->dynamic_model.mfs = atoi(opt.c_str());
+ mod_file->static_dll_model.mfs = atoi(opt.c_str());
+ }
+ else if (name_option == "block_mfs_dll")
+ mod_file->static_dll_model.mode = (StaticDllModel::mode_t)DynamicModel::eSparseDLLMode;
options_list.num_options[name_option] = opt;
}
diff --git a/Statement.cc b/Statement.cc
index 4350af575e931385a51ece1f1d12cb4cd6add0f3..527b9af3f502009f4779d3743653155973949ed0 100644
--- a/Statement.cc
+++ b/Statement.cc
@@ -31,7 +31,8 @@ ModFileStructure::ModFileStructure() :
bvar_density_present(false),
bvar_forecast_present(false),
identification_present(false),
- steady_block_mfs_option(false)
+ steady_block_mfs_option(false),
+ steady_block_mfs_dll_option(false)
{
}
diff --git a/Statement.hh b/Statement.hh
index 85f7686e9b8ad6e9859b32d27c175e382b8759ea..fd519c52cf125e5d343fb4a6c8cf9ca8e97b538f 100644
--- a/Statement.hh
+++ b/Statement.hh
@@ -62,6 +62,8 @@ public:
bool identification_present;
//! Whether the option "block_mfs" is used on steady statement
bool steady_block_mfs_option;
+ //! Whether the option "block_mfs_dll" is used on steady statement
+ bool steady_block_mfs_dll_option;
};
class Statement