diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 30fceaca6dbfdf3b47ec435b62cf33c4efa541b5..ea49debcdb563c084f68535297b5209fd028da3c 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 44004f0e42e92290903114f02df6bd3f8e3c9f1d..698b22e0b5144d7536501d9e05f7f72245ec6a3d 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 7e70c5320a795a8ce28ed84d3e442cf7674c0b0a..cbcfefcd224c3a4cc70729074c1dba1b897ab3b1 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 663525de482fcc58a4f9315879ba992af1323c29..29080f9b2f1e930f0bba17e93d01f6b101dd7deb 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 0aba63f5f08b4fe586a12595366259c62433f81d..777dabadb7d5edc41dc8cbeb9bf462bb28b3d9f0 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 0cd0db8bba44841a1a5a1237a63f3364a06934c1..f63f91c249b058dcefcee94a375941b9af54790e 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 1f5ef4fa92ab8c001200fbdcee01cfbf6cc87717..5ac211b62f548136be76ffaa083f9cecc9259e88 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: