2008-02-03 11:28:36 +01:00
|
|
|
|
/*
|
2023-01-04 16:18:53 +01:00
|
|
|
|
* Copyright © 2003-2023 Dynare Team
|
2008-02-03 11:28:36 +01:00
|
|
|
|
*
|
|
|
|
|
* 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
|
2021-06-09 16:52:20 +02:00
|
|
|
|
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
2008-02-03 11:28:36 +01:00
|
|
|
|
*/
|
|
|
|
|
|
|
|
|
|
#include "ModelTree.hh"
|
2020-04-09 17:44:17 +02:00
|
|
|
|
#include "VariableDependencyGraph.hh"
|
2022-06-16 17:52:14 +02:00
|
|
|
|
|
2019-04-23 14:51:14 +02:00
|
|
|
|
#pragma GCC diagnostic push
|
|
|
|
|
#pragma GCC diagnostic ignored "-Wold-style-cast"
|
2019-09-11 16:24:09 +02:00
|
|
|
|
#pragma GCC diagnostic ignored "-Wsign-compare"
|
|
|
|
|
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
|
2009-12-16 14:21:31 +01:00
|
|
|
|
#include <boost/graph/adjacency_list.hpp>
|
|
|
|
|
#include <boost/graph/max_cardinality_matching.hpp>
|
|
|
|
|
#include <boost/graph/topological_sort.hpp>
|
2019-04-23 14:51:14 +02:00
|
|
|
|
#pragma GCC diagnostic pop
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2019-11-04 15:50:26 +01:00
|
|
|
|
#ifdef __APPLE__
|
2019-12-20 16:59:30 +01:00
|
|
|
|
# include <mach-o/dyld.h>
|
2019-11-04 15:50:26 +01:00
|
|
|
|
#endif
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
#include <algorithm>
|
2020-01-07 12:41:27 +01:00
|
|
|
|
#include <regex>
|
2020-03-24 16:57:45 +01:00
|
|
|
|
#include <utility>
|
2020-01-07 12:41:27 +01:00
|
|
|
|
|
2022-10-17 13:57:34 +02:00
|
|
|
|
/* NB: The workers must be listed *after* all the other static variables
|
|
|
|
|
related to MEX compilation, so that when the preprocessor exits, the workers
|
|
|
|
|
are destroyed *before* those variables (since the former rely on the latter
|
|
|
|
|
for their functioning). */
|
|
|
|
|
condition_variable_any ModelTree::mex_compilation_cv;
|
2022-10-05 16:38:16 +02:00
|
|
|
|
mutex ModelTree::mex_compilation_mut;
|
2022-10-14 14:52:11 +02:00
|
|
|
|
vector<tuple<filesystem::path, set<filesystem::path>, string>> ModelTree::mex_compilation_queue;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
set<filesystem::path> ModelTree::mex_compilation_ongoing, ModelTree::mex_compilation_done,
|
|
|
|
|
ModelTree::mex_compilation_failed;
|
2022-10-17 13:57:34 +02:00
|
|
|
|
vector<jthread> ModelTree::mex_compilation_workers;
|
2022-09-30 15:39:41 +02:00
|
|
|
|
|
2018-10-09 18:27:19 +02:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::copyHelper(const ModelTree& m)
|
2018-10-09 18:27:19 +02:00
|
|
|
|
{
|
2018-10-10 13:07:25 +02:00
|
|
|
|
auto f = [this](expr_t e) { return e->clone(*this); };
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
|
|
|
|
// Equations
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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](const 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;
|
|
|
|
|
};
|
2018-11-15 16:39:53 +01:00
|
|
|
|
|
2018-10-09 18:27:19 +02:00
|
|
|
|
// Derivatives
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.derivatives)
|
2018-11-15 16:39:53 +01:00
|
|
|
|
derivatives.push_back(convert_deriv_map(it));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.params_derivatives)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
params_derivatives.emplace(it.first, convert_deriv_map(it.second));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.jacobian_sparse_column_major_order)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
jacobian_sparse_column_major_order.emplace(it.first, f(it.second));
|
2018-11-15 16:39:53 +01:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto convert_temporary_terms_t = [f](const temporary_terms_t& tt) {
|
|
|
|
|
temporary_terms_t tt2;
|
|
|
|
|
for (const auto& it : tt)
|
|
|
|
|
tt2.insert(f(it));
|
|
|
|
|
return tt2;
|
|
|
|
|
};
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
|
|
|
|
// Temporary terms
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.temporary_terms_derivatives)
|
2018-11-15 16:39:53 +01:00
|
|
|
|
temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.temporary_terms_idxs)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
temporary_terms_idxs.emplace(f(it.first), it.second);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.params_derivs_temporary_terms)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
params_derivs_temporary_terms.emplace(it.first, convert_temporary_terms_t(it.second));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.params_derivs_temporary_terms_idxs)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
params_derivs_temporary_terms_idxs.emplace(f(it.first), it.second);
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
|
|
|
|
// Other stuff
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.trend_symbols_map)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
trend_symbols_map.emplace(it.first, f(it.second));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.nonstationary_symbols_map)
|
|
|
|
|
nonstationary_symbols_map.emplace(it.first, pair {it.second.first, f(it.second.second)});
|
2020-03-17 18:58:34 +01:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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)));
|
2020-03-17 18:58:34 +01:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.blocks_derivatives)
|
2020-03-17 18:58:34 +01:00
|
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
|
map<tuple<int, int, int>, expr_t> v;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it2 : it)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
v.emplace(it2.first, f(it2.second));
|
2020-03-17 18:58:34 +01:00
|
|
|
|
blocks_derivatives.push_back(v);
|
|
|
|
|
}
|
2020-05-06 17:13:47 +02:00
|
|
|
|
|
2023-12-13 10:28:54 +01:00
|
|
|
|
auto convert_vector_tt = [f](const vector<temporary_terms_t>& vtt) {
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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)
|
2020-05-13 16:58:19 +02:00
|
|
|
|
blocks_temporary_terms.push_back(convert_vector_tt(it));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.blocks_temporary_terms_idxs)
|
2022-09-21 15:13:41 +02:00
|
|
|
|
blocks_temporary_terms_idxs.emplace(f(it.first), it.second);
|
2022-09-28 19:18:15 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : m.blocks_jacobian_sparse_column_major_order)
|
2022-09-28 19:18:15 +02:00
|
|
|
|
{
|
|
|
|
|
map<pair<int, int>, expr_t, columnMajorOrderLess> v;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it2 : it)
|
2022-09-28 19:18:15 +02:00
|
|
|
|
v.emplace(it2.first, f(it2.second));
|
|
|
|
|
blocks_jacobian_sparse_column_major_order.push_back(v);
|
|
|
|
|
}
|
2018-10-09 18:27:19 +02:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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)
|
2018-11-22 14:36:03 +01:00
|
|
|
|
{
|
2022-09-26 17:16:22 +02:00
|
|
|
|
// Ensure that elements accessed by writeParamsDerivativesFileHelper() exist
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& ord :
|
|
|
|
|
{pair {0, 1}, pair {1, 1}, pair {0, 2}, pair {1, 2}, pair {2, 1}, pair {3, 1}})
|
2023-03-23 12:58:42 +01:00
|
|
|
|
params_derivatives.try_emplace(ord);
|
2018-11-22 14:36:03 +01:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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},
|
|
|
|
|
jacobian_sparse_colptr {m.jacobian_sparse_colptr},
|
|
|
|
|
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},
|
|
|
|
|
block_decomposed {m.block_decomposed},
|
|
|
|
|
time_recursive_block_decomposition {m.time_recursive_block_decomposition},
|
|
|
|
|
blocks {m.blocks},
|
|
|
|
|
endo2block {m.endo2block},
|
|
|
|
|
eq2block {m.eq2block},
|
|
|
|
|
blocks_jacobian_sparse_colptr {m.blocks_jacobian_sparse_colptr},
|
|
|
|
|
endo2eq {m.endo2eq},
|
|
|
|
|
cutoff {m.cutoff}
|
2018-10-09 18:27:19 +02:00
|
|
|
|
{
|
|
|
|
|
copyHelper(m);
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree&
|
|
|
|
|
ModelTree::operator=(const ModelTree& m)
|
2018-10-09 18:27:19 +02:00
|
|
|
|
{
|
|
|
|
|
DataTree::operator=(m);
|
|
|
|
|
|
|
|
|
|
equations.clear();
|
|
|
|
|
equations_lineno = m.equations_lineno;
|
|
|
|
|
aux_equations.clear();
|
|
|
|
|
equation_tags = m.equation_tags;
|
2020-01-20 17:22:32 +01:00
|
|
|
|
computed_derivs_order = m.computed_derivs_order;
|
2018-10-09 18:27:19 +02:00
|
|
|
|
NNZDerivatives = m.NNZDerivatives;
|
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
derivatives.clear();
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
|
|
|
|
jacobian_sparse_column_major_order.clear();
|
|
|
|
|
jacobian_sparse_colptr = m.jacobian_sparse_colptr;
|
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
params_derivatives.clear();
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
temporary_terms_derivatives.clear();
|
2018-10-09 18:27:19 +02:00
|
|
|
|
params_derivs_temporary_terms.clear();
|
|
|
|
|
params_derivs_temporary_terms_idxs.clear();
|
|
|
|
|
|
|
|
|
|
trend_symbols_map.clear();
|
|
|
|
|
nonstationary_symbols_map.clear();
|
|
|
|
|
|
2020-04-17 14:55:55 +02:00
|
|
|
|
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;
|
2020-03-17 18:58:34 +01:00
|
|
|
|
equation_type_and_normalized_equation.clear();
|
|
|
|
|
blocks_derivatives.clear();
|
2022-11-29 12:38:11 +01:00
|
|
|
|
block_decomposed = m.block_decomposed;
|
2022-11-30 14:43:44 +01:00
|
|
|
|
time_recursive_block_decomposition = m.time_recursive_block_decomposition;
|
2020-04-21 18:10:46 +02:00
|
|
|
|
blocks = m.blocks;
|
2020-04-29 16:41:33 +02:00
|
|
|
|
endo2block = m.endo2block;
|
|
|
|
|
eq2block = m.eq2block;
|
2020-05-06 17:13:47 +02:00
|
|
|
|
blocks_temporary_terms.clear();
|
|
|
|
|
blocks_temporary_terms_idxs.clear();
|
2022-09-28 19:18:15 +02:00
|
|
|
|
blocks_jacobian_sparse_column_major_order.clear();
|
|
|
|
|
blocks_jacobian_sparse_colptr = m.blocks_jacobian_sparse_colptr;
|
2018-10-09 18:27:19 +02:00
|
|
|
|
endo2eq = m.endo2eq;
|
|
|
|
|
cutoff = m.cutoff;
|
|
|
|
|
|
2019-12-02 19:21:14 +01:00
|
|
|
|
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;
|
|
|
|
|
|
2018-10-09 18:27:19 +02:00
|
|
|
|
copyHelper(m);
|
|
|
|
|
|
|
|
|
|
return *this;
|
|
|
|
|
}
|
|
|
|
|
|
2023-04-24 17:47:29 +02:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::computeNormalization(const jacob_map_t& contemporaneous_jacobian)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2023-04-24 15:05:22 +02:00
|
|
|
|
const int n {static_cast<int>(equations.size())};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
|
|
assert(n == symbol_table.endo_nbr());
|
|
|
|
|
|
2018-06-05 14:56:34 +02:00
|
|
|
|
using BipartiteGraph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
|
2023-04-24 15:05:22 +02:00
|
|
|
|
using vertex_descriptor_t = boost::graph_traits<BipartiteGraph>::vertex_descriptor;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
|
|
/*
|
|
|
|
|
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
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [eq_and_endo, val] : contemporaneous_jacobian)
|
2020-03-20 18:42:59 +01:00
|
|
|
|
add_edge(eq_and_endo.first + n, eq_and_endo.second, g);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
|
|
// Compute maximum cardinality matching
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<vertex_descriptor_t> mate_map(2 * n);
|
2023-04-24 15:05:22 +02:00
|
|
|
|
edmonds_maximum_cardinality_matching(g, &mate_map[0]);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
|
|
// Check if all variables are normalized
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (auto it = find(mate_map.begin(), mate_map.begin() + n,
|
|
|
|
|
boost::graph_traits<BipartiteGraph>::null_vertex());
|
2019-12-16 19:42:59 +01:00
|
|
|
|
it != mate_map.begin() + n)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
throw ModelNormalizationFailed {
|
|
|
|
|
symbol_table.getName(symbol_table.getID(SymbolType::endogenous, it - mate_map.begin()))};
|
2023-04-24 15:05:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to
|
|
|
|
|
// them
|
2023-04-24 15:05:22 +02:00
|
|
|
|
endo2eq.resize(equations.size());
|
2023-11-30 15:28:57 +01:00
|
|
|
|
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(),
|
|
|
|
|
[=](vertex_descriptor_t i) { return i - n; });
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2022-09-13 15:43:13 +02:00
|
|
|
|
bool
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::computeNonSingularNormalization(const eval_context_t& eval_context)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2023-04-24 15:05:22 +02:00
|
|
|
|
const int n {static_cast<int>(equations.size())};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2022-09-13 15:43:13 +02:00
|
|
|
|
/* Optimal policy models (discretionary, or Ramsey before computing FOCs) do
|
|
|
|
|
not have as many equations as variables. */
|
|
|
|
|
if (n != symbol_table.endo_nbr())
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
cout << "The " << modelClassName()
|
|
|
|
|
<< " cannot be normalized, since it does not have as many equations as variables."
|
|
|
|
|
<< endl;
|
2022-09-13 15:43:13 +02:00
|
|
|
|
return false;
|
|
|
|
|
}
|
|
|
|
|
|
2022-09-14 17:48:33 +02:00
|
|
|
|
cout << "Normalizing the " << modelClassName() << "..." << endl;
|
2022-09-13 15:43:13 +02:00
|
|
|
|
|
2023-04-24 17:49:54 +02:00
|
|
|
|
/* If the model is purely backward, determine whether all original equations
|
|
|
|
|
have a single contemporaneous endogenous on the LHS. If this is the case,
|
|
|
|
|
then first try a normalization by enforcing that each original equation is
|
|
|
|
|
matched with the endogenous on the LHS. */
|
|
|
|
|
if (time_recursive_block_decomposition)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [normalize_by_lhs, lhs_symbolic_jacobian] {computeLeftHandSideSymbolicJacobian()};
|
2023-04-24 17:49:54 +02:00
|
|
|
|
if (normalize_by_lhs)
|
|
|
|
|
try
|
|
|
|
|
{
|
|
|
|
|
computeNormalization(lhs_symbolic_jacobian);
|
|
|
|
|
return true;
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ModelNormalizationFailed& e)
|
2023-04-24 17:49:54 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
cerr << "WARNING: All equations are written so that a single contemporaneous "
|
|
|
|
|
"endogenous variable appears on the left-hand side. This suggests a natural "
|
|
|
|
|
"normalization of the model. However, variable "
|
|
|
|
|
<< e.unmatched_endo
|
|
|
|
|
<< " could not be matched with an equation. Check whether this is desired."
|
|
|
|
|
<< endl;
|
2023-04-24 17:49:54 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2023-04-24 15:05:22 +02:00
|
|
|
|
auto contemporaneous_jacobian {evaluateAndReduceJacobian(eval_context)};
|
|
|
|
|
|
2020-03-20 18:42:59 +01:00
|
|
|
|
// Compute the maximum value of each row of the contemporaneous Jacobian matrix
|
2022-06-03 16:21:30 +02:00
|
|
|
|
vector max_val(n, 0.0);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [eq_and_endo, val] : contemporaneous_jacobian)
|
2020-03-20 18:42:59 +01:00
|
|
|
|
max_val[eq_and_endo.first] = max(max_val[eq_and_endo.first], fabs(val));
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2020-03-24 16:57:45 +01:00
|
|
|
|
// Compute normalized contemporaneous Jacobian
|
|
|
|
|
jacob_map_t normalized_contemporaneous_jacobian(contemporaneous_jacobian);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& [eq_and_endo, val] : normalized_contemporaneous_jacobian)
|
2020-03-20 18:42:59 +01:00
|
|
|
|
val /= max_val[eq_and_endo.first];
|
2009-12-21 11:29:21 +01:00
|
|
|
|
|
2020-03-20 18:42:59 +01:00
|
|
|
|
// We start with the highest value of the cutoff and try to normalize the model
|
2009-12-21 11:29:21 +01:00
|
|
|
|
double current_cutoff = 0.99999999;
|
2020-03-20 18:42:59 +01:00
|
|
|
|
const double cutoff_lower_limit = 1e-19;
|
2009-12-21 11:29:21 +01:00
|
|
|
|
|
2020-03-24 16:57:45 +01:00
|
|
|
|
int last_suppressed = 0;
|
2023-04-24 17:47:29 +02:00
|
|
|
|
while (current_cutoff > cutoff_lower_limit)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-03-24 16:57:45 +01:00
|
|
|
|
// Drop elements below cutoff from normalized contemporaneous Jacobian
|
|
|
|
|
jacob_map_t normalized_contemporaneous_jacobian_above_cutoff;
|
|
|
|
|
int suppressed = 0;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [eq_and_endo, val] : normalized_contemporaneous_jacobian)
|
2020-03-20 18:42:59 +01:00
|
|
|
|
if (fabs(val) > max(current_cutoff, cutoff))
|
2020-03-24 16:57:45 +01:00
|
|
|
|
normalized_contemporaneous_jacobian_above_cutoff[eq_and_endo] = val;
|
2009-12-21 11:29:21 +01:00
|
|
|
|
else
|
2020-03-24 16:57:45 +01:00
|
|
|
|
suppressed++;
|
2009-12-21 11:29:21 +01:00
|
|
|
|
|
2020-03-24 16:57:45 +01:00
|
|
|
|
if (suppressed != last_suppressed)
|
2023-04-24 17:47:29 +02:00
|
|
|
|
try
|
|
|
|
|
{
|
|
|
|
|
computeNormalization(normalized_contemporaneous_jacobian_above_cutoff);
|
|
|
|
|
return true;
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ModelNormalizationFailed& e)
|
2023-04-24 17:47:29 +02:00
|
|
|
|
{
|
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
last_suppressed = suppressed;
|
2023-04-24 17:47:29 +02:00
|
|
|
|
current_cutoff /= 2;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2023-04-24 17:47:29 +02:00
|
|
|
|
// Try to normalize with the complete numerical jacobian
|
|
|
|
|
try
|
|
|
|
|
{
|
|
|
|
|
computeNormalization(normalized_contemporaneous_jacobian);
|
|
|
|
|
return true;
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ModelNormalizationFailed& e)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
|
|
|
|
}
|
2009-12-21 11:29:21 +01:00
|
|
|
|
|
2023-04-24 17:47:29 +02:00
|
|
|
|
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. */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto symbolic_jacobian {computeSymbolicJacobian(true)};
|
2023-04-24 17:47:29 +02:00
|
|
|
|
try
|
|
|
|
|
{
|
|
|
|
|
computeNormalization(symbolic_jacobian);
|
|
|
|
|
return true;
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ModelNormalizationFailed& e)
|
2023-04-24 17:47:29 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
cerr << "Could not normalize the " << modelClassName() << ". Variable " << e.unmatched_endo
|
|
|
|
|
<< " is not in the maximum cardinality matching." << endl;
|
2023-04-24 17:47:29 +02:00
|
|
|
|
}
|
2022-10-25 16:38:53 +02:00
|
|
|
|
|
2023-04-24 17:47:29 +02:00
|
|
|
|
return false;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-05-06 18:09:44 +02:00
|
|
|
|
ModelTree::jacob_map_t
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::evaluateAndReduceJacobian(const eval_context_t& eval_context) const
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-05-06 18:09:44 +02:00
|
|
|
|
jacob_map_t contemporaneous_jacobian;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : derivatives[1])
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2019-12-16 19:42:59 +01:00
|
|
|
|
int deriv_id = indices[1];
|
2018-07-17 18:34:07 +02:00
|
|
|
|
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2019-12-16 19:42:59 +01:00
|
|
|
|
int eq = indices[0];
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var {getTypeSpecificIDByDerivID(deriv_id)};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
int lag = getLagByDerivID(deriv_id);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
double val {[&] {
|
2022-09-13 15:43:13 +02:00
|
|
|
|
try
|
|
|
|
|
{
|
2023-12-04 16:37:00 +01:00
|
|
|
|
return d1->eval(eval_context); // NOLINT(clang-analyzer-core.NullDereference)
|
2022-09-13 15:43:13 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ExprNode::EvalExternalFunctionException& e)
|
2022-09-13 15:43:13 +02:00
|
|
|
|
{
|
|
|
|
|
return 1.0;
|
|
|
|
|
}
|
|
|
|
|
/* Other types of EvalException should not happen (all symbols should
|
|
|
|
|
have a value; we don’t evaluate an equal sign) */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
}()};
|
2022-09-13 15:43:13 +02:00
|
|
|
|
|
2020-05-07 16:51:49 +02:00
|
|
|
|
if ((isnan(val) || fabs(val) >= cutoff) && lag == 0)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
contemporaneous_jacobian[{eq, var}] = val;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2020-05-06 18:09:44 +02:00
|
|
|
|
return contemporaneous_jacobian;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-04-30 12:48:16 +02:00
|
|
|
|
pair<int, int>
|
2020-04-28 18:00:26 +02:00
|
|
|
|
ModelTree::computePrologueAndEpilogue()
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-03-24 16:57:45 +01:00
|
|
|
|
const int n = equations.size();
|
|
|
|
|
|
2020-04-17 19:21:37 +02:00
|
|
|
|
/* Initialize “eq_idx_block2orig” and “endo_idx_block2orig” to the identity
|
2020-03-24 16:57:45 +01:00
|
|
|
|
permutation. */
|
2020-04-17 14:55:55 +02:00
|
|
|
|
eq_idx_block2orig.resize(n);
|
|
|
|
|
endo_idx_block2orig.resize(n);
|
2020-03-24 16:57:45 +01:00
|
|
|
|
for (int i = 0; i < n; i++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-17 14:55:55 +02:00
|
|
|
|
eq_idx_block2orig[i] = i;
|
2020-04-17 19:21:37 +02:00
|
|
|
|
endo_idx_block2orig[endo2eq[i]] = i;
|
2017-06-14 07:01:31 +02:00
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
|
|
|
|
|
/* Compute incidence matrix, equations in rows, variables in columns. Row
|
|
|
|
|
(resp. column) indices are to be interpreted according to
|
2020-04-17 14:55:55 +02:00
|
|
|
|
“eq_idx_block2orig” (resp. “endo_idx_block2orig”). Stored in row-major
|
2020-03-24 16:57:45 +01:00
|
|
|
|
order. */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector IM(n * n, false);
|
2020-04-28 18:00:26 +02:00
|
|
|
|
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)
|
2022-11-30 14:43:44 +01:00
|
|
|
|
if (!time_recursive_block_decomposition || lag == 0)
|
|
|
|
|
IM[i * n + endo2eq[endo]] = true;
|
2020-04-28 18:00:26 +02:00
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
|
|
|
|
|
bool something_has_been_done;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
// Find the prologue equations and place first the AR(1) shock equations first
|
2020-04-30 12:48:16 +02:00
|
|
|
|
int prologue = 0;
|
2020-03-24 16:57:45 +01:00
|
|
|
|
do
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
|
|
|
|
something_has_been_done = false;
|
2020-03-24 16:57:45 +01:00
|
|
|
|
int new_prologue = prologue;
|
2009-12-16 18:13:23 +01:00
|
|
|
|
for (int i = prologue; i < n; i++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
|
|
|
|
int nze = 0;
|
2020-03-24 16:57:45 +01:00
|
|
|
|
int k = 0;
|
|
|
|
|
for (int j = new_prologue; j < n; j++)
|
2009-12-16 18:13:23 +01:00
|
|
|
|
if (IM[i * n + j])
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2009-12-16 18:13:23 +01:00
|
|
|
|
nze++;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
k = j;
|
|
|
|
|
}
|
2009-12-16 18:13:23 +01:00
|
|
|
|
if (nze == 1)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-03-24 16:57:45 +01:00
|
|
|
|
// Swap equations indexed by “new_prologue” and i
|
2009-12-16 18:13:23 +01:00
|
|
|
|
for (int j = 0; j < n; j++)
|
2020-03-24 16:57:45 +01:00
|
|
|
|
swap(IM[new_prologue * n + j], IM[i * n + j]);
|
2020-04-17 14:55:55 +02:00
|
|
|
|
swap(eq_idx_block2orig[new_prologue], eq_idx_block2orig[i]);
|
2020-03-24 16:57:45 +01:00
|
|
|
|
|
|
|
|
|
// Swap variables indexed by “new_prologue” and k (in the matching)
|
2009-12-16 18:13:23 +01:00
|
|
|
|
for (int j = 0; j < n; j++)
|
2020-03-24 16:57:45 +01:00
|
|
|
|
swap(IM[j * n + new_prologue], IM[j * n + k]);
|
2020-04-17 14:55:55 +02:00
|
|
|
|
swap(endo_idx_block2orig[new_prologue], endo_idx_block2orig[k]);
|
2020-03-24 16:57:45 +01:00
|
|
|
|
|
|
|
|
|
new_prologue++;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
something_has_been_done = true;
|
|
|
|
|
}
|
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
prologue = new_prologue;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
while (something_has_been_done);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
|
// Find the epilogue equations
|
2020-04-30 12:48:16 +02:00
|
|
|
|
int epilogue = 0;
|
2020-03-24 16:57:45 +01:00
|
|
|
|
do
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
|
|
|
|
something_has_been_done = false;
|
2020-03-24 16:57:45 +01:00
|
|
|
|
int new_epilogue = epilogue;
|
2020-03-24 18:26:06 +01:00
|
|
|
|
for (int i = prologue; i < n - epilogue; i++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
|
|
|
|
int nze = 0;
|
2020-03-24 16:57:45 +01:00
|
|
|
|
int k = 0;
|
|
|
|
|
for (int j = prologue; j < n - new_epilogue; j++)
|
2009-12-16 18:13:23 +01:00
|
|
|
|
if (IM[j * n + i])
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2009-12-16 18:13:23 +01:00
|
|
|
|
nze++;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
k = j;
|
|
|
|
|
}
|
2009-12-16 18:13:23 +01:00
|
|
|
|
if (nze == 1)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2009-12-16 18:13:23 +01:00
|
|
|
|
for (int j = 0; j < n; j++)
|
2020-03-24 16:57:45 +01:00
|
|
|
|
swap(IM[(n - 1 - new_epilogue) * n + j], IM[k * n + j]);
|
2020-04-17 14:55:55 +02:00
|
|
|
|
swap(eq_idx_block2orig[n - 1 - new_epilogue], eq_idx_block2orig[k]);
|
2020-03-24 16:57:45 +01:00
|
|
|
|
|
2009-12-16 18:13:23 +01:00
|
|
|
|
for (int j = 0; j < n; j++)
|
2020-03-24 16:57:45 +01:00
|
|
|
|
swap(IM[j * n + n - 1 - new_epilogue], IM[j * n + i]);
|
2020-04-17 14:55:55 +02:00
|
|
|
|
swap(endo_idx_block2orig[n - 1 - new_epilogue], endo_idx_block2orig[i]);
|
2020-03-24 16:57:45 +01:00
|
|
|
|
|
|
|
|
|
new_epilogue++;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
something_has_been_done = true;
|
|
|
|
|
}
|
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
epilogue = new_epilogue;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
2020-03-24 16:57:45 +01:00
|
|
|
|
while (something_has_been_done);
|
2020-04-15 17:56:28 +02:00
|
|
|
|
|
|
|
|
|
updateReverseVariableEquationOrderings();
|
2020-04-30 12:48:16 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {prologue, epilogue};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-03-20 15:23:23 +01:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::equationTypeDetermination(
|
|
|
|
|
const map<tuple<int, int, int>, expr_t>& first_order_endo_derivatives)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-03-20 15:23:23 +01:00
|
|
|
|
equation_type_and_normalized_equation.clear();
|
|
|
|
|
equation_type_and_normalized_equation.resize(equations.size());
|
2020-03-24 18:26:06 +01:00
|
|
|
|
for (int i = 0; i < static_cast<int>(equations.size()); i++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-17 14:55:55 +02:00
|
|
|
|
int eq = eq_idx_block2orig[i];
|
|
|
|
|
int var = endo_idx_block2orig[i];
|
2020-03-30 17:06:48 +02:00
|
|
|
|
expr_t lhs = equations[eq]->arg1;
|
|
|
|
|
EquationType Equation_Simulation_Type = EquationType::solve;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
BinaryOpNode* normalized_eq = nullptr;
|
|
|
|
|
if (auto it = first_order_endo_derivatives.find({eq, var, 0});
|
2020-03-30 17:06:48 +02:00
|
|
|
|
it != first_order_endo_derivatives.end())
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-03-30 17:06:48 +02:00
|
|
|
|
expr_t derivative = it->second;
|
|
|
|
|
// Determine whether the equation can be evaluated rather than solved
|
2020-04-17 14:55:55 +02:00
|
|
|
|
if (lhs->isVariableNodeEqualTo(SymbolType::endogenous, endo_idx_block2orig[i], 0)
|
2020-03-30 17:06:48 +02:00
|
|
|
|
&& derivative->isNumConstNodeEqualTo(1))
|
2020-03-20 18:00:56 +01:00
|
|
|
|
Equation_Simulation_Type = EquationType::evaluate;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
else
|
|
|
|
|
{
|
2020-03-30 17:06:48 +02:00
|
|
|
|
set<pair<int, int>> result;
|
|
|
|
|
derivative->collectEndogenous(result);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bool variable_not_in_derivative = !result.contains({var, 0});
|
2020-03-30 17:06:48 +02:00
|
|
|
|
|
2020-04-02 14:36:26 +02:00
|
|
|
|
try
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
normalized_eq = equations[eq]->normalizeEquation(
|
|
|
|
|
symbol_table.getID(SymbolType::endogenous, var), 0);
|
2023-10-16 16:24:31 +02:00
|
|
|
|
if ((getMFS() == 2 && variable_not_in_derivative) || getMFS() == 3)
|
2020-05-20 11:49:32 +02:00
|
|
|
|
Equation_Simulation_Type = EquationType::evaluateRenormalized;
|
2020-04-02 14:36:26 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ExprNode::NormalizationFailed& e)
|
2020-04-02 14:36:26 +02:00
|
|
|
|
{
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
equation_type_and_normalized_equation[eq] = {Equation_Simulation_Type, normalized_eq};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2020-04-30 11:15:55 +02:00
|
|
|
|
void
|
|
|
|
|
ModelTree::computeDynamicStructureOfBlock(int blk)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector max_endo_lag_lead(blocks[blk].size, pair {0, 0});
|
2020-04-30 11:15:55 +02:00
|
|
|
|
blocks[blk].max_endo_lag = blocks[blk].max_endo_lead = 0;
|
|
|
|
|
for (int eq = 0; eq < blocks[blk].size; eq++)
|
|
|
|
|
{
|
|
|
|
|
set<pair<int, int>> endos_and_lags;
|
|
|
|
|
expr_t e = getBlockEquationExpr(blk, eq);
|
|
|
|
|
|
2023-01-17 14:41:02 +01:00
|
|
|
|
/* Compute max lags/leads and per-variable structure for endos belonging
|
|
|
|
|
to this block */
|
2020-04-30 11:15:55 +02:00
|
|
|
|
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);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto& [max_endo_lag, max_endo_lead]
|
|
|
|
|
= max_endo_lag_lead[getBlockInitialVariableID(blk, endo)];
|
2020-04-30 11:15:55 +02:00
|
|
|
|
max_endo_lag = max(max_endo_lag, -lag);
|
|
|
|
|
max_endo_lead = max(max_endo_lead, lag);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2020-04-30 12:48:16 +02:00
|
|
|
|
void
|
|
|
|
|
ModelTree::computeSimulationTypeOfBlock(int blk)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto& type = blocks[blk].simulation_type;
|
2020-04-30 12:48:16 +02:00
|
|
|
|
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
|
2020-05-20 11:49:32 +02:00
|
|
|
|
|| getBlockEquationType(blk, 0) == EquationType::evaluateRenormalized);
|
2020-04-30 12:48:16 +02:00
|
|
|
|
if (blocks[blk].max_endo_lead > 0)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
type = can_eval ? BlockSimulationType::evaluateBackward
|
|
|
|
|
: BlockSimulationType::solveBackwardSimple;
|
2020-04-30 12:48:16 +02:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
type = can_eval ? BlockSimulationType::evaluateForward
|
|
|
|
|
: BlockSimulationType::solveForwardSimple;
|
2020-04-30 12:48:16 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2020-03-19 17:46:10 +01:00
|
|
|
|
pair<lag_lead_vector_t, lag_lead_vector_t>
|
2020-04-29 16:41:33 +02:00
|
|
|
|
ModelTree::getVariableLeadLagByBlock() const
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
|
|
|
|
int nb_endo = symbol_table.endo_nbr();
|
2020-04-17 19:21:37 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
lag_lead_vector_t variable_lag_lead(nb_endo, {0, 0}), equation_lag_lead(nb_endo, {0, 0});
|
2020-04-28 18:00:26 +02:00
|
|
|
|
for (int eq = 0; eq < nb_endo; eq++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-28 18:00:26 +02:00
|
|
|
|
set<pair<int, int>> endos_and_lags;
|
|
|
|
|
equations[eq]->collectEndogenous(endos_and_lags);
|
|
|
|
|
for (auto [endo, lag] : endos_and_lags)
|
2020-04-29 16:41:33 +02:00
|
|
|
|
if (endo2block[endo] == eq2block[eq])
|
2020-04-28 18:00:26 +02:00
|
|
|
|
{
|
|
|
|
|
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);
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {equation_lag_lead, variable_lag_lead};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-04-29 17:44:57 +02:00
|
|
|
|
void
|
2020-04-30 12:48:16 +02:00
|
|
|
|
ModelTree::computeBlockDecomposition(int prologue, int epilogue)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-17 19:21:37 +02:00
|
|
|
|
int nb_var = symbol_table.endo_nbr();
|
|
|
|
|
int nb_simvars = nb_var - prologue - epilogue;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2020-04-09 17:44:17 +02:00
|
|
|
|
/* Construct the graph representing the dependencies between all
|
2020-04-15 17:56:28 +02:00
|
|
|
|
variables that do not belong to the prologue or the epilogue.
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2020-04-28 18:00:26 +02:00
|
|
|
|
For detecting dependencies between variables, use the symbolic adjacency
|
2020-04-15 17:56:28 +02:00
|
|
|
|
matrix */
|
2020-04-17 19:21:37 +02:00
|
|
|
|
VariableDependencyGraph G(nb_simvars);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [key, value] : computeSymbolicJacobian(time_recursive_block_decomposition))
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-15 17:56:28 +02:00
|
|
|
|
auto [eq, endo] = key;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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
|
2020-04-15 17:56:28 +02:00
|
|
|
|
&& eq != endo2eq[endo])
|
2023-11-30 15:28:57 +01:00
|
|
|
|
add_edge(vertex(eq_idx_orig2block[endo2eq[endo]] - prologue, G),
|
|
|
|
|
vertex(eq_idx_orig2block[eq] - prologue, G), G);
|
2017-06-14 07:01:31 +02:00
|
|
|
|
}
|
2020-04-09 17:44:17 +02:00
|
|
|
|
|
2020-04-15 17:56:28 +02:00
|
|
|
|
/* Identify the simultaneous blocks. Each simultaneous block is given an
|
|
|
|
|
index, starting from 0, in recursive order */
|
2020-04-29 18:48:42 +02:00
|
|
|
|
auto [num_simblocks, simvar2simblock] = G.sortedStronglyConnectedComponents();
|
2020-04-15 17:56:28 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int num_blocks = prologue + num_simblocks + epilogue;
|
2020-04-21 18:10:46 +02:00
|
|
|
|
|
|
|
|
|
blocks.clear();
|
|
|
|
|
blocks.resize(num_blocks);
|
2020-04-29 16:41:33 +02:00
|
|
|
|
endo2block.resize(nb_var);
|
|
|
|
|
eq2block.resize(nb_var);
|
2020-04-21 18:10:46 +02:00
|
|
|
|
|
2020-04-29 16:41:33 +02:00
|
|
|
|
// Initialize size and mfs_size for prologue and epilogue, plus eq/endo→block mappings
|
2020-04-29 18:48:42 +02:00
|
|
|
|
for (int blk = 0; blk < num_blocks; blk++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (blk < prologue || blk >= num_blocks - epilogue)
|
2020-04-29 18:48:42 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var_eq = (blk < prologue ? blk : blk - num_simblocks + nb_simvars);
|
2020-04-29 18:48:42 +02:00
|
|
|
|
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;
|
|
|
|
|
}
|
2020-04-15 17:56:28 +02:00
|
|
|
|
|
2020-04-29 18:48:42 +02:00
|
|
|
|
// 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++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-29 18:48:42 +02:00
|
|
|
|
simblock2simvars[simvar2simblock[i]].push_back(i);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int blk = prologue + simvar2simblock[i];
|
2020-04-29 16:41:33 +02:00
|
|
|
|
blocks[blk].size++;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
endo2block[endo_idx_block2orig[prologue + i]] = blk;
|
|
|
|
|
eq2block[eq_idx_block2orig[prologue + i]] = blk;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-04-24 11:36:05 +02:00
|
|
|
|
// Determine the dynamic structure of each block
|
2020-04-29 16:41:33 +02:00
|
|
|
|
auto [equation_lag_lead, variable_lag_lead] = getVariableLeadLagByBlock();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2022-11-04 14:22:20 +01:00
|
|
|
|
/* For each simultaneous block, the minimum set of feedback variables is
|
|
|
|
|
computed. Then, the variables within the blocks are reordered so that
|
|
|
|
|
recursive (non-feedback) appear first, in recursive order. They are
|
|
|
|
|
followed by feedback variables, which are reordered according to their
|
|
|
|
|
dynamic status (static first, then backward, mixed and forward). */
|
2020-04-24 11:36:05 +02:00
|
|
|
|
|
|
|
|
|
/* 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 */
|
2020-04-17 19:21:37 +02:00
|
|
|
|
for (int i = 0; i < nb_simvars; i++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (equation_type_and_normalized_equation[eq_idx_block2orig[i + prologue]].first
|
|
|
|
|
== EquationType::solve
|
|
|
|
|
|| (!time_recursive_block_decomposition
|
|
|
|
|
&& (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))
|
2023-10-16 16:24:31 +02:00
|
|
|
|
|| getMFS() == 0)
|
2020-04-17 19:21:37 +02:00
|
|
|
|
add_edge(vertex(i, G), vertex(i, G), G);
|
2019-12-16 19:42:59 +01:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const vector<int> old_eq_idx_block2orig(eq_idx_block2orig),
|
|
|
|
|
old_endo_idx_block2orig(endo_idx_block2orig);
|
2020-04-15 17:56:28 +02:00
|
|
|
|
int ordidx = prologue;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int blk = prologue; blk < prologue + num_simblocks; blk++)
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
blocks[blk].first_equation
|
|
|
|
|
= (blk == 0 ? 0 : blocks[blk - 1].first_equation + blocks[blk - 1].size);
|
|
|
|
|
auto subG = G.extractSubgraph(simblock2simvars[blk - prologue]);
|
2020-04-17 19:21:37 +02:00
|
|
|
|
auto feed_back_vertices = subG.minimalSetOfFeedbackVertices();
|
2020-04-29 18:48:42 +02:00
|
|
|
|
blocks[blk].mfs_size = feed_back_vertices.size();
|
2020-04-29 18:46:15 +02:00
|
|
|
|
auto recursive_vertices = subG.reorderRecursiveVariables(feed_back_vertices);
|
|
|
|
|
auto v_index1 = get(boost::vertex_index1, subG);
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2022-11-04 14:22:20 +01:00
|
|
|
|
/* First the recursive variables conditional on feedback variables, in
|
|
|
|
|
recursive order */
|
|
|
|
|
for (int vtx : recursive_vertices)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int simvar {v_index1[vertex(vtx, subG)]};
|
|
|
|
|
eq_idx_block2orig[ordidx] = old_eq_idx_block2orig[simvar + prologue];
|
|
|
|
|
endo_idx_block2orig[ordidx] = old_endo_idx_block2orig[simvar + prologue];
|
2022-11-04 14:22:20 +01:00
|
|
|
|
ordidx++;
|
|
|
|
|
}
|
2019-12-16 19:42:59 +01:00
|
|
|
|
|
2022-11-04 14:22:20 +01:00
|
|
|
|
// Then the feedback variables, reordered by dynamic status
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto max_lag_lead : {pair {0, 0}, pair {1, 0}, pair {1, 1}, pair {0, 1}})
|
2020-04-29 18:46:15 +02:00
|
|
|
|
for (int vtx : feed_back_vertices)
|
|
|
|
|
if (int simvar = v_index1[vertex(vtx, subG)];
|
2023-11-30 15:28:57 +01:00
|
|
|
|
variable_lag_lead[old_endo_idx_block2orig[simvar + prologue]] == max_lag_lead)
|
2020-04-24 11:36:05 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
eq_idx_block2orig[ordidx] = old_eq_idx_block2orig[simvar + prologue];
|
|
|
|
|
endo_idx_block2orig[ordidx] = old_endo_idx_block2orig[simvar + prologue];
|
2020-04-24 11:36:05 +02:00
|
|
|
|
ordidx++;
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
2020-03-19 17:46:10 +01:00
|
|
|
|
|
2020-04-15 17:56:28 +02:00
|
|
|
|
updateReverseVariableEquationOrderings();
|
2020-04-30 11:15:55 +02:00
|
|
|
|
|
|
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
2020-04-30 12:48:16 +02:00
|
|
|
|
{
|
|
|
|
|
computeDynamicStructureOfBlock(blk);
|
|
|
|
|
computeSimulationTypeOfBlock(blk);
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2009-12-16 18:13:23 +01:00
|
|
|
|
void
|
2020-04-15 17:56:28 +02:00
|
|
|
|
ModelTree::printBlockDecomposition() const
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-21 23:06:51 +02:00
|
|
|
|
int largest_block = 0, Nb_SimulBlocks = 0, Nb_feedback_variable = 0;
|
2020-04-23 16:04:02 +02:00
|
|
|
|
int Nb_TotalBlocks = blocks.size();
|
2020-03-24 18:26:06 +01:00
|
|
|
|
for (int block = 0; block < Nb_TotalBlocks; block++)
|
2020-04-23 16:04:02 +02:00
|
|
|
|
if (BlockSimulationType simulation_type = blocks[block].simulation_type;
|
2020-04-21 23:06:51 +02:00
|
|
|
|
simulation_type == BlockSimulationType::solveForwardComplete
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete)
|
|
|
|
|
{
|
|
|
|
|
Nb_SimulBlocks++;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (int size = blocks[block].size; size > largest_block)
|
2020-04-21 23:06:51 +02:00
|
|
|
|
{
|
|
|
|
|
largest_block = size;
|
2020-04-23 16:04:02 +02:00
|
|
|
|
Nb_feedback_variable = blocks[block].mfs_size;
|
2020-04-21 23:06:51 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
|
|
int Nb_RecursBlocks = Nb_TotalBlocks - Nb_SimulBlocks;
|
|
|
|
|
cout << Nb_TotalBlocks << " block(s) found:" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " " << Nb_RecursBlocks << " recursive block(s) and " << Nb_SimulBlocks
|
|
|
|
|
<< " simultaneous block(s)." << endl
|
2009-12-16 14:21:31 +01:00
|
|
|
|
<< " the largest simultaneous block has " << largest_block << " equation(s)" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " and " << Nb_feedback_variable
|
|
|
|
|
<< " feedback variable(s)." << endl;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-03-20 15:23:23 +01:00
|
|
|
|
void
|
2020-04-30 12:48:16 +02:00
|
|
|
|
ModelTree::reduceBlockDecomposition()
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-30 12:48:16 +02:00
|
|
|
|
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;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int var = 0; var < blocks[blk - 1].size; var++)
|
2020-04-30 12:48:16 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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});
|
2020-04-30 12:48:16 +02:00
|
|
|
|
}
|
2020-04-21 18:10:46 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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))
|
2020-04-30 12:48:16 +02:00
|
|
|
|
{
|
|
|
|
|
// Merge the current block into the previous one
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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)
|
2020-04-30 12:48:16 +02:00
|
|
|
|
if (b >= blk)
|
|
|
|
|
b--;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& b : eq2block)
|
2020-04-30 12:48:16 +02:00
|
|
|
|
if (b >= blk)
|
|
|
|
|
b--;
|
|
|
|
|
blk--;
|
|
|
|
|
continue;
|
|
|
|
|
}
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2020-03-20 15:23:23 +01:00
|
|
|
|
void
|
|
|
|
|
ModelTree::determineLinearBlocks()
|
2009-12-16 14:21:31 +01:00
|
|
|
|
{
|
2020-04-21 18:10:46 +02:00
|
|
|
|
// Note that field “linear” in class BlockInfo defaults to true
|
2020-04-28 12:35:18 +02:00
|
|
|
|
for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
|
switch (blocks[blk].simulation_type)
|
|
|
|
|
{
|
2022-03-15 16:47:03 +01:00
|
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
2020-04-28 12:35:18 +02:00
|
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
2022-03-15 16:47:03 +01:00
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
2020-04-28 12:35:18 +02:00
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : blocks_derivatives[blk])
|
2019-12-16 19:42:59 +01:00
|
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
|
int lag = get<2>(indices);
|
2019-12-16 19:42:59 +01:00
|
|
|
|
if (lag == 0)
|
|
|
|
|
{
|
|
|
|
|
set<pair<int, int>> endogenous;
|
|
|
|
|
d1->collectEndogenous(endogenous);
|
2020-04-28 12:35:18 +02:00
|
|
|
|
for (int l = 0; l < blocks[blk].size; l++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (endogenous.contains({endo_idx_block2orig[blocks[blk].first_equation + l], 0}))
|
2020-04-28 12:35:18 +02:00
|
|
|
|
{
|
|
|
|
|
blocks[blk].linear = false;
|
|
|
|
|
goto the_end;
|
|
|
|
|
}
|
2019-12-16 19:42:59 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
2020-04-28 12:35:18 +02:00
|
|
|
|
the_end:
|
|
|
|
|
break;
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesSimple:
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : blocks_derivatives[blk])
|
2019-12-16 19:42:59 +01:00
|
|
|
|
{
|
2020-04-24 12:29:02 +02:00
|
|
|
|
int lag = get<2>(indices);
|
2019-12-16 19:42:59 +01:00
|
|
|
|
set<pair<int, int>> endogenous;
|
|
|
|
|
d1->collectEndogenous(endogenous);
|
2020-04-28 12:35:18 +02:00
|
|
|
|
for (int l = 0; l < blocks[blk].size; l++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (endogenous.contains({endo_idx_block2orig[blocks[blk].first_equation + l], lag}))
|
2020-04-28 12:35:18 +02:00
|
|
|
|
{
|
|
|
|
|
blocks[blk].linear = false;
|
|
|
|
|
goto the_end2;
|
|
|
|
|
}
|
2019-12-16 19:42:59 +01:00
|
|
|
|
}
|
2020-04-28 12:35:18 +02:00
|
|
|
|
the_end2:
|
|
|
|
|
break;
|
|
|
|
|
default:
|
|
|
|
|
break;
|
|
|
|
|
}
|
2009-12-16 14:21:31 +01:00
|
|
|
|
}
|
|
|
|
|
|
2008-02-03 11:28:36 +01:00
|
|
|
|
int
|
|
|
|
|
ModelTree::equation_number() const
|
2009-01-23 11:59:37 +01:00
|
|
|
|
{
|
2009-12-16 18:13:23 +01:00
|
|
|
|
return (equations.size());
|
2009-01-23 11:59:37 +01:00
|
|
|
|
}
|
2008-02-03 11:28:36 +01:00
|
|
|
|
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::computeDerivatives(int order, const set<int>& vars)
|
2008-02-03 11:28:36 +01:00
|
|
|
|
{
|
2019-12-20 16:59:30 +01:00
|
|
|
|
assert(order >= 1);
|
2008-02-03 11:28:36 +01:00
|
|
|
|
|
2020-01-20 17:22:32 +01:00
|
|
|
|
computed_derivs_order = order;
|
|
|
|
|
|
2018-11-22 14:32:40 +01:00
|
|
|
|
// Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
derivatives.resize(max(static_cast<size_t>(order + 1), derivatives.size()));
|
|
|
|
|
NNZDerivatives.resize(max(static_cast<size_t>(order + 1), NNZDerivatives.size()), 0);
|
2009-04-20 12:48:54 +02:00
|
|
|
|
|
2018-11-22 14:32:40 +01:00
|
|
|
|
// First-order derivatives
|
|
|
|
|
for (int var : vars)
|
2019-04-12 15:41:52 +02:00
|
|
|
|
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
2018-11-22 14:32:40 +01:00
|
|
|
|
{
|
|
|
|
|
expr_t d1 = equations[eq]->getDerivative(var);
|
|
|
|
|
if (d1 == Zero)
|
|
|
|
|
continue;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
derivatives[1][{eq, var}] = d1;
|
2018-11-22 14:32:40 +01:00
|
|
|
|
++NNZDerivatives[1];
|
|
|
|
|
}
|
2009-04-20 12:48:54 +02:00
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
// Compute the sparse representation of the Jacobian
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : derivatives[1])
|
|
|
|
|
jacobian_sparse_column_major_order.try_emplace({indices[0], getJacobianCol(indices[1], true)},
|
|
|
|
|
d1);
|
|
|
|
|
jacobian_sparse_colptr
|
|
|
|
|
= computeCSCColPtr(jacobian_sparse_column_major_order, getJacobianColsNbr(true));
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
2018-11-22 14:32:40 +01:00
|
|
|
|
// Higher-order derivatives
|
|
|
|
|
for (int o = 2; o <= order; o++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [lower_indices, lower_d] : derivatives[o - 1])
|
2018-11-22 14:32:40 +01:00
|
|
|
|
for (int var : vars)
|
2008-02-03 11:28:36 +01:00
|
|
|
|
{
|
2022-01-21 11:31:35 +01:00
|
|
|
|
if (lower_indices.back() > var)
|
2009-04-20 12:48:54 +02:00
|
|
|
|
continue;
|
|
|
|
|
|
2022-01-21 11:31:35 +01:00
|
|
|
|
expr_t d = lower_d->getDerivative(var);
|
2018-11-22 14:32:40 +01:00
|
|
|
|
if (d == Zero)
|
2009-04-20 12:48:54 +02:00
|
|
|
|
continue;
|
2018-11-22 14:32:40 +01:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<int> indices {lower_indices};
|
2018-11-22 14:32:40 +01:00
|
|
|
|
indices.push_back(var);
|
|
|
|
|
// At this point, indices of endogenous variables are sorted in non-decreasing order
|
|
|
|
|
derivatives[o][indices] = d;
|
2019-06-17 15:28:33 +02:00
|
|
|
|
// We output symmetric elements at order = 2
|
|
|
|
|
if (o == 2 && indices[1] != indices[2])
|
|
|
|
|
NNZDerivatives[o] += 2;
|
2019-04-12 15:41:52 +02:00
|
|
|
|
else
|
2018-11-22 14:32:40 +01:00
|
|
|
|
NNZDerivatives[o]++;
|
2008-02-03 11:28:36 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void
|
2018-09-25 15:56:52 +02:00
|
|
|
|
ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
|
2008-02-03 11:28:36 +01:00
|
|
|
|
{
|
2022-09-26 14:53:36 +02:00
|
|
|
|
/* Ensure that we don’t have any model-local variable in the model at this
|
|
|
|
|
point (we used to treat them as temporary terms) */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
assert([&] {
|
2022-09-26 14:53:36 +02:00
|
|
|
|
set<int> used_local_vars;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& equation : equations)
|
2022-09-26 14:53:36 +02:00
|
|
|
|
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
|
|
|
|
return used_local_vars.empty();
|
|
|
|
|
}());
|
2018-03-27 17:14:30 +02:00
|
|
|
|
|
2018-11-30 12:22:13 +01:00
|
|
|
|
// Compute the temporary terms in equations and derivatives
|
2023-04-05 14:53:47 +02:00
|
|
|
|
map<pair<int, int>, unordered_set<expr_t>> temp_terms_map;
|
|
|
|
|
unordered_map<expr_t, pair<int, pair<int, int>>> reference_count;
|
2019-08-14 15:28:41 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& equation : equations)
|
|
|
|
|
equation->computeTemporaryTerms({0, 0}, temp_terms_map, reference_count, is_matlab);
|
2019-08-14 15:28:41 +02:00
|
|
|
|
|
|
|
|
|
for (int order = 1; order < static_cast<int>(derivatives.size()); order++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : derivatives[order])
|
|
|
|
|
it.second->computeTemporaryTerms({order, 0}, temp_terms_map, reference_count, is_matlab);
|
2019-08-14 15:28:41 +02:00
|
|
|
|
|
|
|
|
|
/* 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)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& it : temp_terms_map)
|
|
|
|
|
erase_if(it.second, [](expr_t e) { return !dynamic_cast<AbstractExternalFunctionNode*>(e); });
|
2015-09-03 13:50:02 +02:00
|
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
|
// Fill the structures
|
2018-11-30 12:22:13 +01:00
|
|
|
|
temporary_terms_derivatives.clear();
|
|
|
|
|
temporary_terms_derivatives.resize(derivatives.size());
|
2019-04-23 11:07:32 +02:00
|
|
|
|
for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
copy(temp_terms_map[{order, 0}].begin(), temp_terms_map[{order, 0}].end(),
|
|
|
|
|
inserter(temporary_terms_derivatives.at(order),
|
|
|
|
|
temporary_terms_derivatives.at(order).begin()));
|
2018-03-27 17:14:30 +02:00
|
|
|
|
|
2018-11-30 12:22:13 +01:00
|
|
|
|
// Compute indices in MATLAB/Julia vector
|
2022-09-26 14:53:36 +02:00
|
|
|
|
for (int order {0}, idx {0}; order < static_cast<int>(derivatives.size()); order++)
|
2020-06-05 14:52:22 +02:00
|
|
|
|
for (auto it : temporary_terms_derivatives[order])
|
2018-11-30 12:22:13 +01:00
|
|
|
|
temporary_terms_idxs[it] = idx++;
|
2018-03-27 17:14:30 +02:00
|
|
|
|
}
|
|
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
|
void
|
2022-11-04 12:00:11 +01:00
|
|
|
|
ModelTree::computeBlockTemporaryTerms(bool no_tmp_terms)
|
2020-05-06 17:13:47 +02:00
|
|
|
|
{
|
|
|
|
|
int nb_blocks = blocks.size();
|
|
|
|
|
|
2023-04-05 14:53:47 +02:00
|
|
|
|
unordered_map<expr_t, tuple<int, int, int>> reference_count;
|
|
|
|
|
vector<vector<unordered_set<expr_t>>> temp_terms(nb_blocks);
|
2020-05-06 17:13:47 +02:00
|
|
|
|
for (int blk = 0; blk < nb_blocks; blk++)
|
|
|
|
|
{
|
2023-04-05 14:53:47 +02:00
|
|
|
|
temp_terms[blk].resize(blocks[blk].size + 1);
|
2020-05-06 17:13:47 +02:00
|
|
|
|
for (int eq = 0; eq < blocks[blk].size; eq++)
|
|
|
|
|
{
|
2022-10-26 17:33:26 +02:00
|
|
|
|
/* It is important to compute temporary terms of the renormalized
|
|
|
|
|
equation if the latter is going to be used in the output files.
|
|
|
|
|
Otherwise, for an equation of the form log(x) = RHS, a temporary
|
|
|
|
|
term could be associated to log(x), and since it would be
|
|
|
|
|
associated to this equation, it would be printed and thus computed
|
|
|
|
|
*before* x is actually evaluated, and thus would be incorrect. */
|
|
|
|
|
if ((blocks[blk].simulation_type == BlockSimulationType::evaluateBackward
|
|
|
|
|
|| blocks[blk].simulation_type == BlockSimulationType::evaluateForward
|
|
|
|
|
|| eq < blocks[blk].getRecursiveSize())
|
|
|
|
|
&& isBlockEquationRenormalized(blk, eq))
|
2023-11-30 15:28:57 +01:00
|
|
|
|
getBlockEquationRenormalizedExpr(blk, eq)->computeBlockTemporaryTerms(
|
|
|
|
|
blk, eq, temp_terms, reference_count);
|
2020-05-06 17:13:47 +02:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
getBlockEquationExpr(blk, eq)->computeBlockTemporaryTerms(blk, eq, temp_terms,
|
|
|
|
|
reference_count);
|
2020-05-06 17:13:47 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [ignore, d] : blocks_derivatives[blk])
|
2023-04-05 14:53:47 +02:00
|
|
|
|
d->computeBlockTemporaryTerms(blk, blocks[blk].size, temp_terms, reference_count);
|
2020-05-06 17:13:47 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-11-04 12:00:11 +01:00
|
|
|
|
/* 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)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& it : temp_terms)
|
|
|
|
|
for (auto& it2 : it)
|
|
|
|
|
erase_if(it2, [](expr_t e) { return !dynamic_cast<AbstractExternalFunctionNode*>(e); });
|
2022-11-04 12:00:11 +01:00
|
|
|
|
|
2023-04-05 14:53:47 +02:00
|
|
|
|
blocks_temporary_terms.resize(nb_blocks);
|
|
|
|
|
for (int blk {0}; blk < nb_blocks; blk++)
|
|
|
|
|
{
|
|
|
|
|
blocks_temporary_terms.at(blk).resize(temp_terms.at(blk).size());
|
|
|
|
|
for (size_t i {0}; i < temp_terms.at(blk).size(); i++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
copy(temp_terms.at(blk).at(i).begin(), temp_terms.at(blk).at(i).end(),
|
|
|
|
|
inserter(blocks_temporary_terms.at(blk).at(i),
|
|
|
|
|
blocks_temporary_terms.at(blk).at(i).begin()));
|
2023-04-05 14:53:47 +02:00
|
|
|
|
}
|
|
|
|
|
|
2020-05-06 17:13:47 +02:00
|
|
|
|
// Compute indices in the temporary terms vector
|
|
|
|
|
blocks_temporary_terms_idxs.clear();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int idx {0}; auto& blk_tt : blocks_temporary_terms)
|
|
|
|
|
for (auto& eq_tt : blk_tt)
|
2020-05-13 16:58:19 +02:00
|
|
|
|
for (auto tt : eq_tt)
|
|
|
|
|
blocks_temporary_terms_idxs[tt] = idx++;
|
2020-05-06 17:13:47 +02:00
|
|
|
|
}
|
|
|
|
|
|
2017-02-20 12:18:11 +01:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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
|
2017-02-20 12:18:11 +01:00
|
|
|
|
{
|
|
|
|
|
// Local var used to keep track of temp nodes already written
|
2019-04-18 14:34:48 +02:00
|
|
|
|
temporary_terms_t tt2 = temp_term_union;
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"("external_functions_temporary_terms_)" << concat << R"(": [)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_term {false}; auto it : tt)
|
2019-04-18 14:34:48 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (dynamic_cast<AbstractExternalFunctionNode*>(it))
|
2019-04-18 14:34:48 +02:00
|
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
|
if (exchange(printed_term, true))
|
2019-04-18 14:34:48 +02:00
|
|
|
|
output << ", ";
|
|
|
|
|
vector<string> efout;
|
|
|
|
|
it->writeJsonExternalFunctionOutput(efout, tt2, tef_terms);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_efout {false}; auto& it : efout)
|
2019-04-18 14:34:48 +02:00
|
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
|
if (exchange(printed_efout, true))
|
2019-04-18 14:34:48 +02:00
|
|
|
|
output << ", ";
|
2022-06-03 16:24:26 +02:00
|
|
|
|
output << it;
|
2019-04-18 14:34:48 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
tt2.insert(it);
|
|
|
|
|
}
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
|
|
|
|
output << "]"
|
2019-04-03 16:32:52 +02:00
|
|
|
|
<< R"(, "temporary_terms_)" << concat << R"(": [)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_term {false}; const auto& it : tt)
|
2019-04-18 14:34:48 +02:00
|
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
|
if (exchange(printed_term, true))
|
2019-04-18 14:34:48 +02:00
|
|
|
|
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);
|
|
|
|
|
}
|
2017-02-20 12:18:11 +01:00
|
|
|
|
output << "]";
|
|
|
|
|
}
|
|
|
|
|
|
2016-12-28 14:02:50 +01:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::fixNestedParenthesis(ostringstream& output, map<string, string>& tmp_paren_vars,
|
|
|
|
|
bool& message_printed) const
|
2016-12-28 14:02:50 +01:00
|
|
|
|
{
|
|
|
|
|
string str = output.str();
|
2016-12-30 18:32:20 +01:00
|
|
|
|
if (!testNestedParenthesis(str))
|
|
|
|
|
return;
|
2016-12-28 14:02:50 +01:00
|
|
|
|
int open = 0;
|
2016-12-30 11:26:56 +01:00
|
|
|
|
int first_open_paren = 0;
|
|
|
|
|
int matching_paren = 0;
|
|
|
|
|
bool hit_limit = false;
|
|
|
|
|
int i1 = 0;
|
2017-01-02 11:37:15 +01:00
|
|
|
|
for (size_t i = 0; i < str.length(); i++)
|
2016-12-28 14:02:50 +01:00
|
|
|
|
{
|
2016-12-30 11:26:56 +01:00
|
|
|
|
if (str.at(i) == '(')
|
|
|
|
|
{
|
|
|
|
|
if (open == 0)
|
|
|
|
|
first_open_paren = i;
|
|
|
|
|
open++;
|
|
|
|
|
}
|
|
|
|
|
else if (str.at(i) == ')')
|
|
|
|
|
{
|
|
|
|
|
open--;
|
|
|
|
|
if (open == 0)
|
|
|
|
|
matching_paren = i;
|
|
|
|
|
}
|
2016-12-28 14:02:50 +01:00
|
|
|
|
if (open > 32)
|
2016-12-30 11:26:56 +01:00
|
|
|
|
hit_limit = true;
|
|
|
|
|
|
|
|
|
|
if (hit_limit && open == 0)
|
2016-12-28 14:02:50 +01:00
|
|
|
|
{
|
2017-01-05 15:19:13 +01:00
|
|
|
|
if (!message_printed)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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;
|
2017-01-05 15:19:13 +01:00
|
|
|
|
message_printed = true;
|
|
|
|
|
}
|
2016-12-30 11:26:56 +01:00
|
|
|
|
string str1 = str.substr(first_open_paren, matching_paren - first_open_paren + 1);
|
2019-12-16 19:42:59 +01:00
|
|
|
|
string repstr, varname;
|
2016-12-30 11:26:56 +01:00
|
|
|
|
while (testNestedParenthesis(str1))
|
|
|
|
|
{
|
2019-12-20 16:59:30 +01:00
|
|
|
|
size_t open_paren_idx = string::npos;
|
2017-01-02 11:37:15 +01:00
|
|
|
|
size_t match_paren_idx = string::npos;
|
|
|
|
|
size_t last_open_paren = string::npos;
|
|
|
|
|
for (size_t j = 0; j < str1.length(); j++)
|
2016-12-30 11:26:56 +01:00
|
|
|
|
{
|
|
|
|
|
if (str1.at(j) == '(')
|
|
|
|
|
{
|
|
|
|
|
// don't match, e.g. y(1)
|
2019-12-16 19:42:59 +01:00
|
|
|
|
if (size_t idx = str1.find_last_of("*/-+", j - 1);
|
|
|
|
|
j == 0 || (idx != string::npos && idx == j - 1))
|
2016-12-30 11:26:56 +01:00
|
|
|
|
open_paren_idx = j;
|
|
|
|
|
last_open_paren = j;
|
|
|
|
|
}
|
|
|
|
|
else if (str1.at(j) == ')')
|
|
|
|
|
{
|
|
|
|
|
// don't match, e.g. y(1)
|
2019-12-16 19:42:59 +01:00
|
|
|
|
if (size_t idx = str1.find_last_not_of("0123456789", j - 1);
|
|
|
|
|
idx != string::npos && idx != last_open_paren)
|
2016-12-30 11:26:56 +01:00
|
|
|
|
match_paren_idx = j;
|
|
|
|
|
}
|
|
|
|
|
|
2017-01-02 11:37:15 +01:00
|
|
|
|
if (open_paren_idx != string::npos && match_paren_idx != string::npos)
|
2016-12-30 11:26:56 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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())
|
2016-12-30 11:26:56 +01:00
|
|
|
|
{
|
2022-06-02 10:50:21 +02:00
|
|
|
|
varname = "paren32_tmp_var_" + to_string(i1++);
|
2023-12-13 16:29:46 +01:00
|
|
|
|
repstr += varname + " = " + val + ";\n";
|
2016-12-30 11:26:56 +01:00
|
|
|
|
tmp_paren_vars[val] = varname;
|
|
|
|
|
}
|
|
|
|
|
else
|
|
|
|
|
varname = it->second;
|
|
|
|
|
str1.replace(open_paren_idx, match_paren_idx - open_paren_idx + 1, varname);
|
|
|
|
|
break;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (auto it = tmp_paren_vars.find(str1); it == tmp_paren_vars.end())
|
2016-12-30 11:26:56 +01:00
|
|
|
|
{
|
2022-06-02 10:50:21 +02:00
|
|
|
|
varname = "paren32_tmp_var_" + to_string(i1++);
|
2023-12-13 16:29:46 +01:00
|
|
|
|
repstr += varname + " = " + str1 + ";\n";
|
2016-12-30 11:26:56 +01:00
|
|
|
|
}
|
|
|
|
|
else
|
|
|
|
|
varname = it->second;
|
|
|
|
|
str.replace(first_open_paren, matching_paren - first_open_paren + 1, varname);
|
2023-12-13 10:50:24 +01:00
|
|
|
|
size_t insertLoc = str.find_last_of('\n', first_open_paren);
|
2016-12-30 11:26:56 +01:00
|
|
|
|
str.insert(insertLoc + 1, repstr);
|
|
|
|
|
hit_limit = false;
|
|
|
|
|
i = -1;
|
|
|
|
|
first_open_paren = matching_paren = open = 0;
|
2016-12-28 14:02:50 +01:00
|
|
|
|
}
|
|
|
|
|
}
|
2016-12-30 11:26:56 +01:00
|
|
|
|
output.str(str);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
bool
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::testNestedParenthesis(const string& str) const
|
2016-12-30 11:26:56 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int open {0}; char i : str)
|
2016-12-30 11:26:56 +01:00
|
|
|
|
{
|
2018-06-04 12:26:16 +02:00
|
|
|
|
if (i == '(')
|
2016-12-30 11:26:56 +01:00
|
|
|
|
open++;
|
2018-06-04 12:26:16 +02:00
|
|
|
|
else if (i == ')')
|
2016-12-30 11:26:56 +01:00
|
|
|
|
open--;
|
|
|
|
|
if (open > 32)
|
|
|
|
|
return true;
|
|
|
|
|
}
|
|
|
|
|
return false;
|
2016-12-28 14:02:50 +01:00
|
|
|
|
}
|
|
|
|
|
|
2017-02-20 12:18:11 +01:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeJsonModelLocalVariables(ostream& output, bool write_tef_terms,
|
|
|
|
|
deriv_node_temp_terms_t& tef_terms) const
|
2017-02-20 12:18:11 +01:00
|
|
|
|
{
|
|
|
|
|
/* 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;
|
|
|
|
|
|
2018-06-04 12:26:16 +02:00
|
|
|
|
for (auto equation : equations)
|
2018-07-17 18:34:07 +02:00
|
|
|
|
equation->collectVariables(SymbolType::modelLocalVariable, used_local_vars);
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"("model_local_variables": [)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something {false}; int id : local_variables_vector)
|
2022-05-04 16:01:34 +02:00
|
|
|
|
if (used_local_vars.contains(id))
|
2017-08-28 15:14:11 +02:00
|
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
|
if (exchange(printed_something, true))
|
2017-08-28 15:14:11 +02:00
|
|
|
|
output << ", ";
|
|
|
|
|
|
2022-09-27 12:51:22 +02:00
|
|
|
|
expr_t value = local_variables_table.at(id);
|
2020-06-05 16:07:52 +02:00
|
|
|
|
if (write_tef_terms)
|
2017-11-06 15:21:11 +01:00
|
|
|
|
{
|
2020-06-05 16:07:52 +02:00
|
|
|
|
vector<string> efout;
|
2020-06-05 14:52:22 +02:00
|
|
|
|
value->writeJsonExternalFunctionOutput(efout, {}, tef_terms);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_efout {false}; auto& it : efout)
|
2020-06-05 16:07:52 +02:00
|
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
|
if (exchange(printed_efout, true))
|
2020-06-05 16:07:52 +02:00
|
|
|
|
output << ", ";
|
2022-06-03 16:24:26 +02:00
|
|
|
|
output << it;
|
2020-06-05 16:07:52 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (!efout.empty())
|
2017-11-06 15:21:11 +01:00
|
|
|
|
output << ", ";
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << R"({"variable": ")" << symbol_table.getName(id) << R"(", "value": ")";
|
2020-06-05 14:52:22 +02:00
|
|
|
|
value->writeJsonOutput(output, {}, tef_terms);
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"("})" << endl;
|
2017-08-28 15:14:11 +02:00
|
|
|
|
}
|
2017-02-20 12:18:11 +01:00
|
|
|
|
output << "]";
|
|
|
|
|
}
|
|
|
|
|
|
2022-07-13 13:04:10 +02:00
|
|
|
|
int
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeBytecodeBinFile(const filesystem::path& filename, bool is_two_boundaries) const
|
2010-01-22 11:03:29 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ofstream SaveCode {filename, ios::out | ios::binary};
|
2010-01-22 11:03:29 +01:00
|
|
|
|
if (!SaveCode.is_open())
|
|
|
|
|
{
|
2023-01-05 16:40:04 +01:00
|
|
|
|
cerr << R"(Error : Can't open file ")" << filename.string() << R"(" for writing)" << endl;
|
2010-01-22 11:03:29 +01:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
2022-07-13 13:04:10 +02:00
|
|
|
|
int u_count {0};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : derivatives[1])
|
|
|
|
|
if (int deriv_id {indices[1]}; getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
|
|
|
|
int eq {indices[0]};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
SaveCode.write(reinterpret_cast<char*>(&eq), sizeof eq);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
int tsid {getTypeSpecificIDByDerivID(deriv_id)};
|
|
|
|
|
int lag {getLagByDerivID(deriv_id)};
|
|
|
|
|
int varr {tsid + lag * symbol_table.endo_nbr()};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
SaveCode.write(reinterpret_cast<char*>(&varr), sizeof varr);
|
|
|
|
|
SaveCode.write(reinterpret_cast<char*>(&lag), sizeof lag);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
int u {u_count + symbol_table.endo_nbr()};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
SaveCode.write(reinterpret_cast<char*>(&u), sizeof u);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
u_count++;
|
|
|
|
|
}
|
2010-01-22 11:03:29 +01:00
|
|
|
|
if (is_two_boundaries)
|
2022-07-13 13:04:10 +02:00
|
|
|
|
u_count += symbol_table.endo_nbr();
|
|
|
|
|
for (int j {0}; j < symbol_table.endo_nbr(); j++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
SaveCode.write(reinterpret_cast<char*>(&j), sizeof j);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
for (int j {0}; j < symbol_table.endo_nbr(); j++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
SaveCode.write(reinterpret_cast<char*>(&j), sizeof j);
|
2010-01-22 11:03:29 +01:00
|
|
|
|
SaveCode.close();
|
2022-07-13 13:04:10 +02:00
|
|
|
|
return u_count;
|
2010-01-22 11:03:29 +01:00
|
|
|
|
}
|
|
|
|
|
|
2022-07-19 18:24:36 +02:00
|
|
|
|
int
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeBlockBytecodeBinFile(ofstream& bin_file, int block) const
|
2022-07-19 18:24:36 +02:00
|
|
|
|
{
|
|
|
|
|
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};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, ignore] : blocks_derivatives[block])
|
2022-07-19 18:24:36 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const auto& [eq, var, lag] {indices};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
if (lag != 0 && !is_two_boundaries)
|
|
|
|
|
continue;
|
|
|
|
|
if (eq >= block_recursive && var >= block_recursive)
|
|
|
|
|
{
|
|
|
|
|
int v {eq - block_recursive};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bin_file.write(reinterpret_cast<char*>(&v), sizeof v);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
int varr {var - block_recursive + lag * block_mfs};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bin_file.write(reinterpret_cast<char*>(&varr), sizeof varr);
|
|
|
|
|
bin_file.write(reinterpret_cast<const char*>(&lag), sizeof lag);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
int u {u_count + block_mfs};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bin_file.write(reinterpret_cast<char*>(&u), sizeof u);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
u_count++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (is_two_boundaries)
|
|
|
|
|
u_count += block_mfs;
|
|
|
|
|
for (int j {block_recursive}; j < block_size; j++)
|
|
|
|
|
{
|
|
|
|
|
int varr {getBlockVariableID(block, j)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bin_file.write(reinterpret_cast<char*>(&varr), sizeof varr);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
for (int j {block_recursive}; j < block_size; j++)
|
|
|
|
|
{
|
|
|
|
|
int eqr {getBlockEquationID(block, j)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bin_file.write(reinterpret_cast<char*>(&eqr), sizeof eqr);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
return u_count;
|
|
|
|
|
}
|
|
|
|
|
|
2009-04-30 15:14:33 +02:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_basename,
|
|
|
|
|
ExprNodeOutputType output_type, bool write_equation_tags) const
|
2009-04-30 15:14:33 +02:00
|
|
|
|
{
|
2019-04-04 17:01:37 +02:00
|
|
|
|
filesystem::create_directories(mod_basename + "/latex");
|
2019-07-11 17:33:53 +02:00
|
|
|
|
|
2023-01-05 16:40:04 +01:00
|
|
|
|
const filesystem::path filename {mod_basename + "/latex/" + latex_basename + ".tex"},
|
2023-11-30 15:28:57 +01:00
|
|
|
|
content_filename {mod_basename + "/latex/" + latex_basename + "_content" + ".tex"};
|
|
|
|
|
ofstream output {filename, ios::out | ios::binary};
|
2009-04-30 15:14:33 +02:00
|
|
|
|
if (!output.is_open())
|
|
|
|
|
{
|
2023-01-05 16:40:04 +01:00
|
|
|
|
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
2009-04-30 15:14:33 +02:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ofstream content_output {content_filename, ios::out | ios::binary};
|
2015-07-15 08:58:15 +02:00
|
|
|
|
if (!content_output.is_open())
|
|
|
|
|
{
|
2023-01-05 16:40:04 +01:00
|
|
|
|
cerr << "ERROR: Can't open file " << content_filename.string() << " for writing" << endl;
|
2015-07-15 08:58:15 +02:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
|
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;
|
2009-04-30 15:14:33 +02:00
|
|
|
|
|
|
|
|
|
// Write model local variables
|
2018-06-04 12:26:16 +02:00
|
|
|
|
for (int id : local_variables_vector)
|
2009-04-30 15:14:33 +02:00
|
|
|
|
{
|
2022-09-27 12:51:22 +02:00
|
|
|
|
expr_t value = local_variables_table.at(id);
|
2009-04-30 15:14:33 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
content_output << R"(\begin{dmath*})" << endl << symbol_table.getTeXName(id) << " = ";
|
2009-04-30 15:14:33 +02:00
|
|
|
|
// Use an empty set for the temporary terms
|
2015-07-15 08:58:15 +02:00
|
|
|
|
value->writeOutput(content_output, output_type);
|
2019-04-03 16:32:52 +02:00
|
|
|
|
content_output << endl << R"(\end{dmath*})" << endl;
|
2009-04-30 15:14:33 +02:00
|
|
|
|
}
|
|
|
|
|
|
2019-04-23 11:07:32 +02:00
|
|
|
|
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
2009-04-30 15:14:33 +02:00
|
|
|
|
{
|
2017-04-05 10:28:37 +02:00
|
|
|
|
content_output << "% Equation " << eq + 1 << endl;
|
2017-04-04 15:28:27 +02:00
|
|
|
|
if (write_equation_tags)
|
2020-02-20 15:29:10 +01:00
|
|
|
|
equation_tags.writeLatexOutput(content_output, eq);
|
2017-04-04 15:28:27 +02:00
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
|
content_output << R"(\begin{dmath})" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// 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);
|
2019-04-03 16:32:52 +02:00
|
|
|
|
content_output << endl << R"(\end{dmath})" << endl;
|
2009-04-30 15:14:33 +02:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << R"(\include{)" << latex_basename + "_content"
|
|
|
|
|
<< "}" << endl
|
2019-04-03 16:32:52 +02:00
|
|
|
|
<< R"(\end{document})" << endl;
|
2009-04-30 15:14:33 +02:00
|
|
|
|
|
|
|
|
|
output.close();
|
2015-07-15 08:58:15 +02:00
|
|
|
|
content_output.close();
|
2009-04-30 15:14:33 +02:00
|
|
|
|
}
|
|
|
|
|
|
2008-02-03 11:28:36 +01:00
|
|
|
|
void
|
2023-12-13 15:37:07 +01:00
|
|
|
|
ModelTree::addEquation(expr_t eq, const optional<int>& lineno)
|
2009-01-23 11:59:37 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto beq = dynamic_cast<BinaryOpNode*>(eq);
|
2019-12-16 19:42:59 +01:00
|
|
|
|
assert(beq && beq->op_code == BinaryOpcode::equal);
|
2008-02-03 11:28:36 +01:00
|
|
|
|
|
|
|
|
|
equations.push_back(beq);
|
2023-12-13 15:37:07 +01:00
|
|
|
|
equations_lineno.push_back(lineno);
|
2008-02-03 11:28:36 +01:00
|
|
|
|
}
|
2009-09-02 16:37:59 +02:00
|
|
|
|
|
2019-01-28 14:57:30 +01:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const
|
2019-01-28 14:57:30 +01:00
|
|
|
|
{
|
2020-01-27 15:58:13 +01:00
|
|
|
|
for (size_t i = 0; i < equations.size(); i++)
|
2023-01-04 16:18:53 +01:00
|
|
|
|
if (!equation_tags.exists(i, "mcp"))
|
2020-01-27 15:58:13 +01:00
|
|
|
|
equations[i]->findConstantEquations(subst_table);
|
2019-01-28 14:57:30 +01:00
|
|
|
|
}
|
|
|
|
|
|
2009-09-02 16:37:59 +02:00
|
|
|
|
void
|
2023-12-13 15:37:07 +01:00
|
|
|
|
ModelTree::addEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags)
|
2009-09-02 16:37:59 +02:00
|
|
|
|
{
|
2023-11-14 15:55:30 +01:00
|
|
|
|
equation_tags.add(equations.size(), move(eq_tags));
|
2023-12-13 15:37:07 +01:00
|
|
|
|
addEquation(eq, lineno);
|
2009-09-02 16:37:59 +02:00
|
|
|
|
}
|
2009-09-30 17:10:31 +02:00
|
|
|
|
|
|
|
|
|
void
|
2010-09-16 19:18:45 +02:00
|
|
|
|
ModelTree::addAuxEquation(expr_t eq)
|
2009-09-30 17:10:31 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto beq = dynamic_cast<BinaryOpNode*>(eq);
|
2019-12-16 19:42:59 +01:00
|
|
|
|
assert(beq && beq->op_code == BinaryOpcode::equal);
|
2009-09-30 17:10:31 +02:00
|
|
|
|
|
|
|
|
|
aux_equations.push_back(beq);
|
|
|
|
|
}
|
2010-10-15 19:05:16 +02:00
|
|
|
|
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::addTrendVariables(const vector<int>& trend_vars, expr_t growth_factor) noexcept(false)
|
2010-10-15 19:05:16 +02:00
|
|
|
|
{
|
2019-07-05 18:36:10 +02:00
|
|
|
|
for (int id : trend_vars)
|
2022-05-04 16:01:34 +02:00
|
|
|
|
if (trend_symbols_map.contains(id))
|
2023-11-30 15:28:57 +01:00
|
|
|
|
throw TrendException {symbol_table.getName(id)};
|
2010-10-15 19:05:16 +02:00
|
|
|
|
else
|
2019-07-05 18:36:10 +02:00
|
|
|
|
trend_symbols_map[id] = growth_factor;
|
2010-10-15 19:05:16 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::addNonstationaryVariables(const vector<int>& nonstationary_vars, bool log_deflator,
|
|
|
|
|
expr_t deflator) noexcept(false)
|
2010-10-15 19:05:16 +02:00
|
|
|
|
{
|
2019-07-05 18:36:10 +02:00
|
|
|
|
for (int id : nonstationary_vars)
|
2022-05-04 16:01:34 +02:00
|
|
|
|
if (nonstationary_symbols_map.contains(id))
|
2023-11-30 15:28:57 +01:00
|
|
|
|
throw TrendException {symbol_table.getName(id)};
|
2010-10-15 19:05:16 +02:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
nonstationary_symbols_map[id] = {log_deflator, deflator};
|
2010-10-15 19:05:16 +02:00
|
|
|
|
}
|
2011-06-22 11:56:07 +02:00
|
|
|
|
|
|
|
|
|
void
|
|
|
|
|
ModelTree::initializeVariablesAndEquations()
|
|
|
|
|
{
|
2017-06-14 23:49:10 +02:00
|
|
|
|
for (size_t j = 0; j < equations.size(); j++)
|
2020-04-17 14:55:55 +02:00
|
|
|
|
eq_idx_block2orig.push_back(j);
|
2018-10-10 17:06:19 +02:00
|
|
|
|
|
|
|
|
|
for (int j = 0; j < symbol_table.endo_nbr(); j++)
|
2020-04-17 14:55:55 +02:00
|
|
|
|
endo_idx_block2orig.push_back(j);
|
2011-06-22 11:56:07 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void
|
|
|
|
|
ModelTree::set_cutoff_to_zero()
|
|
|
|
|
{
|
|
|
|
|
cutoff = 0;
|
|
|
|
|
}
|
2011-08-19 15:05:57 +02:00
|
|
|
|
|
2012-11-29 18:07:48 +01:00
|
|
|
|
void
|
2016-05-18 12:26:19 +02:00
|
|
|
|
ModelTree::computeParamsDerivatives(int paramsDerivsOrder)
|
2012-11-29 18:07:48 +01:00
|
|
|
|
{
|
2018-11-22 15:41:11 +01:00
|
|
|
|
assert(paramsDerivsOrder >= 1);
|
|
|
|
|
|
2012-11-29 18:07:48 +01:00
|
|
|
|
set<int> deriv_id_set;
|
|
|
|
|
addAllParamDerivId(deriv_id_set);
|
2016-03-25 15:38:49 +01:00
|
|
|
|
|
2018-11-22 15:41:11 +01:00
|
|
|
|
// First-order derivatives w.r.t. params
|
2018-06-04 12:26:16 +02:00
|
|
|
|
for (int param : deriv_id_set)
|
2012-11-29 18:07:48 +01:00
|
|
|
|
{
|
2019-04-23 11:07:32 +02:00
|
|
|
|
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
2012-11-29 18:07:48 +01:00
|
|
|
|
{
|
2018-11-22 15:41:11 +01:00
|
|
|
|
expr_t d = equations[eq]->getDerivative(param);
|
|
|
|
|
if (d == Zero)
|
2012-11-29 18:07:48 +01:00
|
|
|
|
continue;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
params_derivatives[{0, 1}][{eq, param}] = d;
|
2012-11-29 18:07:48 +01:00
|
|
|
|
}
|
|
|
|
|
|
2019-04-23 11:07:32 +02:00
|
|
|
|
for (int endoOrd = 1; endoOrd < static_cast<int>(derivatives.size()); endoOrd++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [lower_indices, lower_d] : derivatives[endoOrd])
|
2016-05-12 12:02:34 +02:00
|
|
|
|
{
|
2022-01-21 11:31:35 +01:00
|
|
|
|
expr_t d = lower_d->getDerivative(param);
|
2018-11-22 15:41:11 +01:00
|
|
|
|
if (d == Zero)
|
2016-05-12 12:02:34 +02:00
|
|
|
|
continue;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<int> indices {lower_indices};
|
2022-01-21 11:31:35 +01:00
|
|
|
|
indices.push_back(param);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
params_derivatives[{endoOrd, 1}][indices] = d;
|
2016-05-12 12:02:34 +02:00
|
|
|
|
}
|
2018-11-22 15:41:11 +01:00
|
|
|
|
}
|
2012-11-29 18:07:48 +01:00
|
|
|
|
|
2018-11-22 15:41:11 +01:00
|
|
|
|
// Higher-order derivatives w.r.t. parameters
|
2019-04-23 11:07:32 +02:00
|
|
|
|
for (int endoOrd = 0; endoOrd < static_cast<int>(derivatives.size()); endoOrd++)
|
2018-11-22 15:41:11 +01:00
|
|
|
|
for (int paramOrd = 2; paramOrd <= paramsDerivsOrder; paramOrd++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [lower_indices, lower_d] : params_derivatives[{endoOrd, paramOrd - 1}])
|
2018-11-22 15:41:11 +01:00
|
|
|
|
for (int param : deriv_id_set)
|
|
|
|
|
{
|
2022-01-21 11:31:35 +01:00
|
|
|
|
if (lower_indices.back() > param)
|
2018-11-22 15:41:11 +01:00
|
|
|
|
continue;
|
2016-05-18 10:33:45 +02:00
|
|
|
|
|
2022-01-21 11:31:35 +01:00
|
|
|
|
expr_t d = lower_d->getDerivative(param);
|
2018-11-22 15:41:11 +01:00
|
|
|
|
if (d == Zero)
|
|
|
|
|
continue;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<int> indices {lower_indices};
|
2022-01-21 11:31:35 +01:00
|
|
|
|
indices.push_back(param);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// At this point, indices of both endogenous and parameters are sorted in non-decreasing
|
|
|
|
|
// order
|
|
|
|
|
params_derivatives[{endoOrd, paramOrd}][indices] = d;
|
2018-11-22 15:41:11 +01:00
|
|
|
|
}
|
2012-11-29 18:07:48 +01:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void
|
|
|
|
|
ModelTree::computeParamsDerivativesTemporaryTerms()
|
|
|
|
|
{
|
2023-04-05 14:53:47 +02:00
|
|
|
|
unordered_map<expr_t, pair<int, pair<int, int>>> reference_count;
|
2018-11-16 18:24:06 +01:00
|
|
|
|
|
|
|
|
|
/* The temp terms should be constructed in the same order as the for loops in
|
|
|
|
|
{Static,Dynamic}Model::write{Json,}ParamsDerivativesFile() */
|
2023-04-05 14:53:47 +02:00
|
|
|
|
map<pair<int, int>, unordered_set<expr_t>> temp_terms_map;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [order, derivs] : params_derivatives)
|
|
|
|
|
for (const auto& [indices, d] : derivs)
|
2023-04-05 14:53:47 +02:00
|
|
|
|
d->computeTemporaryTerms(order, temp_terms_map, reference_count, true);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [order, tts] : temp_terms_map)
|
2023-04-05 14:53:47 +02:00
|
|
|
|
copy(temp_terms_map[order].begin(), temp_terms_map[order].end(),
|
2023-11-30 15:28:57 +01:00
|
|
|
|
inserter(params_derivs_temporary_terms[order],
|
|
|
|
|
params_derivs_temporary_terms[order].begin()));
|
2018-05-28 15:23:15 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int idx {0}; const auto& [order, tts] : params_derivs_temporary_terms)
|
|
|
|
|
for (const auto& tt : tts)
|
2018-11-30 12:22:13 +01:00
|
|
|
|
params_derivs_temporary_terms_idxs[tt] = idx++;
|
2012-11-29 18:07:48 +01:00
|
|
|
|
}
|
2013-10-29 11:47:59 +01:00
|
|
|
|
|
2017-06-14 07:01:31 +02:00
|
|
|
|
bool
|
|
|
|
|
ModelTree::isNonstationary(int symb_id) const
|
2013-10-29 11:47:59 +01:00
|
|
|
|
{
|
2022-05-04 16:01:34 +02:00
|
|
|
|
return nonstationary_symbols_map.contains(symb_id);
|
2013-10-29 11:47:59 +01:00
|
|
|
|
}
|
|
|
|
|
|
2017-02-02 15:09:43 +01:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeJsonModelEquations(ostream& output, bool residuals) const
|
2017-02-02 15:09:43 +01:00
|
|
|
|
{
|
2017-02-20 12:18:11 +01:00
|
|
|
|
if (residuals)
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << endl << R"("residuals":[)" << endl;
|
2017-02-20 12:18:11 +01:00
|
|
|
|
else
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << endl << R"("model":[)" << endl;
|
2019-04-23 11:07:32 +02:00
|
|
|
|
for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
|
2017-02-02 15:09:43 +01:00
|
|
|
|
{
|
2017-02-20 12:18:11 +01:00
|
|
|
|
if (eq > 0)
|
|
|
|
|
output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
BinaryOpNode* eq_node = equations[eq];
|
2018-11-28 14:32:26 +01:00
|
|
|
|
expr_t lhs = eq_node->arg1;
|
|
|
|
|
expr_t rhs = eq_node->arg2;
|
2017-03-15 12:52:55 +01:00
|
|
|
|
|
2017-02-20 12:18:11 +01:00
|
|
|
|
if (residuals)
|
2017-02-02 15:09:43 +01:00
|
|
|
|
{
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"({"residual": {)"
|
|
|
|
|
<< R"("lhs": ")";
|
2020-05-06 17:13:47 +02:00
|
|
|
|
lhs->writeJsonOutput(output, {}, {});
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"(")";
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"(, "rhs": ")";
|
2020-05-06 17:13:47 +02:00
|
|
|
|
rhs->writeJsonOutput(output, {}, {});
|
2022-07-05 15:54:52 +02:00
|
|
|
|
output << R"("})";
|
2017-02-02 15:09:43 +01:00
|
|
|
|
}
|
2017-02-20 12:18:11 +01:00
|
|
|
|
else
|
|
|
|
|
{
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"({"lhs": ")";
|
2019-12-16 19:42:59 +01:00
|
|
|
|
lhs->writeJsonOutput(output, {}, {});
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"(", "rhs": ")";
|
2019-12-16 19:42:59 +01:00
|
|
|
|
rhs->writeJsonOutput(output, {}, {});
|
2022-05-05 18:39:27 +02:00
|
|
|
|
output << R"(")";
|
|
|
|
|
if (equations_lineno[eq])
|
|
|
|
|
output << R"(, "line": )" << *equations_lineno[eq];
|
2017-02-20 12:18:11 +01:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (auto eqtags = equation_tags.getTagsByEqn(eq); !eqtags.empty())
|
2017-02-20 12:18:11 +01:00
|
|
|
|
{
|
2019-04-03 16:32:52 +02:00
|
|
|
|
output << R"(, "tags": {)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something {false}; const auto& [name, value] : eqtags)
|
2017-02-20 12:18:11 +01:00
|
|
|
|
{
|
2022-06-03 16:24:26 +02:00
|
|
|
|
if (exchange(printed_something, true))
|
2017-02-20 12:18:11 +01:00
|
|
|
|
output << ", ";
|
2019-12-16 19:42:59 +01:00
|
|
|
|
output << R"(")" << name << R"(": ")" << value << R"(")";
|
2017-02-20 12:18:11 +01:00
|
|
|
|
}
|
|
|
|
|
output << "}";
|
|
|
|
|
eqtags.clear();
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
output << "}" << endl;
|
2017-02-02 15:09:43 +01:00
|
|
|
|
}
|
|
|
|
|
output << endl << "]" << endl;
|
|
|
|
|
}
|
2018-10-26 11:44:26 +02:00
|
|
|
|
|
|
|
|
|
string
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::matlab_arch(const string& mexext)
|
2018-10-26 11:44:26 +02:00
|
|
|
|
{
|
|
|
|
|
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")
|
2019-11-04 15:50:26 +01:00
|
|
|
|
{
|
|
|
|
|
cerr << "32-bit MATLAB not supported on macOS" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
2018-10-26 11:44:26 +02:00
|
|
|
|
else if (mexext == "mexmaci64")
|
|
|
|
|
return "maci64";
|
2022-05-23 10:22:57 +02:00
|
|
|
|
else if (mexext == "mexmaca64")
|
|
|
|
|
return "maca64";
|
2018-10-26 11:44:26 +02:00
|
|
|
|
else
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
cerr << "ERROR: 'mexext' option to preprocessor incorrectly set, needed with 'use_dll'"
|
|
|
|
|
<< endl;
|
2018-10-26 11:44:26 +02:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2021-09-13 15:14:10 +02:00
|
|
|
|
#ifdef __APPLE__
|
2023-07-11 22:12:03 +02:00
|
|
|
|
|
|
|
|
|
pair<filesystem::path, bool>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::findCompilerOnMacos(const string& mexext)
|
2021-09-13 15:14:10 +02:00
|
|
|
|
{
|
2023-07-11 22:12:03 +02:00
|
|
|
|
/* Try to find gcc, otherwise use Apple’s clang compiler.
|
|
|
|
|
Homebrew binaries are located in /usr/local/bin/ on x86_64 systems and in
|
|
|
|
|
/opt/homebrew/bin/ on arm64 systems.
|
|
|
|
|
Apple’s clang is located both in /usr/bin/gcc and /usr/bin/clang, it
|
|
|
|
|
automatically selects x86_64 or arm64 depending on the compile-time
|
|
|
|
|
environment. */
|
2023-06-07 10:27:20 +02:00
|
|
|
|
const string macos_gcc_version {"13"};
|
2021-09-13 15:14:10 +02:00
|
|
|
|
|
2023-07-11 22:12:03 +02:00
|
|
|
|
if (filesystem::path global_gcc_path {"/usr/local/bin/gcc-" + macos_gcc_version};
|
|
|
|
|
exists(global_gcc_path) && mexext == "mexmaci64")
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {global_gcc_path, false};
|
2023-01-10 11:48:17 +01:00
|
|
|
|
else if (filesystem::path global_gcc_path {"/opt/homebrew/bin/gcc-" + macos_gcc_version};
|
|
|
|
|
exists(global_gcc_path) && mexext == "mexmaca64")
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {global_gcc_path, false};
|
2023-07-11 22:12:03 +02:00
|
|
|
|
else if (filesystem::path global_clang_path {"/usr/bin/clang"}; exists(global_clang_path))
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {global_clang_path, true};
|
2021-09-13 15:14:10 +02:00
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
cerr << "ERROR: You must install gcc-" << macos_gcc_version
|
|
|
|
|
<< " on your system before using the `use_dll` option of Dynare. "
|
2023-07-11 22:12:03 +02:00
|
|
|
|
<< "You should install Homebrew";
|
|
|
|
|
if (mexext == "mexmaca64")
|
|
|
|
|
cerr << " for arm64";
|
|
|
|
|
else if (mexext == "mexmaci64")
|
|
|
|
|
cerr << " for x86_64";
|
|
|
|
|
cerr << " and run `brew install gcc-" << macos_gcc_version << "` in a terminal." << endl;
|
2021-09-13 15:14:10 +02:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
#endif
|
|
|
|
|
|
2022-10-05 18:34:21 +02:00
|
|
|
|
filesystem::path
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::compileMEX(const filesystem::path& output_dir, const string& output_basename,
|
|
|
|
|
const string& mexext, const vector<filesystem::path>& input_files,
|
|
|
|
|
const filesystem::path& matlabroot, bool link) const
|
2018-10-26 11:44:26 +02:00
|
|
|
|
{
|
2022-10-14 14:52:11 +02:00
|
|
|
|
assert(!mex_compilation_workers.empty());
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string gcc_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"};
|
|
|
|
|
const string clang_opt_flags {
|
|
|
|
|
"-O3 -g0 --param ira-max-conflict-table-size=1 -Wno-unused-command-line-argument"};
|
2018-10-26 11:44:26 +02:00
|
|
|
|
|
2019-04-04 17:01:37 +02:00
|
|
|
|
filesystem::path compiler;
|
2018-10-26 11:44:26 +02:00
|
|
|
|
ostringstream flags;
|
|
|
|
|
string libs;
|
2023-07-11 22:12:03 +02:00
|
|
|
|
bool is_clang {false};
|
2018-10-26 11:44:26 +02:00
|
|
|
|
|
2020-09-18 15:00:47 +02:00
|
|
|
|
if (matlabroot.empty())
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
cerr << "ERROR: 'matlabroot' option to preprocessor is not set, needed with 'use_dll'"
|
|
|
|
|
<< endl;
|
2020-09-18 15:00:47 +02:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
|
2018-10-26 11:44:26 +02:00
|
|
|
|
if (mexext == "mex")
|
|
|
|
|
{
|
|
|
|
|
// Octave
|
|
|
|
|
compiler = matlabroot / "bin" / "mkoctfile";
|
|
|
|
|
flags << "--mex";
|
|
|
|
|
}
|
|
|
|
|
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";
|
2019-07-11 16:11:24 +02:00
|
|
|
|
libs = "-lmex -lmx";
|
2021-05-24 18:21:23 +02:00
|
|
|
|
if (mexext == "mexa64")
|
2018-10-26 11:44:26 +02:00
|
|
|
|
{
|
|
|
|
|
// GNU/Linux
|
|
|
|
|
flags << " -D_GNU_SOURCE -fPIC -pthread"
|
|
|
|
|
<< " -shared -Wl,--no-undefined -Wl,-rpath-link," << bin_dir;
|
2021-05-24 18:21:23 +02:00
|
|
|
|
libs += " -lm";
|
2018-10-26 11:44:26 +02:00
|
|
|
|
}
|
2023-03-20 17:54:32 +01:00
|
|
|
|
else if (mexext == "mexw64") // Windows
|
|
|
|
|
flags << " -static-libgcc -shared";
|
2021-05-24 18:21:23 +02:00
|
|
|
|
#ifdef __APPLE__
|
2022-05-23 10:22:57 +02:00
|
|
|
|
else if (mexext == "mexmaci64" || mexext == "mexmaca64")
|
2018-10-26 11:44:26 +02:00
|
|
|
|
{
|
2023-07-11 22:12:03 +02:00
|
|
|
|
tie(compiler, is_clang) = findCompilerOnMacos(mexext);
|
2021-01-12 17:34:44 +01:00
|
|
|
|
flags << " -fno-common -Wl,-twolevel_namespace -undefined error -bundle";
|
2021-05-24 18:21:23 +02:00
|
|
|
|
libs += " -lm";
|
|
|
|
|
}
|
|
|
|
|
#endif
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
cerr << "ERROR: unsupported value '" << mexext << "' for 'mexext' option" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
2018-10-26 11:44:26 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2022-10-05 18:34:21 +02:00
|
|
|
|
filesystem::path output_filename {output_dir / (output_basename + "." + (link ? mexext : "o"))};
|
2018-10-26 11:44:26 +02:00
|
|
|
|
|
|
|
|
|
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
|
2019-12-02 19:21:14 +01:00
|
|
|
|
|
|
|
|
|
if (user_set_compiler.empty())
|
|
|
|
|
cmd << compiler << " ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
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);
|
|
|
|
|
}
|
2019-12-02 19:21:14 +01:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
cmd << user_set_compiler << " ";
|
2019-12-02 19:21:14 +01:00
|
|
|
|
|
|
|
|
|
if (user_set_subst_flags.empty())
|
2023-07-11 22:12:03 +02:00
|
|
|
|
cmd << (is_clang ? clang_opt_flags : gcc_opt_flags) << " " << flags.str() << " ";
|
2019-12-02 19:21:14 +01:00
|
|
|
|
else
|
|
|
|
|
cmd << user_set_subst_flags << " ";
|
|
|
|
|
|
|
|
|
|
if (!user_set_add_flags.empty())
|
|
|
|
|
cmd << user_set_add_flags << " ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& f : input_files)
|
2022-10-05 18:34:21 +02:00
|
|
|
|
cmd << f << " ";
|
|
|
|
|
cmd << "-o " << output_filename << " ";
|
2019-12-02 19:21:14 +01:00
|
|
|
|
|
2022-10-05 18:34:21 +02:00
|
|
|
|
if (link)
|
|
|
|
|
{
|
|
|
|
|
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;
|
|
|
|
|
}
|
2019-12-02 19:21:14 +01:00
|
|
|
|
else
|
2022-10-05 18:34:21 +02:00
|
|
|
|
cmd << " -c";
|
2018-10-26 11:44:26 +02:00
|
|
|
|
|
|
|
|
|
#ifdef _WIN32
|
|
|
|
|
cmd << '"';
|
|
|
|
|
#endif
|
|
|
|
|
|
2022-10-11 11:00:50 +02:00
|
|
|
|
cout << "Compiling " << output_filename.string() << endl;
|
2022-10-05 18:34:21 +02:00
|
|
|
|
|
|
|
|
|
// The prerequisites are the object files among the input files
|
|
|
|
|
set<filesystem::path> prerequisites;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
copy_if(input_files.begin(), input_files.end(), inserter(prerequisites, prerequisites.end()),
|
|
|
|
|
[](const auto& p) { return p.extension() == ".o"; });
|
2018-10-26 11:44:26 +02:00
|
|
|
|
|
2022-10-14 14:52:11 +02:00
|
|
|
|
unique_lock<mutex> lk {mex_compilation_mut};
|
|
|
|
|
mex_compilation_queue.emplace_back(output_filename, prerequisites, cmd.str());
|
|
|
|
|
lk.unlock();
|
|
|
|
|
mex_compilation_cv.notify_one();
|
2022-10-05 18:34:21 +02:00
|
|
|
|
|
|
|
|
|
return output_filename;
|
2018-10-26 11:44:26 +02:00
|
|
|
|
}
|
2019-12-03 14:19:32 +01:00
|
|
|
|
|
|
|
|
|
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++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto varexpr = dynamic_cast<VariableNode*>(aux_equations[i]->arg1);
|
2019-12-03 14:19:32 +01:00
|
|
|
|
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)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (auto it = auxEndoToEq.find(endo); it != auxEndoToEq.end() && it->second != i)
|
2019-12-16 19:42:59 +01:00
|
|
|
|
add_edge(i, it->second, g);
|
2019-12-03 14:19:32 +01:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// 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]]];
|
|
|
|
|
}
|
2020-03-30 14:51:53 +02:00
|
|
|
|
|
|
|
|
|
map<tuple<int, int, int>, expr_t>
|
|
|
|
|
ModelTree::collectFirstOrderDerivativesEndogenous()
|
|
|
|
|
{
|
|
|
|
|
map<tuple<int, int, int>, expr_t> endo_derivatives;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& [indices, d1] : derivatives[1])
|
2020-03-30 14:51:53 +02:00
|
|
|
|
if (getTypeByDerivID(indices[1]) == SymbolType::endogenous)
|
|
|
|
|
{
|
|
|
|
|
int eq = indices[0];
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var {getTypeSpecificIDByDerivID(indices[1])};
|
2020-03-30 14:51:53 +02:00
|
|
|
|
int lag = getLagByDerivID(indices[1]);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
endo_derivatives[{eq, var, lag}] = d1;
|
2020-03-30 14:51:53 +02:00
|
|
|
|
}
|
|
|
|
|
return endo_derivatives;
|
|
|
|
|
}
|
2020-04-15 17:56:28 +02:00
|
|
|
|
|
|
|
|
|
ModelTree::jacob_map_t
|
2022-11-28 13:40:49 +01:00
|
|
|
|
ModelTree::computeSymbolicJacobian(bool contemporaneous_only) const
|
2020-04-15 17:56:28 +02:00
|
|
|
|
{
|
|
|
|
|
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);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [endo, lag] : endos_and_lags)
|
2022-11-28 13:40:49 +01:00
|
|
|
|
if (!contemporaneous_only || lag == 0)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
symbolic_jacobian.try_emplace({i, endo}, 1);
|
2020-04-15 17:56:28 +02:00
|
|
|
|
}
|
|
|
|
|
return symbolic_jacobian;
|
|
|
|
|
}
|
|
|
|
|
|
2023-04-24 17:49:54 +02:00
|
|
|
|
pair<bool, ModelTree::jacob_map_t>
|
|
|
|
|
ModelTree::computeLeftHandSideSymbolicJacobian() const
|
|
|
|
|
{
|
|
|
|
|
jacob_map_t lhs_symbolic_jacobian;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto not_contemporaneous = [](const pair<int, int>& p) { return p.second != 0; };
|
2023-04-24 17:49:54 +02:00
|
|
|
|
|
|
|
|
|
for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (equations_lineno[eq]) // Hand-written equation: test whether LHS has single contemporaneous
|
|
|
|
|
// endo
|
2023-04-24 17:49:54 +02:00
|
|
|
|
{
|
|
|
|
|
set<pair<int, int>> endos_and_lags;
|
|
|
|
|
equations[eq]->arg1->collectEndogenous(endos_and_lags);
|
|
|
|
|
erase_if(endos_and_lags, not_contemporaneous);
|
|
|
|
|
if (endos_and_lags.size() == 1)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
lhs_symbolic_jacobian.try_emplace({eq, endos_and_lags.begin()->first}, 1);
|
2023-04-24 17:49:54 +02:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {false, {}};
|
2023-04-24 17:49:54 +02:00
|
|
|
|
}
|
|
|
|
|
else // Generated equation: keep all endos on both LHS and RHS
|
|
|
|
|
{
|
|
|
|
|
set<pair<int, int>> endos_and_lags;
|
|
|
|
|
equations[eq]->collectEndogenous(endos_and_lags);
|
|
|
|
|
erase_if(endos_and_lags, not_contemporaneous);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [endo, lag] : endos_and_lags)
|
|
|
|
|
lhs_symbolic_jacobian.try_emplace({eq, endo}, 1);
|
2023-04-24 17:49:54 +02:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {true, lhs_symbolic_jacobian};
|
2023-04-24 17:49:54 +02:00
|
|
|
|
}
|
|
|
|
|
|
2020-04-15 17:56:28 +02:00
|
|
|
|
void
|
|
|
|
|
ModelTree::updateReverseVariableEquationOrderings()
|
|
|
|
|
{
|
|
|
|
|
int n = equations.size();
|
2020-04-17 14:55:55 +02:00
|
|
|
|
eq_idx_orig2block.resize(n);
|
|
|
|
|
endo_idx_orig2block.resize(n);
|
2020-04-15 17:56:28 +02:00
|
|
|
|
for (int i = 0; i < n; i++)
|
|
|
|
|
{
|
2020-04-17 14:55:55 +02:00
|
|
|
|
endo_idx_orig2block[endo_idx_block2orig[i]] = i;
|
|
|
|
|
eq_idx_orig2block[eq_idx_block2orig[i]] = i;
|
2020-04-15 17:56:28 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
2021-10-26 18:06:26 +02:00
|
|
|
|
|
|
|
|
|
expr_t
|
|
|
|
|
ModelTree::getRHSFromLHS(expr_t lhs) const
|
|
|
|
|
{
|
|
|
|
|
for (auto eq : equations)
|
|
|
|
|
if (eq->arg1 == lhs)
|
|
|
|
|
return eq->arg2;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
throw ExprNode::MatchFailureException {"Cannot find an equation with the requested LHS"};
|
2021-10-26 18:06:26 +02:00
|
|
|
|
}
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
2022-09-30 15:39:41 +02:00
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::initializeMEXCompilationWorkers(int numworkers, const filesystem::path& dynareroot,
|
|
|
|
|
const string& mexext)
|
2022-10-14 14:52:11 +02:00
|
|
|
|
{
|
|
|
|
|
assert(numworkers > 0);
|
|
|
|
|
assert(mex_compilation_workers.empty());
|
|
|
|
|
|
|
|
|
|
cout << "Spawning " << numworkers << " threads for compiling MEX files." << endl;
|
|
|
|
|
|
|
|
|
|
for (int i {0}; i < numworkers; i++)
|
2023-12-13 10:33:37 +01:00
|
|
|
|
/* Passing the stop_token by const reference is ok (and makes clang-tidy happier),
|
|
|
|
|
since the std::jthread constructor calls the lambda with the return argument of the
|
|
|
|
|
get_stop_token() method, which returns a stop_token by value; hence there is no lifetime
|
|
|
|
|
issue. See:
|
|
|
|
|
https://stackoverflow.com/questions/72990607/const-stdstop-token-or-just-stdstop-token-as-parameter-for-thread-funct
|
|
|
|
|
*/
|
2023-12-13 10:28:54 +01:00
|
|
|
|
mex_compilation_workers.emplace_back([](const stop_token& stoken) {
|
2022-10-14 14:52:11 +02:00
|
|
|
|
unique_lock<mutex> lk {mex_compilation_mut};
|
2022-10-17 13:57:34 +02:00
|
|
|
|
filesystem::path output;
|
|
|
|
|
string cmd;
|
2022-10-14 14:52:11 +02:00
|
|
|
|
|
2022-10-17 13:57:34 +02:00
|
|
|
|
/* Look for an object to compile, whose prerequisites are already
|
|
|
|
|
compiled. If found, remove it from the queue, save the output path and
|
|
|
|
|
the compilation command, and return true. Must be run under the lock. */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto pick_job = [&cmd, &output] {
|
2022-10-17 13:57:34 +02:00
|
|
|
|
for (auto it {mex_compilation_queue.begin()}; it != mex_compilation_queue.end(); ++it)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (const auto& prerequisites {get<1>(*it)}; // Will become dangling after erase
|
2022-10-17 13:57:34 +02:00
|
|
|
|
includes(mex_compilation_done.begin(), mex_compilation_done.end(),
|
2022-10-14 14:52:11 +02:00
|
|
|
|
prerequisites.begin(), prerequisites.end()))
|
|
|
|
|
{
|
2022-10-17 13:57:34 +02:00
|
|
|
|
output = get<0>(*it);
|
|
|
|
|
cmd = get<2>(*it);
|
2022-10-14 14:52:11 +02:00
|
|
|
|
mex_compilation_queue.erase(it);
|
2022-12-20 14:47:32 +01:00
|
|
|
|
mex_compilation_ongoing.insert(output);
|
2022-10-17 13:57:34 +02:00
|
|
|
|
return true;
|
2022-10-14 14:52:11 +02:00
|
|
|
|
}
|
2022-10-17 13:57:34 +02:00
|
|
|
|
return false;
|
|
|
|
|
};
|
2022-10-14 14:52:11 +02:00
|
|
|
|
|
2022-10-17 13:57:34 +02:00
|
|
|
|
while (!stoken.stop_requested())
|
|
|
|
|
if (mex_compilation_cv.wait(lk, stoken, pick_job))
|
|
|
|
|
{
|
|
|
|
|
lk.unlock();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int r {system(cmd.c_str())};
|
2022-10-17 13:57:34 +02:00
|
|
|
|
lk.lock();
|
2022-12-20 14:47:32 +01:00
|
|
|
|
mex_compilation_ongoing.erase(output);
|
|
|
|
|
if (r)
|
|
|
|
|
mex_compilation_failed.insert(output);
|
|
|
|
|
else
|
|
|
|
|
mex_compilation_done.insert(output);
|
2022-10-17 13:57:34 +02:00
|
|
|
|
/* The object just compiled may be a prerequisite for several
|
|
|
|
|
other objects, so notify all waiting workers. Also needed to
|
|
|
|
|
notify the main thread when in
|
|
|
|
|
ModelTree::waitForMEXCompilationWorkers().*/
|
|
|
|
|
mex_compilation_cv.notify_all();
|
|
|
|
|
}
|
2022-10-14 14:52:11 +02:00
|
|
|
|
});
|
2023-03-20 17:54:32 +01:00
|
|
|
|
|
|
|
|
|
/* Set some environment variables needed for compilation on Windows/MATLAB
|
|
|
|
|
and macOS/Octave.
|
|
|
|
|
For Windows/MATLAB, this should be done only once, because otherwise
|
|
|
|
|
the PATH variable can become too long and GCC will not be found. */
|
|
|
|
|
if (mexext == "mexw64")
|
|
|
|
|
{
|
|
|
|
|
// 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 an internal copy of the string on MinGW,
|
|
|
|
|
contrary to what is done on GNU/Linux and macOS. */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (putenv(const_cast<char*>(newpath.c_str())) != 0)
|
2023-03-20 17:54:32 +01:00
|
|
|
|
{
|
|
|
|
|
cerr << "Can't set PATH" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
#ifdef __APPLE__
|
|
|
|
|
else if (mexext == "mex")
|
|
|
|
|
{
|
2023-07-11 22:12:03 +02:00
|
|
|
|
/* On macOS, with Octave, enforce our compiler. In particular this is
|
|
|
|
|
necessary if we’ve selected GCC; otherwise Clang will be used, and
|
|
|
|
|
it does not accept the same optimization flags (see dynare#1797) */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [compiler_path, is_clang] {findCompilerOnMacos(mexext)};
|
2023-07-11 22:12:03 +02:00
|
|
|
|
if (setenv("CC", compiler_path.c_str(), 1) != 0)
|
2023-03-20 17:54:32 +01:00
|
|
|
|
{
|
|
|
|
|
cerr << "Can't set CC environment variable" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
// We also define CXX, because that is used for linking
|
2023-07-11 22:12:03 +02:00
|
|
|
|
if (setenv("CXX", compiler_path.c_str(), 1) != 0)
|
2023-03-20 17:54:32 +01:00
|
|
|
|
{
|
|
|
|
|
cerr << "Can't set CXX environment variable" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
#endif
|
2022-10-14 14:52:11 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
void
|
2022-10-17 13:57:34 +02:00
|
|
|
|
ModelTree::waitForMEXCompilationWorkers()
|
2022-09-30 15:39:41 +02:00
|
|
|
|
{
|
2022-10-14 14:52:11 +02:00
|
|
|
|
unique_lock<mutex> lk {mex_compilation_mut};
|
2022-12-20 14:47:32 +01:00
|
|
|
|
mex_compilation_cv.wait(lk, [] {
|
|
|
|
|
return (mex_compilation_queue.empty() && mex_compilation_ongoing.empty())
|
2023-11-30 15:28:57 +01:00
|
|
|
|
|| !mex_compilation_failed.empty();
|
|
|
|
|
});
|
2022-12-20 14:47:32 +01:00
|
|
|
|
if (!mex_compilation_failed.empty())
|
|
|
|
|
{
|
|
|
|
|
cerr << "Compilation failed for: ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& p : mex_compilation_failed)
|
2022-12-20 14:47:32 +01:00
|
|
|
|
cerr << p.string() << " ";
|
|
|
|
|
cerr << endl;
|
|
|
|
|
lk.unlock(); // So that threads can process their stoken
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
2022-09-30 15:39:41 +02:00
|
|
|
|
}
|
2022-09-28 16:31:51 +02:00
|
|
|
|
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::computingPassBlock(const eval_context_t& eval_context, bool no_tmp_terms)
|
2022-09-28 16:31:51 +02:00
|
|
|
|
{
|
2023-04-24 15:05:22 +02:00
|
|
|
|
if (!computeNonSingularNormalization(eval_context))
|
2022-09-28 16:31:51 +02:00
|
|
|
|
return;
|
|
|
|
|
auto [prologue, epilogue] = computePrologueAndEpilogue();
|
|
|
|
|
auto first_order_endo_derivatives = collectFirstOrderDerivativesEndogenous();
|
2023-10-13 23:02:36 +02:00
|
|
|
|
equationTypeDetermination(first_order_endo_derivatives);
|
2022-09-28 16:31:51 +02:00
|
|
|
|
cout << "Finding the optimal block decomposition of the " << modelClassName() << "..." << endl;
|
|
|
|
|
computeBlockDecomposition(prologue, epilogue);
|
|
|
|
|
reduceBlockDecomposition();
|
|
|
|
|
printBlockDecomposition();
|
|
|
|
|
computeChainRuleJacobian();
|
|
|
|
|
determineLinearBlocks();
|
2022-11-04 12:00:11 +01:00
|
|
|
|
computeBlockTemporaryTerms(no_tmp_terms);
|
2022-09-28 16:31:51 +02:00
|
|
|
|
block_decomposed = true;
|
|
|
|
|
}
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
|
|
|
|
vector<int>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::computeCSCColPtr(const SparseColumnMajorOrderMatrix& matrix, int ncols)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<int> colptr(ncols + 1, matrix.size());
|
|
|
|
|
for (int k {0}, current_col {0}; const auto& [indices, d1] : matrix)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
while (indices.second >= current_col)
|
|
|
|
|
colptr[current_col++] = k;
|
|
|
|
|
k++;
|
|
|
|
|
}
|
|
|
|
|
return colptr;
|
|
|
|
|
}
|
2023-03-28 16:37:05 +02:00
|
|
|
|
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType output_type) const
|
2023-03-28 16:37:05 +02:00
|
|
|
|
{
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
for (auto aux_equation : aux_equations)
|
|
|
|
|
if (aux_equation->containsExternalFunction())
|
|
|
|
|
aux_equation->writeExternalFunctionOutput(output, output_type, {}, {}, tef_terms);
|
|
|
|
|
for (auto aux_equation : aux_equations)
|
|
|
|
|
{
|
|
|
|
|
aux_equation->writeOutput(output, output_type, {}, {}, tef_terms);
|
|
|
|
|
output << ";" << endl;
|
|
|
|
|
}
|
|
|
|
|
}
|