From a2bea00fee97dd4cdc2db32690e0d492fda2edfb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 13 May 2020 16:58:19 +0200
Subject: [PATCH] Block decomposition: another fix related to temporary terms
refactoring
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Temporary terms need to be computed per equation (as was done previously), and
not simply per block.
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.
---
src/DynamicModel.cc | 129 +++++++++++++++++++++++++-------------------
src/DynamicModel.hh | 4 +-
src/ExprNode.cc | 106 ++++++++++++++++++------------------
src/ExprNode.hh | 42 ++++++++-------
src/ModelTree.cc | 34 ++++++++----
src/ModelTree.hh | 18 +++++--
src/StaticModel.cc | 127 +++++++++++++++++++++++--------------------
7 files changed, 262 insertions(+), 198 deletions(-)
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 30fceaca..ea49debc 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -212,15 +212,15 @@ DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &inst
void
DynamicModel::additionalBlockTemporaryTerms(int blk,
- vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+ vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
for (const auto &[ignore, d] : blocks_derivatives_exo[blk])
- d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
for (const auto &[ignore, d] : blocks_derivatives_exo_det[blk])
- d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
for (const auto &[ignore, d] : blocks_derivatives_other_endo[blk])
- d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
}
void
@@ -322,17 +322,18 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
// Declare global temp terms from this block and the previous ones
bool global_keyword_written = false;
for (int blk2 = 0; blk2 <= block; blk2++)
- for (auto tt : blocks_temporary_terms[blk2])
- {
- if (!global_keyword_written)
- {
- output << " global";
- global_keyword_written = true;
- }
- output << " ";
- // In the following, "Static" is used to avoid getting the "(it_)" subscripting
- tt->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, blocks_temporary_terms[blk2], {});
- }
+ for (auto &eq_tt : blocks_temporary_terms[blk2])
+ for (auto tt : eq_tt)
+ {
+ if (!global_keyword_written)
+ {
+ output << " global";
+ global_keyword_written = true;
+ }
+ output << " ";
+ // In the following, "Static" is used to avoid getting the "(it_)" subscripting
+ tt->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, eq_tt, {});
+ }
if (global_keyword_written)
output << ";" << endl;
@@ -342,13 +343,14 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
{
output << " " << "% //Temporary variables initialization" << endl
<< " " << "T_zeros = zeros(y_kmin+periods, 1);" << endl;
- for (auto it : blocks_temporary_terms[block])
- {
- output << " ";
- // In the following, "Static" is used to avoid getting the "(it_)" subscripting
- it->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, blocks_temporary_terms[block], {});
- output << " = T_zeros;" << endl;
- }
+ for (auto &eq_tt : blocks_temporary_terms[block])
+ for (auto tt : eq_tt)
+ {
+ output << " ";
+ // In the following, "Static" is used to avoid getting the "(it_)" subscripting
+ tt->writeOutput(output, ExprNodeOutputType::matlabStaticModelSparse, eq_tt, {});
+ output << " = T_zeros;" << endl;
+ }
}
if (simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
@@ -381,23 +383,28 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
sps = "";
deriv_node_temp_terms_t tef_terms;
- for (auto it : blocks_temporary_terms[block])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, temporary_terms_idxs, tef_terms);
-
- output << " " << sps;
- it->writeOutput(output, local_output_type, blocks_temporary_terms[block], {}, tef_terms);
- output << " = ";
- it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
- temporary_terms.insert(it);
- output << ";" << endl;
- }
+
+ auto write_eq_tt = [&](int eq)
+ {
+ for (auto it : blocks_temporary_terms[block][eq])
+ {
+ if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+ it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, temporary_terms_idxs, tef_terms);
+
+ output << " " << sps;
+ it->writeOutput(output, local_output_type, blocks_temporary_terms[block][eq], {}, tef_terms);
+ output << " = ";
+ it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+ temporary_terms.insert(it);
+ output << ";" << endl;
+ }
+ };
// The equations
- temporary_terms_idxs_t temporary_terms_idxs;
for (int i = 0; i < block_size; i++)
{
+ write_eq_tt(i);
+
int variable_ID = getBlockVariableID(block, i);
int equation_ID = getBlockEquationID(block, i);
EquationType equ_type = getBlockEquationType(block, i);
@@ -483,6 +490,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
}
}
// The Jacobian if we have to solve the block
+
+ // Write temporary terms for derivatives
+ write_eq_tt(blocks[block].size);
+
if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple
|| simulation_type == BlockSimulationType::solveTwoBoundariesComplete)
output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl;
@@ -995,30 +1006,36 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, bool linear_
//The Temporary terms
deriv_node_temp_terms_t tef_terms;
- if (!linear_decomposition)
- for (auto it : blocks_temporary_terms[block])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
-
- FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
- fnumexpr.write(code_file, instruction_number);
- it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
- FSTPT_ fstpt(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
- fstpt.write(code_file, instruction_number);
- temporary_terms_union.insert(it);
+
+ auto write_eq_tt = [&](int eq)
+ {
+ for (auto it : blocks_temporary_terms[block][eq])
+ {
+ if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+ it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+
+ FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+ fnumexpr.write(code_file, instruction_number);
+ it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, true, false, tef_terms);
+ FSTPT_ fstpt(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+ fstpt.write(code_file, instruction_number);
+ temporary_terms_union.insert(it);
#ifdef DEBUGC
- cout << "FSTPT " << v << endl;
- instruction_number++;
- code_file.write(&FOK, sizeof(FOK));
- code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
- ki++;
+ cout << "FSTPT " << v << endl;
+ instruction_number++;
+ code_file.write(&FOK, sizeof(FOK));
+ code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
+ ki++;
#endif
- }
+ }
+ };
// The equations
for (i = 0; i < block_size; i++)
{
+ if (!linear_decomposition)
+ write_eq_tt(i);
+
int variable_ID, equation_ID;
EquationType equ_type;
@@ -1088,6 +1105,10 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, bool linear_
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
{
+ // Write temporary terms for derivatives
+ if (!linear_decomposition)
+ write_eq_tt(blocks[block].size);
+
switch (simulation_type)
{
case BlockSimulationType::solveBackwardSimple:
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 44004f0e..698b22e0 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -162,8 +162,8 @@ private:
string reform(const string &name) const;
void additionalBlockTemporaryTerms(int blk,
- vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs) const;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 7e70c532..cbcfefcd 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -75,7 +75,7 @@ ExprNode::cost(int cost, bool is_matlab) const
}
int
-ExprNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
+ExprNode::cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const
{
// For a terminal node, the cost is null
return 0;
@@ -161,8 +161,8 @@ ExprNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
}
void
-ExprNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+ExprNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
// Nothing to do for a terminal node
}
@@ -2192,12 +2192,13 @@ UnaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map,
}
int
-UnaryOpNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
+UnaryOpNode::cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const
{
// For a temporary term, the cost is null
- for (const auto &it : blocks_temporary_terms)
- if (it.find(const_cast<UnaryOpNode *>(this)) != it.end())
- return 0;
+ for (const auto &blk_tt : blocks_temporary_terms)
+ for (const auto &eq_tt : blk_tt)
+ if (eq_tt.find(const_cast<UnaryOpNode *>(this)) != eq_tt.end())
+ return 0;
return cost(arg->cost(blocks_temporary_terms, is_matlab), is_matlab);
}
@@ -2332,22 +2333,22 @@ UnaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
}
void
-UnaryOpNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+UnaryOpNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
expr_t this2 = const_cast<UnaryOpNode *>(this);
if (auto it = reference_count.find(this2);
it == reference_count.end())
{
- reference_count[this2] = { 1, blk };
- arg->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ reference_count[this2] = { 1, blk, eq };
+ arg->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
}
else
{
- auto &[nref, first_blk] = it->second;
+ auto &[nref, first_blk, first_eq] = it->second;
nref++;
if (nref * cost(blocks_temporary_terms, false) > min_cost_c)
- blocks_temporary_terms[first_blk].insert(this2);
+ blocks_temporary_terms[first_blk][first_eq].insert(this2);
}
}
@@ -4042,12 +4043,13 @@ BinaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map,
}
int
-BinaryOpNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
+BinaryOpNode::cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const
{
// For a temporary term, the cost is null
- for (const auto &it : blocks_temporary_terms)
- if (it.find(const_cast<BinaryOpNode *>(this)) != it.end())
- return 0;
+ for (const auto &blk_tt : blocks_temporary_terms)
+ for (const auto &eq_tt : blk_tt)
+ if (eq_tt.find(const_cast<BinaryOpNode *>(this)) != eq_tt.end())
+ return 0;
int arg_cost = arg1->cost(blocks_temporary_terms, is_matlab) + arg2->cost(blocks_temporary_terms, is_matlab);
@@ -4144,24 +4146,24 @@ BinaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
}
void
-BinaryOpNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+BinaryOpNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
expr_t this2 = const_cast<BinaryOpNode *>(this);
if (auto it = reference_count.find(this2);
it == reference_count.end())
{
- reference_count[this2] = { 1, blk };
- arg1->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
- arg2->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ reference_count[this2] = { 1, blk, eq };
+ arg1->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
+ arg2->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
}
else
{
- auto &[nref, first_blk] = it->second;
+ auto &[nref, first_blk, first_eq] = it->second;
nref++;
if (nref * cost(blocks_temporary_terms, false) > min_cost_c
&& op_code != BinaryOpcode::equal)
- blocks_temporary_terms[first_blk].insert(this2);
+ blocks_temporary_terms[first_blk][first_eq].insert(this2);
}
}
@@ -5840,12 +5842,13 @@ TrinaryOpNode::cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map
}
int
-TrinaryOpNode::cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const
+TrinaryOpNode::cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const
{
// For a temporary term, the cost is null
- for (const auto &it : blocks_temporary_terms)
- if (it.find(const_cast<TrinaryOpNode *>(this)) != it.end())
- return 0;
+ for (const auto &blk_tt : blocks_temporary_terms)
+ for (const auto &eq_tt : blk_tt)
+ if (eq_tt.find(const_cast<TrinaryOpNode *>(this)) != eq_tt.end())
+ return 0;
int arg_cost = arg1->cost(blocks_temporary_terms, is_matlab)
+ arg2->cost(blocks_temporary_terms, is_matlab)
@@ -5906,24 +5909,24 @@ TrinaryOpNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
}
void
-TrinaryOpNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+TrinaryOpNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
expr_t this2 = const_cast<TrinaryOpNode *>(this);
if (auto it = reference_count.find(this2);
it == reference_count.end())
{
- reference_count[this2] = { 1, blk };
- arg1->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
- arg2->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
- arg3->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ reference_count[this2] = { 1, blk, eq };
+ arg1->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
+ arg2->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
+ arg3->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
}
else
{
- auto &[nref, first_blk] = it->second;
+ auto &[nref, first_blk, first_eq] = it->second;
nref++;
if (nref * cost(blocks_temporary_terms, false) > min_cost_c)
- blocks_temporary_terms[first_blk].insert(this2);
+ blocks_temporary_terms[first_blk][first_eq].insert(this2);
}
}
@@ -6976,20 +6979,21 @@ AbstractExternalFunctionNode::computeTemporaryTerms(const pair<int, int> &derivO
}
void
-AbstractExternalFunctionNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+AbstractExternalFunctionNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
// See comments in computeTemporaryTerms() for the logic
expr_t this2 = const_cast<AbstractExternalFunctionNode *>(this);
- for (auto &tt : blocks_temporary_terms)
- if (auto it = find_if(tt.cbegin(), tt.cend(), sameTefTermPredicate());
- it != tt.cend())
- {
- tt.insert(this2);
- return;
- }
+ for (auto &btt : blocks_temporary_terms)
+ for (auto &tt : btt)
+ if (auto it = find_if(tt.cbegin(), tt.cend(), sameTefTermPredicate());
+ it != tt.cend())
+ {
+ tt.insert(this2);
+ return;
+ }
- blocks_temporary_terms[blk].insert(this2);
+ blocks_temporary_terms[blk][eq].insert(this2);
}
bool
@@ -8236,8 +8240,8 @@ VarExpectationNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
}
void
-VarExpectationNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+VarExpectationNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
cerr << "VarExpectationNode::computeBlocksTemporaryTerms not implemented." << endl;
exit(EXIT_FAILURE);
@@ -8682,10 +8686,10 @@ PacExpectationNode::computeTemporaryTerms(const pair<int, int> &derivOrder,
}
void
-PacExpectationNode::computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+PacExpectationNode::computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
- blocks_temporary_terms[blk].insert(const_cast<PacExpectationNode *>(this));
+ blocks_temporary_terms[blk][eq].insert(const_cast<PacExpectationNode *>(this));
}
expr_t
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 663525de..29080f9b 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -225,7 +225,7 @@ protected:
//! Cost of computing current node
/*! Nodes included in temporary_terms are considered having a null cost */
virtual int cost(int cost, bool is_matlab) const;
- virtual int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const;
+ virtual int cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const;
virtual int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const;
//! For creating equation cross references
@@ -304,17 +304,19 @@ public:
//! Compute temporary terms in this expression for block decomposed model
/*!
\param[in] blk the block number
+ \param[in] eq the equation number (within the block)
\param[out] blocks_temporary_terms the computed temporary terms, per block
+ and per equation in the block
\param[out] reference_count a temporary structure, used to count
references to each node (first integer is the
reference count, second integer is the number of the block in which the
- expression first appears)
+ expression first appears, third integer is the equation number within the block)
Same rules as computeTemporaryTerms() for computing cost.
*/
- virtual void computeBlockTemporaryTerms(int blk,
- vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const;
+ virtual void computeBlockTemporaryTerms(int blk, int eq,
+ vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const;
//! Writes output of node, using a Txxx notation for nodes in temporary_terms, and specifiying the set of already written external functions
/*!
@@ -899,7 +901,7 @@ public:
private:
expr_t computeDerivative(int deriv_id) override;
int cost(int cost, bool is_matlab) const override;
- int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const override;
+ int cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const override;
int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const override;
//! Returns the derivative of this node if darg is the derivative of the argument
expr_t composeDerivatives(expr_t darg, int deriv_id);
@@ -910,8 +912,8 @@ public:
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const override;
- void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ void computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
void writeJsonAST(ostream &output) const override;
void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
@@ -1003,7 +1005,7 @@ public:
private:
expr_t computeDerivative(int deriv_id) override;
int cost(int cost, bool is_matlab) const override;
- int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const override;
+ int cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const override;
int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const override;
//! Returns the derivative of this node if darg1 and darg2 are the derivatives of the arguments
expr_t composeDerivatives(expr_t darg1, expr_t darg2);
@@ -1017,8 +1019,8 @@ public:
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const override;
- void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ void computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
void writeJsonAST(ostream &output) const override;
void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
@@ -1135,7 +1137,7 @@ public:
private:
expr_t computeDerivative(int deriv_id) override;
int cost(int cost, bool is_matlab) const override;
- int cost(const vector<temporary_terms_t> &blocks_temporary_terms, bool is_matlab) const override;
+ int cost(const vector<vector<temporary_terms_t>> &blocks_temporary_terms, bool is_matlab) const override;
int cost(const map<pair<int, int>, temporary_terms_t> &temp_terms_map, bool is_matlab) const override;
//! Returns the derivative of this node if darg1, darg2 and darg3 are the derivatives of the arguments
expr_t composeDerivatives(expr_t darg1, expr_t darg2, expr_t darg3);
@@ -1148,8 +1150,8 @@ public:
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const override;
- void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ void computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
void writeJsonAST(ostream &output) const override;
void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic) const override;
@@ -1260,8 +1262,8 @@ public:
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const override;
- void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ void computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override = 0;
void writeJsonAST(ostream &output) const override = 0;
void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms, bool isdynamic = true) const override = 0;
@@ -1462,8 +1464,8 @@ public:
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const override;
- void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ void computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
expr_t toStatic(DataTree &static_datatree) const override;
expr_t clone(DataTree &datatree) const override;
@@ -1539,8 +1541,8 @@ public:
map<pair<int, int>, temporary_terms_t> &temp_terms_map,
map<expr_t, pair<int, pair<int, int>>> &reference_count,
bool is_matlab) const override;
- void computeBlockTemporaryTerms(int blk, vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const override;
+ void computeBlockTemporaryTerms(int blk, int eq, vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const override;
void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const override;
expr_t toStatic(DataTree &static_datatree) const override;
expr_t clone(DataTree &datatree) const override;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 0aba63f5..777dabad 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -97,8 +97,20 @@ ModelTree::copyHelper(const ModelTree &m)
blocks_derivatives.push_back(v);
}
+ auto convert_vector_tt = [f](vector<temporary_terms_t> vtt)
+ {
+ vector<temporary_terms_t> vtt2;
+ for (const auto &tt : vtt)
+ {
+ temporary_terms_t tt2;
+ for (const auto &it : tt)
+ tt2.insert(f(it));
+ vtt2.push_back(tt2);
+ }
+ return vtt2;
+ };
for (const auto &it : m.blocks_temporary_terms)
- blocks_temporary_terms.push_back(convert_temporary_terms_t(it));
+ blocks_temporary_terms.push_back(convert_vector_tt(it));
for (const auto &it : m.blocks_temporary_terms_idxs)
blocks_temporary_terms_idxs[f(it.first)] = it.second;
}
@@ -1056,18 +1068,19 @@ ModelTree::computeBlockTemporaryTerms()
int nb_blocks = blocks.size();
blocks_temporary_terms.resize(nb_blocks);
- map<expr_t, pair<int, int>> reference_count;
+ map<expr_t, tuple<int, int, int>> reference_count;
for (int blk = 0; blk < nb_blocks; blk++)
{
+ blocks_temporary_terms[blk].resize(blocks[blk].size + 1);
for (int eq = 0; eq < blocks[blk].size; eq++)
{
if (eq < blocks[blk].getRecursiveSize() && isBlockEquationRenormalized(blk, eq))
- getBlockEquationRenormalizedExpr(blk, eq)->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ getBlockEquationRenormalizedExpr(blk, eq)->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
else
- getBlockEquationExpr(blk, eq)->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ getBlockEquationExpr(blk, eq)->computeBlockTemporaryTerms(blk, eq, blocks_temporary_terms, reference_count);
}
for (const auto &[ignore, d] : blocks_derivatives[blk])
- d->computeBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
+ d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
additionalBlockTemporaryTerms(blk, blocks_temporary_terms, reference_count);
}
@@ -1075,15 +1088,16 @@ ModelTree::computeBlockTemporaryTerms()
// Compute indices in the temporary terms vector
int idx = 0;
blocks_temporary_terms_idxs.clear();
- for (int blk = 0; blk < nb_blocks; blk++)
- for (auto tt : blocks_temporary_terms[blk])
- blocks_temporary_terms_idxs[tt] = idx++;
+ for (auto &blk_tt : blocks_temporary_terms)
+ for (auto &eq_tt : blk_tt)
+ for (auto tt : eq_tt)
+ blocks_temporary_terms_idxs[tt] = idx++;
}
void
ModelTree::additionalBlockTemporaryTerms(int blk,
- vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const
+ vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const
{
}
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 0cd0db8b..f63f91c2 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -193,8 +193,18 @@ protected:
It verifies: ∀i, eq2block[endo2eq[i]] = endo2block[i] */
vector<int> eq2block;
- //! Temporary terms for block decomposed models (one block per element of the vector)
- vector<temporary_terms_t> blocks_temporary_terms;
+ /* 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;
/* Stores, for each temporary term in block decomposed models, its index in
the vector of all temporary terms */
@@ -223,8 +233,8 @@ protected:
be overriden by subclasses (actually by DynamicModel, who needs extra
temporary terms for derivatives w.r.t. exogenous and other endogenous) */
virtual void additionalBlockTemporaryTerms(int blk,
- vector<temporary_terms_t> &blocks_temporary_terms,
- map<expr_t, pair<int, int>> &reference_count) const;
+ vector<vector<temporary_terms_t>> &blocks_temporary_terms,
+ map<expr_t, tuple<int, int, int>> &reference_count) const;
//! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms();
//! Writes temporary terms
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 1f5ef4fa..5ac211b6 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -173,16 +173,17 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
// Declare global temp terms from this block and the previous ones
bool global_keyword_written = false;
for (int blk2 = 0; blk2 <= block; blk2++)
- for (auto tt : blocks_temporary_terms[blk2])
- {
- if (!global_keyword_written)
- {
- output << " global";
- global_keyword_written = true;
- }
- output << " ";
- tt->writeOutput(output, local_output_type, blocks_temporary_terms[blk2], {});
- }
+ for (auto &eq_tt : blocks_temporary_terms[blk2])
+ for (auto tt : eq_tt)
+ {
+ if (!global_keyword_written)
+ {
+ output << " global";
+ global_keyword_written = true;
+ }
+ output << " ";
+ tt->writeOutput(output, local_output_type, eq_tt, {});
+ }
if (global_keyword_written)
output << ";" << endl;
@@ -191,23 +192,28 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
output << " residual=zeros(" << block_mfs << ",1);" << endl;
// The equations
-
deriv_node_temp_terms_t tef_terms;
- for (auto it : blocks_temporary_terms[block])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, {}, tef_terms);
-
- output << " ";
- it->writeOutput(output, local_output_type, blocks_temporary_terms[block], {}, tef_terms);
- output << " = ";
- it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
- temporary_terms.insert(it);
- output << ";" << endl;
- }
+
+ auto write_eq_tt = [&](int eq)
+ {
+ for (auto it : blocks_temporary_terms[block][eq])
+ {
+ if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+ it->writeExternalFunctionOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+
+ output << " ";
+ it->writeOutput(output, local_output_type, blocks_temporary_terms[block][eq], {}, tef_terms);
+ output << " = ";
+ it->writeOutput(output, local_output_type, temporary_terms, {}, tef_terms);
+ temporary_terms.insert(it);
+ output << ";" << endl;
+ }
+ };
for (int i = 0; i < block_size; i++)
{
+ write_eq_tt(i);
+
int variable_ID = getBlockVariableID(block, i);
int equation_ID = getBlockEquationID(block, i);
EquationType equ_type = getBlockEquationType(block, i);
@@ -266,6 +272,9 @@ StaticModel::writeModelEquationsOrdered_M(const string &basename) const
case BlockSimulationType::solveForwardSimple:
case BlockSimulationType::solveBackwardComplete:
case BlockSimulationType::solveForwardComplete:
+ // Write temporary terms for derivatives
+ write_eq_tt(blocks[block].size);
+
for (const auto &[indices, id] : blocks_derivatives[block])
{
auto [eq, var, ignore] = indices;
@@ -528,24 +537,30 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
//The Temporary terms
deriv_node_temp_terms_t tef_terms;
- /* We are force to use a copy of tt_union here, since temp. terms are
+ /* Keep a backup of temporary_terms_union here, since temp. terms are
written a second time below. This is probably unwanted… */
- temporary_terms_t tt2 = temporary_terms_union;
- for (auto it : blocks_temporary_terms[block])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, blocks_temporary_terms_idxs, false, false, tef_terms);
-
- FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
- fnumexpr.write(code_file, instruction_number);
- it->compile(code_file, instruction_number, false, tt2, blocks_temporary_terms_idxs, false, false, tef_terms);
- FSTPST_ fstpst(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
- fstpst.write(code_file, instruction_number);
- tt2.insert(it);
- }
+ temporary_terms_t ttu_old = temporary_terms_union;
+
+ auto write_eq_tt = [&](int eq)
+ {
+ for (auto it : blocks_temporary_terms[block][eq])
+ {
+ if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+ it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+
+ FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+ fnumexpr.write(code_file, instruction_number);
+ it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
+ FSTPST_ fstpst(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
+ fstpst.write(code_file, instruction_number);
+ temporary_terms_union.insert(it);
+ }
+ };
for (i = 0; i < block_size; i++)
{
+ write_eq_tt(i);
+
// The equations
int variable_ID, equation_ID;
EquationType equ_type;
@@ -564,16 +579,16 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
eq_node = getBlockEquationExpr(block, i);
lhs = eq_node->arg1;
rhs = eq_node->arg2;
- rhs->compile(code_file, instruction_number, false, tt2, blocks_temporary_terms_idxs, false, false);
- lhs->compile(code_file, instruction_number, true, tt2, blocks_temporary_terms_idxs, false, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false);
+ lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false);
}
else if (equ_type == EquationType::evaluate_s)
{
eq_node = getBlockEquationRenormalizedExpr(block, i);
lhs = eq_node->arg1;
rhs = eq_node->arg2;
- rhs->compile(code_file, instruction_number, false, tt2, blocks_temporary_terms_idxs, false, false);
- lhs->compile(code_file, instruction_number, true, tt2, blocks_temporary_terms_idxs, false, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false);
+ lhs->compile(code_file, instruction_number, true, temporary_terms_union, blocks_temporary_terms_idxs, false, false);
}
break;
case BlockSimulationType::solveBackwardComplete:
@@ -592,8 +607,8 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
eq_node = getBlockEquationExpr(block, i);
lhs = eq_node->arg1;
rhs = eq_node->arg2;
- lhs->compile(code_file, instruction_number, false, tt2, blocks_temporary_terms_idxs, false, false);
- rhs->compile(code_file, instruction_number, false, tt2, blocks_temporary_terms_idxs, false, false);
+ lhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false);
+ rhs->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false);
FBINARY_ fbinary{static_cast<int>(BinaryOpcode::minus)};
fbinary.write(code_file, instruction_number);
@@ -609,6 +624,9 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
{
+ // Write temporary terms for derivatives
+ write_eq_tt(blocks[block].size);
+
switch (simulation_type)
{
case BlockSimulationType::solveBackwardSimple:
@@ -617,7 +635,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
fnumexpr.write(code_file, instruction_number);
}
- compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), tt2, blocks_temporary_terms_idxs);
+ compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs);
{
FSTPG_ fstpg(0);
fstpg.write(code_file, instruction_number);
@@ -649,7 +667,7 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
Uf[eqr].Ufl->var = varr;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
fnumexpr.write(code_file, instruction_number);
- compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, tt2, blocks_temporary_terms_idxs);
+ compileChainRuleDerivative(code_file, instruction_number, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs);
FSTPSU_ fstpsu(count_u);
fstpsu.write(code_file, instruction_number);
count_u++;
@@ -713,21 +731,12 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
prev_instruction_number = instruction_number;
tef_terms.clear();
- for (auto it : blocks_temporary_terms[block])
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->compileExternalFunctionOutput(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
-
- FNUMEXPR_ fnumexpr(TemporaryTerm, static_cast<int>(blocks_temporary_terms_idxs.at(it)));
- fnumexpr.write(code_file, instruction_number);
- it->compile(code_file, instruction_number, false, temporary_terms_union, blocks_temporary_terms_idxs, false, false, tef_terms);
- FSTPST_ fstpst(static_cast<int>(blocks_temporary_terms_idxs.at(it)));
- fstpst.write(code_file, instruction_number);
- temporary_terms_union.insert(it);
- }
+ temporary_terms_union = ttu_old;
for (i = 0; i < block_size; i++)
{
+ write_eq_tt(i);
+
// The equations
int variable_ID, equation_ID;
EquationType equ_type;
@@ -788,6 +797,10 @@ StaticModel::writeModelEquationsCode_Block(const string &basename) const
fendequ_l.write(code_file, instruction_number);
// The Jacobian if we have to solve the block determinsitic bloc
+
+ // Write temporary terms for derivatives
+ write_eq_tt(blocks[block].size);
+
switch (simulation_type)
{
case BlockSimulationType::solveBackwardSimple:
--
GitLab