1877 lines
66 KiB
C++
1877 lines
66 KiB
C++
/*
|
||
* Copyright © 2003-2022 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 <https://www.gnu.org/licenses/>.
|
||
*/
|
||
|
||
#include "ModelTree.hh"
|
||
#include "VariableDependencyGraph.hh"
|
||
|
||
#pragma GCC diagnostic push
|
||
#pragma GCC diagnostic ignored "-Wold-style-cast"
|
||
#pragma GCC diagnostic ignored "-Wsign-compare"
|
||
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
|
||
#include <boost/graph/adjacency_list.hpp>
|
||
#include <boost/graph/max_cardinality_matching.hpp>
|
||
#include <boost/graph/topological_sort.hpp>
|
||
#pragma GCC diagnostic pop
|
||
|
||
#ifdef __APPLE__
|
||
# include <mach-o/dyld.h>
|
||
#endif
|
||
|
||
#include <regex>
|
||
#include <utility>
|
||
|
||
void
|
||
ModelTree::copyHelper(const ModelTree &m)
|
||
{
|
||
auto f = [this](expr_t e) { return e->clone(*this); };
|
||
|
||
// Equations
|
||
for (const auto &it : m.equations)
|
||
equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
|
||
for (const auto &it : m.aux_equations)
|
||
aux_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
|
||
|
||
auto convert_deriv_map = [f](map<vector<int>, expr_t> dm)
|
||
{
|
||
map<vector<int>, expr_t> dm2;
|
||
for (const auto &it : dm)
|
||
dm2.emplace(it.first, f(it.second));
|
||
return dm2;
|
||
};
|
||
|
||
// Derivatives
|
||
for (const auto &it : m.derivatives)
|
||
derivatives.push_back(convert_deriv_map(it));
|
||
for (const auto &it : m.params_derivatives)
|
||
params_derivatives[it.first] = convert_deriv_map(it.second);
|
||
|
||
auto convert_temporary_terms_t = [f](temporary_terms_t tt)
|
||
{
|
||
temporary_terms_t tt2;
|
||
for (const auto &it : tt)
|
||
tt2.insert(f(it));
|
||
return tt2;
|
||
};
|
||
|
||
// Temporary terms
|
||
for (const auto &it : m.temporary_terms_mlv)
|
||
temporary_terms_mlv[dynamic_cast<VariableNode *>(f(it.first))] = f(it.second);
|
||
for (const auto &it : m.temporary_terms_derivatives)
|
||
temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
|
||
for (const auto &it : m.temporary_terms_idxs)
|
||
temporary_terms_idxs[f(it.first)] = it.second;
|
||
for (const auto &it : m.params_derivs_temporary_terms)
|
||
params_derivs_temporary_terms[it.first] = convert_temporary_terms_t(it.second);
|
||
for (const auto &it : m.params_derivs_temporary_terms_idxs)
|
||
params_derivs_temporary_terms_idxs[f(it.first)] = it.second;
|
||
|
||
// Other stuff
|
||
for (const auto &it : m.trend_symbols_map)
|
||
trend_symbols_map[it.first] = f(it.second);
|
||
for (const auto &it : m.nonstationary_symbols_map)
|
||
nonstationary_symbols_map[it.first] = {it.second.first, f(it.second.second)};
|
||
|
||
for (const auto &it : m.equation_type_and_normalized_equation)
|
||
equation_type_and_normalized_equation.emplace_back(it.first, dynamic_cast<BinaryOpNode *>(f(it.second)));
|
||
|
||
for (const auto &it : m.blocks_derivatives)
|
||
{
|
||
map<tuple<int, int, int>, expr_t> v;
|
||
for (const auto &it2 : it)
|
||
v[it2.first] = f(it2.second);
|
||
blocks_derivatives.push_back(v);
|
||
}
|
||
|
||
auto convert_vector_tt = [f](vector<temporary_terms_t> vtt)
|
||
{
|
||
vector<temporary_terms_t> vtt2;
|
||
for (const auto &tt : vtt)
|
||
{
|
||
temporary_terms_t tt2;
|
||
for (const auto &it : tt)
|
||
tt2.insert(f(it));
|
||
vtt2.push_back(tt2);
|
||
}
|
||
return vtt2;
|
||
};
|
||
for (const auto &it : m.blocks_temporary_terms)
|
||
blocks_temporary_terms.push_back(convert_vector_tt(it));
|
||
for (const auto &it : m.blocks_temporary_terms_idxs)
|
||
blocks_temporary_terms_idxs[f(it.first)] = it.second;
|
||
}
|
||
|
||
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
|
||
NumericalConstants &num_constants_arg,
|
||
ExternalFunctionsTable &external_functions_table_arg,
|
||
bool is_dynamic_arg) :
|
||
DataTree{symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
|
||
derivatives(4),
|
||
NNZDerivatives(4, 0),
|
||
temporary_terms_derivatives(4)
|
||
{
|
||
}
|
||
|
||
ModelTree::ModelTree(const ModelTree &m) :
|
||
DataTree{m},
|
||
user_set_add_flags{m.user_set_add_flags},
|
||
user_set_subst_flags{m.user_set_subst_flags},
|
||
user_set_add_libs{m.user_set_add_libs},
|
||
user_set_subst_libs{m.user_set_subst_libs},
|
||
user_set_compiler{m.user_set_compiler},
|
||
equations_lineno{m.equations_lineno},
|
||
equation_tags{m.equation_tags},
|
||
computed_derivs_order{m.computed_derivs_order},
|
||
NNZDerivatives{m.NNZDerivatives},
|
||
eq_idx_block2orig{m.eq_idx_block2orig},
|
||
endo_idx_block2orig{m.endo_idx_block2orig},
|
||
eq_idx_orig2block{m.eq_idx_orig2block},
|
||
endo_idx_orig2block{m.endo_idx_orig2block},
|
||
blocks{m.blocks},
|
||
endo2block{m.endo2block},
|
||
eq2block{m.eq2block},
|
||
endo2eq{m.endo2eq},
|
||
cutoff{m.cutoff},
|
||
mfs{m.mfs}
|
||
{
|
||
copyHelper(m);
|
||
}
|
||
|
||
ModelTree &
|
||
ModelTree::operator=(const ModelTree &m)
|
||
{
|
||
DataTree::operator=(m);
|
||
|
||
equations.clear();
|
||
equations_lineno = m.equations_lineno;
|
||
aux_equations.clear();
|
||
equation_tags = m.equation_tags;
|
||
computed_derivs_order = m.computed_derivs_order;
|
||
NNZDerivatives = m.NNZDerivatives;
|
||
|
||
derivatives.clear();
|
||
params_derivatives.clear();
|
||
|
||
temporary_terms_mlv.clear();
|
||
temporary_terms_derivatives.clear();
|
||
params_derivs_temporary_terms.clear();
|
||
params_derivs_temporary_terms_idxs.clear();
|
||
|
||
trend_symbols_map.clear();
|
||
nonstationary_symbols_map.clear();
|
||
|
||
eq_idx_block2orig = m.eq_idx_block2orig;
|
||
endo_idx_block2orig = m.endo_idx_block2orig;
|
||
eq_idx_orig2block = m.eq_idx_orig2block;
|
||
endo_idx_orig2block = m.endo_idx_orig2block;
|
||
equation_type_and_normalized_equation.clear();
|
||
blocks_derivatives.clear();
|
||
blocks = m.blocks;
|
||
endo2block = m.endo2block;
|
||
eq2block = m.eq2block;
|
||
blocks_temporary_terms.clear();
|
||
blocks_temporary_terms_idxs.clear();
|
||
endo2eq = m.endo2eq;
|
||
cutoff = m.cutoff;
|
||
mfs = m.mfs;
|
||
|
||
user_set_add_flags = m.user_set_add_flags;
|
||
user_set_subst_flags = m.user_set_subst_flags;
|
||
user_set_add_libs = m.user_set_add_libs;
|
||
user_set_subst_libs = m.user_set_subst_libs;
|
||
user_set_compiler = m.user_set_compiler;
|
||
|
||
copyHelper(m);
|
||
|
||
return *this;
|
||
}
|
||
|
||
bool
|
||
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool dynamic, bool verbose)
|
||
{
|
||
const int n = equations.size();
|
||
|
||
assert(n == symbol_table.endo_nbr());
|
||
|
||
using BipartiteGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
|
||
|
||
/*
|
||
Vertices 0 to n-1 are for endogenous (using type specific ID)
|
||
Vertices n to 2*n-1 are for equations (using equation no.)
|
||
*/
|
||
BipartiteGraph g(2 * n);
|
||
|
||
// Fill in the graph
|
||
for (const auto &[eq_and_endo, val] : contemporaneous_jacobian)
|
||
add_edge(eq_and_endo.first + n, eq_and_endo.second, g);
|
||
|
||
// Compute maximum cardinality matching
|
||
vector<int> mate_map(2*n);
|
||
|
||
bool check = checked_edmonds_maximum_cardinality_matching(g, &mate_map[0]);
|
||
|
||
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(equations.size());
|
||
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), [=](int i) { return i-n; });
|
||
|
||
// Check if all variables are normalized
|
||
if (auto it = find(mate_map.begin(), mate_map.begin() + n, boost::graph_traits<BipartiteGraph>::null_vertex());
|
||
it != mate_map.begin() + n)
|
||
{
|
||
if (verbose)
|
||
cerr << "ERROR: Could not normalize the " << (dynamic ? "dynamic" : "static") << " model. Variable "
|
||
<< symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))
|
||
<< " is not in the maximum cardinality matching." << endl;
|
||
check = false;
|
||
}
|
||
return check;
|
||
}
|
||
|
||
void
|
||
ModelTree::computeNonSingularNormalization(const jacob_map_t &contemporaneous_jacobian, bool dynamic)
|
||
{
|
||
cout << "Normalizing the " << (dynamic ? "dynamic" : "static") << " model..." << endl;
|
||
|
||
int n = equations.size();
|
||
|
||
// Compute the maximum value of each row of the contemporaneous Jacobian matrix
|
||
vector max_val(n, 0.0);
|
||
for (const auto &[eq_and_endo, val] : contemporaneous_jacobian)
|
||
max_val[eq_and_endo.first] = max(max_val[eq_and_endo.first], fabs(val));
|
||
|
||
// Compute normalized contemporaneous Jacobian
|
||
jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
|
||
for (auto &[eq_and_endo, val] : normalized_contemporaneous_jacobian)
|
||
val /= max_val[eq_and_endo.first];
|
||
|
||
// We start with the highest value of the cutoff and try to normalize the model
|
||
double current_cutoff = 0.99999999;
|
||
const double cutoff_lower_limit = 1e-19;
|
||
|
||
bool found_normalization = false;
|
||
int last_suppressed = 0;
|
||
while (!found_normalization && current_cutoff > cutoff_lower_limit)
|
||
{
|
||
// Drop elements below cutoff from normalized contemporaneous Jacobian
|
||
jacob_map_t normalized_contemporaneous_jacobian_above_cutoff;
|
||
int suppressed = 0;
|
||
for (const auto &[eq_and_endo, val] : normalized_contemporaneous_jacobian)
|
||
if (fabs(val) > max(current_cutoff, cutoff))
|
||
normalized_contemporaneous_jacobian_above_cutoff[eq_and_endo] = val;
|
||
else
|
||
suppressed++;
|
||
|
||
if (suppressed != last_suppressed)
|
||
found_normalization = computeNormalization(normalized_contemporaneous_jacobian_above_cutoff, dynamic, false);
|
||
last_suppressed = suppressed;
|
||
if (!found_normalization)
|
||
{
|
||
current_cutoff /= 2;
|
||
// In this last case try to normalize with the complete jacobian
|
||
if (current_cutoff <= cutoff_lower_limit)
|
||
found_normalization = computeNormalization(normalized_contemporaneous_jacobian, dynamic, false);
|
||
}
|
||
}
|
||
|
||
if (!found_normalization)
|
||
{
|
||
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.
|
||
TODO: Explain why symbolic_jacobian is not contemporaneous. */
|
||
auto symbolic_jacobian = computeSymbolicJacobian();
|
||
found_normalization = computeNormalization(symbolic_jacobian, dynamic, true);
|
||
}
|
||
|
||
if (!found_normalization)
|
||
{
|
||
/* Some models don’t have a steady state, and this can cause the
|
||
normalization to fail (e.g. if some variable only appears in a diff(),
|
||
it will disappear from the static model). Suggest the “no_static”
|
||
option as a possible solution. */
|
||
if (!dynamic)
|
||
cerr << "If your model does not have a steady state, you may want to try the 'no_static' option of the 'model' block." << endl;
|
||
/* The last call to computeNormalization(), which was verbose, already
|
||
printed an error message, so we can immediately exit. */
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
ModelTree::jacob_map_t
|
||
ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context) const
|
||
{
|
||
jacob_map_t contemporaneous_jacobian;
|
||
for (const auto &[indices, d1] : derivatives[1])
|
||
{
|
||
int deriv_id = indices[1];
|
||
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
||
{
|
||
int eq = indices[0];
|
||
int var { getTypeSpecificIDByDerivID(deriv_id) };
|
||
int lag = getLagByDerivID(deriv_id);
|
||
double val = 0;
|
||
try
|
||
{
|
||
val = d1->eval(eval_context);
|
||
}
|
||
catch (ExprNode::EvalExternalFunctionException &e)
|
||
{
|
||
val = 1;
|
||
}
|
||
catch (ExprNode::EvalException &e)
|
||
{
|
||
cerr << "ERROR: evaluation of Jacobian failed for equation " << eq+1;
|
||
if (equations_lineno[eq])
|
||
cerr << " (line " << *equations_lineno[eq] << ")";
|
||
cerr << " and variable " << getNameByDerivID(deriv_id) << "(" << lag << ") !" << endl;
|
||
d1->writeOutput(cerr, ExprNodeOutputType::matlabDynamicModel, {}, {});
|
||
cerr << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
if ((isnan(val) || fabs(val) >= cutoff) && lag == 0)
|
||
contemporaneous_jacobian[{ eq, var }] = val;
|
||
}
|
||
}
|
||
|
||
return contemporaneous_jacobian;
|
||
}
|
||
|
||
pair<int, int>
|
||
ModelTree::computePrologueAndEpilogue()
|
||
{
|
||
const int n = equations.size();
|
||
|
||
/* Initialize “eq_idx_block2orig” and “endo_idx_block2orig” to the identity
|
||
permutation. */
|
||
eq_idx_block2orig.resize(n);
|
||
endo_idx_block2orig.resize(n);
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
eq_idx_block2orig[i] = i;
|
||
endo_idx_block2orig[endo2eq[i]] = i;
|
||
}
|
||
|
||
/* Compute incidence matrix, equations in rows, variables in columns. Row
|
||
(resp. column) indices are to be interpreted according to
|
||
“eq_idx_block2orig” (resp. “endo_idx_block2orig”). Stored in row-major
|
||
order. */
|
||
vector IM(n*n, false);
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
set<pair<int, int>> endos_and_lags;
|
||
equations[i]->collectEndogenous(endos_and_lags);
|
||
for (auto [endo, lag] : endos_and_lags)
|
||
IM[i * n + endo2eq[endo]] = true;
|
||
}
|
||
|
||
bool something_has_been_done;
|
||
// Find the prologue equations and place first the AR(1) shock equations first
|
||
int prologue = 0;
|
||
do
|
||
{
|
||
something_has_been_done = false;
|
||
int new_prologue = prologue;
|
||
for (int i = prologue; i < n; i++)
|
||
{
|
||
int nze = 0;
|
||
int k = 0;
|
||
for (int j = new_prologue; j < n; j++)
|
||
if (IM[i * n + j])
|
||
{
|
||
nze++;
|
||
k = j;
|
||
}
|
||
if (nze == 1)
|
||
{
|
||
// Swap equations indexed by “new_prologue” and i
|
||
for (int j = 0; j < n; j++)
|
||
swap(IM[new_prologue * n + j], IM[i * n + j]);
|
||
swap(eq_idx_block2orig[new_prologue], eq_idx_block2orig[i]);
|
||
|
||
// Swap variables indexed by “new_prologue” and k (in the matching)
|
||
for (int j = 0; j < n; j++)
|
||
swap(IM[j * n + new_prologue], IM[j * n + k]);
|
||
swap(endo_idx_block2orig[new_prologue], endo_idx_block2orig[k]);
|
||
|
||
new_prologue++;
|
||
something_has_been_done = true;
|
||
}
|
||
}
|
||
prologue = new_prologue;
|
||
}
|
||
while (something_has_been_done);
|
||
|
||
// Find the epilogue equations
|
||
int epilogue = 0;
|
||
do
|
||
{
|
||
something_has_been_done = false;
|
||
int new_epilogue = epilogue;
|
||
for (int i = prologue; i < n - epilogue; i++)
|
||
{
|
||
int nze = 0;
|
||
int k = 0;
|
||
for (int j = prologue; j < n - new_epilogue; j++)
|
||
if (IM[j * n + i])
|
||
{
|
||
nze++;
|
||
k = j;
|
||
}
|
||
if (nze == 1)
|
||
{
|
||
for (int j = 0; j < n; j++)
|
||
swap(IM[(n - 1 - new_epilogue) * n + j], IM[k * n + j]);
|
||
swap(eq_idx_block2orig[n - 1 - new_epilogue], eq_idx_block2orig[k]);
|
||
|
||
for (int j = 0; j < n; j++)
|
||
swap(IM[j * n + n - 1 - new_epilogue], IM[j * n + i]);
|
||
swap(endo_idx_block2orig[n - 1 - new_epilogue], endo_idx_block2orig[i]);
|
||
|
||
new_epilogue++;
|
||
something_has_been_done = true;
|
||
}
|
||
}
|
||
epilogue = new_epilogue;
|
||
}
|
||
while (something_has_been_done);
|
||
|
||
updateReverseVariableEquationOrderings();
|
||
|
||
return { prologue, epilogue };
|
||
}
|
||
|
||
void
|
||
ModelTree::equationTypeDetermination(const map<tuple<int, int, int>, expr_t> &first_order_endo_derivatives, int mfs)
|
||
{
|
||
equation_type_and_normalized_equation.clear();
|
||
equation_type_and_normalized_equation.resize(equations.size());
|
||
for (int i = 0; i < static_cast<int>(equations.size()); i++)
|
||
{
|
||
int eq = eq_idx_block2orig[i];
|
||
int var = endo_idx_block2orig[i];
|
||
expr_t lhs = equations[eq]->arg1;
|
||
EquationType Equation_Simulation_Type = EquationType::solve;
|
||
BinaryOpNode *normalized_eq = nullptr;
|
||
if (auto it = first_order_endo_derivatives.find({ eq, var, 0 });
|
||
it != first_order_endo_derivatives.end())
|
||
{
|
||
expr_t derivative = it->second;
|
||
// Determine whether the equation can be evaluated rather than solved
|
||
if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, endo_idx_block2orig[i], 0)
|
||
&& derivative->isNumConstNodeEqualTo(1))
|
||
Equation_Simulation_Type = EquationType::evaluate;
|
||
else
|
||
{
|
||
set<pair<int, int>> result;
|
||
derivative->collectEndogenous(result);
|
||
bool variable_not_in_derivative = !result.contains({ var, 0 });
|
||
|
||
try
|
||
{
|
||
normalized_eq = equations[eq]->normalizeEquation(symbol_table.getID(SymbolType::endogenous, var), 0);
|
||
if ((mfs == 2 && variable_not_in_derivative) || mfs == 3)
|
||
Equation_Simulation_Type = EquationType::evaluateRenormalized;
|
||
}
|
||
catch (ExprNode::NormalizationFailed &e)
|
||
{
|
||
}
|
||
}
|
||
}
|
||
equation_type_and_normalized_equation[eq] = { Equation_Simulation_Type, normalized_eq };
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::computeDynamicStructureOfBlock(int blk)
|
||
{
|
||
vector max_endo_lag_lead(blocks[blk].size, pair{0, 0});
|
||
blocks[blk].max_endo_lag = blocks[blk].max_endo_lead = 0;
|
||
blocks[blk].max_other_endo_lag = blocks[blk].max_other_endo_lead = 0;
|
||
blocks[blk].max_exo_lag = blocks[blk].max_exo_lead = 0;
|
||
blocks[blk].max_exo_det_lag = blocks[blk].max_exo_det_lead = 0;
|
||
for (int eq = 0; eq < blocks[blk].size; eq++)
|
||
{
|
||
set<pair<int, int>> endos_and_lags;
|
||
expr_t e = getBlockEquationExpr(blk, eq);
|
||
|
||
/* Compute max lags/leads for endogenous. Also fill per-variable structure
|
||
for endos belonging to this block */
|
||
e->collectEndogenous(endos_and_lags);
|
||
for (auto [endo, lag] : endos_and_lags)
|
||
if (endo2block[endo] == blk)
|
||
{
|
||
blocks[blk].max_endo_lag = max(blocks[blk].max_endo_lag, -lag);
|
||
blocks[blk].max_endo_lead = max(blocks[blk].max_endo_lead, lag);
|
||
auto &[max_endo_lag, max_endo_lead] = max_endo_lag_lead[getBlockInitialVariableID(blk, endo)];
|
||
max_endo_lag = max(max_endo_lag, -lag);
|
||
max_endo_lead = max(max_endo_lead, lag);
|
||
}
|
||
else
|
||
{
|
||
blocks[blk].max_other_endo_lag = max(blocks[blk].max_other_endo_lag, -lag);
|
||
blocks[blk].max_other_endo_lead = max(blocks[blk].max_other_endo_lead, lag);
|
||
}
|
||
|
||
// Compute max lags/leads for exogenous
|
||
blocks[blk].max_exo_lag = max(e->maxExoLag(), blocks[blk].max_exo_lag);
|
||
blocks[blk].max_exo_lead = max(e->maxExoLead(), blocks[blk].max_exo_lead);
|
||
|
||
// Compute max lags/leads for deterministic exogenous
|
||
set<pair<int, int>> dynvars;
|
||
e->collectDynamicVariables(SymbolType::exogenousDet, dynvars);
|
||
for (auto [symb_id, lag] : dynvars)
|
||
{
|
||
blocks[blk].max_exo_det_lag = max(-lag, blocks[blk].max_exo_det_lag);
|
||
blocks[blk].max_exo_det_lead = max(lag, blocks[blk].max_exo_det_lead);
|
||
}
|
||
}
|
||
|
||
// Compute max lags/leads over all variables
|
||
blocks[blk].max_lag = max(blocks[blk].max_endo_lag, max(blocks[blk].max_other_endo_lag,
|
||
max(blocks[blk].max_exo_lag,
|
||
blocks[blk].max_exo_det_lag)));
|
||
blocks[blk].max_lead = max(blocks[blk].max_endo_lead, max(blocks[blk].max_other_endo_lead,
|
||
max(blocks[blk].max_exo_lead,
|
||
blocks[blk].max_exo_det_lead)));
|
||
|
||
// Categorize endos that belong to the block
|
||
blocks[blk].n_mixed = blocks[blk].n_forward = blocks[blk].n_backward = blocks[blk].n_static = 0;
|
||
for (int var = 0; var < blocks[blk].size; var++)
|
||
{
|
||
auto [max_lag, max_lead] = max_endo_lag_lead[var];
|
||
if (max_lag != 0 && max_lead != 0)
|
||
blocks[blk].n_mixed++;
|
||
else if (max_lag == 0 && max_lead != 0)
|
||
blocks[blk].n_forward++;
|
||
else if (max_lag != 0 && max_lead == 0)
|
||
blocks[blk].n_backward++;
|
||
else
|
||
blocks[blk].n_static++;
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::computeSimulationTypeOfBlock(int blk)
|
||
{
|
||
auto &type = blocks[blk].simulation_type;
|
||
if (blocks[blk].max_endo_lag > 0 && blocks[blk].max_endo_lead > 0)
|
||
{
|
||
if (blocks[blk].size == 1)
|
||
type = BlockSimulationType::solveTwoBoundariesSimple;
|
||
else
|
||
type = BlockSimulationType::solveTwoBoundariesComplete;
|
||
}
|
||
else if (blocks[blk].size > 1)
|
||
{
|
||
if (blocks[blk].max_endo_lead > 0)
|
||
type = BlockSimulationType::solveBackwardComplete;
|
||
else
|
||
type = BlockSimulationType::solveForwardComplete;
|
||
}
|
||
else
|
||
{
|
||
bool can_eval = (getBlockEquationType(blk, 0) == EquationType::evaluate
|
||
|| getBlockEquationType(blk, 0) == EquationType::evaluateRenormalized);
|
||
if (blocks[blk].max_endo_lead > 0)
|
||
type = can_eval ? BlockSimulationType::evaluateBackward :
|
||
BlockSimulationType::solveBackwardSimple;
|
||
else
|
||
type = can_eval ? BlockSimulationType::evaluateForward :
|
||
BlockSimulationType::solveForwardSimple;
|
||
}
|
||
}
|
||
|
||
pair<lag_lead_vector_t, lag_lead_vector_t>
|
||
ModelTree::getVariableLeadLagByBlock() const
|
||
{
|
||
int nb_endo = symbol_table.endo_nbr();
|
||
|
||
lag_lead_vector_t variable_lag_lead(nb_endo, { 0, 0 }), equation_lag_lead(nb_endo, { 0, 0 });
|
||
for (int eq = 0; eq < nb_endo; eq++)
|
||
{
|
||
set<pair<int, int>> endos_and_lags;
|
||
equations[eq]->collectEndogenous(endos_and_lags);
|
||
for (auto [endo, lag] : endos_and_lags)
|
||
if (endo2block[endo] == eq2block[eq])
|
||
{
|
||
variable_lag_lead[endo].first = max(variable_lag_lead[endo].first, -lag);
|
||
variable_lag_lead[endo].second = max(variable_lag_lead[endo].second, lag);
|
||
equation_lag_lead[eq].first = max(equation_lag_lead[eq].first, -lag);
|
||
equation_lag_lead[eq].second = max(equation_lag_lead[eq].second, lag);
|
||
}
|
||
}
|
||
return { equation_lag_lead, variable_lag_lead };
|
||
}
|
||
|
||
void
|
||
ModelTree::computeBlockDecomposition(int prologue, int epilogue)
|
||
{
|
||
int nb_var = symbol_table.endo_nbr();
|
||
int nb_simvars = nb_var - prologue - epilogue;
|
||
|
||
/* Construct the graph representing the dependencies between all
|
||
variables that do not belong to the prologue or the epilogue.
|
||
|
||
For detecting dependencies between variables, use the symbolic adjacency
|
||
matrix */
|
||
VariableDependencyGraph G(nb_simvars);
|
||
for (const auto &[key, value] : computeSymbolicJacobian())
|
||
{
|
||
auto [eq, endo] = key;
|
||
if (eq_idx_orig2block[eq] >= prologue
|
||
&& eq_idx_orig2block[eq] < nb_var - epilogue
|
||
&& endo_idx_orig2block[endo] >= prologue
|
||
&& endo_idx_orig2block[endo] < nb_var - epilogue
|
||
&& eq != endo2eq[endo])
|
||
add_edge(vertex(eq_idx_orig2block[endo2eq[endo]]-prologue, G),
|
||
vertex(eq_idx_orig2block[eq]-prologue, G), G);
|
||
}
|
||
|
||
/* Identify the simultaneous blocks. Each simultaneous block is given an
|
||
index, starting from 0, in recursive order */
|
||
auto [num_simblocks, simvar2simblock] = G.sortedStronglyConnectedComponents();
|
||
|
||
int num_blocks = prologue+num_simblocks+epilogue;
|
||
|
||
blocks.clear();
|
||
blocks.resize(num_blocks);
|
||
endo2block.resize(nb_var);
|
||
eq2block.resize(nb_var);
|
||
|
||
// Initialize size and mfs_size for prologue and epilogue, plus eq/endo→block mappings
|
||
for (int blk = 0; blk < num_blocks; blk++)
|
||
if (blk < prologue || blk >= num_blocks-epilogue)
|
||
{
|
||
int var_eq = (blk < prologue ? blk : blk-num_simblocks+nb_simvars);
|
||
blocks[blk].size = 1;
|
||
blocks[blk].mfs_size = 1;
|
||
blocks[blk].first_equation = var_eq;
|
||
endo2block[endo_idx_block2orig[var_eq]] = blk;
|
||
eq2block[eq_idx_block2orig[var_eq]] = blk;
|
||
}
|
||
|
||
// Initialize size for simultaneous blocks, plus eq/endo→block mappings
|
||
vector<vector<int>> simblock2simvars(num_simblocks);
|
||
for (int i = 0; i < static_cast<int>(simvar2simblock.size()); i++)
|
||
{
|
||
simblock2simvars[simvar2simblock[i]].push_back(i);
|
||
int blk = prologue+simvar2simblock[i];
|
||
blocks[blk].size++;
|
||
endo2block[endo_idx_block2orig[prologue+i]] = blk;
|
||
eq2block[eq_idx_block2orig[prologue+i]] = blk;
|
||
}
|
||
|
||
// Determine the dynamic structure of each block
|
||
auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock();
|
||
|
||
/* For each simultaneous block, the minimum set of feedback variable is computed.
|
||
Then, the variables within the blocks are reordered so that recursive
|
||
(non-feedback) appear first, to get a sub-recursive block without feedback variables.
|
||
Within each of the two sub-blocks, variables are reordered depending
|
||
on their dynamic status: static first, then backward, mixed and forward. */
|
||
|
||
/* First, add a loop on vertices which could not be normalized or vertices
|
||
related to lead/lag variables. This forces those vertices to belong to the
|
||
feedback set */
|
||
for (int i = 0; i < nb_simvars; i++)
|
||
if (equation_type_and_normalized_equation[eq_idx_block2orig[i+prologue]].first == EquationType::solve
|
||
|| variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0
|
||
|| variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
|
||
|| equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
|
||
|| equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0
|
||
|| mfs == 0)
|
||
add_edge(vertex(i, G), vertex(i, G), G);
|
||
|
||
const vector<int> old_eq_idx_block2orig(eq_idx_block2orig), old_endo_idx_block2orig(endo_idx_block2orig);
|
||
int ordidx = prologue;
|
||
for (int blk = prologue; blk < prologue+num_simblocks; blk++)
|
||
{
|
||
blocks[blk].first_equation = (blk == 0 ? 0 : blocks[blk-1].first_equation + blocks[blk-1].size);
|
||
auto subG = G.extractSubgraph(simblock2simvars[blk-prologue]);
|
||
auto feed_back_vertices = subG.minimalSetOfFeedbackVertices();
|
||
blocks[blk].mfs_size = feed_back_vertices.size();
|
||
auto recursive_vertices = subG.reorderRecursiveVariables(feed_back_vertices);
|
||
auto v_index1 = get(boost::vertex_index1, subG);
|
||
|
||
const vector dynamic_order{pair{0, 0}, pair{1, 0}, pair{1, 1}, pair{0, 1}};
|
||
|
||
// First the recursive equations conditional on feedback variables
|
||
for (auto max_lag_lead : dynamic_order)
|
||
for (int vtx : recursive_vertices)
|
||
if (int simvar = v_index1[vertex(vtx, subG)];
|
||
variable_lag_lead[old_endo_idx_block2orig[simvar+prologue]] == max_lag_lead)
|
||
{
|
||
eq_idx_block2orig[ordidx] = old_eq_idx_block2orig[simvar+prologue];
|
||
endo_idx_block2orig[ordidx] = old_endo_idx_block2orig[simvar+prologue];
|
||
ordidx++;
|
||
}
|
||
|
||
// Then the equations related to the feedback variables
|
||
for (auto max_lag_lead : dynamic_order)
|
||
for (int vtx : feed_back_vertices)
|
||
if (int simvar = v_index1[vertex(vtx, subG)];
|
||
variable_lag_lead[old_endo_idx_block2orig[simvar+prologue]] == max_lag_lead)
|
||
{
|
||
eq_idx_block2orig[ordidx] = old_eq_idx_block2orig[simvar+prologue];
|
||
endo_idx_block2orig[ordidx] = old_endo_idx_block2orig[simvar+prologue];
|
||
ordidx++;
|
||
}
|
||
}
|
||
|
||
updateReverseVariableEquationOrderings();
|
||
|
||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||
{
|
||
computeDynamicStructureOfBlock(blk);
|
||
computeSimulationTypeOfBlock(blk);
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::printBlockDecomposition() const
|
||
{
|
||
int largest_block = 0, Nb_SimulBlocks = 0, Nb_feedback_variable = 0;
|
||
int Nb_TotalBlocks = blocks.size();
|
||
for (int block = 0; block < Nb_TotalBlocks; block++)
|
||
if (BlockSimulationType simulation_type = blocks[block].simulation_type;
|
||
simulation_type == BlockSimulationType::solveForwardComplete
|
||
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
||
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete)
|
||
{
|
||
Nb_SimulBlocks++;
|
||
if (int size = blocks[block].size;
|
||
size > largest_block)
|
||
{
|
||
largest_block = size;
|
||
Nb_feedback_variable = blocks[block].mfs_size;
|
||
}
|
||
}
|
||
|
||
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;
|
||
}
|
||
|
||
void
|
||
ModelTree::reduceBlockDecomposition()
|
||
{
|
||
for (int blk = 1; blk < static_cast<int>(blocks.size()); blk++)
|
||
if (blocks[blk].size == 1)
|
||
{
|
||
/* Try to merge this block with the previous one.
|
||
This is only possible if the two blocks can simply be evaluated
|
||
(in the same direction), and if the merge does not break the
|
||
restrictions on leads/lags. */
|
||
set<pair<int, int>> endos_and_lags;
|
||
getBlockEquationExpr(blk, 0)->collectEndogenous(endos_and_lags);
|
||
bool is_lead = false, is_lag = false;
|
||
for (int var = 0; var < blocks[blk-1].size; var++)
|
||
{
|
||
is_lag = is_lag || endos_and_lags.contains({ getBlockVariableID(blk-1, var), -1 });
|
||
is_lead = is_lead || endos_and_lags.contains({ getBlockVariableID(blk-1, var), 1 });
|
||
}
|
||
|
||
if ((blocks[blk-1].simulation_type == BlockSimulationType::evaluateForward
|
||
&& blocks[blk].simulation_type == BlockSimulationType::evaluateForward
|
||
&& !is_lead)
|
||
|| (blocks[blk-1].simulation_type == BlockSimulationType::evaluateBackward
|
||
&& blocks[blk].simulation_type == BlockSimulationType::evaluateBackward
|
||
&& !is_lag))
|
||
{
|
||
// Merge the current block into the previous one
|
||
blocks[blk-1].size++;
|
||
blocks[blk-1].mfs_size = blocks[blk-1].size;
|
||
computeDynamicStructureOfBlock(blk-1);
|
||
blocks.erase(blocks.begin()+blk);
|
||
for (auto &b : endo2block)
|
||
if (b >= blk)
|
||
b--;
|
||
for (auto &b : eq2block)
|
||
if (b >= blk)
|
||
b--;
|
||
blk--;
|
||
continue;
|
||
}
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::determineLinearBlocks()
|
||
{
|
||
// Note that field “linear” in class BlockInfo defaults to true
|
||
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
||
switch (blocks[blk].simulation_type)
|
||
{
|
||
case BlockSimulationType::solveBackwardSimple:
|
||
case BlockSimulationType::solveBackwardComplete:
|
||
case BlockSimulationType::solveForwardSimple:
|
||
case BlockSimulationType::solveForwardComplete:
|
||
for (const auto &[indices, d1] : blocks_derivatives[blk])
|
||
{
|
||
int lag = get<2>(indices);
|
||
if (lag == 0)
|
||
{
|
||
set<pair<int, int>> endogenous;
|
||
d1->collectEndogenous(endogenous);
|
||
for (int l = 0; l < blocks[blk].size; l++)
|
||
if (endogenous.contains({ endo_idx_block2orig[blocks[blk].first_equation+l], 0 }))
|
||
{
|
||
blocks[blk].linear = false;
|
||
goto the_end;
|
||
}
|
||
}
|
||
}
|
||
the_end:
|
||
break;
|
||
case BlockSimulationType::solveTwoBoundariesComplete:
|
||
case BlockSimulationType::solveTwoBoundariesSimple:
|
||
for (const auto &[indices, d1] : blocks_derivatives[blk])
|
||
{
|
||
int lag = get<2>(indices);
|
||
set<pair<int, int>> endogenous;
|
||
d1->collectEndogenous(endogenous);
|
||
for (int l = 0; l < blocks[blk].size; l++)
|
||
if (endogenous.contains({ endo_idx_block2orig[blocks[blk].first_equation+l], lag }))
|
||
{
|
||
blocks[blk].linear = false;
|
||
goto the_end2;
|
||
}
|
||
}
|
||
the_end2:
|
||
break;
|
||
default:
|
||
break;
|
||
}
|
||
}
|
||
|
||
int
|
||
ModelTree::equation_number() const
|
||
{
|
||
return (equations.size());
|
||
}
|
||
|
||
void
|
||
ModelTree::computeDerivatives(int order, const set<int> &vars)
|
||
{
|
||
assert(order >= 1);
|
||
|
||
computed_derivs_order = order;
|
||
|
||
// Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
|
||
derivatives.resize(max(static_cast<size_t>(order+1), derivatives.size()));
|
||
NNZDerivatives.resize(max(static_cast<size_t>(order+1), NNZDerivatives.size()), 0);
|
||
|
||
// First-order derivatives
|
||
for (int var : vars)
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
{
|
||
expr_t d1 = equations[eq]->getDerivative(var);
|
||
if (d1 == Zero)
|
||
continue;
|
||
derivatives[1][{ eq, var }] = d1;
|
||
++NNZDerivatives[1];
|
||
}
|
||
|
||
// Higher-order derivatives
|
||
for (int o = 2; o <= order; o++)
|
||
for (const auto &[lower_indices, lower_d] : derivatives[o-1])
|
||
for (int var : vars)
|
||
{
|
||
if (lower_indices.back() > var)
|
||
continue;
|
||
|
||
expr_t d = lower_d->getDerivative(var);
|
||
if (d == Zero)
|
||
continue;
|
||
|
||
vector<int> indices{lower_indices};
|
||
indices.push_back(var);
|
||
// At this point, indices of endogenous variables are sorted in non-decreasing order
|
||
derivatives[o][indices] = d;
|
||
// We output symmetric elements at order = 2
|
||
if (o == 2 && indices[1] != indices[2])
|
||
NNZDerivatives[o] += 2;
|
||
else
|
||
NNZDerivatives[o]++;
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
|
||
{
|
||
/* Collect all model local variables appearing in equations (and only those,
|
||
because printing unused model local variables can lead to a crash,
|
||
see Dynare/dynare#101).
|
||
Then store them in a dedicated structure (temporary_terms_mlv), that will
|
||
be treated as the rest of temporary terms. */
|
||
temporary_terms_mlv.clear();
|
||
set<int> used_local_vars;
|
||
for (auto &equation : equations)
|
||
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
||
for (int used_local_var : used_local_vars)
|
||
{
|
||
VariableNode *v = AddVariable(used_local_var);
|
||
temporary_terms_mlv[v] = local_variables_table.find(used_local_var)->second;
|
||
}
|
||
|
||
// Compute the temporary terms in equations and derivatives
|
||
map<pair<int, int>, temporary_terms_t> temp_terms_map;
|
||
map<expr_t, pair<int, pair<int, int>>> reference_count;
|
||
|
||
for (auto &equation : equations)
|
||
equation->computeTemporaryTerms({ 0, 0 },
|
||
temp_terms_map,
|
||
reference_count,
|
||
is_matlab);
|
||
|
||
for (int order = 1; order < static_cast<int>(derivatives.size()); order++)
|
||
for (const auto &it : derivatives[order])
|
||
it.second->computeTemporaryTerms({ 0, order },
|
||
temp_terms_map,
|
||
reference_count,
|
||
is_matlab);
|
||
|
||
/* If the user has specified the notmpterms option, clear all temporary
|
||
terms, except those that correspond to external functions (since they are
|
||
not optional) */
|
||
if (no_tmp_terms)
|
||
for (auto &it : temp_terms_map)
|
||
erase_if(it.second,
|
||
[](expr_t e) { return !dynamic_cast<AbstractExternalFunctionNode *>(e); });
|
||
|
||
// Fill the structures
|
||
temporary_terms_derivatives.clear();
|
||
temporary_terms_derivatives.resize(derivatives.size());
|
||
for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
|
||
temporary_terms_derivatives[order] = move(temp_terms_map[{ 0, order }]);
|
||
|
||
// Compute indices in MATLAB/Julia vector
|
||
int idx = 0;
|
||
for (auto [mlv, value] : temporary_terms_mlv)
|
||
temporary_terms_idxs[mlv] = idx++;
|
||
for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
|
||
for (auto it : temporary_terms_derivatives[order])
|
||
temporary_terms_idxs[it] = idx++;
|
||
}
|
||
|
||
void
|
||
ModelTree::computeBlockTemporaryTerms()
|
||
{
|
||
int nb_blocks = blocks.size();
|
||
blocks_temporary_terms.resize(nb_blocks);
|
||
|
||
map<expr_t, tuple<int, int, int>> reference_count;
|
||
for (int blk = 0; blk < nb_blocks; blk++)
|
||
{
|
||
blocks_temporary_terms[blk].resize(blocks[blk].size + 1);
|
||
for (int eq = 0; eq < blocks[blk].size; eq++)
|
||
{
|
||
if (eq < blocks[blk].getRecursiveSize() && isBlockEquationRenormalized(blk, eq))
|
||
getBlockEquationRenormalizedExpr(blk, eq)->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
|
||
else
|
||
getBlockEquationExpr(blk, eq)->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
|
||
}
|
||
for (const auto &[ignore, d] : blocks_derivatives[blk])
|
||
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
|
||
|
||
additionalBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
|
||
}
|
||
|
||
// Compute indices in the temporary terms vector
|
||
blocks_temporary_terms_idxs.clear();
|
||
for (int idx{0};
|
||
auto &blk_tt : blocks_temporary_terms)
|
||
for (auto &eq_tt : blk_tt)
|
||
for (auto tt : eq_tt)
|
||
blocks_temporary_terms_idxs[tt] = idx++;
|
||
}
|
||
|
||
void
|
||
ModelTree::additionalBlockTemporaryTerms([[maybe_unused]] int blk,
|
||
[[maybe_unused]] vector<vector<temporary_terms_t>> &blocks_temporary_terms,
|
||
[[maybe_unused]] map<expr_t, tuple<int, int, int>> &reference_count) const
|
||
{
|
||
}
|
||
|
||
void
|
||
ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt,
|
||
temporary_terms_t &temp_term_union,
|
||
ostream &output,
|
||
deriv_node_temp_terms_t &tef_terms, const string &concat) const
|
||
{
|
||
// Local var used to keep track of temp nodes already written
|
||
temporary_terms_t tt2 = temp_term_union;
|
||
|
||
output << R"("external_functions_temporary_terms_)" << concat << R"(": [)";
|
||
for (bool printed_term{false};
|
||
auto it : tt)
|
||
{
|
||
if (dynamic_cast<AbstractExternalFunctionNode *>(it))
|
||
{
|
||
if (exchange(printed_term, true))
|
||
output << ", ";
|
||
vector<string> efout;
|
||
it->writeJsonExternalFunctionOutput(efout, tt2, tef_terms);
|
||
for (bool printed_efout{false};
|
||
auto &it : efout)
|
||
{
|
||
if (exchange(printed_efout, true))
|
||
output << ", ";
|
||
output << it;
|
||
}
|
||
}
|
||
tt2.insert(it);
|
||
}
|
||
|
||
output << "]"
|
||
<< R"(, "temporary_terms_)" << concat << R"(": [)";
|
||
for (bool printed_term{false};
|
||
const auto &it : tt)
|
||
{
|
||
if (exchange(printed_term, true))
|
||
output << ", ";
|
||
output << R"({"temporary_term": ")";
|
||
it->writeJsonOutput(output, tt, tef_terms);
|
||
output << R"(")"
|
||
<< R"(, "value": ")";
|
||
it->writeJsonOutput(output, temp_term_union, tef_terms);
|
||
output << R"("})" << endl;
|
||
|
||
temp_term_union.insert(it);
|
||
}
|
||
output << "]";
|
||
}
|
||
|
||
void
|
||
ModelTree::fixNestedParenthesis(ostringstream &output, map<string, string> &tmp_paren_vars, bool &message_printed) const
|
||
{
|
||
string str = output.str();
|
||
if (!testNestedParenthesis(str))
|
||
return;
|
||
int open = 0;
|
||
int first_open_paren = 0;
|
||
int matching_paren = 0;
|
||
bool hit_limit = false;
|
||
int i1 = 0;
|
||
for (size_t i = 0; i < str.length(); i++)
|
||
{
|
||
if (str.at(i) == '(')
|
||
{
|
||
if (open == 0)
|
||
first_open_paren = i;
|
||
open++;
|
||
}
|
||
else if (str.at(i) == ')')
|
||
{
|
||
open--;
|
||
if (open == 0)
|
||
matching_paren = i;
|
||
}
|
||
if (open > 32)
|
||
hit_limit = true;
|
||
|
||
if (hit_limit && open == 0)
|
||
{
|
||
if (!message_printed)
|
||
{
|
||
cerr << "Warning: A .m file created by Dynare will have more than 32 nested parenthesis. MATLAB cannot support this. " << endl
|
||
<< " We are going to modify, albeit inefficiently, this output to have fewer than 32 nested parenthesis. " << endl
|
||
<< " It would hence behoove you to use the use_dll option of the model block to circumnavigate this problem." << endl
|
||
<< " If you have not yet set up a compiler on your system, see the MATLAB documentation for doing so." << endl
|
||
<< " For Windows, see: https://www.mathworks.com/help/matlab/matlab_external/install-mingw-support-package.html" << endl << endl;
|
||
message_printed = true;
|
||
}
|
||
string str1 = str.substr(first_open_paren, matching_paren - first_open_paren + 1);
|
||
string repstr, varname;
|
||
while (testNestedParenthesis(str1))
|
||
{
|
||
size_t open_paren_idx = string::npos;
|
||
size_t match_paren_idx = string::npos;
|
||
size_t last_open_paren = string::npos;
|
||
for (size_t j = 0; j < str1.length(); j++)
|
||
{
|
||
if (str1.at(j) == '(')
|
||
{
|
||
// don't match, e.g. y(1)
|
||
if (size_t idx = str1.find_last_of("*/-+", j - 1);
|
||
j == 0 || (idx != string::npos && idx == j - 1))
|
||
open_paren_idx = j;
|
||
last_open_paren = j;
|
||
}
|
||
else if (str1.at(j) == ')')
|
||
{
|
||
// don't match, e.g. y(1)
|
||
if (size_t idx = str1.find_last_not_of("0123456789", j - 1);
|
||
idx != string::npos && idx != last_open_paren)
|
||
match_paren_idx = j;
|
||
}
|
||
|
||
if (open_paren_idx != string::npos && match_paren_idx != string::npos)
|
||
{
|
||
string val = str1.substr(open_paren_idx, match_paren_idx - open_paren_idx + 1);
|
||
if (auto it = tmp_paren_vars.find(val);
|
||
it == tmp_paren_vars.end())
|
||
{
|
||
varname = "paren32_tmp_var_" + to_string(i1++);
|
||
repstr = repstr + varname + " = " + val + ";\n";
|
||
tmp_paren_vars[val] = varname;
|
||
}
|
||
else
|
||
varname = it->second;
|
||
str1.replace(open_paren_idx, match_paren_idx - open_paren_idx + 1, varname);
|
||
break;
|
||
}
|
||
}
|
||
}
|
||
if (auto it = tmp_paren_vars.find(str1);
|
||
it == tmp_paren_vars.end())
|
||
{
|
||
varname = "paren32_tmp_var_" + to_string(i1++);
|
||
repstr = repstr + varname + " = " + str1 + ";\n";
|
||
}
|
||
else
|
||
varname = it->second;
|
||
str.replace(first_open_paren, matching_paren - first_open_paren + 1, varname);
|
||
size_t insertLoc = str.find_last_of("\n", first_open_paren);
|
||
str.insert(insertLoc + 1, repstr);
|
||
hit_limit = false;
|
||
i = -1;
|
||
first_open_paren = matching_paren = open = 0;
|
||
}
|
||
}
|
||
output.str(str);
|
||
}
|
||
|
||
bool
|
||
ModelTree::testNestedParenthesis(const string &str) const
|
||
{
|
||
for (int open{0};
|
||
char i : str)
|
||
{
|
||
if (i == '(')
|
||
open++;
|
||
else if (i == ')')
|
||
open--;
|
||
if (open > 32)
|
||
return true;
|
||
}
|
||
return false;
|
||
}
|
||
|
||
void
|
||
ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, 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;
|
||
|
||
for (auto equation : equations)
|
||
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
||
|
||
output << R"("model_local_variables": [)";
|
||
for (bool printed_something{false};
|
||
int id : local_variables_vector)
|
||
if (used_local_vars.contains(id))
|
||
{
|
||
if (exchange(printed_something, true))
|
||
output << ", ";
|
||
|
||
expr_t value = local_variables_table.find(id)->second;
|
||
if (write_tef_terms)
|
||
{
|
||
vector<string> efout;
|
||
value->writeJsonExternalFunctionOutput(efout, {}, tef_terms);
|
||
for (bool printed_efout{false};
|
||
auto &it : efout)
|
||
{
|
||
if (exchange(printed_efout, true))
|
||
output << ", ";
|
||
output << it;
|
||
}
|
||
|
||
if (!efout.empty())
|
||
output << ", ";
|
||
}
|
||
|
||
output << R"({"variable": ")" << symbol_table.getName(id)
|
||
<< R"(", "value": ")";
|
||
value->writeJsonOutput(output, {}, tef_terms);
|
||
output << R"("})" << endl;
|
||
}
|
||
output << "]";
|
||
}
|
||
|
||
int
|
||
ModelTree::writeBytecodeBinFile(const string &filename, bool is_two_boundaries) const
|
||
{
|
||
ofstream SaveCode { filename, ios::out | ios::binary };
|
||
if (!SaveCode.is_open())
|
||
{
|
||
cerr << R"(Error : Can't open file ")" << filename << R"(" for writing)" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
int u_count {0};
|
||
for (const auto &[indices, d1] : derivatives[1])
|
||
if (int deriv_id {indices[1]};
|
||
getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
||
{
|
||
int eq {indices[0]};
|
||
SaveCode.write(reinterpret_cast<char *>(&eq), sizeof eq);
|
||
int tsid {getTypeSpecificIDByDerivID(deriv_id)};
|
||
int lag {getLagByDerivID(deriv_id)};
|
||
int varr {tsid + lag * symbol_table.endo_nbr()};
|
||
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof varr);
|
||
SaveCode.write(reinterpret_cast<char *>(&lag), sizeof lag);
|
||
int u {u_count + symbol_table.endo_nbr()};
|
||
SaveCode.write(reinterpret_cast<char *>(&u), sizeof u);
|
||
u_count++;
|
||
}
|
||
if (is_two_boundaries)
|
||
u_count += symbol_table.endo_nbr();
|
||
for (int j {0}; j < symbol_table.endo_nbr(); j++)
|
||
SaveCode.write(reinterpret_cast<char *>(&j), sizeof j);
|
||
for (int j {0}; j < symbol_table.endo_nbr(); j++)
|
||
SaveCode.write(reinterpret_cast<char *>(&j), sizeof j);
|
||
SaveCode.close();
|
||
return u_count;
|
||
}
|
||
|
||
int
|
||
ModelTree::writeBlockBytecodeBinFile(ofstream &bin_file, int block) const
|
||
{
|
||
int u_count {0};
|
||
int block_size {blocks[block].size};
|
||
int block_mfs {blocks[block].mfs_size};
|
||
int block_recursive {blocks[block].getRecursiveSize()};
|
||
BlockSimulationType simulation_type {blocks[block].simulation_type};
|
||
bool is_two_boundaries {simulation_type == BlockSimulationType::solveTwoBoundariesComplete
|
||
|| simulation_type == BlockSimulationType::solveTwoBoundariesSimple};
|
||
for (const auto &[indices, ignore] : blocks_derivatives[block])
|
||
{
|
||
const auto &[eq, var, lag] {indices};
|
||
if (lag != 0 && !is_two_boundaries)
|
||
continue;
|
||
if (eq >= block_recursive && var >= block_recursive)
|
||
{
|
||
int v {eq - block_recursive};
|
||
bin_file.write(reinterpret_cast<char *>(&v), sizeof v);
|
||
int varr {var - block_recursive + lag * block_mfs};
|
||
bin_file.write(reinterpret_cast<char *>(&varr), sizeof varr);
|
||
bin_file.write(reinterpret_cast<const char *>(&lag), sizeof lag);
|
||
int u {u_count + block_mfs};
|
||
bin_file.write(reinterpret_cast<char *>(&u), sizeof u);
|
||
u_count++;
|
||
}
|
||
}
|
||
|
||
if (is_two_boundaries)
|
||
u_count += block_mfs;
|
||
for (int j {block_recursive}; j < block_size; j++)
|
||
{
|
||
int varr {getBlockVariableID(block, j)};
|
||
bin_file.write(reinterpret_cast<char *>(&varr), sizeof varr);
|
||
}
|
||
for (int j {block_recursive}; j < block_size; j++)
|
||
{
|
||
int eqr {getBlockEquationID(block, j)};
|
||
bin_file.write(reinterpret_cast<char *>(&eqr), sizeof eqr);
|
||
}
|
||
return u_count;
|
||
}
|
||
|
||
void
|
||
ModelTree::writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const
|
||
{
|
||
filesystem::create_directories(mod_basename + "/latex");
|
||
|
||
string filename = mod_basename + "/latex/" + latex_basename + ".tex";
|
||
string content_filename = mod_basename + "/latex/" + latex_basename + "_content" + ".tex";
|
||
ofstream output{filename, ios::out | ios::binary};
|
||
if (!output.is_open())
|
||
{
|
||
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
ofstream content_output{content_filename, ios::out | ios::binary};
|
||
if (!content_output.is_open())
|
||
{
|
||
cerr << "ERROR: Can't open file " << content_filename << " for writing" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
output << R"(\documentclass[10pt,a4paper]{article})" << endl
|
||
<< R"(\usepackage[landscape]{geometry})" << endl
|
||
<< R"(\usepackage{fullpage})" << endl
|
||
<< R"(\usepackage{amsfonts})" << endl
|
||
<< R"(\usepackage{breqn})" << endl
|
||
<< R"(\begin{document})" << endl
|
||
<< R"(\footnotesize)" << endl;
|
||
|
||
// Write model local variables
|
||
for (int id : local_variables_vector)
|
||
{
|
||
expr_t value = local_variables_table.find(id)->second;
|
||
|
||
content_output << R"(\begin{dmath*})" << endl
|
||
<< symbol_table.getTeXName(id) << " = ";
|
||
// Use an empty set for the temporary terms
|
||
value->writeOutput(content_output, output_type);
|
||
content_output << endl << R"(\end{dmath*})" << endl;
|
||
}
|
||
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
{
|
||
content_output << "% Equation " << eq + 1 << endl;
|
||
if (write_equation_tags)
|
||
equation_tags.writeLatexOutput(content_output, eq);
|
||
|
||
content_output << R"(\begin{dmath})" << endl;
|
||
// Here it is necessary to cast to superclass ExprNode, otherwise the overloaded writeOutput() method is not found
|
||
dynamic_cast<ExprNode *>(equations[eq])->writeOutput(content_output, output_type);
|
||
content_output << endl << R"(\end{dmath})" << endl;
|
||
}
|
||
|
||
output << R"(\include{)" << latex_basename + "_content" << "}" << endl
|
||
<< R"(\end{document})" << endl;
|
||
|
||
output.close();
|
||
content_output.close();
|
||
}
|
||
|
||
void
|
||
ModelTree::addEquation(expr_t eq, optional<int> lineno)
|
||
{
|
||
auto beq = dynamic_cast<BinaryOpNode *>(eq);
|
||
assert(beq && beq->op_code == BinaryOpcode::equal);
|
||
|
||
equations.push_back(beq);
|
||
equations_lineno.push_back(move(lineno));
|
||
}
|
||
|
||
void
|
||
ModelTree::findConstantEquationsWithoutMcpTag(map<VariableNode *, NumConstNode *> &subst_table) const
|
||
{
|
||
for (size_t i = 0; i < equations.size(); i++)
|
||
if (auto tags = getEquationTags(i);
|
||
tags.find("mcp") == tags.end())
|
||
equations[i]->findConstantEquations(subst_table);
|
||
}
|
||
|
||
void
|
||
ModelTree::addEquation(expr_t eq, optional<int> lineno, const map<string, string> &eq_tags)
|
||
{
|
||
equation_tags.add(equations.size(), eq_tags);
|
||
addEquation(eq, move(lineno));
|
||
}
|
||
|
||
void
|
||
ModelTree::addAuxEquation(expr_t eq)
|
||
{
|
||
auto beq = dynamic_cast<BinaryOpNode *>(eq);
|
||
assert(beq && beq->op_code == BinaryOpcode::equal);
|
||
|
||
aux_equations.push_back(beq);
|
||
}
|
||
|
||
void
|
||
ModelTree::addTrendVariables(const vector<int> &trend_vars, expr_t growth_factor) noexcept(false)
|
||
{
|
||
for (int id : trend_vars)
|
||
if (trend_symbols_map.contains(id))
|
||
throw TrendException(symbol_table.getName(id));
|
||
else
|
||
trend_symbols_map[id] = growth_factor;
|
||
}
|
||
|
||
void
|
||
ModelTree::addNonstationaryVariables(const vector<int> &nonstationary_vars, bool log_deflator, expr_t deflator) noexcept(false)
|
||
{
|
||
for (int id : nonstationary_vars)
|
||
if (nonstationary_symbols_map.contains(id))
|
||
throw TrendException(symbol_table.getName(id));
|
||
else
|
||
nonstationary_symbols_map[id] = { log_deflator, deflator };
|
||
}
|
||
|
||
void
|
||
ModelTree::initializeVariablesAndEquations()
|
||
{
|
||
for (size_t j = 0; j < equations.size(); j++)
|
||
eq_idx_block2orig.push_back(j);
|
||
|
||
for (int j = 0; j < symbol_table.endo_nbr(); j++)
|
||
endo_idx_block2orig.push_back(j);
|
||
}
|
||
|
||
void
|
||
ModelTree::set_cutoff_to_zero()
|
||
{
|
||
cutoff = 0;
|
||
}
|
||
|
||
void
|
||
ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
|
||
{
|
||
assert(paramsDerivsOrder >= 1);
|
||
|
||
set<int> deriv_id_set;
|
||
addAllParamDerivId(deriv_id_set);
|
||
|
||
// First-order derivatives w.r.t. params
|
||
for (int param : deriv_id_set)
|
||
{
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
{
|
||
expr_t d = equations[eq]->getDerivative(param);
|
||
if (d == Zero)
|
||
continue;
|
||
params_derivatives[{ 0, 1 }][{ eq, param }] = d;
|
||
}
|
||
|
||
for (int endoOrd = 1; endoOrd < static_cast<int>(derivatives.size()); endoOrd++)
|
||
for (const auto &[lower_indices, lower_d] : derivatives[endoOrd])
|
||
{
|
||
expr_t d = lower_d->getDerivative(param);
|
||
if (d == Zero)
|
||
continue;
|
||
vector<int> indices{lower_indices};
|
||
indices.push_back(param);
|
||
params_derivatives[{ endoOrd, 1 }][indices] = d;
|
||
}
|
||
}
|
||
|
||
// Higher-order derivatives w.r.t. parameters
|
||
for (int endoOrd = 0; endoOrd < static_cast<int>(derivatives.size()); endoOrd++)
|
||
for (int paramOrd = 2; paramOrd <= paramsDerivsOrder; paramOrd++)
|
||
for (const auto &[lower_indices, lower_d] : params_derivatives[{ endoOrd, paramOrd-1 }])
|
||
for (int param : deriv_id_set)
|
||
{
|
||
if (lower_indices.back() > param)
|
||
continue;
|
||
|
||
expr_t d = lower_d->getDerivative(param);
|
||
if (d == Zero)
|
||
continue;
|
||
vector<int> indices{lower_indices};
|
||
indices.push_back(param);
|
||
// At this point, indices of both endogenous and parameters are sorted in non-decreasing order
|
||
params_derivatives[{ endoOrd, paramOrd }][indices] = d;
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::computeParamsDerivativesTemporaryTerms()
|
||
{
|
||
map<expr_t, pair<int, pair<int, int>>> reference_count;
|
||
|
||
/* The temp terms should be constructed in the same order as the for loops in
|
||
{Static,Dynamic}Model::write{Json,}ParamsDerivativesFile() */
|
||
params_derivs_temporary_terms.clear();
|
||
for (const auto &[order, derivs] : params_derivatives)
|
||
for (const auto &[indices, d] : derivs)
|
||
d->computeTemporaryTerms(order, params_derivs_temporary_terms,
|
||
reference_count, true);
|
||
|
||
int idx = 0;
|
||
for (auto &[mlv, value] : temporary_terms_mlv)
|
||
params_derivs_temporary_terms_idxs[mlv] = idx++;
|
||
for (const auto &[order, tts] : params_derivs_temporary_terms)
|
||
for (const auto &tt : tts)
|
||
params_derivs_temporary_terms_idxs[tt] = idx++;
|
||
}
|
||
|
||
bool
|
||
ModelTree::isNonstationary(int symb_id) const
|
||
{
|
||
return nonstationary_symbols_map.contains(symb_id);
|
||
}
|
||
|
||
void
|
||
ModelTree::writeJsonModelEquations(ostream &output, bool residuals) const
|
||
{
|
||
if (residuals)
|
||
output << endl << R"("residuals":[)" << endl;
|
||
else
|
||
output << endl << R"("model":[)" << endl;
|
||
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
||
{
|
||
if (eq > 0)
|
||
output << ", ";
|
||
|
||
BinaryOpNode *eq_node = equations[eq];
|
||
expr_t lhs = eq_node->arg1;
|
||
expr_t rhs = eq_node->arg2;
|
||
|
||
if (residuals)
|
||
{
|
||
output << R"({"residual": {)"
|
||
<< R"("lhs": ")";
|
||
lhs->writeJsonOutput(output, {}, {});
|
||
output << R"(")";
|
||
|
||
output << R"(, "rhs": ")";
|
||
rhs->writeJsonOutput(output, {}, {});
|
||
output << R"("})";
|
||
}
|
||
else
|
||
{
|
||
output << R"({"lhs": ")";
|
||
lhs->writeJsonOutput(output, {}, {});
|
||
output << R"(", "rhs": ")";
|
||
rhs->writeJsonOutput(output, {}, {});
|
||
output << R"(")";
|
||
if (equations_lineno[eq])
|
||
output << R"(, "line": )" << *equations_lineno[eq];
|
||
|
||
if (auto eqtags = getEquationTags(eq);
|
||
!eqtags.empty())
|
||
{
|
||
output << R"(, "tags": {)";
|
||
for (bool printed_something{false};
|
||
const auto &[name, value] : eqtags)
|
||
{
|
||
if (exchange(printed_something, true))
|
||
output << ", ";
|
||
output << R"(")" << name << R"(": ")" << value << R"(")";
|
||
}
|
||
output << "}";
|
||
eqtags.clear();
|
||
}
|
||
}
|
||
output << "}" << endl;
|
||
}
|
||
output << endl << "]" << endl;
|
||
}
|
||
|
||
string
|
||
ModelTree::matlab_arch(const string &mexext)
|
||
{
|
||
if (mexext == "mexglx")
|
||
return "glnx86";
|
||
else if (mexext == "mexa64")
|
||
return "glnxa64";
|
||
if (mexext == "mexw32")
|
||
return "win32";
|
||
else if (mexext == "mexw64")
|
||
return "win64";
|
||
else if (mexext == "mexmaci")
|
||
{
|
||
cerr << "32-bit MATLAB not supported on macOS" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
else if (mexext == "mexmaci64")
|
||
return "maci64";
|
||
else if (mexext == "mexmaca64")
|
||
return "maca64";
|
||
else
|
||
{
|
||
cerr << "ERROR: 'mexext' option to preprocessor incorrectly set, needed with 'use_dll'" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
#ifdef __APPLE__
|
||
string
|
||
ModelTree::findGccOnMacos(const string &mexext)
|
||
{
|
||
const string macos_gcc_version {"12"}; // doc/manual/source/installation-and-configuration.rst
|
||
// should be updated when this is changed
|
||
char dynare_preprocessor_path[PATH_MAX];
|
||
uint32_t size = PATH_MAX;
|
||
string local_gcc_path;
|
||
if (_NSGetExecutablePath(dynare_preprocessor_path, &size) == 0)
|
||
{
|
||
string s = dynare_preprocessor_path;
|
||
local_gcc_path = s.substr(0, s.find_last_of("/")) + "/../.brew/bin/gcc-" + macos_gcc_version;
|
||
}
|
||
|
||
// if user did not choose to install gcc locally via the pkg-installer then we need to find GNU gcc
|
||
// homebrew binaries are located in /usr/local/bin/ on x86_64 systems and in /opt/homebrew/bin/ on arm64 systems
|
||
if (filesystem::exists(local_gcc_path))
|
||
return local_gcc_path;
|
||
else if (string global_gcc_path = "/usr/local/bin/gcc-" + macos_gcc_version;
|
||
filesystem::exists(global_gcc_path) && mexext == "mexmaci64")
|
||
return global_gcc_path;
|
||
else if (string global_gcc_path = "/opt/homebrew/bin/gcc-" + macos_gcc_version;
|
||
filesystem::exists(global_gcc_path) && mexext == "mexmaca64")
|
||
return global_gcc_path;
|
||
else
|
||
{
|
||
cerr << "ERROR: You must install gcc-" << macos_gcc_version
|
||
<< " on your system before using the `use_dll` option of Dynare. "
|
||
<< "If using MATLAB, you can do this via the Dynare installation package. If using Octave, you should run `brew install gcc-" << macos_gcc_version << "` in a terminal." << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
#endif
|
||
|
||
void
|
||
ModelTree::compileMEX(const string &basename, const string &funcname, const string &mexext, const vector<filesystem::path> &src_files, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const
|
||
{
|
||
const string opt_flags = "-O3 -g0 --param ira-max-conflict-table-size=1 -fno-forward-propagate -fno-gcse -fno-dce -fno-dse -fno-tree-fre -fno-tree-pre -fno-tree-cselim -fno-tree-dse -fno-tree-dce -fno-tree-pta -fno-gcse-after-reload";
|
||
|
||
filesystem::path compiler;
|
||
ostringstream flags;
|
||
string libs;
|
||
|
||
if (matlabroot.empty())
|
||
{
|
||
cerr << "ERROR: 'matlabroot' option to preprocessor is not set, needed with 'use_dll'" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
|
||
if (mexext == "mex")
|
||
{
|
||
// Octave
|
||
compiler = matlabroot / "bin" / "mkoctfile";
|
||
flags << "--mex";
|
||
#ifdef __APPLE__
|
||
/* On macOS, enforce GCC, otherwise Clang will be used, and it does not
|
||
accept our custom optimization flags (see dynare#1797) */
|
||
string gcc_path = findGccOnMacos(mexext);
|
||
if (setenv("CC", gcc_path.c_str(), 1) != 0)
|
||
{
|
||
cerr << "Can't set CC environment variable" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
// We also define CXX, because that is used for linking
|
||
if (setenv("CXX", gcc_path.c_str(), 1) != 0)
|
||
{
|
||
cerr << "Can't set CXX environment variable" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
#endif
|
||
}
|
||
else
|
||
{
|
||
// MATLAB
|
||
compiler = "gcc";
|
||
string arch = matlab_arch(mexext);
|
||
auto include_dir = matlabroot / "extern" / "include";
|
||
flags << "-I " << include_dir;
|
||
auto bin_dir = matlabroot / "bin" / arch;
|
||
flags << " -L " << bin_dir;
|
||
flags << " -fexceptions -DNDEBUG";
|
||
libs = "-lmex -lmx";
|
||
if (mexext == "mexa64")
|
||
{
|
||
// GNU/Linux
|
||
flags << " -D_GNU_SOURCE -fPIC -pthread"
|
||
<< " -shared -Wl,--no-undefined -Wl,-rpath-link," << bin_dir;
|
||
libs += " -lm";
|
||
}
|
||
else if (mexext == "mexw64")
|
||
{
|
||
// Windows
|
||
flags << " -static-libgcc -shared";
|
||
// Put the MinGW environment shipped with Dynare in the path
|
||
auto mingwpath = dynareroot / "mingw64" / "bin";
|
||
string newpath = "PATH=" + mingwpath.string() + ';' + getenv("PATH");
|
||
/* We can’t use setenv() since it is not available on MinGW. Note
|
||
that putenv() seems to make a copy of the string on MinGW, contrary
|
||
to what is done on GNU/Linux and macOS. */
|
||
if (putenv(const_cast<char *>(newpath.c_str())) != 0)
|
||
{
|
||
cerr << "Can't set PATH" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
#ifdef __APPLE__
|
||
else if (mexext == "mexmaci64" || mexext == "mexmaca64")
|
||
{
|
||
compiler = findGccOnMacos(mexext);
|
||
flags << " -fno-common -Wl,-twolevel_namespace -undefined error -bundle";
|
||
libs += " -lm";
|
||
}
|
||
#endif
|
||
else
|
||
{
|
||
cerr << "ERROR: unsupported value '" << mexext << "' for 'mexext' option" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
filesystem::path mex_dir{"+" + basename};
|
||
filesystem::path binary{mex_dir / (funcname + "." + mexext)};
|
||
|
||
ostringstream cmd;
|
||
|
||
#ifdef _WIN32
|
||
/* On Windows, system() hands the command over to "cmd.exe /C". We need to
|
||
enclose the whole command line within double quotes if we want the inner
|
||
quotes to be correctly handled. See "cmd /?" for more details. */
|
||
cmd << '"';
|
||
#endif
|
||
|
||
if (user_set_compiler.empty())
|
||
cmd << compiler << " ";
|
||
else
|
||
if (!filesystem::exists(user_set_compiler))
|
||
{
|
||
cerr << "Error: The specified compiler '" << user_set_compiler << "' cannot be found on your system" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
else
|
||
cmd << user_set_compiler << " ";
|
||
|
||
if (user_set_subst_flags.empty())
|
||
cmd << opt_flags << " " << flags.str() << " ";
|
||
else
|
||
cmd << user_set_subst_flags << " ";
|
||
|
||
if (!user_set_add_flags.empty())
|
||
cmd << user_set_add_flags << " ";
|
||
|
||
for (auto &src : src_files)
|
||
cmd << src << " ";
|
||
cmd << "-o " << binary << " ";
|
||
|
||
if (user_set_subst_libs.empty())
|
||
cmd << libs;
|
||
else
|
||
cmd << user_set_subst_libs;
|
||
|
||
if (!user_set_add_libs.empty())
|
||
cmd << " " << user_set_add_libs;
|
||
|
||
#ifdef _WIN32
|
||
cmd << '"';
|
||
#endif
|
||
|
||
cout << "Compiling " << funcname << " MEX..." << endl << cmd.str() << endl;
|
||
|
||
if (system(cmd.str().c_str()))
|
||
{
|
||
cerr << "Compilation failed" << endl;
|
||
exit(EXIT_FAILURE);
|
||
}
|
||
}
|
||
|
||
void
|
||
ModelTree::reorderAuxiliaryEquations()
|
||
{
|
||
using namespace boost;
|
||
|
||
// Create the mapping between auxiliary variables and auxiliary equations
|
||
int n = static_cast<int>(aux_equations.size());
|
||
map<int, int> auxEndoToEq;
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
auto varexpr = dynamic_cast<VariableNode *>(aux_equations[i]->arg1);
|
||
assert(varexpr && symbol_table.getType(varexpr->symb_id) == SymbolType::endogenous);
|
||
auxEndoToEq[varexpr->symb_id] = i;
|
||
}
|
||
assert(static_cast<int>(auxEndoToEq.size()) == n);
|
||
|
||
/* Construct the directed acyclic graph where auxiliary equations are
|
||
vertices and edges represent dependency relationships. */
|
||
using Graph = adjacency_list<vecS, vecS, directedS>;
|
||
Graph g(n);
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
set<int> endos;
|
||
aux_equations[i]->collectVariables(SymbolType::endogenous, endos);
|
||
for (int endo : endos)
|
||
if (auto it = auxEndoToEq.find(endo);
|
||
it != auxEndoToEq.end() && it->second != i)
|
||
add_edge(i, it->second, g);
|
||
}
|
||
|
||
// Topological sort of the graph
|
||
using Vertex = graph_traits<Graph>::vertex_descriptor;
|
||
vector<Vertex> ordered;
|
||
topological_sort(g, back_inserter(ordered));
|
||
|
||
// Reorder auxiliary equations accordingly
|
||
auto aux_equations_old = aux_equations;
|
||
auto index = get(vertex_index, g); // Maps vertex descriptors to their index
|
||
for (int i = 0; i < n; i++)
|
||
aux_equations[i] = aux_equations_old[index[ordered[i]]];
|
||
}
|
||
|
||
map<tuple<int, int, int>, expr_t>
|
||
ModelTree::collectFirstOrderDerivativesEndogenous()
|
||
{
|
||
map<tuple<int, int, int>, expr_t> endo_derivatives;
|
||
for (auto &[indices, d1] : derivatives[1])
|
||
if (getTypeByDerivID(indices[1]) == SymbolType::endogenous)
|
||
{
|
||
int eq = indices[0];
|
||
int var { getTypeSpecificIDByDerivID(indices[1]) };
|
||
int lag = getLagByDerivID(indices[1]);
|
||
endo_derivatives[{ eq, var, lag }] = d1;
|
||
}
|
||
return endo_derivatives;
|
||
}
|
||
|
||
ModelTree::jacob_map_t
|
||
ModelTree::computeSymbolicJacobian() const
|
||
{
|
||
jacob_map_t symbolic_jacobian;
|
||
for (int i = 0; i < static_cast<int>(equations.size()); i++)
|
||
{
|
||
set<pair<int, int>> endos_and_lags;
|
||
equations[i]->collectEndogenous(endos_and_lags);
|
||
for (const auto &[endo, lag] : endos_and_lags)
|
||
symbolic_jacobian[{ i, endo }] = 1;
|
||
}
|
||
return symbolic_jacobian;
|
||
}
|
||
|
||
void
|
||
ModelTree::updateReverseVariableEquationOrderings()
|
||
{
|
||
int n = equations.size();
|
||
eq_idx_orig2block.resize(n);
|
||
endo_idx_orig2block.resize(n);
|
||
for (int i = 0; i < n; i++)
|
||
{
|
||
endo_idx_orig2block[endo_idx_block2orig[i]] = i;
|
||
eq_idx_orig2block[eq_idx_block2orig[i]] = i;
|
||
}
|
||
}
|
||
|
||
expr_t
|
||
ModelTree::getRHSFromLHS(expr_t lhs) const
|
||
{
|
||
for (auto eq : equations)
|
||
if (eq->arg1 == lhs)
|
||
return eq->arg2;
|
||
throw ExprNode::MatchFailureException{"Cannot find an equation with the requested LHS"};
|
||
}
|
||
|
||
void
|
||
ModelTree::writeBlockBytecodeAdditionalDerivatives([[maybe_unused]] BytecodeWriter &code_file,
|
||
[[maybe_unused]] int block,
|
||
[[maybe_unused]] const temporary_terms_t &temporary_terms_union,
|
||
[[maybe_unused]] const deriv_node_temp_terms_t &tef_terms) const
|
||
{
|
||
}
|