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

Various simplifications and modernizations of the block/bytecode part

issue#70
Sébastien Villemot 2018-11-23 17:19:59 +01:00
parent fc9cc2dc50
commit d733f3fb5a
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 702 additions and 796 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

970
src/DynamicModel.cc
View File

File diff suppressed because it is too large Load Diff

44
src/DynamicModel.hh
View File

@ -99,8 +99,8 @@ private:
vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
//! Store the derivatives or the chainrule derivatives:map<pair< equation, pair< variable, lead_lag >, expr_t>
using first_chain_rule_derivatives_t = map< pair< int, pair< int, int>>, expr_t>;
//! Store the derivatives or the chainrule derivatives:map<tuple<equation, variable, lead_lag>, expr_t>
using first_chain_rule_derivatives_t = map<tuple<int, int, int>, expr_t>;
first_chain_rule_derivatives_t first_chain_rule_derivatives;
//! Writes dynamic model file (Matlab version)
@ -135,7 +135,7 @@ private:
//void evaluateJacobian(const eval_context_t &eval_context, jacob_map *j_m, bool dynamic);
//! return a map on the block jacobian
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> get_Derivatives(int block);
map<tuple<int, int, int, int, int>, int> get_Derivatives(int block);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
@ -163,7 +163,7 @@ private:
//! Computes derivatives of the Jacobian w.r. to trend vars and tests that they are equal to zero
void testTrendDerivativesEqualToZero(const eval_context_t &eval_context);
//! Collect only the first derivatives
map<pair<int, pair<int, int>>, expr_t> collect_first_order_derivatives_endogenous();
map<tuple<int, int, int>, expr_t> collect_first_order_derivatives_endogenous();
//! Allocates the derivation IDs for all dynamic variables of the model
/*! Also computes max_{endo,exo}_{lead_lag}, and initializes dynJacobianColsNbr to the number of dynamic endos */
@ -200,8 +200,8 @@ private:
//! Vector indicating if the block is linear in endogenous variable (true) or not (false)
vector<bool> blocks_linear;
//! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
using derivative_t = map<pair< int, pair<int, int>>, expr_t>;
//! Map the derivatives for a block tuple<lag, eq, var>
using derivative_t = map<tuple<int, int, int>, expr_t>;
//! Vector of derivative for each blocks
vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
@ -216,8 +216,8 @@ private:
map< int, map<int, int>> block_exo_index, block_det_exo_index, block_other_endo_index;
//! for each block described the number of static, forward, backward and mixed variables in the block
/*! pair< pair<static, forward>, pair<backward,mixed>> */
vector<pair< pair<int, int>, pair<int, int>>> block_col_type;
/*! tuple<static, forward, backward, mixed> */
vector<tuple<int, int, int, int>> block_col_type;
//! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
void computeRevXref(map<pair<int, int>, set<int>> &xrefset, const set<pair<int, int>> &eiref, int eqn);
@ -226,7 +226,7 @@ private:
void writeRevXrefs(ostream &output, const map<pair<int, int>, set<int>> &xrefmap, const string &type) const;
//! List for each variable its block number and its maximum lag and lead inside the block
vector<pair<int, pair<int, int>>> variable_block_lead_lag;
vector<tuple<int, int, int>> variable_block_lead_lag;
//! List for each equation its block number
vector<int> equation_block;
@ -481,19 +481,19 @@ public:
BlockSimulationType
getBlockSimulationType(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].first.first);
return (get<0>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the first equation number of a block
unsigned int
getBlockFirstEquation(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].first.second);
return (get<1>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the size of the block block_number
unsigned int
getBlockSize(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].second.first);
return (get<2>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the number of exogenous variable in the block block_number
unsigned int
@ -511,7 +511,7 @@ public:
unsigned int
getBlockMfs(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].second.second);
return (get<3>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the maximum lag in a block
unsigned int
@ -529,37 +529,37 @@ public:
EquationType
getBlockEquationType(int block_number, int equation_number) const override
{
return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);
return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first);
};
//! Return true if the equation has been normalized
bool
isBlockEquationRenormalized(int block_number, int equation_number) const override
{
return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == E_EVALUATE_S);
};
//! Return the expr_t of the equation equation_number belonging to the block block_number
expr_t
getBlockEquationExpr(int block_number, int equation_number) const override
{
return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
return (equations[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]]);
};
//! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
expr_t
getBlockEquationRenormalizedExpr(int block_number, int equation_number) const override
{
return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].second);
};
//! Return the original number of equation equation_number belonging to the block block_number
int
getBlockEquationID(int block_number, int equation_number) const override
{
return (equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);
return (equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]);
};
//! Return the original number of variable variable_number belonging to the block block_number
int
getBlockVariableID(int block_number, int variable_number) const override
{
return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
return (variable_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+variable_number]);
};
//! Return the original number of the exogenous variable varexo_number belonging to the block block_number
int
@ -572,19 +572,19 @@ public:
int
getBlockInitialEquationID(int block_number, int equation_number) const override
{
return ((int) inv_equation_reordered[equation_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
return ((int) inv_equation_reordered[equation_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the position of variable_number in the block number belonging to the block block_number
int
getBlockInitialVariableID(int block_number, int variable_number) const override
{
return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
return ((int) inv_variable_reordered[variable_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the block number containing the endogenous variable variable_number
int
getBlockVariableID(int variable_number) const
{
return (variable_block_lead_lag[variable_number].first);
return (get<0>(variable_block_lead_lag[variable_number]));
};
//! Return the position of the exogenous variable_number in the block number belonging to the block block_number
int

191
src/ModelTree.cc
View File

@ -185,7 +185,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
#if 1
bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
#else // Alternative way to compute normalization, by giving an initial matching using natural normalizations
fill(mate_map.begin(), mate_map.end(), graph_traits<BipartiteGraph>::null_vertex());
fill(mate_map.begin(), mate_map.end(), boost::graph_traits<BipartiteGraph>::null_vertex());
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
@ -201,15 +201,14 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
put(&mate_map[0], n+j, i);
}
edmonds_augmenting_path_finder<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type> augmentor(g, &mate_map[0], get(vertex_index, g));
bool not_maximum_yet = true;
while (not_maximum_yet)
boost::edmonds_augmenting_path_finder<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type> augmentor(g, &mate_map[0], get(boost::vertex_index, g));
while (augmentor.augment_matching())
{
not_maximum_yet = augmentor.augment_matching();
}
};
augmentor.get_current_matching(&mate_map[0]);
bool check = maximum_cardinality_matching_verifier<BipartiteGraph, size_t *, property_map<BipartiteGraph, vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(vertex_index, g));
bool check = boost::maximum_cardinality_matching_verifier<BipartiteGraph, int *, boost::property_map<BipartiteGraph, boost::vertex_index_t>::type>::verify_matching(g, &mate_map[0], get(boost::vertex_index, g));
#endif
assert(check);
@ -222,7 +221,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
endo2eq.resize(equations.size());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](int i) { return i-n; });
#ifdef DEBUG
multimap<int, int> natural_endo2eqs;
@ -237,9 +236,9 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
n1++;
pair<multimap<int, int>::const_iterator, multimap<int, int>::const_iterator> x = natural_endo2eqs.equal_range(i);
if (find_if(x.first, x.second, compose1(bind2nd(equal_to<int>(), endo2eq[i]), select2nd<multimap<int, int>::value_type>())) == x.second)
cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
auto x = natural_endo2eqs.equal_range(i);
if (find_if(x.first, x.second, [=](auto y) { return y.second == endo2eq[i]; }) == x.second)
cout << "Natural normalization of variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, i))
<< " not used." << endl;
else
n2++;
@ -274,9 +273,9 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
//jacob_map normalized_contemporaneous_jacobian;
jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
vector<double> max_val(n, 0.0);
for (jacob_map_t::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++)
if (fabs(iter->second) > max_val[iter->first.first])
max_val[iter->first.first] = fabs(iter->second);
for (const auto &it : contemporaneous_jacobian)
if (fabs(it.second) > max_val[it.first.first])
max_val[it.first.first] = fabs(it.second);
for (auto & iter : normalized_contemporaneous_jacobian)
iter.second /= max_val[iter.first.first];
@ -324,22 +323,22 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
if (check)
{
// Update the jacobian matrix
for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
for (const auto &it : tmp_normalized_contemporaneous_jacobian)
{
if (static_jacobian.find({ it->first.first, it->first.second }) == static_jacobian.end())
static_jacobian[{ it->first.first, it->first.second }] = 0;
if (dynamic_jacobian.find({ 0, { it->first.first, it->first.second } }) == dynamic_jacobian.end())
dynamic_jacobian[{ 0, { it->first.first, it->first.second } }] = nullptr;
if (contemporaneous_jacobian.find({ it->first.first, it->first.second }) == contemporaneous_jacobian.end())
contemporaneous_jacobian[{ it->first.first, it->first.second }] = 0;
if (static_jacobian.find({ it.first.first, it.first.second }) == static_jacobian.end())
static_jacobian[{ it.first.first, it.first.second }] = 0;
if (dynamic_jacobian.find({ 0, it.first.first, it.first.second }) == dynamic_jacobian.end())
dynamic_jacobian[{ 0, it.first.first, it.first.second }] = nullptr;
if (contemporaneous_jacobian.find({ it.first.first, it.first.second }) == contemporaneous_jacobian.end())
contemporaneous_jacobian[{ it.first.first, it.first.second }] = 0;
try
{
if (derivatives[1].find({ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }) == derivatives[1].end())
derivatives[1][{ it->first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it->first.second), 0) }] = Zero;
if (derivatives[1].find({ it.first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it.first.second), 0) }) == derivatives[1].end())
derivatives[1][{ it.first.first, getDerivID(symbol_table.getID(SymbolType::endogenous, it.first.second), 0) }] = Zero;
}
catch (DataTree::UnknownDerivIDException &e)
{
cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it->first.second))
cerr << "The variable " << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it.first.second))
<< " does not appear at the current period (i.e. with no lead and no lag); this case is not handled by the 'block' option of the 'model' block." << endl;
exit(EXIT_FAILURE);
}
@ -382,14 +381,13 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
{
int nb_elements_contemparenous_Jacobian = 0;
set<vector<int>> jacobian_elements_to_delete;
for (auto it = derivatives[1].begin();
it != derivatives[1].end(); it++)
for (const auto &it : derivatives[1])
{
int deriv_id = it->first[1];
int deriv_id = it.first[1];
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
{
expr_t Id = it->second;
int eq = it->first[0];
expr_t Id = it.second;
int eq = it.first[0];
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
@ -426,7 +424,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
static_jacobian[{ eq, var }] += val;
else
static_jacobian[{ eq, var }] = val;
dynamic_jacobian[{ lag, { eq, var } }] = Id;
dynamic_jacobian[{ lag, eq, var }] = Id;
}
}
}
@ -452,17 +450,17 @@ ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_li
/*equation_reordered.resize(equations.size());
variable_reordered.resize(equations.size());*/
unsigned int num = 0;
for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++)
if (!is_equation_linear[*it])
for (auto it : endo2eq)
if (!is_equation_linear[it])
num++;
vector<int> endo2block = vector<int>(endo2eq.size(), 1);
vector<pair<set<int>, pair<set<int>, vector<int> > > > components_set(num);
int i = 0, j = 0;
for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, j++)
for (auto it : endo2eq)
{
if (!is_equation_linear[*it])
if (!is_equation_linear[it])
{
equation_reordered[i] = *it;
equation_reordered[i] = it;
variable_reordered[i] = j;
endo2block[j] = 0;
components_set[endo2block[j]].first.insert(i);
@ -477,6 +475,7 @@ ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_li
cout.flush();
*/
i++;
j++;
}
}
/* for (unsigned int j = 0; j < is_equation_linear.size() ; j++)
@ -501,9 +500,9 @@ ModelTree::select_non_linear_equations_and_variables(vector<bool> is_equation_li
}
cout.flush();
int nb_endo = is_equation_linear.size();
vector<pair<int, int> > blocks = vector<pair<int, int> >(1, make_pair(i, i));
inv_equation_reordered = vector<int>(nb_endo);
inv_variable_reordered = vector<int>(nb_endo);
vector<pair<int, int>> blocks(1, make_pair(i, i));
inv_equation_reordered.resize(nb_endo);
inv_variable_reordered.resize(nb_endo);
for (int i = 0; i < nb_endo; i++)
{
inv_variable_reordered[variable_reordered[i]] = i;
@ -528,7 +527,7 @@ ModelTree::computeNaturalNormalization()
if (result.size() == 1 && result.begin()->second == 0)
{
//check if the endogenous variable has not been already used in an other match !
vector<int>::iterator it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
auto it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first);
if (it == endo2eq.end())
endo2eq[result.begin()->first] = eq;
else
@ -547,15 +546,15 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, ve
vector<int> eq2endo(equations.size(), 0);
equation_reordered.resize(equations.size());
variable_reordered.resize(equations.size());
bool *IM;
int n = equations.size();
IM = (bool *) calloc(n*n, sizeof(bool));
vector<bool> IM(n*n);
int i = 0;
for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
for (auto it : endo2eq)
{
eq2endo[*it] = i;
eq2endo[it] = i;
equation_reordered[i] = i;
variable_reordered[*it] = i;
variable_reordered[it] = i;
i++;
}
if (cutoff == 0)
{
@ -657,11 +656,10 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, ve
}
epilogue = tmp_epilogue;
}
free(IM);
}
equation_type_and_normalized_equation_t
ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
{
expr_t lhs;
BinaryOpNode *eq_node;
@ -674,7 +672,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t
eq_node = equations[eq];
lhs = eq_node->get_arg1();
Equation_Simulation_Type = E_SOLVE;
auto derivative = first_order_endo_derivatives.find({ eq, { var, 0 } });
auto derivative = first_order_endo_derivatives.find({ eq, var, 0 });
pair<bool, expr_t> res;
if (derivative != first_order_endo_derivatives.end())
{
@ -683,9 +681,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t
auto d_endo_variable = result.find({ var, 0 });
//Determine whether the equation could be evaluated rather than to be solved
if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
{
Equation_Simulation_Type = E_EVALUATE;
}
Equation_Simulation_Type = E_EVALUATE;
else
{
vector<pair<int, pair<expr_t, expr_t>>> List_of_Op_RHS;
@ -704,7 +700,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t
}
V_Equation_Simulation_Type[eq] = { Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second) };
}
return (V_Equation_Simulation_Type);
return V_Equation_Simulation_Type;
}
void
@ -734,9 +730,8 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
}
for (const auto & it : dynamic_jacobian)
{
int lag = it.first.first;
int j_1 = it.first.second.first;
int i_1 = it.first.second.second;
int lag, j_1, i_1;
tie(lag, j_1, i_1) = it.first;
if (variable_blck[i_1] == equation_blck[j_1])
{
if (lag > variable_lead_lag[i_1].second)
@ -786,12 +781,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
}
else
tmp_normalized_contemporaneous_jacobian = static_jacobian;
for (jacob_map_t::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
if (reverse_equation_reordered[it->first.first] >= (int) prologue && reverse_equation_reordered[it->first.first] < (int) (nb_var - epilogue)
&& reverse_variable_reordered[it->first.second] >= (int) prologue && reverse_variable_reordered[it->first.second] < (int) (nb_var - epilogue)
&& it->first.first != endo2eq[it->first.second])
add_edge(vertex(reverse_equation_reordered[endo2eq[it->first.second]]-prologue, G2),
vertex(reverse_equation_reordered[it->first.first]-prologue, G2),
for (const auto &it : tmp_normalized_contemporaneous_jacobian)
if (reverse_equation_reordered[it.first.first] >= (int) prologue && reverse_equation_reordered[it.first.first] < (int) (nb_var - epilogue)
&& reverse_variable_reordered[it.first.second] >= (int) prologue && reverse_variable_reordered[it.first.second] < (int) (nb_var - epilogue)
&& it.first.first != endo2eq[it.first.second])
add_edge(vertex(reverse_equation_reordered[endo2eq[it.first.second]]-prologue, G2),
vertex(reverse_equation_reordered[it.first.first]-prologue, G2),
G2);
vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
@ -832,12 +827,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
// - first set = equations belonging to the block,
// - second set = the feeback variables,
// - third vector = the reordered non-feedback variables.
vector<pair<set<int>, pair<set<int>, vector<int>>>> components_set(num);
vector<tuple<set<int>, set<int>, vector<int>>> components_set(num);
for (unsigned int i = 0; i < endo2block.size(); i++)
{
endo2block[i] = unordered2ordered[endo2block[i]];
blocks[endo2block[i]].first++;
components_set[endo2block[i]].first.insert(i);
get<0>(components_set[endo2block[i]]).insert(i);
}
getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
@ -885,12 +880,12 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
for (int i = 0; i < num; i++)
{
AdjacencyList_t G = extract_subgraph(G2, components_set[i].first);
AdjacencyList_t G = extract_subgraph(G2, get<0>(components_set[i]));
set<int> feed_back_vertices;
//Print(G);
AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
auto v_index = get(boost::vertex_index, G);
components_set[i].second.first = feed_back_vertices;
get<1>(components_set[i]) = feed_back_vertices;
blocks[i].second = feed_back_vertices.size();
vector<int> Reordered_Vertice;
Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
@ -898,7 +893,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
//First we have the recursive equations conditional on feedback variables
for (int j = 0; j < 4; j++)
{
for (int & its : Reordered_Vertice)
for (int its : Reordered_Vertice)
{
bool something_done = false;
if (j == 2 && variable_lag_lead[tmp_variable_reordered[its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[its +prologue]].second != 0)
@ -929,7 +924,7 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
}
}
}
components_set[i].second.second = Reordered_Vertice;
get<2>(components_set[i]) = Reordered_Vertice;
//Second we have the equations related to the feedback variables
for (int j = 0; j < 4; j++)
{
@ -978,8 +973,8 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
n_static[prologue+num+i]++;
}
inv_equation_reordered = vector<int>(nb_var);
inv_variable_reordered = vector<int>(nb_var);
inv_equation_reordered.resize(nb_var);
inv_variable_reordered.resize(nb_var);
for (int i = 0; i < nb_var; i++)
{
inv_variable_reordered[variable_reordered[i]] = i;
@ -1017,7 +1012,7 @@ ModelTree::printBlockDecomposition(const vector<pair<int, int>> &blocks) const
}
block_type_firstequation_size_mfs_t
ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int>>> &block_col_type, bool linear_decomposition)
ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<tuple<int, int, int, int>> &block_col_type, bool linear_decomposition)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
@ -1063,7 +1058,7 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
auto it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
if (linear_decomposition)
{
if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end())
if (dynamic_jacobian.find({ curr_lag, equation_reordered[count_equ], curr_variable }) != dynamic_jacobian.end())
{
if (curr_lag > Lead)
Lead = curr_lag;
@ -1074,7 +1069,7 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
else
{
if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end())
if (dynamic_jacobian.find({ curr_lag, equation_reordered[count_equ], curr_variable }) != dynamic_jacobian.end())
{
if (curr_lag > Lead)
Lead = curr_lag;
@ -1121,17 +1116,17 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
if (i > 0)
{
bool is_lead = false, is_lag = false;
int c_Size = (block_type_size_mfs[block_type_size_mfs.size()-1]).second.first;
int first_equation = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.second;
int c_Size = get<2>(block_type_size_mfs[block_type_size_mfs.size()-1]);
int first_equation = get<1>(block_type_size_mfs[block_type_size_mfs.size()-1]);
if (c_Size > 0 && ((prev_Type == EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
|| (prev_Type == EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag)))
{
for (int j = first_equation; j < first_equation+c_Size; j++)
{
auto it = dynamic_jacobian.find({ -1, { equation_reordered[eq], variable_reordered[j] } });
auto it = dynamic_jacobian.find({ -1, equation_reordered[eq], variable_reordered[j] });
if (it != dynamic_jacobian.end())
is_lag = true;
it = dynamic_jacobian.find({ +1, { equation_reordered[eq], variable_reordered[j] } });
it = dynamic_jacobian.find({ +1, equation_reordered[eq], variable_reordered[j] });
if (it != dynamic_jacobian.end())
is_lead = true;
}
@ -1140,55 +1135,55 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
|| (prev_Type == EVALUATE_BACKWARD && Simulation_Type == EVALUATE_BACKWARD && !is_lag))
{
//merge the current block with the previous one
BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
BlockSimulationType c_Type = get<0>(block_type_size_mfs[block_type_size_mfs.size()-1]);
c_Size++;
block_type_size_mfs[block_type_size_mfs.size()-1] = { { c_Type, first_equation }, { c_Size, c_Size } };
block_type_size_mfs[block_type_size_mfs.size()-1] = { c_Type, first_equation, c_Size, c_Size };
if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
block_lag_lead[block_type_size_mfs.size()-1] = { Lag, Lead };
pair< pair< unsigned int, unsigned int>, pair<unsigned int, unsigned int>> tmp = block_col_type[block_col_type.size()-1];
block_col_type[block_col_type.size()-1] = { { tmp.first.first+l_n_static, tmp.first.second+l_n_forward }, { tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed } };
auto tmp = block_col_type[block_col_type.size()-1];
block_col_type[block_col_type.size()-1] = { get<0>(tmp)+l_n_static, get<1>(tmp)+l_n_forward, get<2>(tmp)+l_n_backward, get<3>(tmp)+l_n_mixed };
}
else
{
block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
}
else
{
block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
}
else
{
block_type_size_mfs.emplace_back(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size));
block_type_size_mfs.emplace_back(Simulation_Type, eq, Blck_Size, MFS_Size);
block_lag_lead.emplace_back(Lag, Lead);
block_col_type.emplace_back(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed));
block_col_type.emplace_back(l_n_static, l_n_forward, l_n_backward, l_n_mixed);
}
prev_Type = Simulation_Type;
eq += Blck_Size;
}
return (block_type_size_mfs);
return block_type_size_mfs;
}
vector<bool>
ModelTree::equationLinear(map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives) const
ModelTree::equationLinear(map<tuple<int, int, int>, expr_t> first_order_endo_derivatives) const
{
vector<bool> is_linear(symbol_table.endo_nbr(), true);
for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_order_endo_derivatives.begin(); it != first_order_endo_derivatives.end(); it++)
for (const auto &it : first_order_endo_derivatives)
{
expr_t Id = it->second;
set<pair<int, int> > endogenous;
expr_t Id = it.second;
set<pair<int, int>> endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
{
int eq = it->first.first;
int eq = get<0>(it.first);
is_linear[eq] = false;
}
}
@ -1208,12 +1203,12 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
int first_variable_position = getBlockFirstEquation(block);
if (simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
{
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
for (const auto &it : derivatives_block)
{
int lag = it->second.first;
int lag = get<2>(it);
if (lag == 0)
{
expr_t Id = it->second.second;
expr_t Id = get<3>(it);
set<pair<int, int>> endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
@ -1232,10 +1227,10 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
}
else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
{
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
for (const auto &it : derivatives_block)
{
int lag = it->second.first;
expr_t Id = it->second.second; //
int lag = get<2>(it);
expr_t Id = get<3>(it); //
set<pair<int, int>> endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
@ -1254,7 +1249,7 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
the_end:
;
}
return (blocks_linear);
return blocks_linear;
}
int

18
src/ModelTree.hh
View File

@ -46,16 +46,16 @@ auto vectorToTuple(const vector<T>& v) {
}
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t >>;
using equation_type_and_normalized_equation_t = vector<pair<EquationType, expr_t>>;
//! Vector describing variables: max_lag in the block, max_lead in the block
using lag_lead_vector_t = vector<pair< int, int>>;
//! for each block contains pair< pair<Simulation_Type, first_equation>, pair < Block_Size, Recursive_part_Size >>
using block_type_firstequation_size_mfs_t = vector<pair< pair< BlockSimulationType, int>, pair<int, int>>>;
//! for each block contains tuple<Simulation_Type, first_equation, Block_Size, Recursive_part_Size>
using block_type_firstequation_size_mfs_t = vector<tuple<BlockSimulationType, int, int, int>>;
//! for a block contains derivatives pair< pair<block_equation_number, block_variable_number> , pair<lead_lag, expr_t>>
using block_derivatives_equation_variable_laglead_nodeid_t = vector< pair<pair<int, int>, pair< int, expr_t >>>;
//! for a block contains derivatives tuple<block_equation_number, block_variable_number, lead_lag, expr_t>
using block_derivatives_equation_variable_laglead_nodeid_t = vector<tuple<int, int, int, expr_t>>;
//! for all blocks derivatives description
using blocks_derivatives_t = vector<block_derivatives_equation_variable_laglead_nodeid_t>;
@ -194,7 +194,7 @@ protected:
//! Sparse matrix of double to store the values of the Jacobian
/*! First index is lag, second index is equation number, third index is endogenous type specific ID */
using dynamic_jacob_map_t = map<pair<int, pair<int, int>>, expr_t>;
using dynamic_jacob_map_t = map<tuple<int, int, int>, expr_t>;
//! Normalization of equations
/*! Maps endogenous type specific IDs to equation numbers */
@ -234,15 +234,15 @@ protected:
//! Search the equations and variables belonging to the prologue and the epilogue of the model
void computePrologueAndEpilogue(const jacob_map_t &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered);
//! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
equation_type_and_normalized_equation_t equationTypeDetermination(const map<pair<int, pair<int, int>>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
equation_type_and_normalized_equation_t equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
//! Compute the block decomposition and for a non-recusive block find the minimum feedback set
void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const;
//! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int, int>>> &block_col_type, bool linear_decomposition);
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int>> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<tuple<int, int, int, int>> &block_col_type, bool linear_decomposition);
//! Determine the maximum number of lead and lag for the endogenous variable in a bloc
void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const;
//! For each equation determine if it is linear or not
vector<bool> equationLinear(map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives) const;
vector<bool> equationLinear(map<tuple<int, int, int>, expr_t> first_order_endo_derivatives) const;
//! Print an abstract of the block structure of the model
void printBlockDecomposition(const vector<pair<int, int>> &blocks) const;
//! Determine for each block if it is linear or not

