2008-02-03 11:28:36 +01:00
|
|
|
|
/*
|
2023-01-04 14:43:57 +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
|
|
|
|
*/
|
|
|
|
|
|
2023-12-01 15:39:01 +01:00
|
|
|
|
#ifndef MODEL_TREE_HH
|
|
|
|
|
#define MODEL_TREE_HH
|
2023-12-01 15:31:55 +01:00
|
|
|
|
|
|
|
|
|
#include <array>
|
|
|
|
|
#include <cassert>
|
|
|
|
|
#include <condition_variable>
|
|
|
|
|
#include <deque>
|
|
|
|
|
#include <filesystem>
|
|
|
|
|
#include <map>
|
|
|
|
|
#include <mutex>
|
|
|
|
|
#include <optional>
|
|
|
|
|
#include <ostream>
|
|
|
|
|
#include <string>
|
|
|
|
|
#include <thread>
|
|
|
|
|
#include <vector>
|
|
|
|
|
|
|
|
|
|
#include "Bytecode.hh"
|
|
|
|
|
#include "DataTree.hh"
|
|
|
|
|
#include "EquationTags.hh"
|
|
|
|
|
#include "ExtendedPreprocessorTypes.hh"
|
2008-02-03 11:28:36 +01:00
|
|
|
|
|
2022-01-06 14:54:57 +01:00
|
|
|
|
using namespace std;
|
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
// Helper to convert a vector into a tuple
|
2019-12-20 16:59:30 +01:00
|
|
|
|
template<typename T, size_t... Indices>
|
|
|
|
|
auto
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vectorToTupleHelper(const vector<T>& v, index_sequence<Indices...>)
|
2019-12-20 16:59:30 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return tuple {v[Indices]...};
|
2018-11-15 16:39:53 +01:00
|
|
|
|
}
|
2019-12-20 16:59:30 +01:00
|
|
|
|
template<size_t N, typename T>
|
|
|
|
|
auto
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vectorToTuple(const vector<T>& v)
|
2019-12-20 16:59:30 +01:00
|
|
|
|
{
|
2018-11-15 16:39:53 +01:00
|
|
|
|
assert(v.size() >= N);
|
|
|
|
|
return vectorToTupleHelper(v, make_index_sequence<N>());
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a
|
|
|
|
|
//! expr_t on the new normalized equation
|
|
|
|
|
using equation_type_and_normalized_equation_t = vector<pair<EquationType, BinaryOpNode*>>;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
|
|
|
|
//! Vector describing variables: max_lag in the block, max_lead in the block
|
2020-03-19 17:46:10 +01:00
|
|
|
|
using lag_lead_vector_t = vector<pair<int, int>>;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2009-04-14 16:39:53 +02:00
|
|
|
|
//! Shared code for static and dynamic models
|
2008-02-03 11:28:36 +01:00
|
|
|
|
class ModelTree : public DataTree
|
|
|
|
|
{
|
2009-12-16 14:21:31 +01:00
|
|
|
|
friend class DynamicModel;
|
|
|
|
|
friend class StaticModel;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
|
2019-12-02 19:21:14 +01:00
|
|
|
|
public:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
// Set via the `compiler` command
|
2023-11-30 15:28:57 +01:00
|
|
|
|
string user_set_add_flags, user_set_subst_flags, user_set_add_libs, user_set_subst_libs,
|
|
|
|
|
user_set_compiler;
|
|
|
|
|
|
2009-04-14 16:39:53 +02:00
|
|
|
|
protected:
|
2019-10-29 12:56:53 +01:00
|
|
|
|
/*
|
|
|
|
|
* ************** BEGIN **************
|
|
|
|
|
* The following structures keep track of the model equations and must all be updated
|
|
|
|
|
* when adding or removing an equation. Hence, if a new parallel structure is added
|
|
|
|
|
* in the future, it must be maintained whereever these structures are updated
|
2020-01-06 18:26:35 +01:00
|
|
|
|
* See in particular methods clearEquations(), replaceMyEquations() and
|
|
|
|
|
* computeRamseyPolicyFOCs() of DynamicModel class.
|
2019-10-29 12:56:53 +01:00
|
|
|
|
* NB: This message added with the introduction of the `exclude_eqs` option, hence
|
|
|
|
|
* that's a place to update future structures.
|
|
|
|
|
*/
|
2009-09-30 17:10:31 +02:00
|
|
|
|
//! Stores declared and generated auxiliary equations
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<BinaryOpNode*> equations;
|
2022-05-05 18:39:27 +02:00
|
|
|
|
/* Stores line numbers of declared equations; undefined in some cases (e.g.
|
|
|
|
|
auxiliary equations) */
|
|
|
|
|
vector<optional<int>> equations_lineno;
|
2019-10-29 12:56:53 +01:00
|
|
|
|
//! Stores equation tags
|
2020-02-20 15:29:10 +01:00
|
|
|
|
EquationTags equation_tags;
|
2019-10-29 12:56:53 +01:00
|
|
|
|
/*
|
|
|
|
|
* ************** END **************
|
|
|
|
|
*/
|
|
|
|
|
|
2009-09-30 17:10:31 +02:00
|
|
|
|
//! Only stores generated auxiliary equations, in an order meaningful for evaluation
|
2018-07-27 14:19:59 +02:00
|
|
|
|
/*! These equations only contain the definition of auxiliary variables, and
|
|
|
|
|
may diverge from those in the main model (equations), if other model
|
|
|
|
|
transformations applied subsequently. This is not a problem, since
|
|
|
|
|
aux_equations is only used for regenerating the values of auxiliaries
|
|
|
|
|
given the others.
|
|
|
|
|
|
|
|
|
|
For example, such a divergence appears when there is an expectation
|
|
|
|
|
operator in a ramsey model, see
|
|
|
|
|
tests/optimal_policy/nk_ramsey_expectation.mod */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<BinaryOpNode*> aux_equations;
|
2020-01-20 17:22:32 +01:00
|
|
|
|
|
|
|
|
|
//! Maximum order at which (endogenous) derivatives have been computed
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int computed_derivs_order {0};
|
2009-09-30 17:10:31 +02:00
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
//! Stores derivatives
|
|
|
|
|
/*! Index 0 is not used, index 1 contains first derivatives, ...
|
|
|
|
|
For each derivation order, stores a map whose key is a vector of integer: the
|
|
|
|
|
first integer is the equation index, the remaining ones are the derivation
|
2018-11-22 14:32:40 +01:00
|
|
|
|
IDs of variables (in non-decreasing order, to avoid storing symmetric
|
2022-06-13 16:24:33 +02:00
|
|
|
|
elements several times). Only non-zero derivatives are stored. */
|
2018-11-15 16:39:53 +01:00
|
|
|
|
vector<map<vector<int>, expr_t>> derivatives;
|
2012-11-29 18:07:48 +01:00
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
//! Number of non-zero derivatives
|
|
|
|
|
/*! Index 0 is not used, index 1 contains number of non-zero first
|
|
|
|
|
derivatives, ... */
|
|
|
|
|
vector<int> NNZDerivatives;
|
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
// Used to order pairs of indices (row, col) according to column-major order
|
|
|
|
|
struct columnMajorOrderLess
|
|
|
|
|
{
|
|
|
|
|
bool
|
2023-11-30 15:28:57 +01:00
|
|
|
|
operator()(const pair<int, int>& p1, const pair<int, int>& p2) const
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
return p1.second < p2.second || (p1.second == p2.second && p1.first < p2.first);
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
using SparseColumnMajorOrderMatrix = map<pair<int, int>, expr_t, columnMajorOrderLess>;
|
|
|
|
|
/* The nonzero values of the sparse Jacobian in column-major order (which is
|
|
|
|
|
the natural order for Compressed Sparse Column (CSC) storage).
|
|
|
|
|
The pair of indices is (row, column). */
|
|
|
|
|
SparseColumnMajorOrderMatrix jacobian_sparse_column_major_order;
|
|
|
|
|
/* Column indices for the sparse Jacobian in Compressed Sparse Column (CSC)
|
|
|
|
|
storage (corresponds to the “jc” vector in MATLAB terminology) */
|
|
|
|
|
vector<int> jacobian_sparse_colptr;
|
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
//! Derivatives with respect to parameters
|
|
|
|
|
/*! The key of the outer map is a pair (derivation order w.r.t. endogenous,
|
|
|
|
|
derivation order w.r.t. parameters). For e.g., { 1, 2 } corresponds to the jacobian
|
|
|
|
|
differentiated twice w.r.t. to parameters.
|
|
|
|
|
In inner maps, the vector of integers consists of: the equation index, then
|
2018-11-22 15:41:11 +01:00
|
|
|
|
the derivation IDs of endogenous (in non-decreasing order),
|
|
|
|
|
then the IDs of parameters (in non-decreasing order)*/
|
2023-11-30 15:28:57 +01:00
|
|
|
|
map<pair<int, int>, map<vector<int>, expr_t>> params_derivatives;
|
2012-11-29 18:07:48 +01:00
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
//! Temporary terms for residuals and derivatives
|
|
|
|
|
/*! Index 0 is temp. terms of residuals, index 1 for first derivatives, ... */
|
|
|
|
|
vector<temporary_terms_t> temporary_terms_derivatives;
|
2012-11-29 18:07:48 +01:00
|
|
|
|
|
2018-11-30 12:22:13 +01:00
|
|
|
|
//! Stores, for each temporary term, its index in the MATLAB/Julia vector
|
2018-03-27 17:14:30 +02:00
|
|
|
|
temporary_terms_idxs_t temporary_terms_idxs;
|
|
|
|
|
|
2018-11-15 16:39:53 +01:00
|
|
|
|
//! Temporary terms for parameter derivatives, under a disaggregated form
|
|
|
|
|
/*! The pair of integers is to be interpreted as in param_derivatives */
|
2019-12-20 16:59:30 +01:00
|
|
|
|
map<pair<int, int>, temporary_terms_t> params_derivs_temporary_terms;
|
2012-11-29 18:07:48 +01:00
|
|
|
|
|
2018-11-30 12:22:13 +01:00
|
|
|
|
//! Stores, for each temporary term in param. derivs, its index in the MATLAB/Julia vector
|
2018-05-28 15:23:15 +02:00
|
|
|
|
temporary_terms_idxs_t params_derivs_temporary_terms_idxs;
|
|
|
|
|
|
2010-10-15 19:05:16 +02:00
|
|
|
|
//! Trend variables and their growth factors
|
2013-03-26 16:46:18 +01:00
|
|
|
|
map<int, expr_t> trend_symbols_map;
|
|
|
|
|
|
|
|
|
|
//! for all trends; the boolean is true if this is a log-trend, false otherwise
|
2018-06-04 14:31:26 +02:00
|
|
|
|
using nonstationary_symbols_map_t = map<int, pair<bool, expr_t>>;
|
2013-03-26 16:46:18 +01:00
|
|
|
|
|
2010-10-15 19:05:16 +02:00
|
|
|
|
//! Nonstationary variables and their deflators
|
2013-03-26 16:46:18 +01:00
|
|
|
|
nonstationary_symbols_map_t nonstationary_symbols_map;
|
2008-02-03 11:28:36 +01:00
|
|
|
|
|
2020-04-17 14:55:55 +02:00
|
|
|
|
/* Maps indices of equations in the block-decomposition order into original
|
|
|
|
|
equation IDs */
|
|
|
|
|
vector<int> eq_idx_block2orig;
|
|
|
|
|
/* Maps indices of (endogenous) variables in the block-decomposition order into original
|
|
|
|
|
type-specific endogenous IDs */
|
|
|
|
|
vector<int> endo_idx_block2orig;
|
|
|
|
|
/* Maps original variable and equation indices into the block-decomposition order.
|
|
|
|
|
Set by updateReverseVariableEquationOrderings() */
|
|
|
|
|
vector<int> eq_idx_orig2block, endo_idx_orig2block;
|
2011-06-22 11:56:07 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a
|
|
|
|
|
//! expr_t on the new normalized equation
|
2020-03-17 18:58:34 +01:00
|
|
|
|
equation_type_and_normalized_equation_t equation_type_and_normalized_equation;
|
|
|
|
|
|
2020-04-30 16:00:16 +02:00
|
|
|
|
/* Stores derivatives of each block w.r.t. endogenous that belong to it.
|
2020-04-24 12:29:02 +02:00
|
|
|
|
The tuple is: equation number (inside the block), variable number (inside
|
|
|
|
|
the block), lead/lag */
|
|
|
|
|
vector<map<tuple<int, int, int>, expr_t>> blocks_derivatives;
|
2020-03-17 18:58:34 +01:00
|
|
|
|
|
2020-04-21 18:10:46 +02:00
|
|
|
|
class BlockInfo
|
|
|
|
|
{
|
|
|
|
|
public:
|
|
|
|
|
BlockSimulationType simulation_type;
|
|
|
|
|
int first_equation; // Stores a block-ordered equation ID
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int size {0};
|
|
|
|
|
int mfs_size {0}; // Size of the minimal feedback set
|
|
|
|
|
bool linear {true}; // Whether the block is linear in endogenous variable
|
|
|
|
|
int max_endo_lag {0},
|
|
|
|
|
max_endo_lead {
|
|
|
|
|
0}; // Maximum lag/lead on endos that appear in and *that belong to* the block
|
2020-04-23 16:04:02 +02:00
|
|
|
|
|
2023-12-01 14:02:38 +01:00
|
|
|
|
[[nodiscard]] int
|
2022-06-24 15:08:49 +02:00
|
|
|
|
getRecursiveSize() const
|
|
|
|
|
{
|
|
|
|
|
return size - mfs_size;
|
|
|
|
|
};
|
2020-04-21 18:10:46 +02:00
|
|
|
|
};
|
|
|
|
|
|
2022-09-28 16:31:51 +02:00
|
|
|
|
// Whether block decomposition has been successfully computed
|
|
|
|
|
bool block_decomposed {false};
|
|
|
|
|
|
2022-11-30 14:43:44 +01:00
|
|
|
|
/* Whether the block decomposition to compute is time-recursive (i.e. the
|
|
|
|
|
model can be simulated as a whole period-by-period).
|
|
|
|
|
|
|
|
|
|
If true, only contemporaneous occurrences of variables are considered when
|
|
|
|
|
computing the block structure; leads and lags are essentially treated as
|
|
|
|
|
exogenous, i.e. they are ignored. Such a decomposition only makes sense
|
|
|
|
|
for models that are purely backward/forward/static. When using the
|
|
|
|
|
resulting block decomposition to simulate the model, periods must be the
|
|
|
|
|
outer loop, and blocks the inner loop.
|
|
|
|
|
|
|
|
|
|
If false, then the full lead/lag structure is taken into account when
|
|
|
|
|
computing the block structure. This is the only option if there are both
|
|
|
|
|
leads and lags. When using the
|
|
|
|
|
resulting block decomposition to simulate the model, blocks must be the
|
|
|
|
|
outer loop, and periods the inner loop.
|
|
|
|
|
|
|
|
|
|
Of course, this setting does not make any difference on StaticModel. */
|
|
|
|
|
bool time_recursive_block_decomposition {false};
|
|
|
|
|
|
2020-04-21 18:10:46 +02:00
|
|
|
|
// Stores various informations on the blocks
|
|
|
|
|
vector<BlockInfo> blocks;
|
|
|
|
|
|
2020-04-29 16:41:33 +02:00
|
|
|
|
// Maps endogenous type-specific IDs to the block number to which it belongs
|
|
|
|
|
vector<int> endo2block;
|
|
|
|
|
/* Maps (original) equation number to the block number to which it belongs.
|
|
|
|
|
It verifies: ∀i, eq2block[endo2eq[i]] = endo2block[i] */
|
|
|
|
|
vector<int> eq2block;
|
|
|
|
|
|
2020-05-13 16:58:19 +02:00
|
|
|
|
/* Temporary terms for block decomposed models.
|
|
|
|
|
- the outer vector has as many elements as there are blocks in the model
|
|
|
|
|
- the inner vector has as many elements as there are equations in the
|
|
|
|
|
block, plus a last one which contains the temporary terms for
|
|
|
|
|
derivatives
|
|
|
|
|
|
|
|
|
|
It’s necessary to track temporary terms per equation, because some
|
|
|
|
|
equations are evaluated instead of solved, and an equation E1 may depend
|
|
|
|
|
on the value of an endogenous Y computed by a previously evaluated equation
|
|
|
|
|
E2; in this case, if some temporary term TT of equation E2 contains Y,
|
|
|
|
|
then TT needs to be computed after E1, but before E2. */
|
|
|
|
|
vector<vector<temporary_terms_t>> blocks_temporary_terms;
|
2020-05-06 17:13:47 +02:00
|
|
|
|
|
|
|
|
|
/* Stores, for each temporary term in block decomposed models, its index in
|
|
|
|
|
the vector of all temporary terms */
|
|
|
|
|
temporary_terms_idxs_t blocks_temporary_terms_idxs;
|
|
|
|
|
|
2022-09-28 19:18:15 +02:00
|
|
|
|
/* The nonzero values of the sparse block Jacobians in column-major order
|
|
|
|
|
(which is the natural order for Compressed Sparse Column (CSC) storage).
|
|
|
|
|
The pair of indices is (row, column). Nothing is stored for blocks of
|
|
|
|
|
type “evaluate”, since their Jacobian is not used in the sparse
|
|
|
|
|
representation (since the latter does not support the stochastic mode,
|
|
|
|
|
which only exists in the legacy representation). */
|
|
|
|
|
vector<SparseColumnMajorOrderMatrix> blocks_jacobian_sparse_column_major_order;
|
|
|
|
|
/* Column indices for the sparse block Jacobian in Compressed Sparse Column
|
|
|
|
|
(CSC) storage (corresponds to the “jc” vector in MATLAB terminology).
|
|
|
|
|
Same remark as above regarding blocks of type “evaluate”. */
|
|
|
|
|
vector<vector<int>> blocks_jacobian_sparse_colptr;
|
|
|
|
|
|
2018-11-22 14:32:40 +01:00
|
|
|
|
//! Computes derivatives
|
|
|
|
|
/*! \param order the derivation order
|
|
|
|
|
\param vars the derivation IDs w.r.t. which compute the derivatives */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void computeDerivatives(int order, const set<int>& vars);
|
2012-11-29 18:07:48 +01:00
|
|
|
|
//! Computes derivatives of the Jacobian and Hessian w.r. to parameters
|
2016-05-18 12:26:19 +02:00
|
|
|
|
void computeParamsDerivatives(int paramsDerivsOrder);
|
2009-04-20 12:48:54 +02:00
|
|
|
|
//! Computes temporary terms (for all equations and derivatives)
|
2018-09-25 15:56:52 +02:00
|
|
|
|
void computeTemporaryTerms(bool is_matlab, bool no_tmp_terms);
|
2020-05-06 17:13:47 +02:00
|
|
|
|
//! Computes temporary terms per block
|
2022-11-04 12:00:11 +01:00
|
|
|
|
void computeBlockTemporaryTerms(bool no_tmp_terms);
|
2022-09-30 14:44:32 +02:00
|
|
|
|
|
2012-11-29 18:07:48 +01:00
|
|
|
|
//! Computes temporary terms for the file containing parameters derivatives
|
|
|
|
|
void computeParamsDerivativesTemporaryTerms();
|
2017-06-14 07:01:31 +02:00
|
|
|
|
//! Writes temporary terms
|
2022-07-12 14:13:27 +02:00
|
|
|
|
template<ExprNodeOutputType output_type>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union,
|
|
|
|
|
const temporary_terms_idxs_t& tt_idxs, ostream& output,
|
|
|
|
|
deriv_node_temp_terms_t& tef_terms) const;
|
|
|
|
|
void 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;
|
2022-06-23 14:28:13 +02:00
|
|
|
|
//! Writes temporary terms in bytecode
|
2022-07-13 13:04:10 +02:00
|
|
|
|
template<ExprNodeBytecodeOutputType output_type>
|
2023-12-14 14:52:50 +01:00
|
|
|
|
void writeBytecodeTemporaryTerms(const temporary_terms_t& tt,
|
|
|
|
|
temporary_terms_t& temporary_terms_union,
|
|
|
|
|
Bytecode::Writer& code_file,
|
|
|
|
|
deriv_node_temp_terms_t& tef_terms) const;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
/* Adds information for (non-block) bytecode simulation in a separate .bin
|
|
|
|
|
file.
|
|
|
|
|
Returns the number of first derivatives w.r.t. endogenous variables */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int writeBytecodeBinFile(const filesystem::path& filename, bool is_two_boundaries) const;
|
2022-07-19 18:24:36 +02:00
|
|
|
|
//! Adds per-block information for bytecode simulation in a separate .bin file
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int writeBlockBytecodeBinFile(ofstream& bin_file, int block) const;
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
2016-12-30 11:26:56 +01:00
|
|
|
|
//! Fixes output when there are more than 32 nested parens, Issue #1201
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void fixNestedParenthesis(ostringstream& output, map<string, string>& tmp_paren_vars,
|
|
|
|
|
bool& message_printed) const;
|
2016-12-30 11:26:56 +01:00
|
|
|
|
//! Tests if string contains more than 32 nested parens, Issue #1201
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bool testNestedParenthesis(const string& str) const;
|
2022-07-12 14:13:27 +02:00
|
|
|
|
|
2008-02-03 11:28:36 +01:00
|
|
|
|
//! Writes model equations
|
2022-07-12 14:13:27 +02:00
|
|
|
|
template<ExprNodeOutputType output_type>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeModelEquations(ostream& output, const temporary_terms_t& temporary_terms) const;
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
|
|
|
|
// Returns outputs for derivatives and temporary terms at each derivation order
|
|
|
|
|
template<ExprNodeOutputType output_type>
|
|
|
|
|
pair<vector<ostringstream>, vector<ostringstream>> writeModelFileHelper() const;
|
|
|
|
|
|
2022-10-05 17:45:18 +02:00
|
|
|
|
// Writes and compiles dynamic/static file (C version, legacy representation)
|
2022-10-07 17:19:22 +02:00
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeModelCFile(const string& basename, const string& mexext,
|
|
|
|
|
const filesystem::path& matlabroot) const;
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
2022-07-21 18:20:35 +02:00
|
|
|
|
// Writes per-block residuals and temporary terms (incl. for derivatives)
|
|
|
|
|
template<ExprNodeOutputType output_type>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writePerBlockHelper(int blk, ostream& output, temporary_terms_t& temporary_terms) const;
|
2022-07-21 18:20:35 +02:00
|
|
|
|
|
2022-09-28 19:18:15 +02:00
|
|
|
|
// Writes per-block Jacobian (sparse representation)
|
|
|
|
|
// Assumes temporary terms for derivatives are already set.
|
|
|
|
|
template<ExprNodeOutputType output_type>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeSparsePerBlockJacobianHelper(int blk, ostream& output,
|
|
|
|
|
temporary_terms_t& temporary_terms) const;
|
2022-09-28 19:18:15 +02:00
|
|
|
|
|
2022-07-12 12:27:22 +02:00
|
|
|
|
/* Helper for writing derivatives w.r.t. parameters.
|
|
|
|
|
Returns { tt, rp, gp, rpp, gpp, hp, g3p }.
|
|
|
|
|
g3p is empty if requesting a static output type. */
|
|
|
|
|
template<ExprNodeOutputType output_type>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
tuple<ostringstream, ostringstream, ostringstream, ostringstream, ostringstream, ostringstream,
|
|
|
|
|
ostringstream>
|
|
|
|
|
writeParamsDerivativesFileHelper() const;
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2022-07-13 13:04:10 +02:00
|
|
|
|
// Helper for writing bytecode (without block decomposition)
|
|
|
|
|
template<bool dynamic>
|
2023-12-14 14:52:50 +01:00
|
|
|
|
void writeBytecodeHelper(Bytecode::Writer& code_file) const;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
2022-07-19 18:24:36 +02:00
|
|
|
|
// Helper for writing blocks in bytecode
|
|
|
|
|
template<bool dynamic>
|
2023-12-14 14:52:50 +01:00
|
|
|
|
void writeBlockBytecodeHelper(Bytecode::Writer& code_file, int block,
|
2023-11-30 15:28:57 +01:00
|
|
|
|
temporary_terms_t& temporary_terms_union) const;
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
// Helper for writing sparse derivatives indices in MATLAB/Octave driver file
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeDriverSparseIndicesHelper(ostream& output) const;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
|
|
|
|
// Helper for writing sparse derivatives indices in JSON
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeJsonSparseIndicesHelper(ostream& output) const;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
2022-09-28 19:18:15 +02:00
|
|
|
|
// Helper for writing sparse block derivatives indices in MATLAB/Octave driver file
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeBlockDriverSparseIndicesHelper(ostream& output) const;
|
2022-09-28 19:18:15 +02:00
|
|
|
|
|
2022-07-12 17:47:02 +02:00
|
|
|
|
/* Helper for writing JSON output for residuals and derivatives.
|
|
|
|
|
Returns mlv and derivatives output at each derivation order. */
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
pair<ostringstream, vector<ostringstream>>
|
|
|
|
|
writeJsonComputingPassOutputHelper(bool writeDetails) const;
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
/* Helper for writing JSON derivatives w.r.t. parameters.
|
|
|
|
|
Returns { mlv, tt, rp, gp, rpp, gpp, hp, g3p }.
|
|
|
|
|
g3p is empty if requesting a static output type. */
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
tuple<ostringstream, ostringstream, ostringstream, ostringstream, ostringstream, ostringstream,
|
|
|
|
|
ostringstream, ostringstream>
|
|
|
|
|
writeJsonParamsDerivativesHelper(bool writeDetails) const;
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2017-02-02 15:09:43 +01:00
|
|
|
|
//! Writes JSON model equations
|
2017-02-20 12:18:11 +01:00
|
|
|
|
//! if residuals = true, we are writing the dynamic/static model.
|
|
|
|
|
//! Otherwise, just the model equations (with line numbers, no tmp terms)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeJsonModelEquations(ostream& output, bool residuals) const;
|
2020-06-05 16:07:52 +02:00
|
|
|
|
/* Writes JSON model local variables.
|
|
|
|
|
Optionally put the external function variable calls into TEF terms */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeJsonModelLocalVariables(ostream& output, bool write_tef_terms,
|
|
|
|
|
deriv_node_temp_terms_t& tef_terms) const;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
2022-06-23 14:28:13 +02:00
|
|
|
|
//! Writes model equations in bytecode
|
2022-07-13 13:04:10 +02:00
|
|
|
|
template<ExprNodeBytecodeOutputType output_type>
|
2023-12-14 14:52:50 +01:00
|
|
|
|
void writeBytecodeModelEquations(Bytecode::Writer& code_file,
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const temporary_terms_t& temporary_terms,
|
|
|
|
|
const deriv_node_temp_terms_t& tef_terms) const;
|
2008-02-03 11:28:36 +01:00
|
|
|
|
|
2022-09-30 17:26:43 +02:00
|
|
|
|
// Writes the sparse representation of the model in MATLAB/Octave
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeSparseModelMFiles(const string& basename) const;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
|
2022-10-05 17:45:18 +02:00
|
|
|
|
// Writes and compiles the sparse representation of the model in C
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeSparseModelCFiles(const string& basename, const string& mexext,
|
|
|
|
|
const filesystem::path& matlabroot) const;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
2022-09-20 18:28:30 +02:00
|
|
|
|
// Writes the sparse representation of the model in Julia
|
|
|
|
|
// Assumes that the directory <MODFILE>/model/julia/ already exists
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeSparseModelJuliaFiles(const string& basename) const;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
|
2009-04-30 15:14:33 +02:00
|
|
|
|
//! Writes LaTeX model file
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void 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
|
|
|
|
|
2022-11-09 14:46:02 +01:00
|
|
|
|
/* Write files for helping a user to debug their model (MATLAB/Octave,
|
|
|
|
|
sparse representation).
|
|
|
|
|
Creates a dynamic/static files which evaluates separately the LHS and RHS
|
|
|
|
|
of each equation.
|
|
|
|
|
They are not optimized for performance (hence in particular the absence of
|
|
|
|
|
a C version, or the non-reuse of temporary terms). */
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeDebugModelMFiles(const string& basename) const;
|
2022-11-09 14:46:02 +01:00
|
|
|
|
|
2023-03-28 16:37:05 +02:00
|
|
|
|
// Write the file that sets auxiliary variables given the original variables
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
|
2023-03-28 16:37:05 +02:00
|
|
|
|
|
2022-09-28 16:31:51 +02:00
|
|
|
|
private:
|
2020-04-28 18:00:26 +02:00
|
|
|
|
//! Sparse matrix of double to store the values of the static Jacobian
|
2009-12-16 14:21:31 +01:00
|
|
|
|
/*! First index is equation number, second index is endogenous type specific ID */
|
2018-06-04 14:31:26 +02:00
|
|
|
|
using jacob_map_t = map<pair<int, int>, double>;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2020-03-24 16:57:45 +01:00
|
|
|
|
//! Normalization of equations, as computed by computeNonSingularNormalization()
|
2009-12-16 14:21:31 +01:00
|
|
|
|
/*! Maps endogenous type specific IDs to equation numbers */
|
|
|
|
|
vector<int> endo2eq;
|
|
|
|
|
|
2022-10-14 14:52:11 +02:00
|
|
|
|
// Stores workers used for compiling MEX files in parallel
|
|
|
|
|
static vector<jthread> mex_compilation_workers;
|
2022-09-30 15:39:41 +02:00
|
|
|
|
|
2022-10-05 18:34:21 +02:00
|
|
|
|
/* The following variables implement the thread synchronization mechanism for
|
|
|
|
|
limiting the number of concurrent GCC processes and tracking dependencies
|
|
|
|
|
between object files. */
|
2022-10-17 13:57:34 +02:00
|
|
|
|
static condition_variable_any mex_compilation_cv;
|
2022-10-05 16:38:16 +02:00
|
|
|
|
static mutex mex_compilation_mut;
|
2022-10-14 14:52:11 +02:00
|
|
|
|
/* Object/MEX files waiting to be compiled (with their prerequisites as 2nd
|
|
|
|
|
element and compilation command as the 3rd element) */
|
|
|
|
|
static vector<tuple<filesystem::path, set<filesystem::path>, string>> mex_compilation_queue;
|
2022-12-20 14:47:32 +01:00
|
|
|
|
// Object/MEX files in the process of being compiled
|
|
|
|
|
static set<filesystem::path> mex_compilation_ongoing;
|
|
|
|
|
// Object/MEX files already compiled successfully
|
2022-10-14 14:52:11 +02:00
|
|
|
|
static set<filesystem::path> mex_compilation_done;
|
2022-12-20 14:47:32 +01:00
|
|
|
|
// Object/MEX files whose compilation failed
|
|
|
|
|
static set<filesystem::path> mex_compilation_failed;
|
2022-10-05 16:38:16 +02:00
|
|
|
|
|
2020-04-15 17:56:28 +02:00
|
|
|
|
/* Compute a pseudo-Jacobian whose all elements are either zero or one,
|
2022-11-28 13:40:49 +01:00
|
|
|
|
depending on whether the variable symbolically appears in the equation. If
|
|
|
|
|
contemporaneous_only=true, only considers contemporaneous occurences of
|
|
|
|
|
variables; otherwise also considers leads and lags. */
|
|
|
|
|
jacob_map_t computeSymbolicJacobian(bool contemporaneous_only) const;
|
2020-04-15 17:56:28 +02:00
|
|
|
|
|
2023-04-24 17:49:54 +02:00
|
|
|
|
/* Compute a pseudo-Jacobian whose all elements are either zero or one.
|
|
|
|
|
For the equations that were originally written by the user (identified as
|
|
|
|
|
those having an associated line number), checks whether there is a single
|
|
|
|
|
contemporaneous endogenous on the left-hand side; if yes, only this
|
|
|
|
|
endogenous is associated with a one on the line of the corresponding
|
|
|
|
|
equation; otherwise, returns false as the first output argument and
|
|
|
|
|
aborts the computation.
|
|
|
|
|
For the other equations, fills the corresponding lines as is done
|
|
|
|
|
by computeSymbolicJacobian(true). */
|
|
|
|
|
pair<bool, jacob_map_t> computeLeftHandSideSymbolicJacobian() const;
|
|
|
|
|
|
2020-04-17 14:55:55 +02:00
|
|
|
|
// Compute {var,eq}_idx_orig2block from {var,eq}_idx_block2orig
|
2020-04-15 17:56:28 +02:00
|
|
|
|
void updateReverseVariableEquationOrderings();
|
|
|
|
|
|
2023-04-24 17:47:29 +02:00
|
|
|
|
struct ModelNormalizationFailed
|
|
|
|
|
{
|
|
|
|
|
string unmatched_endo; // Name of endogenous not in maximum cardinality matching
|
|
|
|
|
};
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
//! Compute the matching between endogenous and variable using the jacobian
|
|
|
|
|
//! contemporaneous_jacobian
|
2009-12-21 11:29:21 +01:00
|
|
|
|
/*!
|
2023-11-30 15:28:57 +01:00
|
|
|
|
\param contemporaneous_jacobian Jacobian used as an incidence matrix: all elements declared in
|
|
|
|
|
the map (even if they are zero), are used as vertices of the incidence matrix \return True if a
|
|
|
|
|
complete normalization has been achieved
|
2009-12-21 11:29:21 +01:00
|
|
|
|
*/
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void computeNormalization(const jacob_map_t& contemporaneous_jacobian);
|
2009-12-21 11:29:21 +01:00
|
|
|
|
|
2009-12-16 14:21:31 +01:00
|
|
|
|
//! Try to compute the matching between endogenous and variable using a decreasing cutoff
|
2009-12-21 11:29:21 +01:00
|
|
|
|
/*!
|
|
|
|
|
Applied to the jacobian contemporaneous_jacobian and stop when a matching is found.
|
2023-11-30 15:28:57 +01:00
|
|
|
|
If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e.
|
|
|
|
|
use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to
|
|
|
|
|
reflect all the edges in the symbolic incidence matrix. If no matching is found with a zero
|
|
|
|
|
cutoff, an error message is printed. The resulting normalization is stored in endo2eq. Returns a
|
|
|
|
|
boolean indicating success.
|
2009-12-21 11:29:21 +01:00
|
|
|
|
*/
|
2023-11-30 15:28:57 +01:00
|
|
|
|
bool computeNonSingularNormalization(const eval_context_t& eval_context);
|
2020-03-20 18:42:59 +01:00
|
|
|
|
//! Evaluate the jacobian (w.r.t. endogenous) and suppress all the elements below the cutoff
|
2020-05-06 18:09:44 +02:00
|
|
|
|
/*! Returns the contemporaneous_jacobian.
|
2020-04-28 18:00:26 +02:00
|
|
|
|
Elements below the cutoff are discarded. External functions are evaluated to 1. */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
jacob_map_t evaluateAndReduceJacobian(const eval_context_t& eval_context) const;
|
2020-04-30 12:48:16 +02:00
|
|
|
|
/* Search the equations and variables belonging to the prologue and the
|
|
|
|
|
epilogue of the model.
|
|
|
|
|
Initializes “eq_idx_block2orig” and “endo_idx_block2orig”.
|
|
|
|
|
Returns the sizes of the prologue and epilogue. */
|
|
|
|
|
pair<int, int> computePrologueAndEpilogue();
|
2020-03-20 18:47:29 +01:00
|
|
|
|
//! Determine the type of each equation of model and try to normalize the unnormalized equation
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void
|
|
|
|
|
equationTypeDetermination(const map<tuple<int, int, int>, expr_t>& first_order_endo_derivatives);
|
2020-04-30 11:15:55 +02:00
|
|
|
|
/* Fills the max lags/leads and n_{static,mixed,forward,backward} fields of a
|
|
|
|
|
given block.
|
|
|
|
|
Needs the fields size and first_equation. */
|
|
|
|
|
void computeDynamicStructureOfBlock(int blk);
|
2020-04-30 12:48:16 +02:00
|
|
|
|
/* Fills the simulation_type field of a given block.
|
|
|
|
|
Needs the fields size, max_endo_lag and max_endo_lead. */
|
|
|
|
|
void computeSimulationTypeOfBlock(int blk);
|
2020-04-21 18:10:46 +02:00
|
|
|
|
/* Compute the block decomposition and for a non-recusive block find the minimum feedback set
|
|
|
|
|
|
2020-04-30 12:48:16 +02:00
|
|
|
|
Initializes the “blocks”, “endo2block” and “eq2block” structures. */
|
|
|
|
|
void computeBlockDecomposition(int prologue, int epilogue);
|
2020-04-21 18:10:46 +02:00
|
|
|
|
/* Reduce the number of block by merging the same type of equations in the
|
2020-04-30 12:48:16 +02:00
|
|
|
|
prologue and the epilogue */
|
|
|
|
|
void reduceBlockDecomposition();
|
2020-04-17 19:21:37 +02:00
|
|
|
|
/* The 1st output gives, for each equation (in original order) the (max_lag,
|
|
|
|
|
max_lead) across all endogenous that appear in the equation and that
|
|
|
|
|
belong to the same block (i.e. those endogenous are solved in the same
|
|
|
|
|
block).
|
|
|
|
|
|
|
|
|
|
The 2nd output gives, for each type-specific endo IDs, its (max_lag,
|
|
|
|
|
max_lead) across all its occurences inside the equations of the block to
|
|
|
|
|
which it belongs. */
|
2020-04-29 16:41:33 +02:00
|
|
|
|
pair<lag_lead_vector_t, lag_lead_vector_t> getVariableLeadLagByBlock() const;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
//! Print an abstract of the block structure of the model
|
2020-04-15 17:56:28 +02:00
|
|
|
|
void printBlockDecomposition() const;
|
2009-12-16 14:21:31 +01:00
|
|
|
|
//! Determine for each block if it is linear or not
|
2020-03-20 15:23:23 +01:00
|
|
|
|
void determineLinearBlocks();
|
2009-12-16 14:21:31 +01:00
|
|
|
|
|
2022-09-28 16:31:51 +02:00
|
|
|
|
protected:
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the type of equation belonging to the block
|
2020-04-17 16:04:19 +02:00
|
|
|
|
EquationType
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockEquationType(int blk, int eq) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
|
|
|
|
|
.first;
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2009-12-16 14:21:31 +01:00
|
|
|
|
//! Return true if the equation has been normalized
|
2020-04-17 16:04:19 +02:00
|
|
|
|
bool
|
2020-04-23 16:04:02 +02:00
|
|
|
|
isBlockEquationRenormalized(int blk, int eq) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
|
|
|
|
|
.first
|
|
|
|
|
== EquationType::evaluateRenormalized;
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the expr_t of equation belonging to the block
|
2023-11-30 15:28:57 +01:00
|
|
|
|
BinaryOpNode*
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockEquationExpr(int blk, int eq) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2020-04-23 16:04:02 +02:00
|
|
|
|
return equations[eq_idx_block2orig[blocks[blk].first_equation + eq]];
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the expr_t of renormalized equation belonging to the block
|
2023-11-30 15:28:57 +01:00
|
|
|
|
BinaryOpNode*
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockEquationRenormalizedExpr(int blk, int eq) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
|
|
|
|
|
.second;
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the original number of equation belonging to the block
|
2020-04-17 16:04:19 +02:00
|
|
|
|
int
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockEquationID(int blk, int eq) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2020-04-23 16:04:02 +02:00
|
|
|
|
return eq_idx_block2orig[blocks[blk].first_equation + eq];
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the original number of variable belonging to the block
|
2020-04-17 16:04:19 +02:00
|
|
|
|
int
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockVariableID(int blk, int var) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2020-04-23 16:04:02 +02:00
|
|
|
|
return endo_idx_block2orig[blocks[blk].first_equation + var];
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the position of an equation (given by its original index) inside its block
|
2020-04-17 16:04:19 +02:00
|
|
|
|
int
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockInitialEquationID(int blk, int eq) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2020-04-23 16:04:02 +02:00
|
|
|
|
return eq_idx_orig2block[eq] - blocks[blk].first_equation;
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2020-04-23 16:04:02 +02:00
|
|
|
|
//! Return the position of a variable (given by its original index) inside its block
|
2020-04-17 16:04:19 +02:00
|
|
|
|
int
|
2020-04-23 16:04:02 +02:00
|
|
|
|
getBlockInitialVariableID(int blk, int var) const
|
2020-04-17 16:04:19 +02:00
|
|
|
|
{
|
2020-04-23 16:04:02 +02:00
|
|
|
|
return endo_idx_orig2block[var] - blocks[blk].first_equation;
|
2020-04-17 16:04:19 +02:00
|
|
|
|
};
|
2011-06-22 11:56:07 +02:00
|
|
|
|
//! Initialize equation_reordered & variable_reordered
|
|
|
|
|
void initializeVariablesAndEquations();
|
2022-09-28 16:31:51 +02:00
|
|
|
|
|
|
|
|
|
private:
|
2020-03-30 14:51:53 +02:00
|
|
|
|
//! Returns the 1st derivatives w.r.t. endogenous in a different format
|
|
|
|
|
/*! Returns a map (equation, type-specific ID, lag) → derivative.
|
|
|
|
|
Assumes that derivatives have already been computed. */
|
|
|
|
|
map<tuple<int, int, int>, expr_t> collectFirstOrderDerivativesEndogenous();
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
2022-09-28 16:31:51 +02:00
|
|
|
|
protected:
|
|
|
|
|
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
|
|
|
|
|
virtual void computeChainRuleJacobian() = 0;
|
|
|
|
|
|
|
|
|
|
/* Compute block decomposition, its derivatives and temporary terms. Meant to
|
|
|
|
|
be overriden in derived classes which don’t support block decomposition
|
|
|
|
|
(currently Epilogue and PlannerObjective). Sets “block_decomposed” to true
|
|
|
|
|
in case of success. */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
virtual void computingPassBlock(const eval_context_t& eval_context, bool no_tmp_terms);
|
2022-09-28 16:31:51 +02:00
|
|
|
|
|
2022-07-19 18:24:36 +02:00
|
|
|
|
/* Get column number within Jacobian of a given block.
|
|
|
|
|
“var” is the block-specific endogenous variable index. */
|
|
|
|
|
virtual int getBlockJacobianEndoCol(int blk, int var, int lag) const = 0;
|
|
|
|
|
|
2022-09-14 17:48:33 +02:00
|
|
|
|
// Returns a human-readable string describing the model class (e.g. “dynamic model”…)
|
|
|
|
|
virtual string modelClassName() const = 0;
|
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
/* Given a sparse matrix in column major order, returns the colptr pointer for
|
|
|
|
|
the CSC storage */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
static vector<int> computeCSCColPtr(const SparseColumnMajorOrderMatrix& matrix, int ncols);
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
2018-10-09 18:27:19 +02:00
|
|
|
|
private:
|
|
|
|
|
//! Internal helper for the copy constructor and assignment operator
|
|
|
|
|
/*! Copies all the structures that contain ExprNode*, by the converting the
|
|
|
|
|
pointers into their equivalent in the new tree */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void copyHelper(const ModelTree& m);
|
2018-10-26 11:44:26 +02:00
|
|
|
|
//! Returns the name of the MATLAB architecture given the extension used for MEX files
|
2023-11-30 15:28:57 +01:00
|
|
|
|
static string matlab_arch(const string& mexext);
|
2023-12-01 15:31:55 +01:00
|
|
|
|
#ifdef __APPLE__
|
2023-07-11 22:12:03 +02:00
|
|
|
|
/* Finds a suitable compiler on macOS.
|
|
|
|
|
The boolean is false if this is GCC and true if this is Clang */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
static pair<filesystem::path, bool> findCompilerOnMacos(const string& mexext);
|
2023-12-01 15:31:55 +01:00
|
|
|
|
#endif
|
2022-10-05 18:34:21 +02:00
|
|
|
|
/* Compiles a MEX file (if link=true) or an object file to be linked later
|
2022-10-14 14:52:11 +02:00
|
|
|
|
into a MEX file (if link=false). The compilation is done in separate
|
|
|
|
|
worker threads working in parallel, so the call to this function is not
|
|
|
|
|
blocking. The dependency of a linked MEX file upon intermediary objects is
|
|
|
|
|
nicely handled. Returns the name of the output file (to be reused later as
|
|
|
|
|
input file if link=false). */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
filesystem::path 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 = true) const;
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
2008-02-03 11:28:36 +01:00
|
|
|
|
public:
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
|
|
|
|
|
ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg = false);
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
2022-05-18 18:24:13 +02:00
|
|
|
|
protected:
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree(const ModelTree& m);
|
|
|
|
|
ModelTree& operator=(const ModelTree& m);
|
2018-10-09 18:27:19 +02:00
|
|
|
|
|
2022-05-18 18:24:13 +02:00
|
|
|
|
public:
|
2011-06-22 11:56:07 +02:00
|
|
|
|
//! Absolute value under which a number is considered to be zero
|
2023-11-30 15:28:57 +01:00
|
|
|
|
double cutoff {1e-15};
|
2014-01-27 16:41:43 +01:00
|
|
|
|
//! Declare a node as an equation of the model; also give its line number
|
2023-12-13 15:37:07 +01:00
|
|
|
|
void addEquation(expr_t eq, const optional<int>& lineno);
|
2013-04-11 17:07:39 +02:00
|
|
|
|
//! Declare a node as an equation of the model, also giving its tags
|
2023-12-13 15:37:07 +01:00
|
|
|
|
void addEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
//! Declare a node as an auxiliary equation of the model, adding it at the end of the list of
|
|
|
|
|
//! auxiliary equations
|
2010-09-16 19:18:45 +02:00
|
|
|
|
void addAuxEquation(expr_t eq);
|
2008-10-17 19:21:22 +02:00
|
|
|
|
//! Returns the number of equations in the model
|
2008-02-03 11:28:36 +01:00
|
|
|
|
int equation_number() const;
|
2010-10-15 19:05:16 +02:00
|
|
|
|
//! Adds a trend variable with its growth factor
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void addTrendVariables(const vector<int>& trend_vars, expr_t growth_factor) noexcept(false);
|
2013-03-26 16:46:18 +01:00
|
|
|
|
//! Adds a nonstationary variables with their (common) deflator
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void addNonstationaryVariables(const vector<int>& nonstationary_vars, bool log_deflator,
|
|
|
|
|
expr_t deflator) noexcept(false);
|
2013-10-29 11:47:59 +01:00
|
|
|
|
//! Is a given variable non-stationary?
|
|
|
|
|
bool isNonstationary(int symb_id) const;
|
2011-06-22 11:56:07 +02:00
|
|
|
|
void set_cutoff_to_zero();
|
2019-12-03 14:19:32 +01:00
|
|
|
|
/*! Reorder auxiliary variables so that they appear in recursive order in
|
|
|
|
|
set_auxiliary_variables.m and dynamic_set_auxiliary_series.m */
|
|
|
|
|
void reorderAuxiliaryEquations();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
//! Find equations of the form “variable=constant”, excluding equations with “mcp” tag (see
|
|
|
|
|
//! dynare#1697)
|
|
|
|
|
void findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const;
|
2021-10-26 18:06:26 +02:00
|
|
|
|
/* Given an expression, searches for the first equation that has exactly this
|
|
|
|
|
expression on the LHS, and returns the RHS of that equation.
|
|
|
|
|
If no such equation can be found, throws an ExprNode::MatchFailureExpression */
|
|
|
|
|
expr_t getRHSFromLHS(expr_t lhs) const;
|
2019-12-13 17:15:03 +01:00
|
|
|
|
|
2023-03-20 17:54:32 +01:00
|
|
|
|
/* Initialize the MEX compilation workers (and some environment variables
|
|
|
|
|
needed for finding GCC) */
|
2023-11-30 15:28:57 +01:00
|
|
|
|
static void initializeMEXCompilationWorkers(int numworkers, const filesystem::path& dynareroot,
|
|
|
|
|
const string& mexext);
|
2022-10-14 14:52:11 +02:00
|
|
|
|
|
2022-10-17 13:57:34 +02:00
|
|
|
|
// Waits until the MEX compilation queue is empty
|
|
|
|
|
static void waitForMEXCompilationWorkers();
|
2022-09-30 15:39:41 +02:00
|
|
|
|
|
2023-03-28 16:37:05 +02:00
|
|
|
|
// Write the definitions of the auxiliary variables (assumed to be in recursive order)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
void writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType output_type) const;
|
2023-03-28 16:37:05 +02:00
|
|
|
|
|
2021-04-22 14:52:29 +02:00
|
|
|
|
//! Returns the vector of non-zero derivative counts
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const vector<int>&
|
2021-04-22 14:52:29 +02:00
|
|
|
|
getNNZDerivatives() const
|
|
|
|
|
{
|
|
|
|
|
return NNZDerivatives;
|
|
|
|
|
}
|
|
|
|
|
|
2021-06-11 12:37:31 +02:00
|
|
|
|
//! Returns the vector of temporary terms derivatives
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const vector<temporary_terms_t>&
|
2021-06-11 12:37:31 +02:00
|
|
|
|
getTemporaryTermsDerivatives() const
|
|
|
|
|
{
|
|
|
|
|
return temporary_terms_derivatives;
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
//! Returns the maximum order of computed derivatives
|
2022-06-24 15:08:49 +02:00
|
|
|
|
int
|
2021-04-22 14:52:29 +02:00
|
|
|
|
getComputedDerivsOrder() const
|
|
|
|
|
{
|
|
|
|
|
return computed_derivs_order;
|
|
|
|
|
}
|
|
|
|
|
|
2022-06-24 15:08:49 +02:00
|
|
|
|
static string
|
2020-03-20 17:31:14 +01:00
|
|
|
|
BlockSim(BlockSimulationType type)
|
2009-12-16 18:13:23 +01:00
|
|
|
|
{
|
|
|
|
|
switch (type)
|
|
|
|
|
{
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::evaluateForward:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "EVALUATE FORWARD ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::evaluateBackward:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "EVALUATE BACKWARD ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "SOLVE FORWARD SIMPLE ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "SOLVE BACKWARD SIMPLE ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::solveTwoBoundariesSimple:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "SOLVE TWO BOUNDARIES SIMPLE ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "SOLVE FORWARD COMPLETE ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "SOLVE BACKWARD COMPLETE ";
|
2020-03-20 17:31:14 +01:00
|
|
|
|
case BlockSimulationType::solveTwoBoundariesComplete:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "SOLVE TWO BOUNDARIES COMPLETE";
|
2009-12-16 18:13:23 +01:00
|
|
|
|
default:
|
2019-12-16 19:42:59 +01:00
|
|
|
|
return "UNKNOWN ";
|
2009-12-16 18:13:23 +01:00
|
|
|
|
}
|
2020-03-20 18:00:56 +01:00
|
|
|
|
}
|
2023-10-16 16:24:31 +02:00
|
|
|
|
|
|
|
|
|
/* Returns the minimum feedback set value (see the “mfs” option of the
|
|
|
|
|
“model” block in the reference manual for the possible values) */
|
|
|
|
|
virtual int getMFS() const = 0;
|
2008-02-03 11:28:36 +01:00
|
|
|
|
};
|
|
|
|
|
|
2022-07-12 14:13:27 +02:00
|
|
|
|
template<ExprNodeOutputType output_type>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union,
|
|
|
|
|
const temporary_terms_idxs_t& tt_idxs, ostream& output,
|
|
|
|
|
deriv_node_temp_terms_t& tef_terms) const
|
2022-07-12 14:13:27 +02:00
|
|
|
|
{
|
|
|
|
|
for (auto it : tt)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (dynamic_cast<AbstractExternalFunctionNode*>(it))
|
2022-07-12 14:13:27 +02:00
|
|
|
|
it->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
|
|
|
|
|
|
2023-12-04 16:37:00 +01:00
|
|
|
|
// NOLINTNEXTLINE(clang-analyzer-core.CallAndMessage)
|
2022-07-12 14:13:27 +02:00
|
|
|
|
it->writeOutput(output, output_type, tt, tt_idxs, tef_terms);
|
|
|
|
|
output << " = ";
|
|
|
|
|
it->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (isCOutput(output_type) || isMatlabOutput(output_type))
|
2022-07-12 14:13:27 +02:00
|
|
|
|
output << ";";
|
|
|
|
|
output << endl;
|
|
|
|
|
|
|
|
|
|
temp_term_union.insert(it);
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<ExprNodeOutputType output_type>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeModelEquations(ostream& output, const temporary_terms_t& temporary_terms) const
|
2022-07-12 14:13:27 +02:00
|
|
|
|
{
|
|
|
|
|
for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
BinaryOpNode* eq_node {equations[eq]};
|
|
|
|
|
expr_t lhs {eq_node->arg1}, rhs {eq_node->arg2};
|
2022-07-12 14:13:27 +02:00
|
|
|
|
|
|
|
|
|
// Test if the right hand side of the equation is empty.
|
|
|
|
|
double vrhs {1.0};
|
|
|
|
|
try
|
|
|
|
|
{
|
|
|
|
|
vrhs = rhs->eval({});
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ExprNode::EvalException& e)
|
2022-07-12 14:13:27 +02:00
|
|
|
|
{
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
|
2023-01-04 14:43:57 +01:00
|
|
|
|
{
|
|
|
|
|
output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
|
2023-01-04 14:43:57 +01:00
|
|
|
|
<< " = (";
|
|
|
|
|
lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
|
|
|
|
|
output << ") - (";
|
|
|
|
|
rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
|
|
|
|
|
output << ");" << endl;
|
|
|
|
|
}
|
2022-07-12 14:13:27 +02:00
|
|
|
|
else // The right hand side of the equation is empty ==> residual=lhs;
|
|
|
|
|
{
|
|
|
|
|
output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
|
2022-07-12 14:13:27 +02:00
|
|
|
|
<< " = ";
|
|
|
|
|
lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
|
|
|
|
|
output << ";" << endl;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2022-07-11 17:33:09 +02:00
|
|
|
|
template<ExprNodeOutputType output_type>
|
|
|
|
|
pair<vector<ostringstream>, vector<ostringstream>>
|
|
|
|
|
ModelTree::writeModelFileHelper() const
|
|
|
|
|
{
|
2022-09-14 17:07:08 +02:00
|
|
|
|
constexpr bool sparse {isSparseModelOutput(output_type)};
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<ostringstream> d_output(
|
|
|
|
|
derivatives.size()); // Derivatives output (at all orders, including 0=residual)
|
2022-07-11 17:33:09 +02:00
|
|
|
|
vector<ostringstream> tt_output(derivatives.size()); // Temp terms output (at all orders)
|
|
|
|
|
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
temporary_terms_t temp_term_union;
|
|
|
|
|
|
2022-07-12 14:13:27 +02:00
|
|
|
|
writeTemporaryTerms<output_type>(temporary_terms_derivatives[0], temp_term_union,
|
|
|
|
|
temporary_terms_idxs, tt_output[0], tef_terms);
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
2022-07-12 14:13:27 +02:00
|
|
|
|
writeModelEquations<output_type>(d_output[0], temp_term_union);
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
|
|
|
|
// Writing Jacobian
|
|
|
|
|
if (!derivatives[1].empty())
|
|
|
|
|
{
|
2022-07-12 14:13:27 +02:00
|
|
|
|
writeTemporaryTerms<output_type>(temporary_terms_derivatives[1], temp_term_union,
|
|
|
|
|
temporary_terms_idxs, tt_output[1], tef_terms);
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (sparse)
|
2022-07-11 17:33:09 +02:00
|
|
|
|
{
|
2022-09-14 17:07:08 +02:00
|
|
|
|
// NB: we iterate over the Jacobian reordered in column-major order
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// Indices of rows and columns are output in M_ and the JSON file (since they are
|
|
|
|
|
// constant)
|
|
|
|
|
for (int k {0}; const auto& [row_col, d1] : jacobian_sparse_column_major_order)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
d_output[1] << "g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d1->writeOutput(d_output[1], output_type, temp_term_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[1] << ";" << endl;
|
|
|
|
|
k++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
else // Legacy representation (dense matrix)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : derivatives[1])
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
auto [eq, var] = vectorToTuple<2>(indices);
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (isMatlabOutput(output_type) || isJuliaOutput(output_type))
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[1] << eq + 1 << "," << getJacobianCol(var, sparse) + 1;
|
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d_output[1] << eq + getJacobianCol(var, sparse) * equations.size();
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d1->writeOutput(d_output[1], output_type, temp_term_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[1] << ";" << endl;
|
|
|
|
|
}
|
2022-07-11 17:33:09 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// Write derivatives for order ≥ 2
|
|
|
|
|
for (size_t i = 2; i < derivatives.size(); i++)
|
|
|
|
|
if (!derivatives[i].empty())
|
|
|
|
|
{
|
2022-07-12 14:13:27 +02:00
|
|
|
|
writeTemporaryTerms<output_type>(temporary_terms_derivatives[i], temp_term_union,
|
|
|
|
|
temporary_terms_idxs, tt_output[i], tef_terms);
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (sparse)
|
2022-07-11 17:33:09 +02:00
|
|
|
|
{
|
2022-09-14 17:07:08 +02:00
|
|
|
|
/* List non-zero elements of the tensor in row-major order (this is
|
|
|
|
|
suitable for the k-order solver according to Normann). */
|
|
|
|
|
// Tensor indices are output in M_ and the JSON file (since they are constant)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int k {0}; const auto& [vidx, d] : derivatives[i])
|
2022-07-11 17:33:09 +02:00
|
|
|
|
{
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[i] << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[i] << ";" << endl;
|
2022-07-11 17:33:09 +02:00
|
|
|
|
k++;
|
|
|
|
|
}
|
2022-09-14 17:07:08 +02:00
|
|
|
|
}
|
|
|
|
|
else // Legacy representation
|
|
|
|
|
{
|
|
|
|
|
/* When creating the sparse matrix (in MATLAB or C mode), since storage
|
|
|
|
|
is in column-major order, output the first column, then the second,
|
|
|
|
|
then the third. This gives a significant performance boost in use_dll
|
|
|
|
|
mode (at both compilation and runtime), because it facilitates memory
|
|
|
|
|
accesses and expression reusage. */
|
|
|
|
|
ostringstream i_output, j_output, v_output;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int k {0}; // Current line index in the 3-column matrix
|
|
|
|
|
const auto& [vidx, d] : derivatives[i])
|
2022-07-11 17:33:09 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int eq {vidx[0]};
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int col_idx {0};
|
2022-09-14 17:07:08 +02:00
|
|
|
|
for (size_t j = 1; j < vidx.size(); j++)
|
|
|
|
|
{
|
|
|
|
|
col_idx *= getJacobianColsNbr(sparse);
|
|
|
|
|
col_idx += getJacobianCol(vidx[j], sparse);
|
|
|
|
|
}
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (isJuliaOutput(output_type))
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
d_output[i] << " g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[i] << endl;
|
|
|
|
|
}
|
2022-07-11 17:33:09 +02:00
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl;
|
2022-07-11 17:33:09 +02:00
|
|
|
|
j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << col_idx + 1 << ";"
|
|
|
|
|
<< endl;
|
2022-07-11 17:33:09 +02:00
|
|
|
|
v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-09-14 17:07:08 +02:00
|
|
|
|
v_output << ";" << endl;
|
2022-07-11 17:33:09 +02:00
|
|
|
|
|
|
|
|
|
k++;
|
|
|
|
|
}
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
|
|
|
|
// Output symetric elements at order 2
|
|
|
|
|
if (i == 2 && vidx[1] != vidx[2])
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int col_idx_sym {getJacobianCol(vidx[2], sparse) * getJacobianColsNbr(sparse)
|
|
|
|
|
+ getJacobianCol(vidx[1], sparse)};
|
2022-09-14 17:07:08 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (isJuliaOutput(output_type))
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[2] << " g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
|
|
|
|
|
<< "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";"
|
|
|
|
|
<< endl;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << col_idx_sym + 1
|
|
|
|
|
<< ";" << endl;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "="
|
|
|
|
|
<< "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< k - 1 + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
|
|
|
|
|
|
|
|
|
|
k++;
|
|
|
|
|
}
|
|
|
|
|
}
|
2022-07-11 17:33:09 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (!isJuliaOutput(output_type))
|
2022-09-14 17:07:08 +02:00
|
|
|
|
d_output[i] << i_output.str() << j_output.str() << v_output.str();
|
2022-07-11 17:33:09 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (isMatlabOutput(output_type))
|
2022-07-12 12:19:17 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// Check that we don't have more than 32 nested parenthesis because MATLAB does not suppor
|
|
|
|
|
// this. See Issue #1201
|
2022-07-12 12:19:17 +02:00
|
|
|
|
map<string, string> tmp_paren_vars;
|
|
|
|
|
bool message_printed {false};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& it : tt_output)
|
2022-07-12 12:19:17 +02:00
|
|
|
|
fixNestedParenthesis(it, tmp_paren_vars, message_printed);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (auto& it : d_output)
|
2022-07-12 12:19:17 +02:00
|
|
|
|
fixNestedParenthesis(it, tmp_paren_vars, message_printed);
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {move(d_output), move(tt_output)};
|
2022-07-11 17:33:09 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-10-07 17:19:22 +02:00
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeModelCFile(const string& basename, const string& mexext,
|
|
|
|
|
const filesystem::path& matlabroot) const
|
2022-10-07 17:19:22 +02:00
|
|
|
|
{
|
|
|
|
|
ofstream output;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto open_file = [&output](const filesystem::path& p) {
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output.open(p, ios::out | ios::binary);
|
|
|
|
|
if (!output.is_open())
|
|
|
|
|
{
|
2022-10-11 11:00:50 +02:00
|
|
|
|
cerr << "ERROR: Can't open file " << p.string() << " for writing" << endl;
|
2022-10-07 17:19:22 +02:00
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path model_src_dir {filesystem::path {basename} / "model" / "src"};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [d_output, tt_output] = writeModelFileHelper < dynamic
|
|
|
|
|
? ExprNodeOutputType::CDynamicModel
|
|
|
|
|
: ExprNodeOutputType::CStaticModel > ();
|
2022-10-07 17:19:22 +02:00
|
|
|
|
vector<filesystem::path> header_files, object_files;
|
|
|
|
|
|
|
|
|
|
// TODO: when C++20 support is complete, mark the following strings constexpr
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prefix {dynamic ? "dynamic_" : "static_"};
|
|
|
|
|
const string ss_it_argin {dynamic ? ", const double *restrict steady_state, int it_" : ""};
|
|
|
|
|
const string ss_it_argout {dynamic ? ", steady_state, it_" : ""};
|
|
|
|
|
const string nb_row_x_argin {dynamic ? ", int nb_row_x" : ""};
|
|
|
|
|
const string nb_row_x_argout {dynamic ? ", nb_row_x" : ""};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
|
|
|
|
for (size_t i {0}; i < d_output.size(); i++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string funcname {prefix + (i == 0 ? "resid" : "g" + to_string(i))};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prototype_tt {"void " + funcname
|
|
|
|
|
+ "_tt(const double *restrict y, const double *restrict x"
|
|
|
|
|
+ nb_row_x_argin + ", const double *restrict params" + ss_it_argin
|
|
|
|
|
+ ", double *restrict T)"};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path header_tt {model_src_dir / (funcname + "_tt.h")};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
open_file(header_tt);
|
|
|
|
|
output << prototype_tt << ";" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
header_files.push_back(header_tt);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path source_tt {model_src_dir / (funcname + "_tt.c")};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
open_file(source_tt);
|
|
|
|
|
output << "#include <math.h>" << endl
|
|
|
|
|
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
|
|
|
|
|
<< endl;
|
2023-04-11 15:29:12 +02:00
|
|
|
|
writeCHelpersDefinition(output);
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output << endl
|
|
|
|
|
<< prototype_tt << endl
|
|
|
|
|
<< "{" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_output[i].str() << "}" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< endl;
|
|
|
|
|
output.close();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
object_files.push_back(
|
|
|
|
|
compileMEX(model_src_dir, funcname + "_tt", mexext, {source_tt}, matlabroot, false));
|
|
|
|
|
|
|
|
|
|
const string prototype_main {[&funcname, &ss_it_argin, &nb_row_x_argin, i] {
|
|
|
|
|
string p = "void " + funcname + "(const double *restrict y, const double *restrict x"
|
|
|
|
|
+ nb_row_x_argin + ", const double *restrict params" + ss_it_argin
|
|
|
|
|
+ ", const double *restrict T, ";
|
|
|
|
|
if (i == 0)
|
|
|
|
|
p += "double *restrict residual";
|
|
|
|
|
else if (i == 1)
|
|
|
|
|
p += "double *restrict g1";
|
|
|
|
|
else
|
|
|
|
|
p += "double *restrict g" + to_string(i) + "_i, double *restrict g" + to_string(i)
|
|
|
|
|
+ "_j, double *restrict g" + to_string(i) + "_v";
|
|
|
|
|
p += ")";
|
|
|
|
|
return p;
|
|
|
|
|
}()};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path header_main {model_src_dir / (funcname + ".h")};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
open_file(header_main);
|
|
|
|
|
output << prototype_main << ";" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
header_files.push_back(header_main);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path source_main {model_src_dir / (funcname + ".c")};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
open_file(source_main);
|
|
|
|
|
output << "#include <math.h>" << endl
|
|
|
|
|
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
|
|
|
|
|
<< endl;
|
2023-04-11 15:29:12 +02:00
|
|
|
|
writeCHelpersDefinition(output);
|
|
|
|
|
if (i == 0)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
writeCHelpersDeclaration(
|
|
|
|
|
output); // Provide external definition of helpers in resid main file
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output << endl
|
|
|
|
|
<< prototype_main << endl
|
2023-01-04 14:43:57 +01:00
|
|
|
|
<< "{" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< d_output[i].str() << "}" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< endl;
|
|
|
|
|
output.close();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
object_files.push_back(
|
|
|
|
|
compileMEX(model_src_dir, funcname, mexext, {source_main}, matlabroot, false));
|
2022-10-07 17:19:22 +02:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path filename {model_src_dir / (dynamic ? "dynamic.c" : "static.c")};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const int ntt {static_cast<int>(
|
|
|
|
|
temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()
|
|
|
|
|
+ temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size())};
|
2022-10-07 17:19:22 +02:00
|
|
|
|
|
|
|
|
|
open_file(filename);
|
|
|
|
|
output << "/*" << endl
|
|
|
|
|
<< " * " << filename << " : Computes " << modelClassName() << " for Dynare" << endl
|
|
|
|
|
<< " *" << endl
|
|
|
|
|
<< " * Warning : this file is generated automatically by Dynare" << endl
|
|
|
|
|
<< " * from model file (.mod)" << endl
|
|
|
|
|
<< " */" << endl
|
|
|
|
|
<< endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "#include <math.h>" << endl // Needed for getPowerDeriv()
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< "#include <stdlib.h>" << endl // Needed for malloc() and free()
|
|
|
|
|
<< R"(#include "mex.h")" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : header_files)
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output << "#include " << it.filename() << endl;
|
|
|
|
|
output << endl
|
|
|
|
|
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
|
|
|
|
|
<< "{" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
|
|
|
|
output
|
|
|
|
|
<< " if (nlhs > " << min(computed_derivs_order + 1, 4) << ")" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Derivatives of higher order than computed have been requested");)"
|
|
|
|
|
<< endl
|
|
|
|
|
<< " if (nrhs != 5)" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Requires exactly 5 input arguments");)" << endl;
|
2022-10-07 17:19:22 +02:00
|
|
|
|
else
|
|
|
|
|
output << " if (nrhs > 3)" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Accepts at most 3 output arguments");)" << endl
|
|
|
|
|
<< " if (nrhs != 3)" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Requires exactly 3 input arguments");)" << endl;
|
|
|
|
|
output << endl
|
|
|
|
|
<< " double *y = mxGetPr(prhs[0]);" << endl
|
|
|
|
|
<< " double *x = mxGetPr(prhs[1]);" << endl
|
|
|
|
|
<< " double *params = mxGetPr(prhs[2]);" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output << " double *steady_state = mxGetPr(prhs[3]);" << endl
|
|
|
|
|
<< " int it_ = (int) mxGetScalar(prhs[4]) - 1;" << endl
|
|
|
|
|
<< " int nb_row_x = mxGetM(prhs[1]);" << endl;
|
|
|
|
|
output << endl
|
|
|
|
|
<< " double *T = (double *) malloc(sizeof(double)*" << ntt << ");" << endl
|
|
|
|
|
<< endl
|
|
|
|
|
<< " if (nlhs >= 1)" << endl
|
|
|
|
|
<< " {" << endl
|
|
|
|
|
<< " plhs[0] = mxCreateDoubleMatrix(" << equations.size() << ",1, mxREAL);" << endl
|
|
|
|
|
<< " double *residual = mxGetPr(plhs[0]);" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " " << prefix << "resid_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T);" << endl
|
|
|
|
|
<< " " << prefix << "resid(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T, residual);" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< endl
|
|
|
|
|
<< " if (nlhs >= 2)" << endl
|
|
|
|
|
<< " {" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", "
|
|
|
|
|
<< getJacobianColsNbr(false) << ", mxREAL);" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " double *g1 = mxGetPr(plhs[1]);" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " " << prefix << "g1_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T);" << endl
|
|
|
|
|
<< " " << prefix << "g1(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T, g1);" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< endl
|
|
|
|
|
<< " if (nlhs >= 3)" << endl
|
|
|
|
|
<< " {" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " mxArray *g2_i = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1
|
|
|
|
|
<< ", mxREAL);" << endl
|
|
|
|
|
<< " mxArray *g2_j = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1
|
|
|
|
|
<< ", mxREAL);" << endl
|
|
|
|
|
<< " mxArray *g2_v = mxCreateDoubleMatrix(" << NNZDerivatives[2] << ", " << 1
|
|
|
|
|
<< ", mxREAL);" << endl
|
|
|
|
|
<< " " << prefix << "g2_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T);" << endl
|
|
|
|
|
<< " " << prefix << "g2(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " mxArray *n = mxCreateDoubleScalar("
|
|
|
|
|
<< getJacobianColsNbr(false) * getJacobianColsNbr(false) << ");" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl
|
|
|
|
|
<< R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
|
|
|
|
|
<< " plhs[2] = plhs_sparse[0];" << endl
|
|
|
|
|
<< " mxDestroyArray(g2_i);" << endl
|
|
|
|
|
<< " mxDestroyArray(g2_j);" << endl
|
|
|
|
|
<< " mxDestroyArray(g2_v);" << endl
|
|
|
|
|
<< " mxDestroyArray(m);" << endl
|
|
|
|
|
<< " mxDestroyArray(n);" << endl
|
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output << " if (nlhs >= 4)" << endl
|
|
|
|
|
<< " {" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " mxArray *g3_i = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1
|
|
|
|
|
<< ", mxREAL);" << endl
|
|
|
|
|
<< " mxArray *g3_j = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1
|
|
|
|
|
<< ", mxREAL);" << endl
|
|
|
|
|
<< " mxArray *g3_v = mxCreateDoubleMatrix(" << NNZDerivatives[3] << ", " << 1
|
|
|
|
|
<< ", mxREAL);" << endl
|
|
|
|
|
<< " " << prefix << "g3_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T);" << endl
|
|
|
|
|
<< " " << prefix << "g3(y, x" << nb_row_x_argout << ", params" << ss_it_argout
|
|
|
|
|
<< ", T, mxGetPr(g3_i), mxGetPr(g3_j), mxGetPr(g3_v));" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " mxArray *n = mxCreateDoubleScalar("
|
|
|
|
|
<< getJacobianColsNbr(false) * getJacobianColsNbr(false) * getJacobianColsNbr(false)
|
|
|
|
|
<< ");" << endl
|
2022-10-07 17:19:22 +02:00
|
|
|
|
<< " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g3_i, g3_j, g3_v, m, n };" << endl
|
|
|
|
|
<< R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
|
|
|
|
|
<< " plhs[3] = plhs_sparse[0];" << endl
|
|
|
|
|
<< " mxDestroyArray(g3_i);" << endl
|
|
|
|
|
<< " mxDestroyArray(g3_j);" << endl
|
|
|
|
|
<< " mxDestroyArray(g3_v);" << endl
|
|
|
|
|
<< " mxDestroyArray(m);" << endl
|
|
|
|
|
<< " mxDestroyArray(n);" << endl
|
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< endl;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " free(T);" << endl << "}" << endl;
|
2022-10-07 17:19:22 +02:00
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
|
|
object_files.push_back(filename);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
compileMEX(packageDir(basename), dynamic ? "dynamic" : "static", mexext, object_files,
|
|
|
|
|
matlabroot);
|
2022-10-07 17:19:22 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-07-21 18:20:35 +02:00
|
|
|
|
template<ExprNodeOutputType output_type>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writePerBlockHelper(int blk, ostream& output, temporary_terms_t& temporary_terms) const
|
2022-07-21 18:20:35 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int block_recursive_size {blocks[blk].getRecursiveSize()};
|
2022-07-21 18:20:35 +02:00
|
|
|
|
|
|
|
|
|
// The equations
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto write_eq_tt = [&](int eq) {
|
|
|
|
|
for (auto it : blocks_temporary_terms[blk][eq])
|
|
|
|
|
{
|
|
|
|
|
if (dynamic_cast<AbstractExternalFunctionNode*>(it))
|
|
|
|
|
it->writeExternalFunctionOutput(output, output_type, temporary_terms,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
|
|
|
|
|
output << " ";
|
|
|
|
|
it->writeOutput(output, output_type, blocks_temporary_terms[blk][eq],
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
output << '=';
|
|
|
|
|
it->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
|
|
|
|
temporary_terms.insert(it);
|
|
|
|
|
output << ';' << endl;
|
|
|
|
|
}
|
|
|
|
|
};
|
2022-07-21 18:20:35 +02:00
|
|
|
|
|
|
|
|
|
for (int eq {0}; eq < blocks[blk].size; eq++)
|
|
|
|
|
{
|
|
|
|
|
write_eq_tt(eq);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
EquationType equ_type {getBlockEquationType(blk, eq)};
|
|
|
|
|
BinaryOpNode* e {getBlockEquationExpr(blk, eq)};
|
|
|
|
|
expr_t lhs {e->arg1}, rhs {e->arg2};
|
2022-07-21 18:20:35 +02:00
|
|
|
|
switch (blocks[blk].simulation_type)
|
|
|
|
|
{
|
|
|
|
|
case BlockSimulationType::evaluateBackward:
|
|
|
|
|
case BlockSimulationType::evaluateForward:
|
2023-11-30 15:28:57 +01:00
|
|
|
|
evaluation:
|
2022-07-21 18:20:35 +02:00
|
|
|
|
if (equ_type == EquationType::evaluateRenormalized)
|
|
|
|
|
{
|
|
|
|
|
e = getBlockEquationRenormalizedExpr(blk, eq);
|
|
|
|
|
lhs = e->arg1;
|
|
|
|
|
rhs = e->arg2;
|
|
|
|
|
}
|
2023-11-03 17:56:06 +01:00
|
|
|
|
else
|
|
|
|
|
assert(equ_type == EquationType::evaluate);
|
2022-07-21 18:20:35 +02:00
|
|
|
|
output << " ";
|
|
|
|
|
lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
|
|
|
|
output << '=';
|
|
|
|
|
rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
|
|
|
|
output << ';' << endl;
|
|
|
|
|
break;
|
|
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesSimple:
|
|
|
|
|
if (eq < block_recursive_size)
|
|
|
|
|
goto evaluation;
|
2023-11-03 17:56:06 +01:00
|
|
|
|
[[fallthrough]];
|
|
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
2022-07-21 18:20:35 +02:00
|
|
|
|
output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< eq - block_recursive_size + ARRAY_SUBSCRIPT_OFFSET(output_type)
|
2022-07-21 18:20:35 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
|
|
|
|
|
lhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
|
|
|
|
output << ")-(";
|
|
|
|
|
rhs->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
|
|
|
|
output << ");" << endl;
|
|
|
|
|
break;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
/* Write temporary terms for derivatives.
|
|
|
|
|
This is done even for “evaluate” blocks, whose derivatives are not
|
|
|
|
|
always computed at runtime; still those temporary terms may be needed by
|
|
|
|
|
subsequent blocks. */
|
|
|
|
|
write_eq_tt(blocks[blk].size);
|
|
|
|
|
}
|
|
|
|
|
|
2022-07-12 12:27:22 +02:00
|
|
|
|
template<ExprNodeOutputType output_type>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
tuple<ostringstream, ostringstream, ostringstream, ostringstream, ostringstream, ostringstream,
|
|
|
|
|
ostringstream>
|
2022-07-12 12:27:22 +02:00
|
|
|
|
ModelTree::writeParamsDerivativesFileHelper() const
|
|
|
|
|
{
|
|
|
|
|
static_assert(!isCOutput(output_type), "C output is not implemented");
|
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
constexpr bool sparse {isSparseModelOutput(output_type)};
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ostringstream tt_output; // Used for storing model temp vars and equations
|
|
|
|
|
ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters
|
|
|
|
|
ostringstream gp_output; // 1st deriv. of Jacobian w.r.t. parameters
|
2022-07-12 12:27:22 +02:00
|
|
|
|
ostringstream rpp_output; // 2nd deriv of residuals w.r.t. parameters
|
|
|
|
|
ostringstream gpp_output; // 2nd deriv of Jacobian w.r.t. parameters
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ostringstream hp_output; // 1st deriv. of Hessian w.r.t. parameters
|
|
|
|
|
ostringstream
|
|
|
|
|
g3p_output; // 1st deriv. of 3rd deriv. matrix w.r.t. parameters (only in dynamic case)
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
|
|
|
|
temporary_terms_t temp_term_union;
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [order, tts] : params_derivs_temporary_terms)
|
2022-07-12 14:13:27 +02:00
|
|
|
|
writeTemporaryTerms<output_type>(tts, temp_term_union, params_derivs_temporary_terms_idxs,
|
|
|
|
|
tt_output, tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : params_derivatives.at({0, 1}))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, param] {vectorToTuple<2>(indices)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
rp_output << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq + 1 << ", " << param_col
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d1->writeOutput(rp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
rp_output << ";" << endl;
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d2] : params_derivatives.at({1, 1}))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var, param] {vectorToTuple<3>(indices)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var_col {getJacobianCol(var, sparse) + 1};
|
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
gp_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq + 1 << ", " << var_col << ", "
|
|
|
|
|
<< param_col << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
|
|
|
|
|
d2->writeOutput(gp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
gp_output << ";" << endl;
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int i {1}; const auto& [indices, d2] : params_derivatives.at({0, 2}))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, param1, param2] {vectorToTuple<3>(indices)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int param1_col {getTypeSpecificIDByDerivID(param1) + 1};
|
|
|
|
|
int param2_col {getTypeSpecificIDByDerivID(param2) + 1};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
|
|
|
|
rpp_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
|
|
|
|
|
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
|
|
|
|
|
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d2->writeOutput(rpp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
rpp_output << ";" << endl;
|
|
|
|
|
|
|
|
|
|
i++;
|
|
|
|
|
|
|
|
|
|
if (param1 != param2)
|
|
|
|
|
{
|
|
|
|
|
// Treat symmetric elements
|
|
|
|
|
rpp_output << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
|
|
|
|
|
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
|
|
|
|
|
<< "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=rpp"
|
|
|
|
|
<< LEFT_ARRAY_SUBSCRIPT(output_type) << i - 1 << ",4"
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
|
|
|
|
|
i++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int i {1}; const auto& [indices, d2] : params_derivatives.at({1, 2}))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var, param1, param2] {vectorToTuple<4>(indices)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var_col {getJacobianCol(var, sparse) + 1};
|
|
|
|
|
int param1_col {getTypeSpecificIDByDerivID(param1) + 1};
|
|
|
|
|
int param2_col {getTypeSpecificIDByDerivID(param2) + 1};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
|
|
|
|
gpp_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
|
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
|
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
|
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d2->writeOutput(gpp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
gpp_output << ";" << endl;
|
|
|
|
|
|
|
|
|
|
i++;
|
|
|
|
|
|
|
|
|
|
if (param1 != param2)
|
|
|
|
|
{
|
|
|
|
|
// Treat symmetric elements
|
|
|
|
|
gpp_output << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
|
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
|
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
|
|
|
|
|
<< "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=gpp"
|
|
|
|
|
<< LEFT_ARRAY_SUBSCRIPT(output_type) << i - 1 << ",5"
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
|
|
|
|
|
i++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int i {1}; const auto& [indices, d2] : params_derivatives.at({2, 1}))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var1, var2, param] {vectorToTuple<4>(indices)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var1_col {getJacobianCol(var1, sparse) + 1};
|
|
|
|
|
int var2_col {getJacobianCol(var2, sparse) + 1};
|
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
|
|
|
|
hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
|
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
|
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
|
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d2->writeOutput(hp_output, output_type, temp_term_union, params_derivs_temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
hp_output << ";" << endl;
|
|
|
|
|
|
|
|
|
|
i++;
|
|
|
|
|
|
|
|
|
|
if (var1 != var2)
|
|
|
|
|
{
|
|
|
|
|
// Treat symmetric elements
|
|
|
|
|
hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
|
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
|
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
|
|
|
|
|
<< "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=hp"
|
|
|
|
|
<< LEFT_ARRAY_SUBSCRIPT(output_type) << i - 1 << ",5"
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
|
|
|
|
|
i++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (output_type == ExprNodeOutputType::matlabDynamicModel
|
|
|
|
|
|| output_type == ExprNodeOutputType::juliaDynamicModel)
|
|
|
|
|
for (int i {1}; const auto& [indices, d2] : params_derivatives.at({3, 1}))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var1, var2, var3, param] {vectorToTuple<5>(indices)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var1_col {getJacobianCol(var1, sparse) + 1};
|
|
|
|
|
int var2_col {getJacobianCol(var2, sparse) + 1};
|
|
|
|
|
int var3_col {getJacobianCol(var3, sparse) + 1};
|
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
|
|
|
|
|
g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq + 1 << ";" << endl
|
2022-07-12 12:27:22 +02:00
|
|
|
|
<< "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
|
|
|
|
|
<< "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
|
|
|
|
|
<< "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var3_col << ";" << endl
|
|
|
|
|
<< "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
|
|
|
|
|
<< "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",6"
|
|
|
|
|
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d2->writeOutput(g3p_output, output_type, temp_term_union,
|
|
|
|
|
params_derivs_temporary_terms_idxs, tef_terms);
|
2022-07-12 12:27:22 +02:00
|
|
|
|
g3p_output << ";" << endl;
|
|
|
|
|
|
|
|
|
|
i++;
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (isMatlabOutput(output_type))
|
2022-07-12 12:27:22 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// Check that we don't have more than 32 nested parenthesis because MATLAB does not support
|
|
|
|
|
// this. See Issue #1201
|
2022-07-12 12:27:22 +02:00
|
|
|
|
map<string, string> tmp_paren_vars;
|
|
|
|
|
bool message_printed {false};
|
|
|
|
|
fixNestedParenthesis(tt_output, tmp_paren_vars, message_printed);
|
|
|
|
|
fixNestedParenthesis(rp_output, tmp_paren_vars, message_printed);
|
|
|
|
|
fixNestedParenthesis(gp_output, tmp_paren_vars, message_printed);
|
|
|
|
|
fixNestedParenthesis(rpp_output, tmp_paren_vars, message_printed);
|
|
|
|
|
fixNestedParenthesis(gpp_output, tmp_paren_vars, message_printed);
|
|
|
|
|
fixNestedParenthesis(hp_output, tmp_paren_vars, message_printed);
|
|
|
|
|
fixNestedParenthesis(g3p_output, tmp_paren_vars, message_printed);
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {move(tt_output), move(rp_output), move(gp_output), move(rpp_output),
|
|
|
|
|
move(gpp_output), move(hp_output), move(g3p_output)};
|
2022-07-12 12:27:22 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-07-13 13:04:10 +02:00
|
|
|
|
template<ExprNodeBytecodeOutputType output_type>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeBytecodeTemporaryTerms(const temporary_terms_t& tt,
|
|
|
|
|
temporary_terms_t& temporary_terms_union,
|
2023-12-14 14:52:50 +01:00
|
|
|
|
Bytecode::Writer& code_file,
|
2023-11-30 15:28:57 +01:00
|
|
|
|
deriv_node_temp_terms_t& tef_terms) const
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
|
|
|
|
for (auto it : tt)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (dynamic_cast<AbstractExternalFunctionNode*>(it))
|
|
|
|
|
it->writeBytecodeExternalFunctionOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
temporary_terms_idxs, tef_terms);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
|
|
|
|
int idx {temporary_terms_idxs.at(it)};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::TemporaryTerm, idx};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
it->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-14 09:17:11 +02:00
|
|
|
|
|
|
|
|
|
static_assert(output_type == ExprNodeBytecodeOutputType::dynamicModel
|
|
|
|
|
|| output_type == ExprNodeBytecodeOutputType::staticModel);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (output_type == ExprNodeBytecodeOutputType::dynamicModel)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPT {idx};
|
2022-07-14 09:17:11 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPST {idx};
|
2022-07-14 09:17:11 +02:00
|
|
|
|
|
|
|
|
|
temporary_terms_union.insert(it);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<ExprNodeBytecodeOutputType output_type>
|
|
|
|
|
void
|
2023-12-14 14:52:50 +01:00
|
|
|
|
ModelTree::writeBytecodeModelEquations(Bytecode::Writer& code_file,
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const temporary_terms_t& temporary_terms,
|
|
|
|
|
const deriv_node_temp_terms_t& tef_terms) const
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
|
|
|
|
for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
BinaryOpNode* eq_node {equations[eq]};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
expr_t lhs {eq_node->arg1}, rhs {eq_node->arg2};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::ModelEquation, eq};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
// Test if the right hand side of the equation is empty.
|
|
|
|
|
double vrhs {1.0};
|
|
|
|
|
try
|
|
|
|
|
{
|
|
|
|
|
vrhs = rhs->eval({});
|
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
catch (ExprNode::EvalException& e)
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (vrhs != 0) // The right hand side of the equation is not empty ⇒ residual=lhs-rhs
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
lhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
|
|
|
|
rhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::minus} << Bytecode::FSTPR {eq};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
else // The right hand side of the equation is empty ⇒ residual=lhs
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
lhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPR {eq};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-12-14 14:52:50 +01:00
|
|
|
|
ModelTree::writeBytecodeHelper(Bytecode::Writer& code_file) const
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
constexpr ExprNodeBytecodeOutputType output_type {
|
|
|
|
|
dynamic ? ExprNodeBytecodeOutputType::dynamicModel : ExprNodeBytecodeOutputType::staticModel};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
|
|
|
|
temporary_terms_t temporary_terms_union;
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
writeBytecodeTemporaryTerms<output_type>(temporary_terms_derivatives[0], temporary_terms_union,
|
|
|
|
|
code_file, tef_terms);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
writeBytecodeModelEquations<output_type>(code_file, temporary_terms_union, tef_terms);
|
|
|
|
|
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FENDEQU {};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
|
|
|
|
// Temporary terms for the Jacobian
|
2023-11-30 15:28:57 +01:00
|
|
|
|
writeBytecodeTemporaryTerms<output_type>(temporary_terms_derivatives[1], temporary_terms_union,
|
|
|
|
|
code_file, tef_terms);
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
|
|
|
|
// Get the current code_file position and jump if “evaluate” mode
|
|
|
|
|
int pos_jmpifeval {code_file.getInstructionCounter()};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FJMPIFEVAL {0}; // Use 0 as jump offset for the time being
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
|
|
|
|
// The Jacobian in “simulate” mode
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<vector<tuple<int, int, int>>> my_derivatives(symbol_table.endo_nbr());
|
|
|
|
|
;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
int count_u {symbol_table.endo_nbr()};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : derivatives[1])
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
|
|
|
|
auto [eq, deriv_id] {vectorToTuple<2>(indices)};
|
|
|
|
|
if (getTypeByDerivID(deriv_id) == SymbolType::endogenous)
|
|
|
|
|
{
|
|
|
|
|
int tsid {getTypeSpecificIDByDerivID(deriv_id)};
|
|
|
|
|
int lag {getLagByDerivID(deriv_id)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::FirstEndoDerivative, eq,
|
|
|
|
|
tsid, lag};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::FirstEndoDerivative, eq,
|
|
|
|
|
tsid};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
if (!my_derivatives[eq].size())
|
|
|
|
|
my_derivatives[eq].clear();
|
|
|
|
|
my_derivatives[eq].emplace_back(tsid, lag, count_u);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d1->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
temporary_terms_idxs, tef_terms);
|
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPU {count_u};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPSU {count_u};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
count_u++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
for (int i {0}; i < symbol_table.endo_nbr(); i++)
|
|
|
|
|
{
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDR {i};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
if (my_derivatives[i].size())
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool first_term {true}; const auto& [tsid, lag, uidx] : my_derivatives[i])
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDU {uidx}
|
|
|
|
|
<< Bytecode::FLDV {SymbolType::endogenous, tsid, lag};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDSU {uidx}
|
|
|
|
|
<< Bytecode::FLDSV {SymbolType::endogenous, tsid};
|
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::times};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
if (!exchange(first_term, false))
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::plus};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::minus};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPU {i};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPSU {i};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// Jump unconditionally after the block
|
|
|
|
|
int pos_jmp {code_file.getInstructionCounter()};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FJMP {0}; // Use 0 as jump offset for the time being
|
2022-07-13 13:04:10 +02:00
|
|
|
|
// Update jump offset for previous JMPIFEVAL
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file.overwriteInstruction(pos_jmpifeval, Bytecode::FJMPIFEVAL {pos_jmp - pos_jmpifeval});
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
|
|
|
|
// The Jacobian in “evaluate” mode
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : derivatives[1])
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
|
|
|
|
auto [eq, deriv_id] {vectorToTuple<2>(indices)};
|
|
|
|
|
int tsid {getTypeSpecificIDByDerivID(deriv_id)};
|
|
|
|
|
int lag {getLagByDerivID(deriv_id)};
|
|
|
|
|
SymbolType type {getTypeByDerivID(deriv_id)};
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-13 13:04:10 +02:00
|
|
|
|
{
|
2023-12-14 14:52:50 +01:00
|
|
|
|
Bytecode::ExpressionType expr_type;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
switch (type)
|
|
|
|
|
{
|
|
|
|
|
case SymbolType::endogenous:
|
2023-12-14 14:52:50 +01:00
|
|
|
|
expr_type = Bytecode::ExpressionType::FirstEndoDerivative;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
break;
|
|
|
|
|
case SymbolType::exogenous:
|
2023-12-14 14:52:50 +01:00
|
|
|
|
expr_type = Bytecode::ExpressionType::FirstExoDerivative;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
break;
|
|
|
|
|
case SymbolType::exogenousDet:
|
2023-12-14 14:52:50 +01:00
|
|
|
|
expr_type = Bytecode::ExpressionType::FirstExodetDerivative;
|
2022-07-13 13:04:10 +02:00
|
|
|
|
break;
|
|
|
|
|
default:
|
|
|
|
|
assert(false);
|
|
|
|
|
break;
|
|
|
|
|
}
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {expr_type, eq, tsid, lag};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
assert(type == SymbolType::endogenous);
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::FirstEndoDerivative, eq, tsid};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d1->writeBytecodeOutput(code_file, output_type, temporary_terms_union, temporary_terms_idxs,
|
|
|
|
|
tef_terms);
|
|
|
|
|
if constexpr (dynamic)
|
|
|
|
|
{
|
|
|
|
|
// Bytecode MEX uses a separate matrix for exogenous and exodet Jacobians
|
|
|
|
|
int jacob_col {type == SymbolType::endogenous ? getJacobianCol(deriv_id, false) : tsid};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPG3 {eq, tsid, lag, jacob_col};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
}
|
2022-07-13 13:04:10 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPG2 {eq, tsid};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// Update jump offset for previous JMP
|
|
|
|
|
int pos_end_block {code_file.getInstructionCounter()};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file.overwriteInstruction(pos_jmp, Bytecode::FJMP {pos_end_block - pos_jmp - 1});
|
2022-07-13 13:04:10 +02:00
|
|
|
|
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FENDBLOCK {} << Bytecode::FEND {};
|
2022-07-13 13:04:10 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-07-19 18:24:36 +02:00
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-12-14 14:52:50 +01:00
|
|
|
|
ModelTree::writeBlockBytecodeHelper(Bytecode::Writer& code_file, int block,
|
2023-11-30 15:28:57 +01:00
|
|
|
|
temporary_terms_t& temporary_terms_union) const
|
2022-07-19 18:24:36 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
constexpr ExprNodeBytecodeOutputType output_type {
|
|
|
|
|
dynamic ? ExprNodeBytecodeOutputType::dynamicModel : ExprNodeBytecodeOutputType::staticModel};
|
|
|
|
|
constexpr ExprNodeBytecodeOutputType assignment_lhs_output_type {
|
|
|
|
|
dynamic ? ExprNodeBytecodeOutputType::dynamicAssignmentLHS
|
|
|
|
|
: ExprNodeBytecodeOutputType::staticAssignmentLHS};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
|
|
|
|
const BlockSimulationType simulation_type {blocks[block].simulation_type};
|
|
|
|
|
const int block_size {blocks[block].size};
|
|
|
|
|
const int block_mfs {blocks[block].mfs_size};
|
|
|
|
|
const int block_recursive {blocks[block].getRecursiveSize()};
|
|
|
|
|
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto write_eq_tt = [&](int eq) {
|
2022-07-19 18:24:36 +02:00
|
|
|
|
for (auto it : blocks_temporary_terms[block][eq])
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (dynamic_cast<AbstractExternalFunctionNode*>(it))
|
|
|
|
|
it->writeBytecodeExternalFunctionOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::TemporaryTerm,
|
|
|
|
|
blocks_temporary_terms_idxs.at(it)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
it->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPT {blocks_temporary_terms_idxs.at(it)};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPST {blocks_temporary_terms_idxs.at(it)};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
temporary_terms_union.insert(it);
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
// The equations
|
|
|
|
|
for (int i {0}; i < block_size; i++)
|
|
|
|
|
{
|
|
|
|
|
write_eq_tt(i);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
EquationType equ_type {getBlockEquationType(block, i)};
|
|
|
|
|
BinaryOpNode* e {getBlockEquationExpr(block, i)};
|
|
|
|
|
expr_t lhs {e->arg1}, rhs {e->arg2};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
switch (simulation_type)
|
|
|
|
|
{
|
|
|
|
|
case BlockSimulationType::evaluateBackward:
|
|
|
|
|
case BlockSimulationType::evaluateForward:
|
2023-11-30 15:28:57 +01:00
|
|
|
|
evaluation:
|
2023-11-03 17:56:06 +01:00
|
|
|
|
if (equ_type == EquationType::evaluateRenormalized)
|
2022-07-19 18:24:36 +02:00
|
|
|
|
{
|
2023-11-03 17:56:06 +01:00
|
|
|
|
e = getBlockEquationRenormalizedExpr(block, i);
|
|
|
|
|
lhs = e->arg1;
|
|
|
|
|
rhs = e->arg2;
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
2023-11-03 17:56:06 +01:00
|
|
|
|
else
|
|
|
|
|
assert(equ_type == EquationType::evaluate);
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::ModelEquation,
|
|
|
|
|
getBlockEquationID(block, i)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
lhs->writeBytecodeOutput(code_file, assignment_lhs_output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
break;
|
|
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesSimple:
|
|
|
|
|
if (i < block_recursive)
|
|
|
|
|
goto evaluation;
|
|
|
|
|
[[fallthrough]];
|
2023-11-03 17:56:06 +01:00
|
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::ModelEquation,
|
|
|
|
|
getBlockEquationID(block, i)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
lhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::minus}
|
|
|
|
|
<< Bytecode::FSTPR {i - block_recursive};
|
2023-11-03 17:56:06 +01:00
|
|
|
|
break;
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2022-07-21 17:22:08 +02:00
|
|
|
|
/* Write temporary terms for derivatives. This is done before FENDEQU,
|
|
|
|
|
because residuals of a subsequent block may depend on temporary terms for
|
|
|
|
|
the derivatives of the present block.
|
|
|
|
|
|
|
|
|
|
Also note that in the case of “evaluate” blocks, derivatives are not
|
|
|
|
|
computed in the “evaluate” mode; still their temporary terms must be
|
|
|
|
|
computed even in that mode, because for the same reason as above they may
|
|
|
|
|
be needed in subsequent blocks. */
|
|
|
|
|
write_eq_tt(blocks[block].size);
|
|
|
|
|
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FENDEQU {};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
|
|
|
|
// Get the current code_file position and jump if evaluating
|
|
|
|
|
int pos_jmpifeval {code_file.getInstructionCounter()};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FJMPIFEVAL {0}; // Use 0 as jump offset for the time being
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
|
|
|
|
/* Write the derivatives for the “simulate” mode (not needed if the block
|
|
|
|
|
is of type “evaluate backward/forward”) */
|
|
|
|
|
if (simulation_type != BlockSimulationType::evaluateBackward
|
|
|
|
|
&& simulation_type != BlockSimulationType::evaluateForward)
|
|
|
|
|
{
|
|
|
|
|
switch (simulation_type)
|
|
|
|
|
{
|
|
|
|
|
case BlockSimulationType::solveBackwardSimple:
|
|
|
|
|
case BlockSimulationType::solveForwardSimple:
|
|
|
|
|
{
|
|
|
|
|
int eqr {getBlockEquationID(block, 0)};
|
|
|
|
|
int varr {getBlockVariableID(block, 0)};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::FirstEndoDerivative, eqr,
|
|
|
|
|
varr, 0};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
// Get contemporaneous derivative of the single variable in the block
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if (auto it {blocks_derivatives[block].find({0, 0, 0})};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
it != blocks_derivatives[block].end())
|
2023-11-30 15:28:57 +01:00
|
|
|
|
it->second->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
2022-07-19 18:24:36 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDZ {};
|
|
|
|
|
code_file << Bytecode::FSTPG {0};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
break;
|
|
|
|
|
|
|
|
|
|
case BlockSimulationType::solveBackwardComplete:
|
|
|
|
|
case BlockSimulationType::solveForwardComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesComplete:
|
|
|
|
|
case BlockSimulationType::solveTwoBoundariesSimple:
|
|
|
|
|
{
|
|
|
|
|
// For each equation, stores a list of tuples (index_u, var, lag)
|
|
|
|
|
vector<vector<tuple<int, int, int>>> Uf(symbol_table.endo_nbr());
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (int count_u {block_mfs}; const auto& [indices, d1] : 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
|
|
|
|
int eqr {getBlockEquationID(block, eq)};
|
|
|
|
|
int varr {getBlockVariableID(block, var)};
|
2023-11-03 11:52:54 +01:00
|
|
|
|
assert(eq >= block_recursive);
|
|
|
|
|
if (var >= block_recursive)
|
2022-07-19 18:24:36 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-19 18:24:36 +02:00
|
|
|
|
if (lag != 0
|
|
|
|
|
&& (simulation_type == BlockSimulationType::solveForwardComplete
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveBackwardComplete))
|
|
|
|
|
continue;
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::FirstEndoDerivative,
|
|
|
|
|
eqr, varr, lag};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d1->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPU {count_u};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPSU {count_u};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
Uf[eqr].emplace_back(count_u, varr, lag);
|
|
|
|
|
count_u++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
for (int i {0}; i < block_size; i++)
|
|
|
|
|
if (i >= block_recursive)
|
|
|
|
|
{
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDR {i - block_recursive} << Bytecode::FLDZ {};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
|
|
|
|
int eqr {getBlockEquationID(block, i)};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [index_u, var, lag] : Uf[eqr])
|
2022-07-19 18:24:36 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDU {index_u}
|
|
|
|
|
<< Bytecode::FLDV {SymbolType::endogenous, var, lag};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FLDSU {index_u}
|
|
|
|
|
<< Bytecode::FLDSV {SymbolType::endogenous, var};
|
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::times}
|
|
|
|
|
<< Bytecode::FBINARY {BinaryOpcode::plus};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FBINARY {BinaryOpcode::minus};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPU {i - block_recursive};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPSU {i - block_recursive};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
break;
|
|
|
|
|
default:
|
|
|
|
|
break;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// Jump unconditionally after the block
|
|
|
|
|
int pos_jmp {code_file.getInstructionCounter()};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FJMP {0}; // Use 0 as jump offset for the time being
|
2022-07-19 18:24:36 +02:00
|
|
|
|
// Update jump offset for previous JMPIFEVAL
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file.overwriteInstruction(pos_jmpifeval, Bytecode::FJMPIFEVAL {pos_jmp - pos_jmpifeval});
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
|
|
|
|
// Write the derivatives for the “evaluate” mode
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d] : 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
|
|
|
|
int eqr {getBlockEquationID(block, eq)};
|
|
|
|
|
int varr {getBlockVariableID(block, var)};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::FirstEndoDerivative, eqr, varr,
|
|
|
|
|
lag};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d->writeBytecodeOutput(code_file, output_type, temporary_terms_union,
|
|
|
|
|
blocks_temporary_terms_idxs, tef_terms);
|
2023-11-03 11:52:54 +01:00
|
|
|
|
assert(eq >= block_recursive);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPG3 {eq - block_recursive, var, lag,
|
|
|
|
|
getBlockJacobianEndoCol(block, var, lag)};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
else
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FSTPG2 {eq - block_recursive,
|
|
|
|
|
getBlockJacobianEndoCol(block, var, lag)};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// Update jump offset for previous JMP
|
|
|
|
|
int pos_end_block {code_file.getInstructionCounter()};
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file.overwriteInstruction(pos_jmp, Bytecode::FJMP {pos_end_block - pos_jmp - 1});
|
2022-07-19 18:24:36 +02:00
|
|
|
|
|
2023-12-14 16:17:22 +01:00
|
|
|
|
code_file << Bytecode::FENDBLOCK {};
|
2022-07-19 18:24:36 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-07-12 17:47:02 +02:00
|
|
|
|
template<bool dynamic>
|
|
|
|
|
pair<ostringstream, vector<ostringstream>>
|
|
|
|
|
ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
|
|
|
|
|
{
|
|
|
|
|
ostringstream mlv_output; // Used for storing model local vars
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<ostringstream> d_output(
|
|
|
|
|
derivatives.size()); // Derivatives output (at all orders, including 0=residual)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
temporary_terms_t temp_term_union;
|
|
|
|
|
|
|
|
|
|
writeJsonModelLocalVariables(mlv_output, true, tef_terms);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms,
|
|
|
|
|
"");
|
2022-07-12 17:47:02 +02:00
|
|
|
|
d_output[0] << ", ";
|
|
|
|
|
writeJsonModelEquations(d_output[0], true);
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int ncols {getJacobianColsNbr(false)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
for (size_t i {1}; i < derivatives.size(); i++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
string matrix_name {i == 1 ? "jacobian"
|
|
|
|
|
: i == 2 ? "hessian"
|
|
|
|
|
: i == 3 ? "third_derivative"
|
|
|
|
|
: to_string(i) + "th_derivative"};
|
|
|
|
|
writeJsonTemporaryTerms(temporary_terms_derivatives[i], temp_term_union, d_output[i],
|
|
|
|
|
tef_terms, matrix_name);
|
|
|
|
|
temp_term_union.insert(temporary_terms_derivatives[i].begin(),
|
|
|
|
|
temporary_terms_derivatives[i].end());
|
|
|
|
|
d_output[i] << R"(, ")" << matrix_name << R"(": {)"
|
|
|
|
|
<< R"( "nrows": )" << equations.size() << R"(, "ncols": )" << ncols
|
2022-07-12 17:47:02 +02:00
|
|
|
|
<< R"(, "entries": [)";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : derivatives[i])
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
d_output[i] << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int eq {vidx[0]};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
int col_idx {0};
|
|
|
|
|
for (size_t j {1}; j < vidx.size(); j++)
|
|
|
|
|
{
|
2022-09-14 17:07:08 +02:00
|
|
|
|
col_idx *= getJacobianColsNbr(false);
|
|
|
|
|
col_idx += getJacobianCol(vidx[j], false);
|
2022-07-12 17:47:02 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
d_output[i] << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
d_output[i] << R"({"row": )" << eq + 1;
|
|
|
|
|
|
|
|
|
|
d_output[i] << R"(, "col": )" << (i > 1 ? "[" : "") << col_idx + 1;
|
|
|
|
|
|
|
|
|
|
if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int col_idx_sym {getJacobianCol(vidx[2], false) * getJacobianColsNbr(false)
|
|
|
|
|
+ getJacobianCol(vidx[1], false)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
d_output[i] << ", " << col_idx_sym + 1;
|
|
|
|
|
}
|
|
|
|
|
if (i > 1)
|
|
|
|
|
d_output[i] << "]";
|
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
for (size_t j = 1; j < vidx.size(); j++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
d_output[i] << R"(, "var)" << (i > 1 ? to_string(j) : "") << R"(": ")"
|
|
|
|
|
<< getNameByDerivID(vidx[j]) << R"(")";
|
|
|
|
|
if constexpr (dynamic)
|
|
|
|
|
d_output[i] << R"(, "shift)" << (i > 1 ? to_string(j) : "") << R"(": )"
|
|
|
|
|
<< getLagByDerivID(vidx[j]);
|
2022-07-12 17:47:02 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
d_output[i] << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(d_output[i], temp_term_union, tef_terms);
|
|
|
|
|
d_output[i] << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
d_output[i] << "]}";
|
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
ncols *= getJacobianColsNbr(false);
|
2022-07-12 17:47:02 +02:00
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {move(mlv_output), move(d_output)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
2023-11-30 15:28:57 +01:00
|
|
|
|
tuple<ostringstream, ostringstream, ostringstream, ostringstream, ostringstream, ostringstream,
|
|
|
|
|
ostringstream, ostringstream>
|
2022-07-12 17:47:02 +02:00
|
|
|
|
ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
|
|
|
|
|
{
|
|
|
|
|
ostringstream mlv_output; // Used for storing model local vars
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ostringstream tt_output; // Used for storing model temp vars and equations
|
|
|
|
|
ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters
|
|
|
|
|
ostringstream gp_output; // 1st deriv. of Jacobian w.r.t. parameters
|
2022-07-12 17:47:02 +02:00
|
|
|
|
ostringstream rpp_output; // 2nd deriv of residuals w.r.t. parameters
|
|
|
|
|
ostringstream gpp_output; // 2nd deriv of Jacobian w.r.t. parameters
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ostringstream hp_output; // 1st deriv. of Hessian w.r.t. parameters
|
2022-07-12 17:47:02 +02:00
|
|
|
|
ostringstream g3p_output; // 1st deriv. of 3rd deriv. matrix w.r.t. parameters
|
|
|
|
|
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
writeJsonModelLocalVariables(mlv_output, true, tef_terms);
|
|
|
|
|
|
|
|
|
|
temporary_terms_t temp_term_union;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [order, tts] : params_derivs_temporary_terms)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
writeJsonTemporaryTerms(tts, temp_term_union, tt_output, tef_terms, "all");
|
|
|
|
|
|
|
|
|
|
rp_output << R"("deriv_wrt_params": {)"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( "neqs": )" << equations.size() << R"(, "nparamcols": )"
|
|
|
|
|
<< symbol_table.param_nbr() << R"(, "entries": [)";
|
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : params_derivatives.at({0, 1}))
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
rp_output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, param] {vectorToTuple<2>(vidx)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
rp_output << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
rp_output << R"({"row": )" << eq + 1;
|
|
|
|
|
|
|
|
|
|
rp_output << R"(, "param_col": )" << param_col;
|
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
rp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
|
|
|
|
|
|
|
|
|
|
rp_output << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(rp_output, temp_term_union, tef_terms);
|
|
|
|
|
rp_output << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
rp_output << "]}";
|
|
|
|
|
|
|
|
|
|
gp_output << R"("deriv_jacobian_wrt_params": {)"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( "neqs": )" << equations.size() << R"(, "nvarcols": )"
|
|
|
|
|
<< getJacobianColsNbr(false) << R"(, "nparamcols": )" << symbol_table.param_nbr()
|
2022-07-12 17:47:02 +02:00
|
|
|
|
<< R"(, "entries": [)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : params_derivatives.at({1, 1}))
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
gp_output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var, param] {vectorToTuple<3>(vidx)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var_col {getJacobianCol(var, false) + 1};
|
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
gp_output << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
gp_output << R"({"row": )" << eq + 1;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
gp_output << R"(, "var_col": )" << var_col << R"(, "param_col": )" << param_col;
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
{
|
|
|
|
|
gp_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
gp_output << R"(, "lag": )" << getLagByDerivID(var);
|
|
|
|
|
gp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
gp_output << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(gp_output, temp_term_union, tef_terms);
|
|
|
|
|
gp_output << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
gp_output << "]}";
|
|
|
|
|
|
|
|
|
|
rpp_output << R"("second_deriv_residuals_wrt_params": {)"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( "nrows": )" << equations.size() << R"(, "nparam1cols": )"
|
|
|
|
|
<< symbol_table.param_nbr() << R"(, "nparam2cols": )" << symbol_table.param_nbr()
|
2022-07-12 17:47:02 +02:00
|
|
|
|
<< R"(, "entries": [)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : params_derivatives.at({0, 2}))
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
rpp_output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, param1, param2] {vectorToTuple<3>(vidx)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int param1_col {getTypeSpecificIDByDerivID(param1) + 1};
|
|
|
|
|
int param2_col {getTypeSpecificIDByDerivID(param2) + 1};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
rpp_output << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
rpp_output << R"({"row": )" << eq + 1;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
rpp_output << R"(, "param1_col": )" << param1_col << R"(, "param2_col": )" << param2_col;
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
rpp_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
|
|
|
|
|
<< R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
|
|
|
|
|
|
|
|
|
|
rpp_output << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(rpp_output, temp_term_union, tef_terms);
|
|
|
|
|
rpp_output << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
rpp_output << "]}";
|
|
|
|
|
|
|
|
|
|
gpp_output << R"("second_deriv_jacobian_wrt_params": {)"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( "neqs": )" << equations.size() << R"(, "nvarcols": )"
|
|
|
|
|
<< getJacobianColsNbr(false) << R"(, "nparam1cols": )" << symbol_table.param_nbr()
|
|
|
|
|
<< R"(, "nparam2cols": )" << symbol_table.param_nbr() << R"(, "entries": [)";
|
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : params_derivatives.at({1, 2}))
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
gpp_output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var, param1, param2] {vectorToTuple<4>(vidx)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var_col {getJacobianCol(var, false) + 1};
|
|
|
|
|
int param1_col {getTypeSpecificIDByDerivID(param1) + 1};
|
|
|
|
|
int param2_col {getTypeSpecificIDByDerivID(param2) + 1};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
gpp_output << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
gpp_output << R"({"row": )" << eq + 1;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
gpp_output << R"(, "var_col": )" << var_col << R"(, "param1_col": )" << param1_col
|
2022-07-12 17:47:02 +02:00
|
|
|
|
<< R"(, "param2_col": )" << param2_col;
|
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
{
|
|
|
|
|
gpp_output << R"(, "var": ")" << getNameByDerivID(var) << R"(")";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
gpp_output << R"(, "lag": )" << getLagByDerivID(var);
|
|
|
|
|
gpp_output << R"(, "param1": ")" << getNameByDerivID(param1) << R"(")"
|
|
|
|
|
<< R"(, "param2": ")" << getNameByDerivID(param2) << R"(")";
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
gpp_output << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(gpp_output, temp_term_union, tef_terms);
|
|
|
|
|
gpp_output << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
gpp_output << "]}" << endl;
|
|
|
|
|
|
|
|
|
|
hp_output << R"("derivative_hessian_wrt_params": {)"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( "neqs": )" << equations.size() << R"(, "nvar1cols": )"
|
|
|
|
|
<< getJacobianColsNbr(false) << R"(, "nvar2cols": )" << getJacobianColsNbr(false)
|
|
|
|
|
<< R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)";
|
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : params_derivatives.at({2, 1}))
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
hp_output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var1, var2, param] {vectorToTuple<4>(vidx)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var1_col {getJacobianCol(var1, false) + 1};
|
|
|
|
|
int var2_col {getJacobianCol(var2, false) + 1};
|
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
hp_output << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
hp_output << R"({"row": )" << eq + 1;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
hp_output << R"(, "var1_col": )" << var1_col << R"(, "var2_col": )" << var2_col
|
2022-07-12 17:47:02 +02:00
|
|
|
|
<< R"(, "param_col": )" << param_col;
|
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
{
|
|
|
|
|
hp_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
hp_output << R"(, "lag1": )" << getLagByDerivID(var1);
|
|
|
|
|
hp_output << R"(, "var2": ")" << getNameByDerivID(var2) << R"(")";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
hp_output << R"(, "lag2": )" << getLagByDerivID(var2);
|
|
|
|
|
hp_output << R"(, "param": ")" << getNameByDerivID(param) << R"(")";
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
hp_output << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(hp_output, temp_term_union, tef_terms);
|
|
|
|
|
hp_output << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
hp_output << "]}" << endl;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
g3p_output << R"("derivative_g3_wrt_params": {)"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( "neqs": )" << equations.size() << R"(, "nvar1cols": )"
|
|
|
|
|
<< getJacobianColsNbr(false) << R"(, "nvar2cols": )" << getJacobianColsNbr(false)
|
|
|
|
|
<< R"(, "nvar3cols": )" << getJacobianColsNbr(false) << R"(, "nparamcols": )"
|
|
|
|
|
<< symbol_table.param_nbr() << R"(, "entries": [)";
|
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : params_derivatives.at({3, 1}))
|
2022-07-12 17:47:02 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
g3p_output << ", ";
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [eq, var1, var2, var3, param] {vectorToTuple<5>(vidx)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int var1_col {getJacobianCol(var1, false) + 1};
|
|
|
|
|
int var2_col {getJacobianCol(var2, false) + 1};
|
|
|
|
|
int var3_col {getJacobianCol(var3, false) + 1};
|
|
|
|
|
int param_col {getTypeSpecificIDByDerivID(param) + 1};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
g3p_output << R"({"eq": )" << eq + 1;
|
|
|
|
|
else
|
|
|
|
|
g3p_output << R"({"row": )" << eq + 1;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
g3p_output << R"(, "var1_col": )" << var1_col + 1 << R"(, "var2_col": )" << var2_col + 1
|
|
|
|
|
<< R"(, "var3_col": )" << var3_col + 1 << R"(, "param_col": )"
|
|
|
|
|
<< param_col + 1;
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
if (writeDetails)
|
|
|
|
|
g3p_output << R"(, "var1": ")" << getNameByDerivID(var1) << R"(")"
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"(, "lag1": )" << getLagByDerivID(var1) << R"(, "var2": ")"
|
|
|
|
|
<< getNameByDerivID(var2) << R"(")"
|
|
|
|
|
<< R"(, "lag2": )" << getLagByDerivID(var2) << R"(, "var3": ")"
|
|
|
|
|
<< getNameByDerivID(var3) << R"(")"
|
|
|
|
|
<< R"(, "lag3": )" << getLagByDerivID(var3) << R"(, "param": ")"
|
|
|
|
|
<< getNameByDerivID(param) << R"(")";
|
2022-07-12 17:47:02 +02:00
|
|
|
|
|
|
|
|
|
g3p_output << R"(, "val": ")";
|
|
|
|
|
d->writeJsonOutput(g3p_output, temp_term_union, tef_terms);
|
|
|
|
|
g3p_output << R"("})" << endl;
|
|
|
|
|
}
|
|
|
|
|
g3p_output << "]}" << endl;
|
|
|
|
|
}
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
return {move(mlv_output), move(tt_output), move(rp_output), move(gp_output),
|
|
|
|
|
move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output)};
|
2022-07-12 17:47:02 +02:00
|
|
|
|
}
|
|
|
|
|
|
2022-09-14 17:07:08 +02:00
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeDriverSparseIndicesHelper(ostream& output) const
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
// TODO: when C++20 support is complete, mark this constexpr
|
|
|
|
|
const string model_name {dynamic ? "dynamic" : "static"};
|
|
|
|
|
|
|
|
|
|
// Write indices for the sparse Jacobian (both naive and CSC storage)
|
|
|
|
|
output << "M_." << model_name << "_g1_sparse_rowval = int32([";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
|
|
|
|
|
output << indices.first + 1 << ' ';
|
|
|
|
|
output << "]);" << endl << "M_." << model_name << "_g1_sparse_colval = int32([";
|
|
|
|
|
for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
|
|
|
|
|
output << indices.second + 1 << ' ';
|
|
|
|
|
output << "]);" << endl << "M_." << model_name << "_g1_sparse_colptr = int32([";
|
2022-09-14 17:07:08 +02:00
|
|
|
|
for (int it : jacobian_sparse_colptr)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << it + 1 << ' ';
|
2022-09-14 17:07:08 +02:00
|
|
|
|
output << "]);" << endl;
|
|
|
|
|
|
|
|
|
|
// Write indices for the sparse higher-order derivatives
|
2022-12-07 16:10:03 +01:00
|
|
|
|
for (int i {2}; i <= computed_derivs_order; i++)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
output << "M_." << model_name << "_g" << i << "_sparse_indices = int32([";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [vidx, d] : derivatives[i])
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
2022-12-12 13:17:16 +01:00
|
|
|
|
for (bool row_number {true}; // First element of vidx is row number
|
|
|
|
|
int it : vidx)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << (exchange(row_number, false) ? it : getJacobianCol(it, true)) + 1 << ' ';
|
2022-09-14 17:07:08 +02:00
|
|
|
|
output << ';' << endl;
|
|
|
|
|
}
|
|
|
|
|
output << "]);" << endl;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeJsonSparseIndicesHelper(ostream& output) const
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
// TODO: when C++20 support is complete, mark this constexpr
|
|
|
|
|
const string model_name {dynamic ? "dynamic" : "static"};
|
|
|
|
|
|
|
|
|
|
// Write indices for the sparse Jacobian (both naive and CSC storage)
|
|
|
|
|
output << '"' << model_name << R"(_g1_sparse_rowval": [)";
|
|
|
|
|
for (bool printed_something {false};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const auto& [indices, d1] : jacobian_sparse_column_major_order)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
output << ", ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << indices.first + 1;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "], " << endl << '"' << model_name << R"(_g1_sparse_colval": [)";
|
2022-09-14 17:07:08 +02:00
|
|
|
|
for (bool printed_something {false};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const auto& [indices, d1] : jacobian_sparse_column_major_order)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
output << ", ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << indices.second + 1;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
}
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "], " << endl << '"' << model_name << R"(_g1_sparse_colptr": [)";
|
|
|
|
|
for (bool printed_something {false}; int it : jacobian_sparse_colptr)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
output << ", ";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << it + 1;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
}
|
|
|
|
|
output << ']' << endl;
|
|
|
|
|
|
|
|
|
|
// Write indices for the sparse higher-order derivatives
|
2022-12-07 16:10:03 +01:00
|
|
|
|
for (int i {2}; i <= computed_derivs_order; i++)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
output << R"(, ")" << model_name << "_g" << i << R"(_sparse_indices": [)";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something {false}; const auto& [vidx, d] : derivatives[i])
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
|
|
|
|
if (exchange(printed_something, true))
|
|
|
|
|
output << ", ";
|
|
|
|
|
output << '[';
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (bool printed_something2 {false}; int it : vidx)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
{
|
2022-12-12 13:17:16 +01:00
|
|
|
|
if (printed_something2)
|
2022-09-14 17:07:08 +02:00
|
|
|
|
output << ", ";
|
2022-12-12 13:17:16 +01:00
|
|
|
|
// First element of vidx is row number
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << (exchange(printed_something2, true) ? getJacobianCol(it, true) : it) + 1;
|
2022-09-14 17:07:08 +02:00
|
|
|
|
}
|
|
|
|
|
output << ']' << endl;
|
|
|
|
|
}
|
|
|
|
|
output << ']' << endl;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2022-09-28 19:18:15 +02:00
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeBlockDriverSparseIndicesHelper(ostream& output) const
|
2022-09-28 19:18:15 +02:00
|
|
|
|
{
|
|
|
|
|
for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string struct_name {"M_.block_structure"s + (dynamic ? "" : "_stat") + ".block("
|
|
|
|
|
+ to_string(blk + 1) + ")."};
|
2022-09-28 19:18:15 +02:00
|
|
|
|
|
|
|
|
|
// Write indices for the sparse Jacobian (both naive and CSC storage)
|
|
|
|
|
output << struct_name << "g1_sparse_rowval = int32([";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [indices, d1] : blocks_jacobian_sparse_column_major_order.at(blk))
|
|
|
|
|
output << indices.first + 1 << ' ';
|
|
|
|
|
output << "]);" << endl << struct_name << "g1_sparse_colval = int32([";
|
|
|
|
|
for (const auto& [indices, d1] : blocks_jacobian_sparse_column_major_order.at(blk))
|
|
|
|
|
output << indices.second + 1 << ' ';
|
|
|
|
|
output << "]);" << endl << struct_name << "g1_sparse_colptr = int32([";
|
2022-09-28 19:18:15 +02:00
|
|
|
|
for (int it : blocks_jacobian_sparse_colptr.at(blk))
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << it + 1 << ' ';
|
2022-09-28 19:18:15 +02:00
|
|
|
|
output << "]);" << endl;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
template<ExprNodeOutputType output_type>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeSparsePerBlockJacobianHelper(int blk, ostream& output,
|
|
|
|
|
temporary_terms_t& temporary_terms) const
|
2022-09-28 19:18:15 +02:00
|
|
|
|
{
|
|
|
|
|
static_assert(isSparseModelOutput(output_type));
|
|
|
|
|
|
|
|
|
|
// NB: stochastic mode is currently unsupported by sparse representation
|
|
|
|
|
/* See also the comment above the definition of
|
|
|
|
|
blocks_jacobian_sparse_column_major_order and
|
|
|
|
|
blocks_jacobian_sparse_column_major_colptr */
|
|
|
|
|
|
|
|
|
|
assert(blocks[blk].simulation_type != BlockSimulationType::evaluateForward
|
|
|
|
|
&& blocks[blk].simulation_type != BlockSimulationType::evaluateBackward);
|
|
|
|
|
|
|
|
|
|
int k {0};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& [row_col, d1] : blocks_jacobian_sparse_column_major_order[blk])
|
2022-09-28 19:18:15 +02:00
|
|
|
|
{
|
|
|
|
|
output << "g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< "=";
|
2022-09-28 19:18:15 +02:00
|
|
|
|
d1->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
|
|
|
|
|
output << ";" << endl;
|
|
|
|
|
k++;
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
|
2022-09-20 18:28:30 +02:00
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeSparseModelJuliaFiles(const string& basename) const
|
2022-09-20 18:28:30 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper < dynamic
|
|
|
|
|
? ExprNodeOutputType::juliaSparseDynamicModel
|
|
|
|
|
: ExprNodeOutputType::juliaSparseStaticModel > ();
|
2022-09-20 18:28:30 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
filesystem::path julia_dir {filesystem::path {basename} / "model" / "julia"};
|
2022-09-20 18:28:30 +02:00
|
|
|
|
// TODO: when C++20 support is complete, mark the following strings constexpr
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prefix {dynamic ? "SparseDynamic" : "SparseStatic"};
|
|
|
|
|
const string ss_argin {dynamic ? ", steady_state::Vector{<: Real}" : ""};
|
|
|
|
|
const string ss_argout {dynamic ? ", steady_state" : ""};
|
|
|
|
|
const int ylen {(dynamic ? 3 : 1) * symbol_table.endo_nbr()};
|
|
|
|
|
const int xlen {symbol_table.exo_nbr() + symbol_table.exo_det_nbr()};
|
2022-09-20 18:28:30 +02:00
|
|
|
|
|
|
|
|
|
size_t ttlen {0};
|
|
|
|
|
|
|
|
|
|
stringstream output;
|
|
|
|
|
|
|
|
|
|
// ResidTT!
|
|
|
|
|
output << "function " << prefix << "ResidTT!(T::Vector{<: Real}, "
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}" << ss_argin << ")"
|
|
|
|
|
<< endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< "@inbounds begin" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[0].str() << "end" << endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " return nothing" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "end" << endl
|
|
|
|
|
<< endl;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
writeToFileIfModified(output, julia_dir / (prefix + "ResidTT!.jl"));
|
|
|
|
|
ttlen += temporary_terms_derivatives[0].size();
|
|
|
|
|
|
|
|
|
|
// Resid!
|
|
|
|
|
output.str("");
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function " << prefix
|
|
|
|
|
<< "Resid!(T::Vector{<: Real}, residual::AbstractVector{<: Real}, "
|
|
|
|
|
<< "y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}" << ss_argin << ")"
|
|
|
|
|
<< endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " @assert length(T) >= " << ttlen << endl
|
|
|
|
|
<< " @assert length(residual) == " << equations.size() << endl
|
|
|
|
|
<< " @assert length(y) == " << ylen << endl
|
|
|
|
|
<< " @assert length(x) == " << xlen << endl
|
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
|
<< "@inbounds begin" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< d_sparse_output[0].str() << "end" << endl;
|
|
|
|
|
output << " return nothing" << endl << "end" << endl << endl;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
writeToFileIfModified(output, julia_dir / (prefix + "Resid!.jl"));
|
|
|
|
|
|
|
|
|
|
// G1TT!
|
|
|
|
|
output.str("");
|
|
|
|
|
output << "function " << prefix << "G1TT!(T::Vector{<: Real}, y::Vector{<: Real}, "
|
|
|
|
|
<< "x::Vector{<: Real}, params::Vector{<: Real}" << ss_argin << ")" << endl
|
|
|
|
|
<< " " << prefix << "ResidTT!(T, y, x, params" << ss_argout << ")" << endl
|
|
|
|
|
<< "@inbounds begin" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[1].str() << "end" << endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " return nothing" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "end" << endl
|
|
|
|
|
<< endl;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
writeToFileIfModified(output, julia_dir / (prefix + "G1TT!.jl"));
|
|
|
|
|
ttlen += temporary_terms_derivatives[1].size();
|
|
|
|
|
|
|
|
|
|
// G1!
|
|
|
|
|
output.str("");
|
|
|
|
|
output << "function " << prefix << "G1!(T::Vector{<: Real}, g1_v::Vector{<: Real}, "
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}" << ss_argin << ")"
|
|
|
|
|
<< endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " @assert length(T) >= " << ttlen << endl
|
|
|
|
|
<< " @assert length(g1_v) == " << derivatives[1].size() << endl
|
|
|
|
|
<< " @assert length(y) == " << ylen << endl
|
|
|
|
|
<< " @assert length(x) == " << xlen << endl
|
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
|
<< "@inbounds begin" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< d_sparse_output[1].str() << "end" << endl;
|
|
|
|
|
output << " return nothing" << endl << "end" << endl << endl;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
writeToFileIfModified(output, julia_dir / (prefix + "G1!.jl"));
|
|
|
|
|
|
|
|
|
|
for (int i {2}; i <= computed_derivs_order; i++)
|
|
|
|
|
{
|
|
|
|
|
// G<i>TT!
|
|
|
|
|
output.str("");
|
|
|
|
|
output << "function " << prefix << "G" << i << "TT!(T::Vector{<: Real}, y::Vector{<: Real}, "
|
|
|
|
|
<< "x::Vector{<: Real}, params::Vector{<: Real}" << ss_argin << ")" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " " << prefix << "G" << to_string(i - 1) << "TT!(T, y, x, params" << ss_argout
|
|
|
|
|
<< ")" << endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< "@inbounds begin" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[i].str() << "end" << endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " return nothing" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "end" << endl
|
|
|
|
|
<< endl;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
writeToFileIfModified(output, julia_dir / (prefix + "G" + to_string(i) + "TT!.jl"));
|
|
|
|
|
ttlen += temporary_terms_derivatives[i].size();
|
|
|
|
|
|
|
|
|
|
// G<i>!
|
|
|
|
|
output.str("");
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function " << prefix << "G" << i << "!(T::Vector{<: Real}, g" << i
|
|
|
|
|
<< "_v::Vector{<: Real}, "
|
|
|
|
|
<< "y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}" << ss_argin << ")"
|
|
|
|
|
<< endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " @assert length(T) >= " << ttlen << endl
|
|
|
|
|
<< " @assert length(g" << i << "_v) == " << derivatives[i].size() << endl
|
|
|
|
|
<< " @assert length(y) == " << ylen << endl
|
|
|
|
|
<< " @assert length(x) == " << xlen << endl
|
|
|
|
|
<< " @assert length(params) == " << symbol_table.param_nbr() << endl
|
|
|
|
|
<< "@inbounds begin" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< d_sparse_output[i].str() << "end" << endl
|
2022-09-20 18:28:30 +02:00
|
|
|
|
<< " return nothing" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "end" << endl
|
|
|
|
|
<< endl;
|
2022-09-20 18:28:30 +02:00
|
|
|
|
writeToFileIfModified(output, julia_dir / (prefix + "G" + to_string(i) + "!.jl"));
|
|
|
|
|
}
|
|
|
|
|
}
|
2022-09-30 17:26:43 +02:00
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeSparseModelMFiles(const string& basename) const
|
2022-09-30 17:26:43 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
constexpr ExprNodeOutputType output_type {dynamic ? ExprNodeOutputType::matlabSparseDynamicModel
|
|
|
|
|
: ExprNodeOutputType::matlabSparseStaticModel};
|
2022-09-30 17:26:43 +02:00
|
|
|
|
auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper<output_type>();
|
|
|
|
|
|
|
|
|
|
const filesystem::path m_dir {packageDir(basename) / "+sparse"};
|
|
|
|
|
// TODO: when C++20 support is complete, mark the following strings constexpr
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prefix {dynamic ? "dynamic_" : "static_"};
|
|
|
|
|
const string full_prefix {basename + ".sparse." + prefix};
|
|
|
|
|
const string ss_arg {dynamic ? ", steady_state" : ""};
|
2022-09-30 17:26:43 +02:00
|
|
|
|
|
|
|
|
|
size_t ttlen {temporary_terms_derivatives[0].size()};
|
|
|
|
|
|
|
|
|
|
ofstream output;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto open_file = [&output](const filesystem::path& p) {
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output.open(p, ios::out | ios::binary);
|
|
|
|
|
if (!output.is_open())
|
|
|
|
|
{
|
|
|
|
|
cerr << "ERROR: Can't open file " << p.string() << " for writing" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
// Residuals (non-block)
|
2023-01-16 16:55:14 +01:00
|
|
|
|
open_file(m_dir / (prefix + "resid_tt.m"));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T)" << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< "if T_order >= 0" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< " return" << endl
|
|
|
|
|
<< "end" << endl
|
|
|
|
|
<< "T_order = 0;" << endl
|
|
|
|
|
<< "if size(T, 1) < " << ttlen << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< " T = [T; NaN(" << ttlen << " - size(T, 1), 1)];" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[0].str() << "end" << endl;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
|
|
open_file(m_dir / (prefix + "resid.m"));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T)" << endl
|
|
|
|
|
<< "if nargin < " << 5 + static_cast<int>(dynamic) << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< " T_order = -1;" << endl
|
|
|
|
|
<< " T = NaN(" << ttlen << ", 1);" << endl
|
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T);" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "residual = NaN(" << equations.size() << ", 1);" << endl
|
|
|
|
|
<< d_sparse_output[0].str();
|
|
|
|
|
output << "end" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
|
|
// Jacobian (non-block)
|
|
|
|
|
ttlen += temporary_terms_derivatives[1].size();
|
|
|
|
|
|
2023-01-16 16:55:14 +01:00
|
|
|
|
open_file(m_dir / (prefix + "g1_tt.m"));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T)" << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< "if T_order >= 1" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< " return" << endl
|
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T);" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "T_order = 1;" << endl
|
|
|
|
|
<< "if size(T, 1) < " << ttlen << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< " T = [T; NaN(" << ttlen << " - size(T, 1), 1)];" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[1].str() << "end" << endl;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
|
|
open_file(m_dir / (prefix + "g1.m"));
|
|
|
|
|
// NB: At first order, sparse indices are passed as extra arguments
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function [g1, T_order, T] = " << prefix << "g1(y, x, params" << ss_arg
|
|
|
|
|
<< ", sparse_rowval, sparse_colval, sparse_colptr, T_order, T)" << endl
|
|
|
|
|
<< "if nargin < " << 8 + static_cast<int>(dynamic) << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< " T_order = -1;" << endl
|
|
|
|
|
<< " T = NaN(" << ttlen << ", 1);" << endl
|
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T);"
|
|
|
|
|
<< endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "g1_v = NaN(" << jacobian_sparse_column_major_order.size() << ", 1);" << endl
|
|
|
|
|
<< d_sparse_output[1].str();
|
2022-12-08 14:32:26 +01:00
|
|
|
|
// On MATLAB < R2020a, sparse() does not accept int32 indices
|
|
|
|
|
output << "if ~isoctave && matlab_ver_less_than('9.8')" << endl
|
|
|
|
|
<< " sparse_rowval = double(sparse_rowval);" << endl
|
|
|
|
|
<< " sparse_colval = double(sparse_colval);" << endl
|
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << equations.size() << ", "
|
|
|
|
|
<< getJacobianColsNbr(true) << ");" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "end" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
|
|
// Higher-order derivatives (non-block)
|
|
|
|
|
for (int i {2}; i <= computed_derivs_order; i++)
|
|
|
|
|
{
|
|
|
|
|
ttlen += temporary_terms_derivatives[i].size();
|
|
|
|
|
|
2023-01-16 16:55:14 +01:00
|
|
|
|
open_file(m_dir / (prefix + "g" + to_string(i) + "_tt.m"));
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output << "function T = " << prefix << "g" << i << "_tt(y, x, params" << ss_arg << ")" << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< "if T_order >= " << i << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< " return" << endl
|
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "[T_order, T] = " << full_prefix << "g" << i - 1 << "_tt(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T);" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "T_order = " << i << ";" << endl
|
|
|
|
|
<< "if size(T, 1) < " << ttlen << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< " T = [T; NaN(" << ttlen << " - size(T, 1), 1)];" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[i].str() << "end" << endl;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output.close();
|
|
|
|
|
|
|
|
|
|
open_file(m_dir / (prefix + "g" + to_string(i) + ".m"));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function [g" << i << "_v, T_order, T] = " << prefix << "g" << i << "(y, x, params"
|
|
|
|
|
<< ss_arg << ", T_order, T)" << endl
|
|
|
|
|
<< "if nargin < " << 5 + static_cast<int>(dynamic) << endl
|
2023-01-10 15:58:45 +01:00
|
|
|
|
<< " T_order = -1;" << endl
|
|
|
|
|
<< " T = NaN(" << ttlen << ", 1);" << endl
|
|
|
|
|
<< "end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << ss_arg
|
|
|
|
|
<< ", T_order, T);" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "g" << i << "_v = NaN(" << derivatives[i].size() << ", 1);" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< d_sparse_output[i].str() << "end" << endl;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output.close();
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
// Block decomposition
|
|
|
|
|
if (block_decomposed)
|
|
|
|
|
{
|
|
|
|
|
const filesystem::path block_dir {m_dir / "+block"};
|
|
|
|
|
temporary_terms_t temporary_terms_written;
|
|
|
|
|
for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string funcname {prefix + to_string(blk + 1)};
|
2022-09-30 17:26:43 +02:00
|
|
|
|
const BlockSimulationType simulation_type {blocks[blk].simulation_type};
|
|
|
|
|
const bool evaluate {simulation_type == BlockSimulationType::evaluateForward
|
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateBackward};
|
|
|
|
|
const string resid_g1_arg {evaluate ? "" : ", residual, g1"};
|
|
|
|
|
open_file(block_dir / (funcname + ".m"));
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "function [y, T" << resid_g1_arg << "] = " << funcname << "(y, x, params"
|
|
|
|
|
<< ss_arg << ", sparse_rowval, sparse_colval, sparse_colptr, T)" << endl;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
if (!evaluate)
|
|
|
|
|
output << "residual=NaN(" << blocks[blk].mfs_size << ", 1);" << endl;
|
|
|
|
|
|
|
|
|
|
// Write residuals and temporary terms (incl. those for derivatives)
|
|
|
|
|
writePerBlockHelper<output_type>(blk, output, temporary_terms_written);
|
|
|
|
|
|
|
|
|
|
// Write Jacobian
|
|
|
|
|
if (!evaluate)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const bool one_boundary {
|
|
|
|
|
simulation_type == BlockSimulationType::solveBackwardSimple
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveForwardSimple
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveForwardComplete};
|
2022-09-30 17:26:43 +02:00
|
|
|
|
output << "if nargout > 3" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " g1_v = NaN(" << blocks_jacobian_sparse_column_major_order[blk].size()
|
|
|
|
|
<< ", 1);" << endl;
|
2022-09-30 17:26:43 +02:00
|
|
|
|
writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
|
2022-12-08 14:32:26 +01:00
|
|
|
|
// On MATLAB < R2020a, sparse() does not accept int32 indices
|
|
|
|
|
output << " if ~isoctave && matlab_ver_less_than('9.8')" << endl
|
|
|
|
|
<< " sparse_rowval = double(sparse_rowval);" << endl
|
|
|
|
|
<< " sparse_colval = double(sparse_colval);" << endl
|
|
|
|
|
<< " end" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " g1 = sparse(sparse_rowval, sparse_colval, g1_v, "
|
|
|
|
|
<< blocks[blk].mfs_size << ", "
|
|
|
|
|
<< (one_boundary ? 1 : 3) * blocks[blk].mfs_size << ");" << endl
|
2022-09-30 17:26:43 +02:00
|
|
|
|
<< "end" << endl;
|
|
|
|
|
}
|
|
|
|
|
output << "end" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
|
|
|
|
|
const filesystem::path& matlabroot) const
|
2022-10-05 17:45:18 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
constexpr ExprNodeOutputType output_type {dynamic ? ExprNodeOutputType::CSparseDynamicModel
|
|
|
|
|
: ExprNodeOutputType::CSparseStaticModel};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper<output_type>();
|
|
|
|
|
|
|
|
|
|
const filesystem::path mex_dir {packageDir(basename) / "+sparse"};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path model_src_dir {filesystem::path {basename} / "model" / "src" / "sparse"};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
// TODO: when C++20 support is complete, mark the following strings constexpr
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prefix {dynamic ? "dynamic_" : "static_"};
|
|
|
|
|
const string ss_argin {dynamic ? ", const double *restrict steady_state" : ""};
|
|
|
|
|
const string ss_argout {dynamic ? ", steady_state" : ""};
|
|
|
|
|
const int ylen {(dynamic ? 3 : 1) * symbol_table.endo_nbr()};
|
|
|
|
|
const int xlen {symbol_table.exo_nbr() + symbol_table.exo_det_nbr()};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
|
|
|
|
vector<filesystem::path> tt_object_files;
|
|
|
|
|
|
|
|
|
|
ofstream output;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto open_file = [&output](const filesystem::path& p) {
|
2022-10-05 17:45:18 +02:00
|
|
|
|
output.open(p, ios::out | ios::binary);
|
|
|
|
|
if (!output.is_open())
|
|
|
|
|
{
|
|
|
|
|
cerr << "ERROR: Can't open file " << p.string() << " for writing" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
size_t ttlen {0};
|
|
|
|
|
|
|
|
|
|
// Helper for dealing with y, x, params and steady_state inputs (shared with block case)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto y_x_params_ss_inputs = [&](bool assign_y) {
|
|
|
|
|
output << " if (!(mxIsDouble(prhs[0]) && !mxIsComplex(prhs[0]) && !mxIsSparse(prhs[0]) && "
|
|
|
|
|
"mxGetNumberOfElements(prhs[0]) == "
|
|
|
|
|
<< ylen << "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("y must be a real dense numeric array with )" << ylen
|
|
|
|
|
<< R"( elements");)" << endl;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
if (assign_y)
|
|
|
|
|
output << " const double *restrict y = mxGetPr(prhs[0]);" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " if (!(mxIsDouble(prhs[1]) && !mxIsComplex(prhs[1]) && !mxIsSparse(prhs[1]) && "
|
|
|
|
|
"mxGetNumberOfElements(prhs[1]) == "
|
|
|
|
|
<< xlen << "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("x must be a real dense numeric array with )" << xlen
|
|
|
|
|
<< R"( elements");)" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " const double *restrict x = mxGetPr(prhs[1]);" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " if (!(mxIsDouble(prhs[2]) && !mxIsComplex(prhs[2]) && !mxIsSparse(prhs[2]) && "
|
|
|
|
|
"mxGetNumberOfElements(prhs[2]) == "
|
|
|
|
|
<< symbol_table.param_nbr() << "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("params must be a real dense numeric array with )"
|
|
|
|
|
<< symbol_table.param_nbr() << R"( elements");)" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " const double *restrict params = mxGetPr(prhs[2]);" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
|
|
|
|
output << " if (!(mxIsDouble(prhs[3]) && !mxIsComplex(prhs[3]) && !mxIsSparse(prhs[3]) && "
|
|
|
|
|
"mxGetNumberOfElements(prhs[3]) == "
|
|
|
|
|
<< symbol_table.endo_nbr() << "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("steady_state must be a real dense numeric array with )"
|
|
|
|
|
<< symbol_table.endo_nbr() << R"( elements");)" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " const double *restrict steady_state = mxGetPr(prhs[3]);" << endl;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
// Helper for dealing with sparse_rowval and sparse_colptr inputs (shared with block case)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto sparse_indices_inputs = [&](int ncols, int nzval) {
|
2022-10-05 17:45:18 +02:00
|
|
|
|
// We use sparse_rowval and sparse_colptr (sparse_colval is unused)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const int row_idx {3 + static_cast<int>(dynamic)}, col_idx {row_idx + 2};
|
|
|
|
|
output << " if (!(mxIsInt32(prhs[" << row_idx << "]) && mxGetNumberOfElements(prhs[" << row_idx
|
|
|
|
|
<< "]) == " << nzval << "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("sparse_rowval must be an int32 array with )" << nzval
|
|
|
|
|
<< R"( elements");)" << endl
|
|
|
|
|
<< " if (!(mxIsInt32(prhs[" << col_idx << "]) && mxGetNumberOfElements(prhs[" << col_idx
|
|
|
|
|
<< "]) == " << ncols + 1 << "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("sparse_colptr must be an int32 array with )" << ncols + 1
|
|
|
|
|
<< R"( elements");)" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< "#if MX_HAS_INTERLEAVED_COMPLEX" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " const int32_T *restrict sparse_rowval = mxGetInt32s(prhs[" << row_idx << "]);"
|
|
|
|
|
<< endl
|
|
|
|
|
<< " const int32_T *restrict sparse_colptr = mxGetInt32s(prhs[" << col_idx << "]);"
|
|
|
|
|
<< endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< "#else" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " const int32_T *restrict sparse_rowval = (int32_T *) mxGetData(prhs[" << row_idx
|
|
|
|
|
<< "]);" << endl
|
|
|
|
|
<< " const int32_T *restrict sparse_colptr = (int32_T *) mxGetData(prhs[" << col_idx
|
|
|
|
|
<< "]);" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< "#endif" << endl;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
// Helper for creating sparse Jacobian (shared with block case)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto sparse_jacobian_create = [&](int argidx, int nrows, int ncols, int nzval) {
|
|
|
|
|
output << " plhs[" << argidx << "] = mxCreateSparse(" << nrows << ", " << ncols << ", "
|
|
|
|
|
<< nzval << ", mxREAL);" << endl
|
|
|
|
|
<< " mwIndex *restrict ir = mxGetIr(plhs[" << argidx
|
|
|
|
|
<< "]), *restrict jc = mxGetJc(plhs[" << argidx << "]);" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " for (mwSize i = 0; i < " << nzval << "; i++)" << endl
|
|
|
|
|
<< " *ir++ = *sparse_rowval++ - 1;" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " for (mwSize i = 0; i < " << ncols + 1 << "; i++)" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " *jc++ = *sparse_colptr++ - 1;" << endl;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
for (int i {0}; i <= computed_derivs_order; i++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string funcname {prefix + (i == 0 ? "resid" : "g" + to_string(i))};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
ttlen += temporary_terms_derivatives[i].size();
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prototype_tt {
|
|
|
|
|
"void " + funcname
|
|
|
|
|
+ "_tt(const double *restrict y, const double *restrict x, const double *restrict params"
|
|
|
|
|
+ ss_argin + ", double *restrict T)"};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path header_tt {model_src_dir / (funcname + "_tt.h")};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
open_file(header_tt);
|
|
|
|
|
output << prototype_tt << ";" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path source_tt {model_src_dir / (funcname + "_tt.c")};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
open_file(source_tt);
|
|
|
|
|
output << "#include <math.h>" << endl
|
|
|
|
|
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
|
2023-04-07 15:45:44 +02:00
|
|
|
|
<< endl;
|
2023-04-11 15:29:12 +02:00
|
|
|
|
writeCHelpersDefinition(output);
|
2023-04-07 15:45:44 +02:00
|
|
|
|
output << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< prototype_tt << endl
|
|
|
|
|
<< "{" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< tt_sparse_output[i].str() << "}" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< endl;
|
|
|
|
|
output.close();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
tt_object_files.push_back(
|
|
|
|
|
compileMEX(model_src_dir, funcname + "_tt", mexext, {source_tt}, matlabroot, false));
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prototype_main {
|
|
|
|
|
"void " + funcname
|
|
|
|
|
+ "(const double *restrict y, const double *restrict x, const double *restrict params"
|
|
|
|
|
+ ss_argin + ", const double *restrict T, double *restrict "
|
|
|
|
|
+ (i == 0 ? "residual" : "g" + to_string(i) + "_v") + ")"};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path header_main {model_src_dir / (funcname + ".h")};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
open_file(header_main);
|
|
|
|
|
output << prototype_main << ";" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path source_main {model_src_dir / (funcname + ".c")};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
open_file(source_main);
|
|
|
|
|
output << "#include <math.h>" << endl
|
|
|
|
|
<< R"(#include "mex.h")" << endl // Needed for calls to external functions
|
|
|
|
|
<< endl;
|
2023-04-11 15:29:12 +02:00
|
|
|
|
writeCHelpersDefinition(output);
|
|
|
|
|
writeCHelpersDeclaration(output); // Provide external definition of helpers in main file
|
2022-10-05 17:45:18 +02:00
|
|
|
|
output << endl
|
|
|
|
|
<< prototype_main << endl
|
2023-01-04 14:43:57 +01:00
|
|
|
|
<< "{" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< d_sparse_output[i].str() << "}" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< endl;
|
|
|
|
|
output.close();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
auto main_object_file {
|
|
|
|
|
compileMEX(model_src_dir, funcname, mexext, {source_main}, matlabroot, false)};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const filesystem::path source_mex {model_src_dir / (funcname + "_mex.c")};
|
|
|
|
|
int nargin {5 + static_cast<int>(dynamic) + 3 * static_cast<int>(i == 1)};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
open_file(source_mex);
|
|
|
|
|
output << "#include <string.h>" << endl // For memcpy()
|
|
|
|
|
<< R"(#include "mex.h")" << endl
|
|
|
|
|
<< R"(#include ")" << funcname << R"(.h")" << endl;
|
|
|
|
|
for (int j {0}; j <= i; j++)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << R"(#include ")" << prefix << (j == 0 ? "resid" : "g" + to_string(j))
|
|
|
|
|
<< R"(_tt.h")" << endl;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
output << endl
|
|
|
|
|
<< "#define max(a, b) ((a > b) ? (a) : (b))" << endl
|
|
|
|
|
<< endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])"
|
|
|
|
|
<< endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< "{" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< " if (nrhs != " << nargin - 2 << " && nrhs != " << nargin << ")" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Accepts exactly )" << nargin - 2 << " or " << nargin
|
|
|
|
|
<< R"( input arguments");)" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " if (nlhs != 1 && nlhs != 3)" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Accepts exactly 1 or 3 output arguments");)" << endl;
|
|
|
|
|
|
|
|
|
|
y_x_params_ss_inputs(true);
|
|
|
|
|
|
|
|
|
|
if (i == 1)
|
|
|
|
|
sparse_indices_inputs(getJacobianColsNbr(true), jacobian_sparse_column_major_order.size());
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output
|
|
|
|
|
<< " mxArray *T_mx, *T_order_mx;" << endl
|
|
|
|
|
<< " int T_order_on_input;" << endl
|
|
|
|
|
<< " if (nrhs > " << nargin - 2 << ")" << endl
|
|
|
|
|
<< " {" << endl
|
|
|
|
|
<< " T_order_mx = (mxArray *) prhs[" << nargin - 2 << "];" << endl
|
|
|
|
|
<< " T_mx = (mxArray *) prhs[" << nargin - 1 << "];" << endl
|
|
|
|
|
<< " if (!(mxIsScalar(T_order_mx) && mxIsNumeric(T_order_mx)))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("T_order should be a numeric scalar");)" << endl
|
|
|
|
|
<< " if (!(mxIsDouble(T_mx) && !mxIsComplex(T_mx) && !mxIsSparse(T_mx) && "
|
|
|
|
|
"mxGetN(T_mx) == 1))"
|
|
|
|
|
<< endl
|
|
|
|
|
<< R"( mexErrMsgTxt("T_mx should be a real dense column vector");)" << endl
|
|
|
|
|
<< " T_order_on_input = mxGetScalar(T_order_mx);" << endl
|
|
|
|
|
<< " if (T_order_on_input < " << i << ")" << endl
|
|
|
|
|
<< " {" << endl
|
|
|
|
|
<< " T_order_mx = mxCreateDoubleScalar(" << i << ");" << endl
|
|
|
|
|
<< " const mxArray *T_old_mx = T_mx;" << endl
|
|
|
|
|
<< " T_mx = mxCreateDoubleMatrix(max(" << ttlen
|
|
|
|
|
<< ", mxGetM(T_old_mx)), 1, mxREAL);" << endl
|
|
|
|
|
<< " memcpy(mxGetPr(T_mx), mxGetPr(T_old_mx), mxGetM(T_old_mx)*sizeof(double));"
|
|
|
|
|
<< endl
|
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< " else if (mxGetM(T_mx) < " << ttlen << ")" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("T_mx should have at least )" << ttlen << R"( elements");)"
|
|
|
|
|
<< endl
|
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< " else" << endl
|
|
|
|
|
<< " {" << endl
|
|
|
|
|
<< " T_order_mx = mxCreateDoubleScalar(" << i << ");" << endl
|
|
|
|
|
<< " T_mx = mxCreateDoubleMatrix(" << ttlen << ", 1, mxREAL);" << endl
|
|
|
|
|
<< " T_order_on_input = -1;" << endl
|
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< " double *restrict T = mxGetPr(T_mx);" << endl
|
|
|
|
|
<< " if (T_order_on_input < " << i << ")" << endl
|
|
|
|
|
<< " switch (T_order_on_input)" << endl
|
|
|
|
|
<< " {" << endl;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
for (int j {0}; j <= i; j++)
|
|
|
|
|
{
|
|
|
|
|
if (j == 0)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " default:" << endl << " " << prefix << "resid";
|
2022-10-05 17:45:18 +02:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " case " << j - 1 << ":" << endl << " " << prefix << "g" << j;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
output << "_tt(y, x, params" << ss_argout << ", T);" << endl;
|
|
|
|
|
}
|
|
|
|
|
output << " }" << endl;
|
|
|
|
|
if (i == 1)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
sparse_jacobian_create(0, equations.size(), getJacobianColsNbr(true),
|
|
|
|
|
jacobian_sparse_column_major_order.size());
|
2022-10-05 17:45:18 +02:00
|
|
|
|
else
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " plhs[0] = mxCreateDoubleMatrix("
|
|
|
|
|
<< (i == 0 ? equations.size() : derivatives[i].size()) << ", 1, mxREAL);" << endl;
|
|
|
|
|
output << " " << prefix << (i == 0 ? "resid" : "g" + to_string(i)) << "(y, x, params"
|
|
|
|
|
<< ss_argout << ", T, mxGetPr(plhs[0]));" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " if (nlhs == 3)" << endl
|
|
|
|
|
<< " {" << endl
|
|
|
|
|
<< " plhs[1] = T_order_mx;" << endl
|
|
|
|
|
<< " plhs[2] = T_mx;" << endl
|
|
|
|
|
<< " }" << endl
|
|
|
|
|
<< "}" << endl;
|
|
|
|
|
output.close();
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
vector<filesystem::path> mex_input_files {main_object_file, source_mex};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
for (int j {0}; j <= i; j++)
|
|
|
|
|
mex_input_files.push_back(tt_object_files[j]);
|
2023-03-20 17:54:32 +01:00
|
|
|
|
compileMEX(mex_dir, funcname, mexext, mex_input_files, matlabroot);
|
2022-10-05 17:45:18 +02:00
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
if (block_decomposed)
|
|
|
|
|
{
|
|
|
|
|
const filesystem::path block_dir {mex_dir / "+block"};
|
|
|
|
|
temporary_terms_t temporary_terms_written;
|
|
|
|
|
size_t total_blk_ttlen {0}; // Temporary terms of all blocks up to the present one
|
|
|
|
|
for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string funcname {prefix + to_string(blk + 1)};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
const filesystem::path source_mex {model_src_dir / "block" / (funcname + ".c")};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int nargin {7 + static_cast<int>(dynamic)};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
const BlockSimulationType simulation_type {blocks[blk].simulation_type};
|
|
|
|
|
const bool evaluate {simulation_type == BlockSimulationType::evaluateForward
|
|
|
|
|
|| simulation_type == BlockSimulationType::evaluateBackward};
|
|
|
|
|
const bool one_boundary {simulation_type == BlockSimulationType::solveBackwardSimple
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveForwardSimple
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveBackwardComplete
|
|
|
|
|
|| simulation_type == BlockSimulationType::solveForwardComplete};
|
|
|
|
|
// The following two variables are not used for “evaluate” blocks
|
2023-11-30 15:28:57 +01:00
|
|
|
|
int g1_ncols {(one_boundary ? 1 : 3) * blocks[blk].mfs_size};
|
2022-10-05 17:45:18 +02:00
|
|
|
|
open_file(source_mex);
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << "#include <math.h>" << endl << R"(#include "mex.h")" << endl << endl;
|
2023-04-11 15:29:12 +02:00
|
|
|
|
writeCHelpersDefinition(output);
|
|
|
|
|
writeCHelpersDeclaration(output); // Provide external definition of helpers
|
2023-04-07 15:45:44 +02:00
|
|
|
|
output << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])"
|
|
|
|
|
<< endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< "{" << endl
|
|
|
|
|
<< " if (nrhs != " << nargin << ")" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( mexErrMsgTxt("Accepts exactly )" << nargin << R"( input arguments");)"
|
|
|
|
|
<< endl;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
if (evaluate)
|
|
|
|
|
output << " if (nlhs != 2)" << endl
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< R"( mexErrMsgTxt("Accepts exactly 2 output arguments");)" << endl;
|
2022-10-05 17:45:18 +02:00
|
|
|
|
else
|
|
|
|
|
/* Only two output arguments make sense if one only wants to
|
|
|
|
|
evaluate the recursive variables. */
|
|
|
|
|
output << " if (nlhs < 2 || nlhs > 4)" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("Accepts 2 to 4 output arguments");)" << endl;
|
|
|
|
|
y_x_params_ss_inputs(false);
|
|
|
|
|
|
|
|
|
|
/* We’d like to avoid copying y if this is a “solve” block without
|
|
|
|
|
recursive variables. Unfortunately “plhs[0]=prhs[0]” leads to
|
|
|
|
|
crashes (at least with MATLAB R2022b), though it seems to work
|
|
|
|
|
under Octave. The undocumented mxCreateSharedDataCopy() could have
|
|
|
|
|
been a solution, but was removed in MATLAB R2018a. The solution
|
|
|
|
|
seems to be to use the C++ MEX API, but it is not supported by
|
|
|
|
|
Octave and older MATLAB (and also C++ does not have the “restrict”
|
|
|
|
|
qualifier, so it may not be a net gain). See:
|
|
|
|
|
https://fr.mathworks.com/matlabcentral/answers/77048-return-large-unchange-mxarray-from-mex
|
2023-11-30 15:28:57 +01:00
|
|
|
|
https://fr.mathworks.com/matlabcentral/answers/422751-how-to-output-a-mexfunction-input-without-copy
|
|
|
|
|
*/
|
2022-10-05 17:45:18 +02:00
|
|
|
|
output << " plhs[0] = mxDuplicateArray(prhs[0]);" << endl
|
|
|
|
|
<< " double *restrict y = mxGetPr(plhs[0]);" << endl;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
// NB: For “evaluate” blocks, sparse_{rowval,colval,colptr} arguments are present but
|
|
|
|
|
// ignored
|
2022-10-05 17:45:18 +02:00
|
|
|
|
if (!evaluate)
|
|
|
|
|
sparse_indices_inputs(g1_ncols, blocks_jacobian_sparse_column_major_order[blk].size());
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
for (const auto& it : blocks_temporary_terms[blk])
|
2022-10-05 17:45:18 +02:00
|
|
|
|
total_blk_ttlen += it.size();
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " if (!(mxIsDouble(prhs[" << nargin - 1 << "]) && !mxIsComplex(prhs["
|
|
|
|
|
<< nargin - 1 << "]) && !mxIsSparse(prhs[" << nargin - 1
|
|
|
|
|
<< "]) && mxGetNumberOfElements(prhs[" << nargin - 1 << "]) >= " << total_blk_ttlen
|
|
|
|
|
<< "))" << endl
|
|
|
|
|
<< R"( mexErrMsgTxt("T must be a real dense numeric array with at least )"
|
|
|
|
|
<< total_blk_ttlen << R"( elements");)"
|
|
|
|
|
<< endl
|
|
|
|
|
/* We’d like to avoid copying T when the block has no temporary
|
|
|
|
|
terms, but the same remark as above applies. */
|
|
|
|
|
<< " plhs[1] = mxDuplicateArray(prhs[" << nargin - 1 << "]);" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " double *restrict T = mxGetPr(plhs[1]);" << endl;
|
|
|
|
|
|
|
|
|
|
if (!evaluate)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " mxArray *residual_mx = mxCreateDoubleMatrix(" << blocks[blk].mfs_size
|
|
|
|
|
<< ", 1, mxREAL);" << endl
|
2022-10-05 17:45:18 +02:00
|
|
|
|
<< " double *restrict residual = mxGetPr(residual_mx);" << endl
|
|
|
|
|
<< " if (nlhs > 2)" << endl
|
|
|
|
|
<< " plhs[2] = residual_mx;" << endl;
|
|
|
|
|
|
|
|
|
|
// Write residuals and temporary terms (incl. those for derivatives)
|
|
|
|
|
writePerBlockHelper<output_type>(blk, output, temporary_terms_written);
|
|
|
|
|
|
|
|
|
|
if (!evaluate)
|
|
|
|
|
{
|
|
|
|
|
// Write Jacobian
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << " if (nlhs > 3)" << endl << " {" << endl;
|
|
|
|
|
sparse_jacobian_create(3, blocks[blk].mfs_size, g1_ncols,
|
|
|
|
|
blocks_jacobian_sparse_column_major_order[blk].size());
|
2022-10-05 17:45:18 +02:00
|
|
|
|
output << " double *restrict g1_v = mxGetPr(plhs[3]);" << endl;
|
|
|
|
|
writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
|
|
|
|
|
output << " }" << endl;
|
|
|
|
|
}
|
|
|
|
|
output << "}" << endl;
|
|
|
|
|
output.close();
|
2023-11-30 15:28:57 +01:00
|
|
|
|
compileMEX(block_dir, funcname, mexext, {source_mex}, matlabroot);
|
2022-10-05 17:45:18 +02:00
|
|
|
|
}
|
|
|
|
|
}
|
|
|
|
|
}
|
2022-11-09 14:46:02 +01:00
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeDebugModelMFiles(const string& basename) const
|
2022-11-09 14:46:02 +01:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
constexpr ExprNodeOutputType output_type {dynamic ? ExprNodeOutputType::matlabSparseDynamicModel
|
|
|
|
|
: ExprNodeOutputType::matlabSparseStaticModel};
|
2022-11-09 14:46:02 +01:00
|
|
|
|
|
|
|
|
|
const filesystem::path m_dir {packageDir(basename) / "+debug"};
|
|
|
|
|
// TODO: when C++20 support is complete, mark the following strings constexpr
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const string prefix {dynamic ? "dynamic_" : "static_"};
|
|
|
|
|
const string ss_arg {dynamic ? ", steady_state" : ""};
|
2022-11-09 14:46:02 +01:00
|
|
|
|
|
|
|
|
|
const filesystem::path resid_filename {m_dir / (prefix + "resid.m")};
|
|
|
|
|
ofstream output {resid_filename, ios::out | ios::binary};
|
|
|
|
|
if (!output.is_open())
|
|
|
|
|
{
|
|
|
|
|
cerr << "ERROR: Can't open file " << resid_filename.string() << " for writing" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
output << "function [lhs, rhs] = " << prefix << "resid(y, x, params" << ss_arg << ")" << endl
|
|
|
|
|
<< "T = NaN(" << temporary_terms_derivatives[0].size() << ", 1);" << endl
|
|
|
|
|
<< "lhs = NaN(" << equations.size() << ", 1);" << endl
|
|
|
|
|
<< "rhs = NaN(" << equations.size() << ", 1);" << endl;
|
|
|
|
|
deriv_node_temp_terms_t tef_terms;
|
|
|
|
|
temporary_terms_t temporary_terms;
|
|
|
|
|
writeTemporaryTerms<output_type>(temporary_terms_derivatives[0], temporary_terms,
|
|
|
|
|
temporary_terms_idxs, output, tef_terms);
|
|
|
|
|
for (size_t eq {0}; eq < equations.size(); eq++)
|
|
|
|
|
{
|
|
|
|
|
output << "lhs" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< " = ";
|
2022-11-09 14:46:02 +01:00
|
|
|
|
equations[eq]->arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
|
|
|
|
|
output << ";" << endl
|
|
|
|
|
<< "rhs" << LEFT_ARRAY_SUBSCRIPT(output_type)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
|
|
|
|
|
<< " = ";
|
2022-11-09 14:46:02 +01:00
|
|
|
|
equations[eq]->arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
|
|
|
|
|
output << ";" << endl;
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
output << "end" << endl;
|
|
|
|
|
|
|
|
|
|
output.close();
|
|
|
|
|
}
|
2023-03-28 16:37:05 +02:00
|
|
|
|
|
|
|
|
|
template<bool dynamic>
|
|
|
|
|
void
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const
|
2023-03-28 16:37:05 +02:00
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
const auto output_type {julia ? (dynamic ? ExprNodeOutputType::juliaTimeDataFrame
|
|
|
|
|
: ExprNodeOutputType::juliaStaticModel)
|
|
|
|
|
: (dynamic ? ExprNodeOutputType::matlabDseries
|
|
|
|
|
: ExprNodeOutputType::matlabStaticModel)};
|
2023-03-28 16:37:05 +02:00
|
|
|
|
|
|
|
|
|
ostringstream output_func_body;
|
|
|
|
|
writeAuxVarRecursiveDefinitions(output_func_body, output_type);
|
|
|
|
|
|
|
|
|
|
if (output_func_body.str().empty())
|
|
|
|
|
return;
|
|
|
|
|
|
2023-11-30 15:28:57 +01:00
|
|
|
|
string func_name {dynamic ? "dynamic_set_auxiliary_series" : "set_auxiliary_variables"};
|
2023-03-28 16:37:05 +02:00
|
|
|
|
if (julia)
|
|
|
|
|
func_name.push_back('!');
|
|
|
|
|
const char comment {julia ? '#' : '%'};
|
|
|
|
|
|
|
|
|
|
stringstream output;
|
|
|
|
|
output << "function ";
|
|
|
|
|
if (!julia)
|
|
|
|
|
{
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-03-28 16:37:05 +02:00
|
|
|
|
output << "ds = ";
|
|
|
|
|
else
|
|
|
|
|
output << "y = ";
|
|
|
|
|
}
|
|
|
|
|
output << func_name + "(";
|
2023-11-30 15:28:57 +01:00
|
|
|
|
if constexpr (dynamic)
|
2023-03-28 16:37:05 +02:00
|
|
|
|
output << "ds";
|
|
|
|
|
else
|
|
|
|
|
output << "y, x";
|
|
|
|
|
output << ", params)" << endl
|
|
|
|
|
<< comment << endl
|
|
|
|
|
<< comment << " Computes auxiliary variables of the " << modelClassName() << endl
|
|
|
|
|
<< comment << endl;
|
|
|
|
|
if (julia)
|
|
|
|
|
output << "@inbounds begin" << endl;
|
2023-11-30 15:28:57 +01:00
|
|
|
|
output << output_func_body.str() << "end" << endl;
|
2023-03-28 16:37:05 +02:00
|
|
|
|
if (julia)
|
|
|
|
|
output << "end" << endl;
|
|
|
|
|
|
|
|
|
|
if (julia)
|
2023-11-30 15:28:57 +01:00
|
|
|
|
writeToFileIfModified(
|
|
|
|
|
output, filesystem::path {basename} / "model" / "julia"
|
|
|
|
|
/ (dynamic ? "DynamicSetAuxiliarySeries.jl" : "SetAuxiliaryVariables.jl"));
|
2023-03-28 16:37:05 +02:00
|
|
|
|
else
|
|
|
|
|
{
|
|
|
|
|
/* Calling writeToFileIfModified() is useless here since we write inside
|
|
|
|
|
a subdirectory deleted at each preprocessor run. */
|
|
|
|
|
filesystem::path filename {packageDir(basename) / (func_name + ".m")};
|
2023-11-30 15:28:57 +01:00
|
|
|
|
ofstream output_file {filename, ios::out | ios::binary};
|
2023-03-28 16:37:05 +02:00
|
|
|
|
if (!output_file.is_open())
|
|
|
|
|
{
|
|
|
|
|
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
|
|
|
|
|
exit(EXIT_FAILURE);
|
|
|
|
|
}
|
|
|
|
|
output_file << output.str();
|
|
|
|
|
output_file.close();
|
|
|
|
|
}
|
|
|
|
|
}
|
2023-12-01 15:31:55 +01:00
|
|
|
|
|
|
|
|
|
#endif
|