Select Git revision
ConfigFile.hh
Forked from
Dynare / preprocessor
Source project has a limited visibility.
-
Houtan Bastani authoredHoutan Bastani authored
ModelTree.cc NaN GiB
/*
* Copyright (C) 2003-2013 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cstdlib>
#include <cassert>
#include <cmath>
#include <iostream>
#include <fstream>
#include "ModelTree.hh"
#include "MinimumFeedbackSet.hh"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <boost/graph/strong_components.hpp>
#include <boost/graph/topological_sort.hpp>
using namespace boost;
using namespace MFS;
bool
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
{
const int n = equation_number();
assert(n == symbol_table.endo_nbr());
typedef adjacency_list<vecS, vecS, undirectedS> BipartiteGraph;
/*
Vertices 0 to n-1 are for endogenous (using type specific ID)
Vertices n to 2*n-1 are for equations (using equation no.)
*/
BipartiteGraph g(2 * n);
// Fill in the graph
set<pair<int, int> > endo;
for (jacob_map_t::const_iterator it = contemporaneous_jacobian.begin(); it != contemporaneous_jacobian.end(); it++)
add_edge(it->first.first + n, it->first.second, g);
// Compute maximum cardinality matching
vector<int> mate_map(2*n);
#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());
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
int j = natural_endo2eqs.find(i)->second;
put(&mate_map[0], i, n+j);
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)
{
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));
#endif
assert(check);
#ifdef DEBUG
for (int i = 0; i < n; i++)
cout << "Endogenous " << symbol_table.getName(symbol_table.getID(eEndogenous, i))
<< " matched with equation " << (mate_map[i]-n+1) << endl;
#endif
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
endo2eq.resize(equation_number());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus<int>(), n));
#ifdef DEBUG
multimap<int, int> natural_endo2eqs;
computeNormalizedEquations(natural_endo2eqs);
int n1 = 0, n2 = 0;
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
if (natural_endo2eqs.count(i) == 0)
continue;
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))
<< " not used." << endl;
else
n2++;
}
cout << "Used " << n2 << " natural normalizations out of " << n1 << ", for a total of " << n << " equations." << endl;
#endif
// Check if all variables are normalized
vector<int>::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
if (it != mate_map.begin() + n)
{
if (verbose)
cerr << "ERROR: Could not normalize the model. Variable "
<< symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
<< " is not in the maximum cardinality matching." << endl;
check = false;
}
return check;
}
void
ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian)
{
bool check = false;
cout << "Normalizing the model..." << endl;
int n = equation_number();
// compute the maximum value of each row of the contemporaneous Jacobian matrix
//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 (jacob_map_t::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
iter->second /= max_val[iter->first.first];
//We start with the highest value of the cutoff and try to normalize the model
double current_cutoff = 0.99999999;
int suppressed = 0;
while (!check && current_cutoff > 1e-19)
{
jacob_map_t tmp_normalized_contemporaneous_jacobian;
int suppress = 0;
for (jacob_map_t::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
if (fabs(iter->second) > max(current_cutoff, cutoff))
tmp_normalized_contemporaneous_jacobian[make_pair(iter->first.first, iter->first.second)] = iter->second;
else
suppress++;
if (suppress != suppressed)
check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
suppressed = suppress;
if (!check)
{
current_cutoff /= 2;
// In this last case try to normalize with the complete jacobian
if (current_cutoff <= 1e-19)
check = computeNormalization(normalized_contemporaneous_jacobian, false);
}
}
if (!check)
{
cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
//if no non-singular normalization can be found, try to find a normalization even with a potential singularity
jacob_map_t tmp_normalized_contemporaneous_jacobian;
set<pair<int, int> > endo;
for (int i = 0; i < n; i++)
{
endo.clear();
equations[i]->collectEndogenous(endo);
for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1;
}
check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
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++)
{
if (static_jacobian.find(make_pair(it->first.first, it->first.second)) == static_jacobian.end())
static_jacobian[make_pair(it->first.first, it->first.second)] = 0;
if (dynamic_jacobian.find(make_pair(0, make_pair(it->first.first, it->first.second))) == dynamic_jacobian.end())
dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
if (first_derivatives.find(make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))) == first_derivatives.end())
first_derivatives[make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))] = Zero;
}
}
}
if (!check)
{
cerr << "No normalization could be computed. Aborting." << endl;
exit(EXIT_FAILURE);
}
}
void
ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
{
for (int i = 0; i < equation_number(); i++)
{
VariableNode *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
if (lhs == NULL)
continue;
int symb_id = lhs->get_symb_id();
if (symbol_table.getType(symb_id) != eEndogenous)
continue;
set<pair<int, int> > endo;
equations[i]->get_arg2()->collectEndogenous(endo);
if (endo.find(make_pair(symbol_table.getTypeSpecificID(symb_id), 0)) != endo.end())
continue;
endo2eqs.insert(make_pair(symbol_table.getTypeSpecificID(symb_id), i));
cout << "Endogenous " << symbol_table.getName(symb_id) << " normalized in equation " << (i+1) << endl;
}
}
void
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose)
{
int nb_elements_contemparenous_Jacobian = 0;
set<pair<int, int> > jacobian_elements_to_delete;
for (first_derivatives_t::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int deriv_id = it->first.second;
if (getTypeByDerivID(deriv_id) == eEndogenous)
{
expr_t Id = it->second;
int eq = it->first.first;
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
double val = 0;
try
{
val = Id->eval(eval_context);
}
catch (ExprNode::EvalExternalFunctionException &e)
{
val = 1;
}
catch (ExprNode::EvalException &e)
{
cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1 << " and variable " << symbol_table.getName(symb) << "(" << lag << ") [" << symb << "] !" << endl;
Id->writeOutput(cerr, oMatlabDynamicModelSparse, temporary_terms);
cerr << endl;
exit(EXIT_FAILURE);
}
if (fabs(val) < cutoff)
{
if (verbose)
cout << "the coefficient related to variable " << var << " with lag " << lag << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix (size=" << symbol_table.endo_nbr() << ")" << endl;
jacobian_elements_to_delete.insert(make_pair(eq, deriv_id));
}
else
{
if (lag == 0)
{
nb_elements_contemparenous_Jacobian++;
contemporaneous_jacobian[make_pair(eq, var)] = val;
}
if (static_jacobian.find(make_pair(eq, var)) != static_jacobian.end())
static_jacobian[make_pair(eq, var)] += val;
else
static_jacobian[make_pair(eq, var)] = val;
dynamic_jacobian[make_pair(lag, make_pair(eq, var))] = Id;
}
}
}
// Get rid of the elements of the Jacobian matrix below the cutoff
for (set<pair<int, int> >::const_iterator it = jacobian_elements_to_delete.begin(); it != jacobian_elements_to_delete.end(); it++)
first_derivatives.erase(*it);
if (jacobian_elements_to_delete.size() > 0)
{
cout << jacobian_elements_to_delete.size() << " elements among " << first_derivatives.size() << " in the incidence matrices are below the cutoff (" << cutoff << ") and are discarded" << endl
<< "The contemporaneous incidence matrix has " << nb_elements_contemparenous_Jacobian << " elements" << endl;
}
}
void
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector<int> &equation_reordered, vector<int> &variable_reordered)
{
vector<int> eq2endo(equation_number(), 0);
equation_reordered.resize(equation_number());
variable_reordered.resize(equation_number());
bool *IM;
int n = equation_number();
IM = (bool *) calloc(n*n, sizeof(bool));
int i = 0;
for (vector<int>::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
{
eq2endo[*it] = i;
equation_reordered[i] = i;
variable_reordered[*it] = i;
}
for (jacob_map_t::const_iterator it = static_jacobian_arg.begin(); it != static_jacobian_arg.end(); it++)
IM[it->first.first * n + endo2eq[it->first.second]] = true;
bool something_has_been_done = true;
prologue = 0;
int k = 0;
// Find the prologue equations and place first the AR(1) shock equations first
while (something_has_been_done)
{
int tmp_prologue = prologue;
something_has_been_done = false;
for (int i = prologue; i < n; i++)
{
int nze = 0;
for (int j = tmp_prologue; j < n; j++)
if (IM[i * n + j])
{
nze++;
k = j;
}
if (nze == 1)
{
for (int j = 0; j < n; j++)
{
bool tmp_bool = IM[tmp_prologue * n + j];
IM[tmp_prologue * n + j] = IM[i * n + j];
IM[i * n + j] = tmp_bool;
}
int tmp = equation_reordered[tmp_prologue];
equation_reordered[tmp_prologue] = equation_reordered[i];
equation_reordered[i] = tmp;
for (int j = 0; j < n; j++)
{
bool tmp_bool = IM[j * n + tmp_prologue];
IM[j * n + tmp_prologue] = IM[j * n + k];
IM[j * n + k] = tmp_bool;
}
tmp = variable_reordered[tmp_prologue];
variable_reordered[tmp_prologue] = variable_reordered[k];
variable_reordered[k] = tmp;
tmp_prologue++;
something_has_been_done = true;
}
}
prologue = tmp_prologue;
}
something_has_been_done = true;
epilogue = 0;
// Find the epilogue equations
while (something_has_been_done)
{
int tmp_epilogue = epilogue;
something_has_been_done = false;
for (int i = prologue; i < n - (int) epilogue; i++)
{
int nze = 0;
for (int j = prologue; j < n - tmp_epilogue; j++)
if (IM[j * n + i])
{
nze++;
k = j;
}
if (nze == 1)
{
for (int j = 0; j < n; j++)
{
bool tmp_bool = IM[(n - 1 - tmp_epilogue) * n + j];
IM[(n - 1 - tmp_epilogue) * n + j] = IM[k * n + j];
IM[k * n + j] = tmp_bool;
}
int tmp = equation_reordered[n - 1 - tmp_epilogue];
equation_reordered[n - 1 - tmp_epilogue] = equation_reordered[k];
equation_reordered[k] = tmp;
for (int j = 0; j < n; j++)
{
bool tmp_bool = IM[j * n + n - 1 - tmp_epilogue];
IM[j * n + n - 1 - tmp_epilogue] = IM[j * n + i];
IM[j * n + i] = tmp_bool;
}
tmp = variable_reordered[n - 1 - tmp_epilogue];
variable_reordered[n - 1 - tmp_epilogue] = variable_reordered[i];
variable_reordered[i] = tmp;
tmp_epilogue++;
something_has_been_done = true;
}
}
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
{
expr_t lhs;
BinaryOpNode *eq_node;
EquationType Equation_Simulation_Type;
equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
for (unsigned int i = 0; i < equations.size(); i++)
{
int eq = Index_Equ_IM[i];
int var = Index_Var_IM[i];
eq_node = equations[eq];
lhs = eq_node->get_arg1();
Equation_Simulation_Type = E_SOLVE;
map<pair<int, pair<int, int> >, expr_t>::const_iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
pair<bool, expr_t> res;
if (derivative != first_order_endo_derivatives.end())
{
set<pair<int, int> > result;
derivative->second->collectEndogenous(result);
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 (lhs->isVariableNodeEqualTo(eEndogenous, Index_Var_IM[i], 0) && derivative->second->isNumConstNodeEqualTo(1))
{
Equation_Simulation_Type = E_EVALUATE;
}
else
{
vector<pair<int, pair<expr_t, expr_t> > > List_of_Op_RHS;
res = equations[eq]->normalizeEquation(var, List_of_Op_RHS);
if (mfs == 2)
{
if (d_endo_variable == result.end() && res.second)
Equation_Simulation_Type = E_EVALUATE_S;
}
else if (mfs == 3)
{
if (res.second) // The equation could be solved analytically
Equation_Simulation_Type = E_EVALUATE_S;
}
}
}
V_Equation_Simulation_Type[eq] = make_pair(Equation_Simulation_Type, dynamic_cast<BinaryOpNode *>(res.second));
}
return (V_Equation_Simulation_Type);
}
void
ModelTree::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
{
int nb_endo = symbol_table.endo_nbr();
variable_lead_lag = lag_lead_vector_t(nb_endo, make_pair(0, 0));
equation_lead_lag = lag_lead_vector_t(nb_endo, make_pair(0, 0));
vector<int> variable_blck(nb_endo), equation_blck(nb_endo);
for (int i = 0; i < nb_endo; i++)
{
if (i < (int) prologue)
{
variable_blck[variable_reordered[i]] = i;
equation_blck[equation_reordered[i]] = i;
}
else if (i < (int) (components_set.size() + prologue))
{
variable_blck[variable_reordered[i]] = components_set[i-prologue] + prologue;
equation_blck[equation_reordered[i]] = components_set[i-prologue] + prologue;
}
else
{
variable_blck[variable_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
equation_blck[equation_reordered[i]] = i- (nb_endo - nb_blck_sim - prologue - epilogue);
}
}
for (dynamic_jacob_map_t::const_iterator it = dynamic_jacobian.begin(); it != dynamic_jacobian.end(); it++)
{
int lag = it->first.first;
int j_1 = it->first.second.first;
int i_1 = it->first.second.second;
if (variable_blck[i_1] == equation_blck[j_1])
{
if (lag > variable_lead_lag[i_1].second)
variable_lead_lag[i_1] = make_pair(variable_lead_lag[i_1].first, lag);
if (lag < -variable_lead_lag[i_1].first)
variable_lead_lag[i_1] = make_pair(-lag, variable_lead_lag[i_1].second);
if (lag > equation_lead_lag[j_1].second)
equation_lead_lag[j_1] = make_pair(equation_lead_lag[j_1].first, lag);
if (lag < -equation_lead_lag[j_1].first)
equation_lead_lag[j_1] = make_pair(-lag, equation_lead_lag[j_1].second);
}
}
}
void
ModelTree::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, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const
{
int nb_var = variable_reordered.size();
int n = nb_var - prologue - epilogue;
AdjacencyList_t G2(n);
// It is necessary to manually initialize vertex_index property since this graph uses listS and not vecS as underlying vertex container
property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G2);
for (int i = 0; i < n; i++)
put(v_index, vertex(i, G2), i);
vector<int> reverse_equation_reordered(nb_var), reverse_variable_reordered(nb_var);
for (int i = 0; i < nb_var; i++)
{
reverse_equation_reordered[equation_reordered[i]] = i;
reverse_variable_reordered[variable_reordered[i]] = i;
}
for (jacob_map_t::const_iterator it = static_jacobian.begin(); it != static_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),
G2);
vector<int> endo2block(num_vertices(G2)), discover_time(num_vertices(G2));
iterator_property_map<int *, property_map<AdjacencyList_t, vertex_index_t>::type, int, int &> endo2block_map(&endo2block[0], get(vertex_index, G2));
// Compute strongly connected components
int num = strong_components(G2, endo2block_map);
blocks = vector<pair<int, int> >(num, make_pair(0, 0));
// Create directed acyclic graph associated to the strongly connected components
typedef adjacency_list<vecS, vecS, directedS> DirectedGraph;
DirectedGraph dag(num);
for (unsigned int i = 0; i < num_vertices(G2); i++)
{
AdjacencyList_t::out_edge_iterator it_out, out_end;
AdjacencyList_t::vertex_descriptor vi = vertex(i, G2);
for (tie(it_out, out_end) = out_edges(vi, G2); it_out != out_end; ++it_out)
{
int t_b = endo2block_map[target(*it_out, G2)];
int s_b = endo2block_map[source(*it_out, G2)];
if (s_b != t_b)
add_edge(s_b, t_b, dag);
}
}
// Compute topological sort of DAG (ordered list of unordered SCC)
deque<int> ordered2unordered;
topological_sort(dag, front_inserter(ordered2unordered)); // We use a front inserter because topological_sort returns the inverse order
// Construct mapping from unordered SCC to ordered SCC
vector<int> unordered2ordered(num);
for (int i = 0; i < num; i++)
unordered2ordered[ordered2unordered[i]] = i;
//This vector contains for each block:
// - first set = equations belonging to the block,
// - second set = the feeback variables,
// - third vector = the reordered non-feedback variables.
vector<pair<set<int>, pair<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);
}
getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
int order = prologue;
//Add a loop on vertices which could not be normalized or vertices related to lead variables => force those vertices to belong to the feedback set
if (select_feedback_variable)
{
for (int i = 0; i < n; i++)
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE
|| variable_lag_lead[variable_reordered[i+prologue]].second > 0
|| variable_lag_lead[variable_reordered[i+prologue]].first > 0
|| equation_lag_lead[equation_reordered[i+prologue]].second > 0
|| equation_lag_lead[equation_reordered[i+prologue]].first > 0
|| mfs == 0)
add_edge(vertex(i, G2), vertex(i, G2), G2);
}
else
{
for (int i = 0; i < n; i++)
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
add_edge(vertex(i, G2), vertex(i, G2), G2);
}
//Determines the dynamic structure of each equation
n_static = vector<unsigned int>(prologue+num+epilogue, 0);
n_forward = vector<unsigned int>(prologue+num+epilogue, 0);
n_backward = vector<unsigned int>(prologue+num+epilogue, 0);
n_mixed = vector<unsigned int>(prologue+num+epilogue, 0);
for (int i = 0; i < (int) prologue; i++)
{
if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
n_mixed[i]++;
else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
n_forward[i]++;
else if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
n_backward[i]++;
else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
n_static[i]++;
}
//For each block, the minimum set of feedback variable is computed
// and the non-feedback variables are reordered to get
// a sub-recursive block without feedback variables
for (int i = 0; i < num; i++)
{
AdjacencyList_t G = extract_subgraph(G2, components_set[i].first);
set<int> feed_back_vertices;
//Print(G);
AdjacencyList_t G1 = Minimal_set_of_feedback_vertex(feed_back_vertices, G);
property_map<AdjacencyList_t, vertex_index_t>::type v_index = get(vertex_index, G);
components_set[i].second.first = feed_back_vertices;
blocks[i].second = feed_back_vertices.size();
vector<int> Reordered_Vertice;
Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
//First we have the recursive equations conditional on feedback variables
for (int j = 0; j < 4; j++)
{
for (vector<int>::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
{
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)
{
n_mixed[prologue+i]++;
something_done = true;
}
else if (j == 3 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
{
n_forward[prologue+i]++;
something_done = true;
}
else if (j == 1 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
{
n_backward[prologue+i]++;
something_done = true;
}
else if (j == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
{
n_static[prologue+i]++;
something_done = true;
}
if (something_done)
{
equation_reordered[order] = tmp_equation_reordered[*its+prologue];
variable_reordered[order] = tmp_variable_reordered[*its+prologue];
order++;
}
}
}
components_set[i].second.second = Reordered_Vertice;
//Second we have the equations related to the feedback variables
for (int j = 0; j < 4; j++)
{
for (set<int>::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
{
bool something_done = false;
if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
{
n_mixed[prologue+i]++;
something_done = true;
}
else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
{
n_forward[prologue+i]++;
something_done = true;
}
else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
{
n_backward[prologue+i]++;
something_done = true;
}
else if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
{
n_static[prologue+i]++;
something_done = true;
}
if (something_done)
{
equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
order++;
}
}
}
}
for (int i = 0; i < (int) epilogue; i++)
{
if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
n_mixed[prologue+num+i]++;
else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second != 0)
n_forward[prologue+num+i]++;
else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
n_backward[prologue+num+i]++;
else if (variable_lag_lead[tmp_variable_reordered[prologue+n+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+n+i]].second == 0)
n_static[prologue+num+i]++;
}
inv_equation_reordered = vector<int>(nb_var);
inv_variable_reordered = vector<int>(nb_var);
for (int i = 0; i < nb_var; i++)
{
inv_variable_reordered[variable_reordered[i]] = i;
inv_equation_reordered[equation_reordered[i]] = i;
}
}
void
ModelTree::printBlockDecomposition(const vector<pair<int, int> > &blocks) const
{
int largest_block = 0;
int Nb_SimulBlocks = 0;
int Nb_feedback_variable = 0;
unsigned int Nb_TotalBlocks = getNbBlocks();
for (unsigned int block = 0; block < Nb_TotalBlocks; block++)
{
BlockSimulationType simulation_type = getBlockSimulationType(block);
if (simulation_type == SOLVE_FORWARD_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE)
{
Nb_SimulBlocks++;
int size = getBlockSize(block);
if (size > largest_block)
{
largest_block = size;
Nb_feedback_variable = getBlockMfs(block);
}
}
}
int Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
cout << Nb_TotalBlocks << " block(s) found:" << endl
<< " " << Nb_RecursBlocks << " recursive block(s) and " << Nb_SimulBlocks << " simultaneous block(s)." << endl
<< " the largest simultaneous block has " << largest_block << " equation(s)" << endl
<< " and " << Nb_feedback_variable << " feedback variable(s)." << endl;
}
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)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
int Blck_Size, MFS_Size;
int Lead, Lag;
block_type_firstequation_size_mfs_t block_type_size_mfs;
BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
int eq = 0;
unsigned int l_n_static = 0;
unsigned int l_n_forward = 0;
unsigned int l_n_backward = 0;
unsigned int l_n_mixed = 0;
for (i = 0; i < (int) (prologue+blocks.size()+epilogue); i++)
{
int first_count_equ = count_equ;
if (i < (int) prologue)
{
Blck_Size = 1;
MFS_Size = 1;
}
else if (i < (int) (prologue+blocks.size()))
{
Blck_Size = blocks[blck_count_simult].first;
MFS_Size = blocks[blck_count_simult].second;
blck_count_simult++;
}
else if (i < (int) (prologue+blocks.size()+epilogue))
{
Blck_Size = 1;
MFS_Size = 1;
}
Lag = Lead = 0;
set<pair<int, int> > endo;
for (count_equ = first_count_equ; count_equ < Blck_Size+first_count_equ; count_equ++)
{
endo.clear();
equations[equation_reordered[count_equ]]->collectEndogenous(endo);
for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
{
int curr_variable = it->first;
int curr_lag = it->second;
vector<int>::const_iterator it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable);
if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size))
if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end())
{
if (curr_lag > Lead)
Lead = curr_lag;
else if (-curr_lag > Lag)
Lag = -curr_lag;
}
}
}
if ((Lag > 0) && (Lead > 0))
{
if (Blck_Size == 1)
Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE;
else
Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE;
}
else if (Blck_Size > 1)
{
if (Lead > 0)
Simulation_Type = SOLVE_BACKWARD_COMPLETE;
else
Simulation_Type = SOLVE_FORWARD_COMPLETE;
}
else
{
if (Lead > 0)
Simulation_Type = SOLVE_BACKWARD_SIMPLE;
else
Simulation_Type = SOLVE_FORWARD_SIMPLE;
}
l_n_static = n_static[i];
l_n_forward = n_forward[i];
l_n_backward = n_backward[i];
l_n_mixed = n_mixed[i];
if (Blck_Size == 1)
{
if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
{
if (Simulation_Type == SOLVE_BACKWARD_SIMPLE)
Simulation_Type = EVALUATE_BACKWARD;
else if (Simulation_Type == SOLVE_FORWARD_SIMPLE)
Simulation_Type = EVALUATE_FORWARD;
}
if (i > 0)
{
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;
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++)
{
dynamic_jacob_map_t::const_iterator it = dynamic_jacobian.find(make_pair(-1, make_pair(equation_reordered[eq], variable_reordered[j])));
if (it != dynamic_jacobian.end())
is_lag = true;
it = dynamic_jacobian.find(make_pair(+1, make_pair(equation_reordered[eq], variable_reordered[j])));
if (it != dynamic_jacobian.end())
is_lead = true;
}
}
if ((prev_Type == EVALUATE_FORWARD && Simulation_Type == EVALUATE_FORWARD && !is_lead)
|| (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;
c_Size++;
block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(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] = make_pair(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] = make_pair(make_pair(tmp.first.first+l_n_static, tmp.first.second+l_n_forward), make_pair(tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed));
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
block_col_type.push_back(make_pair(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed)));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
block_col_type.push_back(make_pair(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed)));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
block_col_type.push_back(make_pair(make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed)));
}
prev_Type = Simulation_Type;
eq += Blck_Size;
}
return (block_type_size_mfs);
}
vector<bool>
ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector<int> &variable_reordered) const
{
unsigned int nb_blocks = getNbBlocks();
vector<bool> blocks_linear(nb_blocks, true);
for (unsigned int block = 0; block < nb_blocks; block++)
{
BlockSimulationType simulation_type = getBlockSimulationType(block);
int block_size = getBlockSize(block);
block_derivatives_equation_variable_laglead_nodeid_t derivatives_block = blocks_derivatives[block];
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++)
{
int lag = it->second.first;
if (lag == 0)
{
expr_t Id = it->second.second;
set<pair<int, int> > endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
{
for (int l = 0; l < block_size; l++)
{
if (endogenous.find(make_pair(variable_reordered[first_variable_position+l], 0)) != endogenous.end())
{
blocks_linear[block] = false;
goto the_end;
}
}
}
}
}
}
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++)
{
int lag = it->second.first;
expr_t Id = it->second.second; //
set<pair<int, int> > endogenous;
Id->collectEndogenous(endogenous);
if (endogenous.size() > 0)
{
for (int l = 0; l < block_size; l++)
{
if (endogenous.find(make_pair(variable_reordered[first_variable_position+l], lag)) != endogenous.end())
{
blocks_linear[block] = false;
goto the_end;
}
}
}
}
}
the_end:
;
}
return (blocks_linear);
}
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
NumericalConstants &num_constants_arg,
ExternalFunctionsTable &external_functions_table_arg) :
DataTree(symbol_table_arg, num_constants_arg, external_functions_table_arg),
cutoff(1e-15),
mfs(0)
{
for (int i = 0; i < 3; i++)
NNZDerivatives[i] = 0;
}
int
ModelTree::equation_number() const
{
return (equations.size());
}
void
ModelTree::writeDerivative(ostream &output, int eq, int symb_id, int lag,
ExprNodeOutputType output_type,
const temporary_terms_t &temporary_terms) const
{
first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symb_id, lag)));
if (it != first_derivatives.end())
(it->second)->writeOutput(output, output_type, temporary_terms);
else
output << 0;
}
void
ModelTree::computeJacobian(const set<int> &vars)
{
for (set<int>::const_iterator it = vars.begin();
it != vars.end(); it++)
{
for (int eq = 0; eq < (int) equations.size(); eq++)
{
expr_t d1 = equations[eq]->getDerivative(*it);
if (d1 == Zero)
continue;
first_derivatives[make_pair(eq, *it)] = d1;
++NNZDerivatives[0];
}
}
}
void
ModelTree::computeHessian(const set<int> &vars)
{
for (first_derivatives_t::const_iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
int eq = it->first.first;
int var1 = it->first.second;
expr_t d1 = it->second;
// Store only second derivatives with var2 <= var1
for (set<int>::const_iterator it2 = vars.begin();
it2 != vars.end(); it2++)
{
int var2 = *it2;
if (var2 > var1)
continue;
expr_t d2 = d1->getDerivative(var2);
if (d2 == Zero)
continue;
second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
if (var2 == var1)
++NNZDerivatives[1];
else
NNZDerivatives[1] += 2;
}
}
}
void
ModelTree::computeThirdDerivatives(const set<int> &vars)
{
for (second_derivatives_t::const_iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
{
int eq = it->first.first;
int var1 = it->first.second.first;
int var2 = it->first.second.second;
// By construction, var2 <= var1
expr_t d2 = it->second;
// Store only third derivatives such that var3 <= var2 <= var1
for (set<int>::const_iterator it2 = vars.begin();
it2 != vars.end(); it2++)
{
int var3 = *it2;
if (var3 > var2)
continue;
expr_t d3 = d2->getDerivative(var3);
if (d3 == Zero)
continue;
third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
if (var3 == var2 && var2 == var1)
++NNZDerivatives[2];
else if (var3 == var2 || var2 == var1)
NNZDerivatives[2] += 3;
else
NNZDerivatives[2] += 6;
}
}
}
void
ModelTree::computeTemporaryTerms(bool is_matlab)
{
map<expr_t, int> reference_count;
temporary_terms.clear();
for (vector<BinaryOpNode *>::iterator it = equations.begin();
it != equations.end(); it++)
(*it)->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
for (first_derivatives_t::iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
for (second_derivatives_t::iterator it = second_derivatives.begin();
it != second_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
for (third_derivatives_t::iterator it = third_derivatives.begin();
it != third_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab);
}
void
ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, ostream &output,
ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const
{
// Local var used to keep track of temp nodes already written
temporary_terms_t tt2;
for (temporary_terms_t::const_iterator it = tt.begin();
it != tt.end(); it++)
{
if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
(*it)->writeExternalFunctionOutput(output, output_type, tt2, tef_terms);
if (IS_C(output_type))
output << "double ";
(*it)->writeOutput(output, output_type, tt, tef_terms);
output << " = ";
(*it)->writeOutput(output, output_type, tt2, tef_terms);
if (IS_C(output_type))
output << ";" << endl;
// Insert current node into tt2
tt2.insert(*it);
if (IS_MATLAB(output_type))
output << ";" << endl;
}
}
void
ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
{
// Local var used to keep track of temp nodes already written
temporary_terms_t tt2;
// To store the functions that have already been written in the form TEF* = ext_fun();
deriv_node_temp_terms_t tef_terms;
for (temporary_terms_t::const_iterator it = tt.begin();
it != tt.end(); it++)
{
if (dynamic_cast<ExternalFunctionNode *>(*it) != NULL)
{
(*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
}
FNUMEXPR_ fnumexpr(TemporaryTerm, (int) (map_idx.find((*it)->idx)->second));
fnumexpr.write(code_file, instruction_number);
(*it)->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
if (dynamic)
{
FSTPT_ fstpt((int) (map_idx.find((*it)->idx)->second));
fstpt.write(code_file, instruction_number);
}
else
{
FSTPST_ fstpst((int) (map_idx.find((*it)->idx)->second));
fstpst.write(code_file, instruction_number);
}
// Insert current node into tt2
tt2.insert(*it);
}
}
void
ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const
{
/* Collect all model local variables appearing in equations, and print only
them. Printing unused model local variables can lead to a crash (see
ticket #101). */
set<int> used_local_vars;
// Use an empty set for the temporary terms
const temporary_terms_t tt;
for (size_t i = 0; i < equations.size(); i++)
equations[i]->collectModelLocalVariables(used_local_vars);
for (set<int>::const_iterator it = used_local_vars.begin();
it != used_local_vars.end(); ++it)
{
int id = *it;
expr_t value = local_variables_table.find(id)->second;
value->writeExternalFunctionOutput(output, output_type, tt, tef_terms);
if (IS_C(output_type))
output << "double ";
/* We append underscores to avoid name clashes with "g1" or "oo_" (see
also VariableNode::writeOutput) */
output << symbol_table.getName(id) << "__ = ";
value->writeOutput(output, output_type, tt, tef_terms);
output << ";" << endl;
}
}
void
ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const
{
for (int eq = 0; eq < (int) equations.size(); eq++)
{
BinaryOpNode *eq_node = equations[eq];
expr_t lhs = eq_node->get_arg1();
expr_t rhs = eq_node->get_arg2();
// Test if the right hand side of the equation is empty.
double vrhs = 1.0;
try
{
vrhs = rhs->eval(eval_context_t());
}
catch (ExprNode::EvalException &e)
{
}
if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
{
output << "lhs =";
lhs->writeOutput(output, output_type, temporary_terms);
output << ";" << endl;
output << "rhs =";
rhs->writeOutput(output, output_type, temporary_terms);
output << ";" << endl;
output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(output_type)
<< "= lhs-rhs;" << endl;
}
else // The right hand side of the equation is empty ==> residual=lhs;
{
output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(output_type)
<< " = ";
lhs->writeOutput(output, output_type, temporary_terms);
output << ";" << endl;
}
}
}
void
ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
for (int eq = 0; eq < (int) equations.size(); eq++)
{
BinaryOpNode *eq_node = equations[eq];
expr_t lhs = eq_node->get_arg1();
expr_t rhs = eq_node->get_arg2();
FNUMEXPR_ fnumexpr(ModelEquation, eq);
fnumexpr.write(code_file, instruction_number);
// Test if the right hand side of the equation is empty.
double vrhs = 1.0;
try
{
vrhs = rhs->eval(eval_context_t());
}
catch (ExprNode::EvalException &e)
{
}
if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
{
lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
FBINARY_ fbinary(oMinus);
fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(eq);
fstpr.write(code_file, instruction_number);
}
else // The right hand side of the equation is empty ==> residual=lhs;
{
lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
FSTPR_ fstpr(eq);
fstpr.write(code_file, instruction_number);
}
}
}
void
ModelTree::Write_Inf_To_Bin_File(const string &basename,
int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const
{
int j;
std::ofstream SaveCode;
const string bin_basename = basename + ".bin";
if (file_open)
SaveCode.open(bin_basename.c_str(), ios::out | ios::in | ios::binary | ios::ate);
else
SaveCode.open(bin_basename.c_str(), ios::out | ios::binary);
if (!SaveCode.is_open())
{
cout << "Error : Can't open file \"" << bin_basename << "\" for writing\n";
exit(EXIT_FAILURE);
}
u_count_int = 0;
for (first_derivatives_t::const_iterator it = first_derivatives.begin(); it != first_derivatives.end(); it++)
{
int deriv_id = it->first.second;
if (getTypeByDerivID(deriv_id) == eEndogenous)
{
int eq = it->first.first;
int symb = getSymbIDByDerivID(deriv_id);
int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
SaveCode.write(reinterpret_cast<char *>(&eq), sizeof(eq));
int varr = var + lag * block_mfs;
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&lag), sizeof(lag));
int u = u_count_int + block_mfs;
SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
u_count_int++;
}
}
if (is_two_boundaries)
u_count_int += symbol_table.endo_nbr();
for (j = 0; j < (int) symbol_table.endo_nbr(); j++)
SaveCode.write(reinterpret_cast<char *>(&j), sizeof(j));
for (j = 0; j < (int) symbol_table.endo_nbr(); j++)
SaveCode.write(reinterpret_cast<char *>(&j), sizeof(j));
SaveCode.close();
}
void
ModelTree::writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const
{
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "\\documentclass[10pt,a4paper]{article}" << endl
<< "\\usepackage[landscape]{geometry}" << endl
<< "\\usepackage{fullpage}" << endl
<< "\\usepackage{breqn}" << endl
<< "\\begin{document}" << endl
<< "\\footnotesize" << endl;
// Write model local variables
for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
it != local_variables_table.end(); it++)
{
int id = it->first;
expr_t value = it->second;
output << "\\begin{dmath*}" << endl
<< symbol_table.getName(id) << " = ";
// Use an empty set for the temporary terms
value->writeOutput(output, output_type);
output << endl << "\\end{dmath*}" << endl;
}
for (int eq = 0; eq < (int) equations.size(); eq++)
{
output << "\\begin{dmath}" << endl
<< "% Equation " << eq+1 << endl;
// Here it is necessary to cast to superclass ExprNode, otherwise the overloaded writeOutput() method is not found
dynamic_cast<ExprNode *>(equations[eq])->writeOutput(output, output_type);
output << endl << "\\end{dmath}" << endl;
}
output << "\\end{document}" << endl;
output.close();
}
void
ModelTree::addEquation(expr_t eq)
{
BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
assert(beq != NULL && beq->get_op_code() == oEqual);
equations.push_back(beq);
}
void
ModelTree::addEquation(expr_t eq, vector<pair<string, string> > &eq_tags)
{
int n = equation_number();
for (size_t i = 0; i < eq_tags.size(); i++)
equation_tags.push_back(make_pair(n, eq_tags[i]));
addEquation(eq);
}
void
ModelTree::addAuxEquation(expr_t eq)
{
BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
assert(beq != NULL && beq->get_op_code() == oEqual);
aux_equations.push_back(beq);
}
void
ModelTree::addTrendVariables(vector<int> trend_vars, expr_t growth_factor) throw (TrendException)
{
while (!trend_vars.empty())
if (trend_symbols_map.find(trend_vars.back()) != trend_symbols_map.end())
throw TrendException(symbol_table.getName(trend_vars.back()));
else
{
trend_symbols_map[trend_vars.back()] = growth_factor;
trend_vars.pop_back();
}
}
void
ModelTree::addNonstationaryVariables(vector<int> nonstationary_vars, bool log_deflator, expr_t deflator) throw (TrendException)
{
while (!nonstationary_vars.empty())
if (nonstationary_symbols_map.find(nonstationary_vars.back()) != nonstationary_symbols_map.end())
throw TrendException(symbol_table.getName(nonstationary_vars.back()));
else
{
nonstationary_symbols_map[nonstationary_vars.back()] = make_pair(log_deflator, deflator);
nonstationary_vars.pop_back();
}
}
void
ModelTree::initializeVariablesAndEquations()
{
for (int j = 0; j < equation_number(); j++)
{
equation_reordered.push_back(j);
variable_reordered.push_back(j);
}
}
void
ModelTree::set_cutoff_to_zero()
{
cutoff = 0;
}
void
ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
{
output << " g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
if (IS_MATLAB(output_type))
output << eq_nb + 1 << "," << col_nb + 1;
else
output << eq_nb + col_nb *equations.size();
output << RIGHT_ARRAY_SUBSCRIPT(output_type);
}
void
ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
{
output << " v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
if (IS_MATLAB(output_type))
output << row_nb + 1 << "," << col_nb + 1;
else
output << row_nb + col_nb * NNZDerivatives[order-1];
output << RIGHT_ARRAY_SUBSCRIPT(output_type);
}
void
ModelTree::computeParamsDerivatives()
{
set<int> deriv_id_set;
addAllParamDerivId(deriv_id_set);
for (set<int>::const_iterator it = deriv_id_set.begin();
it != deriv_id_set.end(); it++)
{
const int param = *it;
for (int eq = 0; eq < (int) equations.size(); eq++)
{
expr_t d1 = equations[eq]->getDerivative(param);
if (d1 == Zero)
continue;
residuals_params_derivatives[make_pair(eq, param)] = d1;
}
for (first_derivatives_t::const_iterator it2 = residuals_params_derivatives.begin();
it2 != residuals_params_derivatives.end(); it2++)
{
int eq = it2->first.first;
int param1 = it2->first.second;
expr_t d1 = it2->second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
residuals_params_second_derivatives[make_pair(eq, make_pair(param1, param))] = d2;
}
for (first_derivatives_t::const_iterator it2 = first_derivatives.begin();
it2 != first_derivatives.end(); it2++)
{
int eq = it2->first.first;
int var = it2->first.second;
expr_t d1 = it2->second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
jacobian_params_derivatives[make_pair(eq, make_pair(var, param))] = d2;
}
for (second_derivatives_t::const_iterator it2 = jacobian_params_derivatives.begin();
it2 != jacobian_params_derivatives.end(); it2++)
{
int eq = it2->first.first;
int var = it2->first.second.first;
int param1 = it2->first.second.second;
expr_t d1 = it2->second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
jacobian_params_second_derivatives[make_pair(eq, make_pair(var, make_pair(param1, param)))] = d2;
}
for (second_derivatives_t::const_iterator it2 = second_derivatives.begin();
it2 != second_derivatives.end(); it2++)
{
int eq = it2->first.first;
int var1 = it2->first.second.first;
int var2 = it2->first.second.second;
expr_t d1 = it2->second;
expr_t d2 = d1->getDerivative(param);
if (d2 == Zero)
continue;
hessian_params_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, param)))] = d2;
}
}
}
void
ModelTree::computeParamsDerivativesTemporaryTerms()
{
map<expr_t, int> reference_count;
params_derivs_temporary_terms.clear();
for (first_derivatives_t::iterator it = residuals_params_derivatives.begin();
it != residuals_params_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
for (second_derivatives_t::iterator it = jacobian_params_derivatives.begin();
it != jacobian_params_derivatives.end(); it++)
it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
for (second_derivatives_t::const_iterator it = residuals_params_second_derivatives.begin();
it != residuals_params_second_derivatives.end(); ++it)
it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
for (third_derivatives_t::const_iterator it = jacobian_params_second_derivatives.begin();
it != jacobian_params_second_derivatives.end(); ++it)
it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
for (third_derivatives_t::const_iterator it = hessian_params_derivatives.begin();
it != hessian_params_derivatives.end(); ++it)
it->second->computeTemporaryTerms(reference_count, params_derivs_temporary_terms, true);
}