231
src/StaticModel.cc
View File

@ -61,7 +61,7 @@ StaticModel::copyHelper(const StaticModel &m)
{
block_derivatives_equation_variable_laglead_nodeid_t v;
for (const auto &it2 : it)
v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
v.emplace_back(get<0>(it2), get<1>(it2), get<2>(it2), f(get<3>(it2)));
blocks_derivatives.push_back(v);
}
@ -100,12 +100,7 @@ StaticModel::StaticModel(const StaticModel &m) :
global_temporary_terms {m.global_temporary_terms},
block_type_firstequation_size_mfs {m.block_type_firstequation_size_mfs},
blocks_linear {m.blocks_linear},
other_endo_block {m.other_endo_block},
exo_block {m.exo_block},
exo_det_block {m.exo_det_block},
block_col_type {m.block_col_type},
variable_block_lead_lag {m.variable_block_lead_lag},
equation_block {m.equation_block},
endo_max_leadlag_block {m.endo_max_leadlag_block},
other_endo_max_leadlag_block {m.other_endo_max_leadlag_block},
exo_max_leadlag_block {m.exo_max_leadlag_block},
@ -145,12 +140,7 @@ StaticModel::operator=(const StaticModel &m)
derivative_exo.clear();
derivative_exo_det.clear();
other_endo_block = m.other_endo_block;
exo_block = m.exo_block;
exo_det_block = m.exo_det_block;
block_col_type = m.block_col_type;
variable_block_lead_lag = m.variable_block_lead_lag;
equation_block = m.equation_block;
endo_max_leadlag_block = m.endo_max_leadlag_block;
other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
exo_max_leadlag_block = m.exo_max_leadlag_block;
@ -231,7 +221,7 @@ StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_nu
void
StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
{
auto it = first_chain_rule_derivatives.find({ eqr, { varr, lag } });
auto it = first_chain_rule_derivatives.find({ eqr, varr, lag });
if (it != first_chain_rule_derivatives.end())
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
@ -249,7 +239,6 @@ StaticModel::computeTemporaryTermsOrdered()
BinaryOpNode *eq_node;
first_chain_rule_derivatives_t::const_iterator it_chr;
ostringstream tmp_s;
v_temporary_terms.clear();
map_idx.clear();
unsigned int nb_blocks = getNbBlocks();
@ -266,11 +255,8 @@ StaticModel::computeTemporaryTermsOrdered()
for (unsigned int block = 0; block < nb_blocks; block++)
{
map<expr_t, int> reference_count_local;
reference_count_local.clear();
map<expr_t, pair<int, int>> first_occurence_local;
first_occurence_local.clear();
temporary_terms_t temporary_terms_l;
temporary_terms_l.clear();
unsigned int block_size = getBlockSize(block);
unsigned int block_nb_mfs = getBlockMfs(block);
@ -287,14 +273,12 @@ StaticModel::computeTemporaryTermsOrdered()
eq_node->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, i);
}
}
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
for (const auto &it : blocks_derivatives[block])
{
expr_t id = it->second.second;
expr_t id = get<3>(it);
id->computeTemporaryTerms(reference_count_local, temporary_terms_l, first_occurence_local, block, v_temporary_terms_local, block_size-1);
}
set<int> temporary_terms_in_use;
temporary_terms_in_use.clear();
v_temporary_terms_inuse[block] = temporary_terms_in_use;
v_temporary_terms_inuse[block] = {};
computeTemporaryTermsMapping(temporary_terms_l, map_idx2[block]);
}
@ -316,9 +300,9 @@ StaticModel::computeTemporaryTermsOrdered()
eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
}
}
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
for (const auto &it : blocks_derivatives[block])
{
expr_t id = it->second.second;
expr_t id = get<3>(it);
id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
}
}
@ -340,15 +324,14 @@ StaticModel::computeTemporaryTermsOrdered()
eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
}
}
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
for (const auto &it : blocks_derivatives[block])
{
expr_t id = it->second.second;
expr_t id = get<3>(it);
id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
}
for (int i = 0; i < (int) getBlockSize(block); i++)
for (auto it = v_temporary_terms[block][i].begin();
it != v_temporary_terms[block][i].end(); it++)
(*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
for (const auto &it : v_temporary_terms[block][i])
it->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
v_temporary_terms_inuse[block] = temporary_terms_in_use;
}
computeTemporaryTermsMapping(temporary_terms, map_idx);
@ -396,16 +379,16 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
tmp1_output.str("");
tmp1_output << packageDir(basename + ".block") << "/static_" << block+1 << ".m";
output.open(tmp1_output.str(), ios::out | ios::binary);
output << "%\n";
output << "% " << tmp1_output.str() << " : Computes static model for Dynare\n";
output << "%\n";
output << "% Warning : this file is generated automatically by Dynare\n";
output << "% from model file (.mod)\n\n";
output << "%/\n";
output << "%" << endl
<< "% " << tmp1_output.str() << " : Computes static model for Dynare" << endl
<< "%" << endl
<< "% Warning : this file is generated automatically by Dynare" << endl
<< "% from model file (.mod)" << endl << endl
<< "%/" << endl;
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
output << "function y = static_" << block+1 << "(y, x, params)\n";
output << "function y = static_" << block+1 << "(y, x, params)" << endl;
else
output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)\n";
output << "function [residual, y, g1] = static_" << block+1 << "(y, x, params)" << endl;
BlockType block_type;
if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE)
@ -436,11 +419,11 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
tmp_output.str("");
for (int it : v_temporary_terms_inuse[block])
tmp_output << " T" << it;
output << " global" << tmp_output.str() << ";\n";
output << " global" << tmp_output.str() << ";" << endl;
}
if (simulation_type != EVALUATE_BACKWARD && simulation_type != EVALUATE_FORWARD)
output << " residual=zeros(" << block_mfs << ",1);\n";
output << " residual=zeros(" << block_mfs << ",1);" << endl;
// The equations
temporary_terms_idxs_t temporary_terms_idxs;
@ -449,7 +432,6 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
if (!global_temporary_terms)
local_temporary_terms = v_temporary_terms[block][i];
temporary_terms_t tt2;
tt2.clear();
if (v_temporary_terms[block].size())
{
output << " " << "% //Temporary variables" << endl;
@ -498,7 +480,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
if (isBlockEquationRenormalized(block, i))
{
rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
output << "\n ";
output << endl << " ";
tmp_output.str("");
eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
lhs = eq_node->get_arg1();
@ -510,10 +492,10 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
}
else
{
cerr << "Type mismatch for equation " << equation_ID+1 << "\n";
cerr << "Type mismatch for equation " << equation_ID+1 << endl;
exit(EXIT_FAILURE);
}
output << ";\n";
output << ";" << endl;
break;
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
@ -531,7 +513,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
output << tmp_output.str();
output << ") - (";
rhs->writeOutput(output, local_output_type, local_temporary_terms, {});
output << ");\n";
output << ");" << endl;
}
}
// The Jacobian if we have to solve the block
@ -544,13 +526,13 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
case SOLVE_FORWARD_SIMPLE:
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
for (const auto &it : blocks_derivatives[block])
{
unsigned int eq = it->first.first;
unsigned int var = it->first.second;
unsigned int eq, var;
expr_t id;
tie(eq, var, ignore, id) = it;
unsigned int eqr = getBlockEquationID(block, eq);
unsigned int varr = getBlockVariableID(block, var);
expr_t id = it->second.second;
output << " g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms, {});
output << "; % variable=" << symbol_table.getName(symbol_table.getID(SymbolType::endogenous, varr))
@ -626,8 +608,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
fjmp_if_eval.write(code_file, instruction_number);
int prev_instruction_number = instruction_number;
vector<vector<pair<int, int>>> my_derivatives;
my_derivatives.resize(symbol_table.endo_nbr());
vector<vector<pair<int, int>>> my_derivatives(symbol_table.endo_nbr());
count_u = symbol_table.endo_nbr();
for (const auto & first_derivative : derivatives[1])
{
@ -657,8 +638,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
fldr.write(code_file, instruction_number);
if (my_derivatives[i].size())
{
for (vector<pair<int, int>>::const_iterator it = my_derivatives[i].begin();
it != my_derivatives[i].end(); it++)
for (auto it = my_derivatives[i].begin(); it != my_derivatives[i].end(); it++)
{
FLDSU_ fldsu(it->second);
fldsu.write(code_file, instruction_number);
@ -690,10 +670,7 @@ StaticModel::writeModelEquationsCode(const string &basename, map_idx_t map_idx)
code_file.seekp(pos3);
prev_instruction_number = instruction_number;
temporary_terms_t tt2;
tt2.clear();
temporary_terms_t tt3;
tt3.clear();
temporary_terms_t tt2, tt3;
// The Jacobian if we have to solve the block determinsitic bloc
for (const auto & first_derivative : derivatives[1])
@ -820,7 +797,6 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
{
//The Temporary terms
temporary_terms_t tt2;
tt2.clear();
if (v_temporary_terms[block].size())
{
for (auto it : v_temporary_terms[block][i])
@ -919,10 +895,10 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
count_u = feedback_variables.size();
for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
for (const auto &it : blocks_derivatives[block])
{
unsigned int eq = it->first.first;
unsigned int var = it->first.second;
unsigned int eq, var;
tie(eq, var, ignore, ignore) = it;
unsigned int eqr = getBlockEquationID(block, eq);
unsigned int varr = getBlockVariableID(block, var);
if (eq >= block_recursive && var >= block_recursive)
@ -1005,32 +981,28 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
code_file.seekp(pos3);
prev_instruction_number = instruction_number;
temporary_terms_t tt2;
tt2.clear();
temporary_terms_t tt3;
tt3.clear();
temporary_terms_t tt2, tt3;
deriv_node_temp_terms_t tef_terms2;
for (i = 0; i < (int) block_size; i++)
{
if (v_temporary_terms_local[block].size())
{
for (auto it = v_temporary_terms_local[block][i].begin();
it != v_temporary_terms_local[block][i].end(); it++)
for (const auto &it : v_temporary_terms_local[block][i])
{
if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != nullptr)
(*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
it->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find((*it)->idx)->second));
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find(it->idx)->second));
fnumexpr.write(code_file, instruction_number);
(*it)->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
it->compile(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms);
FSTPST_ fstpst((int)(map_idx2[block].find((*it)->idx)->second));
FSTPST_ fstpst((int)(map_idx2[block].find(it->idx)->second));
fstpst.write(code_file, instruction_number);
// Insert current node into tt2
tt3.insert(*it);
tt2.insert(*it);
tt3.insert(it);
tt2.insert(it);
}
}
@ -1113,10 +1085,10 @@ StaticModel::writeModelEquationsCode_Block(const string &basename, map_idx_t map
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
count_u = feedback_variables.size();
for (auto it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
for (const auto &it : blocks_derivatives[block])
{
unsigned int eq = it->first.first;
unsigned int var = it->first.second;
unsigned int eq, var;
tie(eq, var, ignore, ignore) = it;
unsigned int eqr = getBlockEquationID(block, eq);
unsigned int varr = getBlockVariableID(block, var);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, 0);
@ -1165,10 +1137,10 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num,
unsigned int block_size = getBlockSize(num);
unsigned int block_mfs = getBlockMfs(num);
unsigned int block_recursive = block_size - block_mfs;
for (auto it = blocks_derivatives[num].begin(); it != (blocks_derivatives[num]).end(); it++)
for (const auto &it : blocks_derivatives[num])
{
unsigned int eq = it->first.first;
unsigned int var = it->first.second;
unsigned int eq, var;
tie(eq, var, ignore, ignore) = it;
int lag = 0;
if (eq >= block_recursive && var >= block_recursive)
{
@ -1196,10 +1168,10 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num,
SaveCode.close();
}
map<pair<int, pair<int, int >>, expr_t>
map<tuple<int, int, int>, expr_t>
StaticModel::collect_first_order_derivatives_endogenous()
{
map<pair<int, pair<int, int >>, expr_t> endo_derivatives;
map<tuple<int, int, int>, expr_t> endo_derivatives;
for (auto & first_derivative : derivatives[1])
{
if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous)
@ -1207,7 +1179,7 @@ StaticModel::collect_first_order_derivatives_endogenous()
int eq = first_derivative.first[0];
int var = symbol_table.getTypeSpecificID(getSymbIDByDerivID(first_derivative.first[1]));
int lag = 0;
endo_derivatives[{ eq, { var, lag } }] = first_derivative.second;
endo_derivatives[{ eq, var, lag }] = first_derivative.second;
}
}
return endo_derivatives;
@ -1270,12 +1242,12 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co
computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
map<pair<int, pair<int, int>>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
map<tuple<int, int, int>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
if (!nopreprocessoroutput)
cout << "Finding the optimal block decomposition of the model ...\n";
cout << "Finding the optimal block decomposition of the model ..." << endl;
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
@ -2209,7 +2181,7 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
output << "function [residual, g1, y, var_index] = static(nblock, y, x, params)" << endl
<< " residual = [];" << endl
<< " g1 = [];" << endl
<< " var_index = [];\n" << endl
<< " var_index = [];" << endl << endl
<< " switch nblock" << endl;
unsigned int nb_blocks = getNbBlocks();
@ -2225,16 +2197,16 @@ StaticModel::writeStaticBlockMFSFile(const string &basename) const
if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD)
{
output << " y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);\n";
output << " y_tmp = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
ostringstream tmp;
for (int i = 0; i < (int) getBlockSize(b); i++)
tmp << " " << getBlockVariableID(b, i)+1;
output << " var_index = [" << tmp.str() << "];\n";
output << " residual = y(var_index) - y_tmp(var_index);\n";
output << " y = y_tmp;\n";
output << " var_index = [" << tmp.str() << "];" << endl
<< " residual = y(var_index) - y_tmp(var_index);" << endl
<< " y = y_tmp;" << endl;
}
else
output << " [residual, y, g1] = " << basename << ".block.static_" << b+1 << "(y, x, params);\n";
output << " [residual, y, g1] = " << basename << ".block.static_" << b+1 << "(y, x, params);" << endl;
}
output << " end" << endl
@ -2267,23 +2239,23 @@ StaticModel::writeOutput(ostream &output, bool block) const
tmp_s << " " << getBlockVariableID(b, i)+1;
tmp_s_eq << " " << getBlockEquationID(b, i)+1;
}
output << "block_structure_stat.block(" << b+1 << ").Simulation_Type = " << simulation_type << ";\n";
output << "block_structure_stat.block(" << b+1 << ").endo_nbr = " << block_size << ";\n";
output << "block_structure_stat.block(" << b+1 << ").mfs = " << getBlockMfs(block) << ";\n";
output << "block_structure_stat.block(" << b+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
output << "block_structure_stat.block(" << b+1 << ").variable = [" << tmp_s.str() << "];\n";
output << "block_structure_stat.block(" << b+1 << ").Simulation_Type = " << simulation_type << ";" << endl
<< "block_structure_stat.block(" << b+1 << ").endo_nbr = " << block_size << ";" << endl
<< "block_structure_stat.block(" << b+1 << ").mfs = " << getBlockMfs(block) << ";" << endl
<< "block_structure_stat.block(" << b+1 << ").equation = [" << tmp_s_eq.str() << "];" << endl
<< "block_structure_stat.block(" << b+1 << ").variable = [" << tmp_s.str() << "];" << endl;
}
output << "M_.block_structure_stat.block = block_structure_stat.block;\n";
output << "M_.block_structure_stat.block = block_structure_stat.block;" << endl;
string cst_s;
int nb_endo = symbol_table.endo_nbr();
output << "M_.block_structure_stat.variable_reordered = [";
for (int i = 0; i < nb_endo; i++)
output << " " << variable_reordered[i]+1;
output << "];\n";
output << "M_.block_structure_stat.equation_reordered = [";
output << "];" << endl
<< "M_.block_structure_stat.equation_reordered = [";
for (int i = 0; i < nb_endo; i++)
output << " " << equation_reordered[i]+1;
output << "];\n";
output << "];" << endl;
map<pair<int, int>, int> row_incidence;
for (const auto & first_derivative : derivatives[1])
@ -2299,11 +2271,9 @@ StaticModel::writeOutput(ostream &output, bool block) const
}
}
output << "M_.block_structure_stat.incidence.sparse_IM = [";
for (map<pair< int, int >, int>::const_iterator it = row_incidence.begin(); it != row_incidence.end(); it++)
{
output << it->first.first+1 << " " << it->first.second+1 << ";\n";
}
output << "];\n";
for (const auto &it : row_incidence)
output << it.first.first+1 << " " << it.first.second+1 << ";" << endl;
output << "];" << endl;
}
SymbolType
@ -2352,11 +2322,10 @@ StaticModel::addAllParamDerivId(set<int> &deriv_id_set)
deriv_id_set.insert(i + symbol_table.endo_nbr());
}
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>
map<tuple<int, int, int, int, int>, int>
StaticModel::get_Derivatives(int block)
{
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> Derivatives;
Derivatives.clear();
map<tuple<int, int, int, int, int>, int> Derivatives;
int block_size = getBlockSize(block);
int block_nb_recursive = block_size - getBlockMfs(block);
int lag = 0;
@ -2366,10 +2335,10 @@ StaticModel::get_Derivatives(int block)
for (int var = 0; var < block_size; var++)
{
int varr = getBlockVariableID(block, var);
if (dynamic_jacobian.find({ lag, { eqr, varr } }) != dynamic_jacobian.end())
if (dynamic_jacobian.find({ lag, eqr, varr }) != dynamic_jacobian.end())
{
bool OK = true;
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>::const_iterator its = Derivatives.find({ { lag, { eq, var } }, { eqr, varr } });
auto its = Derivatives.find({ lag, eq, var, eqr, varr });
if (its != Derivatives.end())
{
if (its->second == 2)
@ -2380,10 +2349,10 @@ StaticModel::get_Derivatives(int block)
{
if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursive)
//It's a normalized equation, we have to recompute the derivative using chain rule derivative function
Derivatives[{ { lag, { eq, var } }, { eqr, varr } }] = 1;
Derivatives[{ lag, eq, var, eqr, varr }] = 1;
else
//It's a feedback equation we can use the derivatives
Derivatives[{ { lag, { eq, var } }, { eqr, varr } }] = 0;
Derivatives[{ lag, eq, var, eqr, varr }] = 0;
}
if (var < block_nb_recursive)
{
@ -2392,15 +2361,15 @@ StaticModel::get_Derivatives(int block)
{
int varrs = getBlockVariableID(block, vars);
//A new derivative needs to be computed using the chain rule derivative function (a feedback variable appears in a recursive equation)
if (Derivatives.find({ { lag, { var, vars } }, { eqs, varrs } }) != Derivatives.end())
Derivatives[{ { lag, { eq, vars } }, { eqr, varrs } }] = 2;
if (Derivatives.find({ lag, var, vars, eqs, varrs }) != Derivatives.end())
Derivatives[{ lag, eq, vars, eqr, varrs }] = 2;
}
}
}
}
}
return (Derivatives);
return Derivatives;
}
void
@ -2427,30 +2396,24 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
else
recursive_variables[getDerivID(symbol_table.getID(SymbolType::endogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
}
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> Derivatives = get_Derivatives(block);
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int>::const_iterator it = Derivatives.begin();
for (int i = 0; i < (int) Derivatives.size(); i++)
auto Derivatives = get_Derivatives(block);
for (const auto &it : Derivatives)
{
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;
int Deriv_type = it.second;
int lag, eq, var, eqr, varr;
tie(lag, eq, var, eqr, varr) = it.first;
if (Deriv_type == 0)
first_chain_rule_derivatives[{ eqr, { varr, lag } }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
first_chain_rule_derivatives[{ eqr, varr, lag }] = derivatives[1][{ eqr, getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag) }];
else if (Deriv_type == 1)
first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
{
if (getBlockEquationType(block, eq) == E_EVALUATE_S && eq < block_nb_recursives)
first_chain_rule_derivatives[{ eqr, { varr, lag } }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
first_chain_rule_derivatives[{ eqr, varr, lag }] = (equation_type_and_normalized_equation[eqr].second)->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
else
first_chain_rule_derivatives[{ eqr, { varr, lag } }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
first_chain_rule_derivatives[{ eqr, varr, lag }] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), lag), recursive_variables);
}
tmp_derivatives.emplace_back(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))]));
tmp_derivatives.emplace_back(eq, var, lag, first_chain_rule_derivatives[{ eqr, varr, lag }]);
}
}
else
@ -2472,8 +2435,8 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(SymbolType::endogenous, varr), 0), recursive_variables);
if (d1 == Zero)
continue;
first_chain_rule_derivatives[{ eqr, { varr, 0 } }] = d1;
tmp_derivatives.emplace_back(make_pair(eq, var), make_pair(0, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))]));
first_chain_rule_derivatives[{ eqr, varr, 0 }] = d1;
tmp_derivatives.emplace_back(eq, var, 0, first_chain_rule_derivatives[{ eqr, varr, 0 }]);
}
}
}
@ -2485,10 +2448,8 @@ void
StaticModel::collect_block_first_order_derivatives()
{
//! vector for an equation or a variable indicates the block number
vector<int> equation_2_block, variable_2_block;
vector<int> equation_2_block(equation_reordered.size()), variable_2_block(variable_reordered.size());
unsigned int nb_blocks = getNbBlocks();
equation_2_block = vector<int>(equation_reordered.size());
variable_2_block = vector<int>(variable_reordered.size());
for (unsigned int block = 0; block < nb_blocks; block++)
{
unsigned int block_size = getBlockSize(block);
@ -2517,7 +2478,7 @@ StaticModel::collect_block_first_order_derivatives()
if (getTypeByDerivID(first_derivative.first[1]) == SymbolType::endogenous && block_eq == block_var)
{
tmp_derivative = derivative_endo[block_eq];
tmp_derivative[{ lag, { eq, var } }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
tmp_derivative[{ lag, eq, var }] = derivatives[1][{ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, var), lag) }];
derivative_endo[block_eq] = tmp_derivative;
}
}

44
src/StaticModel.hh
View File

@ -42,7 +42,7 @@ private:
vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
using first_chain_rule_derivatives_t = map< pair< int, pair< int, int>>, expr_t>;
using first_chain_rule_derivatives_t = map<tuple<int, int, int>, expr_t>;
first_chain_rule_derivatives_t first_chain_rule_derivatives;
//! Writes static model file (standard Matlab version)
@ -99,11 +99,11 @@ private:
//! Compute the column indices of the static Jacobian
void computeStatJacobianCols();
//! return a map on the block jacobian
map<pair<pair<int, pair<int, int>>, pair<int, int>>, int> get_Derivatives(int block);
map<tuple<int, int, int, int, int>, int> get_Derivatives(int block);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
//! Collect only the first derivatives
map<pair<int, pair<int, int>>, expr_t> collect_first_order_derivatives_endogenous();
map<tuple<int, int, int>, expr_t> collect_first_order_derivatives_endogenous();
//! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
void collect_block_first_order_derivatives();
@ -126,24 +126,18 @@ private:
//! Vector indicating if the block is linear in endogenous variable (true) or not (false)
vector<bool> blocks_linear;
//! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
using derivative_t = map<pair< int, pair<int, int>>, expr_t>;
//! Map the derivatives for a block tuple<lag, eq, var>
using derivative_t = map<tuple<int, int, int>, expr_t>;
//! Vector of derivative for each blocks
vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
//!List for each block and for each lag-leag all the other endogenous variables and exogenous variables
using var_t = set<int>;
using lag_var_t = map<int, var_t>;
vector<lag_var_t> other_endo_block, exo_block, exo_det_block;
//! for each block described the number of static, forward, backward and mixed variables in the block
/*! pair< pair<static, forward>, pair<backward,mixed>> */
vector<pair< pair<int, int>, pair<int, int>>> block_col_type;
//! List for each variable its block number and its maximum lag and lead inside the block
vector<pair<int, pair<int, int>>> variable_block_lead_lag;
//! List for each equation its block number
vector<int> equation_block;
/*! tuple<static, forward, backward, mixed> */
vector<tuple<int, int, int, int>> block_col_type;
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector<pair<int, int>> endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
@ -246,19 +240,19 @@ public:
BlockSimulationType
getBlockSimulationType(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].first.first);
return (get<0>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the first equation number of a block
unsigned int
getBlockFirstEquation(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].first.second);
return (get<1>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the size of the block block_number
unsigned int
getBlockSize(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].second.first);
return (get<2>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the number of exogenous variable in the block block_number
unsigned int
@ -276,7 +270,7 @@ public:
unsigned int
getBlockMfs(int block_number) const override
{
return (block_type_firstequation_size_mfs[block_number].second.second);
return (get<3>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the maximum lag in a block
unsigned int
@ -294,37 +288,37 @@ public:
EquationType
getBlockEquationType(int block_number, int equation_number) const override
{
return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);
return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first);
};
//! Return true if the equation has been normalized
bool
isBlockEquationRenormalized(int block_number, int equation_number) const override
{
return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].first == E_EVALUATE_S);
};
//! Return the expr_t of the equation equation_number belonging to the block block_number
expr_t
getBlockEquationExpr(int block_number, int equation_number) const override
{
return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
return (equations[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]]);
};
//! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
expr_t
getBlockEquationRenormalizedExpr(int block_number, int equation_number) const override
{
return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
return (equation_type_and_normalized_equation[equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]].second);
};
//! Return the original number of equation equation_number belonging to the block block_number
int
getBlockEquationID(int block_number, int equation_number) const override
{
return (equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);
return (equation_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+equation_number]);
};
//! Return the original number of variable variable_number belonging to the block block_number
int
getBlockVariableID(int block_number, int variable_number) const override
{
return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
return (variable_reordered[get<1>(block_type_firstequation_size_mfs[block_number])+variable_number]);
};
//! Return the original number of the exogenous variable varexo_number belonging to the block block_number
int
@ -336,13 +330,13 @@ public:
int
getBlockInitialEquationID(int block_number, int equation_number) const override
{
return ((int) inv_equation_reordered[equation_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
return ((int) inv_equation_reordered[equation_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the position of variable_number in the block number belonging to the block block_number
int
getBlockInitialVariableID(int block_number, int variable_number) const override
{
return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
return ((int) inv_variable_reordered[variable_number] - (int) get<1>(block_type_firstequation_size_mfs[block_number]));
};
//! Return the position of variable_number in the block number belonging to the block block_number
int