diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index e22bd2cc4c106e8f84fd4302b77e914cfa11388f..9a10ad93fb7cc2b3e7de2aacf4f0e890db276ff5 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -198,238 +198,6 @@ DynamicModel::additionalBlockTemporaryTerms(int blk,
d->computeBlockTemporaryTerms(blk, blocks[blk].size, blocks_temporary_terms, reference_count);
}
-void
-DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const
-{
- BlockSimulationType simulation_type = blocks[blk].simulation_type;
- int block_size = blocks[blk].size;
- int block_mfs_size = blocks[blk].mfs_size;
- int block_recursive_size = blocks[blk].getRecursiveSize();
-
- deriv_node_temp_terms_t tef_terms;
-
- 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;
- }
- };
-
- // The equations
- for (int eq = 0; eq < block_size; eq++)
- {
- write_eq_tt(eq);
-
- EquationType equ_type = getBlockEquationType(blk, eq);
- BinaryOpNode *e = getBlockEquationExpr(blk, eq);
- expr_t lhs = e->arg1, rhs = e->arg2;
- switch (simulation_type)
- {
- case BlockSimulationType::evaluateBackward:
- case BlockSimulationType::evaluateForward:
- evaluation:
- if (equ_type == EquationType::evaluateRenormalized)
- {
- e = getBlockEquationRenormalizedExpr(blk, eq);
- lhs = e->arg1;
- rhs = e->arg2;
- }
- else if (equ_type != EquationType::evaluate)
- {
- cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
- exit(EXIT_FAILURE);
- }
- 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::solveBackwardSimple:
- case BlockSimulationType::solveForwardSimple:
- case BlockSimulationType::solveBackwardComplete:
- case BlockSimulationType::solveForwardComplete:
- case BlockSimulationType::solveTwoBoundariesComplete:
- case BlockSimulationType::solveTwoBoundariesSimple:
- if (eq < block_recursive_size)
- goto evaluation;
- output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
- goto end;
- default:
- end:
- 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;
- }
- }
-
- // The Jacobian if we have to solve the block
-
- // Write temporary terms for derivatives
- write_eq_tt(blocks[blk].size);
-
- if (isCOutput(output_type))
- output << " if (stochastic_mode) {" << endl;
- else
- output << " if stochastic_mode" << endl;
-
- ostringstream i_output, j_output, v_output;
- int line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
- for (const auto &[indices, d] : blocks_derivatives[blk])
- {
- auto [eq, var, lag] = indices;
- int jacob_col = blocks_jacob_cols_endo[blk].at({ var, lag });
- i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
- j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
- v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- assert(line_counter == nze_stochastic+ARRAY_SUBSCRIPT_OFFSET(output_type));
- output << i_output.str() << j_output.str() << v_output.str();
-
- i_output.str("");
- j_output.str("");
- v_output.str("");
- line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
- for (const auto &[indices, d] : blocks_derivatives_exo[blk])
- {
- auto [eq, var, lag] = indices;
- int jacob_col = blocks_jacob_cols_exo[blk].at({ var, lag });
- i_output << " g1_x_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
- j_output << " g1_x_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
- v_output << " g1_x_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- assert(line_counter == nze_exo+ARRAY_SUBSCRIPT_OFFSET(output_type));
- output << i_output.str() << j_output.str() << v_output.str();
-
- i_output.str("");
- j_output.str("");
- v_output.str("");
- line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
- for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
- {
- auto [eq, var, lag] = indices;
- int jacob_col = blocks_jacob_cols_exo_det[blk].at({ var, lag });
- i_output << " g1_xd_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
- j_output << " g1_xd_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
- v_output << " g1_xd_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- assert(line_counter == nze_exo_det+ARRAY_SUBSCRIPT_OFFSET(output_type));
- output << i_output.str() << j_output.str() << v_output.str();
-
- i_output.str("");
- j_output.str("");
- v_output.str("");
- line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
- for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
- {
- auto [eq, var, lag] = indices;
- int jacob_col = blocks_jacob_cols_other_endo[blk].at({ var, lag });
- i_output << " g1_o_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
- j_output << " g1_o_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
- v_output << " g1_o_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- assert(line_counter == nze_other_endo+ARRAY_SUBSCRIPT_OFFSET(output_type));
- output << i_output.str() << j_output.str() << v_output.str();
-
- // Deterministic mode
- if (simulation_type != BlockSimulationType::evaluateForward
- && simulation_type != BlockSimulationType::evaluateBackward)
- {
- if (isCOutput(output_type))
- output << " } else {" << endl;
- else
- output << " else" << endl;
- i_output.str("");
- j_output.str("");
- v_output.str("");
- line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
- if (simulation_type == BlockSimulationType::solveBackwardSimple
- || simulation_type == BlockSimulationType::solveForwardSimple
- || simulation_type == BlockSimulationType::solveBackwardComplete
- || simulation_type == BlockSimulationType::solveForwardComplete)
- for (const auto &[indices, d] : blocks_derivatives[blk])
- {
- auto [eq, var, lag] = indices;
- if (lag == 0 && eq >= block_recursive_size && var >= block_recursive_size)
- {
- i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
- << eq+1-block_recursive_size << ';' << endl;
- j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
- << var+1-block_recursive_size << ';' << endl;
- v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- }
- else // solveTwoBoundariesSimple || solveTwoBoundariesComplete
- for (const auto &[indices, d] : blocks_derivatives[blk])
- {
- auto [eq, var, lag] = indices;
- assert(lag >= -1 && lag <= 1);
- if (eq >= block_recursive_size && var >= block_recursive_size)
- {
- i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
- << eq+1-block_recursive_size << ';' << endl;
- j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
- << var+1-block_recursive_size+block_mfs_size*(lag+1) << ';' << endl;
- v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- }
- assert(line_counter == nze_deterministic+ARRAY_SUBSCRIPT_OFFSET(output_type));
- output << i_output.str() << j_output.str() << v_output.str();
- }
- if (isCOutput(output_type))
- output << " }" << endl;
- else
- output << " end" << endl;
-}
-
int
DynamicModel::nzeDeterministicJacobianForBlock(int blk) const
{
@@ -527,7 +295,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const
output << " end" << endl
<< endl;
- writeDynamicPerBlockHelper(blk, output, ExprNodeOutputType::matlabDynamicModel, temporary_terms,
+ writeDynamicPerBlockHelper<ExprNodeOutputType::matlabDynamicModel>(blk, output, temporary_terms,
nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);
output << endl
@@ -611,7 +379,7 @@ DynamicModel::writeDynamicPerBlockCFiles(const string &basename) const
output << "void dynamic_" << blk+1 << "(double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, double *restrict T, int it_, bool stochastic_mode, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v, double *restrict g1_x_i, double *restrict g1_x_j, double *restrict g1_x_v, double *restrict g1_xd_i, double *restrict g1_xd_j, double *restrict g1_xd_v, double *restrict g1_o_i, double *restrict g1_o_j, double *restrict g1_o_v)" << endl;
output << '{' << endl;
- writeDynamicPerBlockHelper(blk, output, ExprNodeOutputType::CDynamicModel, temporary_terms,
+ writeDynamicPerBlockHelper<ExprNodeOutputType::CDynamicModel>(blk, output, temporary_terms,
nze_stochastic, nze_deterministic, nze_exo, nze_exo_det, nze_other_endo);
output << '}' << endl
@@ -5369,18 +5137,18 @@ DynamicModel::isChecksumMatching(const string &basename) const
// Write equation tags
equation_tags.writeCheckSumInfo(buffer);
- ExprNodeOutputType buffer_type = ExprNodeOutputType::CDynamicModel;
+ constexpr ExprNodeOutputType buffer_type{ExprNodeOutputType::CDynamicModel};
deriv_node_temp_terms_t tef_terms;
temporary_terms_t temp_term_union;
- writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs,
- buffer, buffer_type, tef_terms);
+ writeModelLocalVariableTemporaryTerms<buffer_type>(temp_term_union, temporary_terms_idxs,
+ buffer, tef_terms);
- writeTemporaryTerms(temporary_terms_derivatives[0],
- temp_term_union, temporary_terms_idxs,
- buffer, buffer_type, tef_terms);
+ writeTemporaryTerms<buffer_type>(temporary_terms_derivatives[0],
+ temp_term_union, temporary_terms_idxs,
+ buffer, tef_terms);
- writeModelEquations(buffer, buffer_type, temp_term_union);
+ writeModelEquations<buffer_type>(buffer, temp_term_union);
size_t result = hash<string>{}(buffer.str());
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 15adc46bc8517d6e138406608bb299be5a7eee53..ce051ec493024cdead844c7d2f30e276a292265f 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -132,7 +132,8 @@ private:
block-decomposed model */
int nzeDeterministicJacobianForBlock(int blk) const;
//! Helper for writing the per-block dynamic files of block decomposed models
- void writeDynamicPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const;
+ template<ExprNodeOutputType output_type>
+ void writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const;
//! Writes the per-block dynamic files of block decomposed model (MATLAB version)
void writeDynamicPerBlockMFiles(const string &basename) const;
//! Writes the per-block dynamic files of block decomposed model (C version)
@@ -659,6 +660,241 @@ public:
set<int> findPacExpectationEquationNumbers() const;
};
+template<ExprNodeOutputType output_type>
+void
+DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms,
+ int nze_stochastic, int nze_deterministic, int nze_exo,
+ int nze_exo_det, int nze_other_endo) const
+{
+ BlockSimulationType simulation_type { blocks[blk].simulation_type };
+ int block_size { blocks[blk].size };
+ int block_mfs_size { blocks[blk].mfs_size };
+ int block_recursive_size { blocks[blk].getRecursiveSize() };
+
+ deriv_node_temp_terms_t tef_terms;
+
+ 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;
+ }
+ };
+
+ // The equations
+ for (int eq {0}; eq < block_size; eq++)
+ {
+ write_eq_tt(eq);
+
+ EquationType equ_type { getBlockEquationType(blk, eq) };
+ BinaryOpNode *e { getBlockEquationExpr(blk, eq) };
+ expr_t lhs { e->arg1 }, rhs { e->arg2 };
+ switch (simulation_type)
+ {
+ case BlockSimulationType::evaluateBackward:
+ case BlockSimulationType::evaluateForward:
+ evaluation:
+ if (equ_type == EquationType::evaluateRenormalized)
+ {
+ e = getBlockEquationRenormalizedExpr(blk, eq);
+ lhs = e->arg1;
+ rhs = e->arg2;
+ }
+ else if (equ_type != EquationType::evaluate)
+ {
+ cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
+ exit(EXIT_FAILURE);
+ }
+ 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::solveBackwardSimple:
+ case BlockSimulationType::solveForwardSimple:
+ case BlockSimulationType::solveBackwardComplete:
+ case BlockSimulationType::solveForwardComplete:
+ case BlockSimulationType::solveTwoBoundariesComplete:
+ case BlockSimulationType::solveTwoBoundariesSimple:
+ if (eq < block_recursive_size)
+ goto evaluation;
+ output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=(";
+ goto end;
+ default:
+ end:
+ 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;
+ }
+ }
+
+ // The Jacobian if we have to solve the block
+
+ // Write temporary terms for derivatives
+ write_eq_tt(blocks[blk].size);
+
+ if constexpr(isCOutput(output_type))
+ output << " if (stochastic_mode) {" << endl;
+ else
+ output << " if stochastic_mode" << endl;
+
+ ostringstream i_output, j_output, v_output;
+ int line_counter { ARRAY_SUBSCRIPT_OFFSET(output_type) };
+ for (const auto &[indices, d] : blocks_derivatives[blk])
+ {
+ const auto &[eq, var, lag] {indices};
+ int jacob_col { blocks_jacob_cols_endo[blk].at({ var, lag }) };
+ i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
+ j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
+ v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ assert(line_counter == nze_stochastic+ARRAY_SUBSCRIPT_OFFSET(output_type));
+ output << i_output.str() << j_output.str() << v_output.str();
+
+ i_output.str("");
+ j_output.str("");
+ v_output.str("");
+ line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
+ for (const auto &[indices, d] : blocks_derivatives_exo[blk])
+ {
+ const auto &[eq, var, lag] {indices};
+ int jacob_col { blocks_jacob_cols_exo[blk].at({ var, lag }) };
+ i_output << " g1_x_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
+ j_output << " g1_x_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
+ v_output << " g1_x_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ assert(line_counter == nze_exo+ARRAY_SUBSCRIPT_OFFSET(output_type));
+ output << i_output.str() << j_output.str() << v_output.str();
+
+ i_output.str("");
+ j_output.str("");
+ v_output.str("");
+ line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
+ for (const auto &[indices, d] : blocks_derivatives_exo_det[blk])
+ {
+ const auto &[eq, var, lag] {indices};
+ int jacob_col { blocks_jacob_cols_exo_det[blk].at({ var, lag }) };
+ i_output << " g1_xd_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
+ j_output << " g1_xd_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
+ v_output << " g1_xd_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ assert(line_counter == nze_exo_det+ARRAY_SUBSCRIPT_OFFSET(output_type));
+ output << i_output.str() << j_output.str() << v_output.str();
+
+ i_output.str("");
+ j_output.str("");
+ v_output.str("");
+ line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
+ for (const auto &[indices, d] : blocks_derivatives_other_endo[blk])
+ {
+ const auto &[eq, var, lag] {indices};
+ int jacob_col { blocks_jacob_cols_other_endo[blk].at({ var, lag }) };
+ i_output << " g1_o_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1 << ';' << endl;
+ j_output << " g1_o_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << jacob_col+1 << ';' << endl;
+ v_output << " g1_o_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ assert(line_counter == nze_other_endo+ARRAY_SUBSCRIPT_OFFSET(output_type));
+ output << i_output.str() << j_output.str() << v_output.str();
+
+ // Deterministic mode
+ if (simulation_type != BlockSimulationType::evaluateForward
+ && simulation_type != BlockSimulationType::evaluateBackward)
+ {
+ if constexpr(isCOutput(output_type))
+ output << " } else {" << endl;
+ else
+ output << " else" << endl;
+ i_output.str("");
+ j_output.str("");
+ v_output.str("");
+ line_counter = ARRAY_SUBSCRIPT_OFFSET(output_type);
+ if (simulation_type == BlockSimulationType::solveBackwardSimple
+ || simulation_type == BlockSimulationType::solveForwardSimple
+ || simulation_type == BlockSimulationType::solveBackwardComplete
+ || simulation_type == BlockSimulationType::solveForwardComplete)
+ for (const auto &[indices, d] : blocks_derivatives[blk])
+ {
+ const auto &[eq, var, lag] {indices};
+ if (lag == 0 && eq >= block_recursive_size && var >= block_recursive_size)
+ {
+ i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
+ << eq+1-block_recursive_size << ';' << endl;
+ j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
+ << var+1-block_recursive_size << ';' << endl;
+ v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ }
+ else // solveTwoBoundariesSimple || solveTwoBoundariesComplete
+ for (const auto &[indices, d] : blocks_derivatives[blk])
+ {
+ const auto &[eq, var, lag] {indices};
+ assert(lag >= -1 && lag <= 1);
+ if (eq >= block_recursive_size && var >= block_recursive_size)
+ {
+ i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
+ << eq+1-block_recursive_size << ';' << endl;
+ j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '='
+ << var+1-block_recursive_size+block_mfs_size*(lag+1) << ';' << endl;
+ v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ }
+ assert(line_counter == nze_deterministic+ARRAY_SUBSCRIPT_OFFSET(output_type));
+ output << i_output.str() << j_output.str() << v_output.str();
+ }
+ if constexpr(isCOutput(output_type))
+ output << " }" << endl;
+ else
+ output << " end" << endl;
+}
+
template<bool julia>
void
DynamicModel::writeParamsDerivativesFile(const string &basename) const
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 43483486c9403301d213f62e620abb4495b603e7..f96efb26429b8c67eaa27094e5fc4be17fc90a03 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1019,61 +1019,6 @@ ModelTree::additionalBlockTemporaryTerms([[maybe_unused]] int blk,
{
}
-void
-ModelTree::writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_union,
- const temporary_terms_idxs_t &tt_idxs,
- ostream &output, ExprNodeOutputType output_type,
- deriv_node_temp_terms_t &tef_terms) const
-{
- temporary_terms_t tto;
- for (auto [mlv, value] : temporary_terms_mlv)
- tto.insert(mlv);
-
- for (auto [mlv, value] : temporary_terms_mlv)
- {
- value->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
-
- if (isJuliaOutput(output_type))
- output << " const ";
-
- mlv->writeOutput(output, output_type, tto, tt_idxs, tef_terms);
- output << " = ";
- value->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
-
- if (isCOutput(output_type) || isMatlabOutput(output_type))
- output << ";";
- output << endl;
-
- /* We put in temp_term_union the VariableNode corresponding to the MLV,
- not its definition, so that when equations use the MLV,
- T(XXX) is printed instead of the MLV name */
- temp_term_union.insert(mlv);
- }
-}
-
-void
-ModelTree::writeTemporaryTerms(const temporary_terms_t &tt,
- temporary_terms_t &temp_term_union,
- const temporary_terms_idxs_t &tt_idxs,
- ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const
-{
- for (auto it : tt)
- {
- if (dynamic_cast<AbstractExternalFunctionNode *>(it))
- it->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
-
- it->writeOutput(output, output_type, tt, tt_idxs, tef_terms);
- output << " = ";
- it->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
-
- if (isCOutput(output_type) || isMatlabOutput(output_type))
- output << ";";
- output << endl;
-
- temp_term_union.insert(it);
- }
-}
-
void
ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt,
temporary_terms_t &temp_term_union,
@@ -1311,62 +1256,6 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, d
output << "]";
}
-void
-ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
- const temporary_terms_t &temporary_terms) const
-{
- for (int eq = 0; eq < static_cast<int>(equations.size()); eq++)
- {
- BinaryOpNode *eq_node = equations[eq];
- expr_t lhs = eq_node->arg1, rhs = eq_node->arg2;
-
- // Test if the right hand side of the equation is empty.
- double vrhs = 1.0;
- try
- {
- vrhs = rhs->eval({});
- }
- catch (ExprNode::EvalException &e)
- {
- }
-
- if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
- if (isJuliaOutput(output_type))
- {
- output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
- << RIGHT_ARRAY_SUBSCRIPT(output_type)
- << " = (";
- lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
- output << ") - (";
- rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
- output << ")" << endl;
- }
- else
- {
- output << "lhs = ";
- lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
- output << ";" << endl
- << "rhs = ";
- rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
- output << ";" << endl
- << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
- << RIGHT_ARRAY_SUBSCRIPT(output_type)
- << " = lhs - rhs;" << endl;
- }
- else // The right hand side of the equation is empty ==> residual=lhs;
- {
- output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
- << RIGHT_ARRAY_SUBSCRIPT(output_type)
- << " = ";
- lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
- output << ";" << endl;
- }
- }
-}
-
void
ModelTree::writeBytecodeModelEquations(BytecodeWriter &code_file, ExprNodeBytecodeOutputType output_type, const temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const
{
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index 9ca733c4048ce5013dfcddd2b96f204b430cae35..7061eb6d5f49a8745f41df9d9a4c532d047e4f07 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -238,7 +238,8 @@ protected:
//! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms();
//! Writes temporary terms
- void writeTemporaryTerms(const temporary_terms_t &tt, temporary_terms_t &temp_term_union, const temporary_terms_idxs_t &tt_idxs, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const;
+ template<ExprNodeOutputType output_type>
+ 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;
//! Writes temporary terms in bytecode
void writeBytecodeTemporaryTerms(BytecodeWriter &code_file, ExprNodeBytecodeOutputType output_type, temporary_terms_t &temporary_terms_union, const temporary_terms_idxs_t &temporary_terms_idxs, deriv_node_temp_terms_t &tef_terms) const;
@@ -248,13 +249,14 @@ protected:
void fixNestedParenthesis(ostringstream &output, map<string, string> &tmp_paren_vars, bool &message_printed) const;
//! Tests if string contains more than 32 nested parens, Issue #1201
bool testNestedParenthesis(const string &str) const;
+
+ template<ExprNodeOutputType output_type>
void writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_union,
const temporary_terms_idxs_t &tt_idxs,
- ostream &output, ExprNodeOutputType output_type,
- deriv_node_temp_terms_t &tef_terms) const;
+ ostream &output, deriv_node_temp_terms_t &tef_terms) const;
//! Writes model equations
- void writeModelEquations(ostream &output, ExprNodeOutputType output_type,
- const temporary_terms_t &temporary_terms) const;
+ template<ExprNodeOutputType output_type>
+ void writeModelEquations(ostream &output, const temporary_terms_t &temporary_terms) const;
// Returns outputs for derivatives and temporary terms at each derivation order
template<ExprNodeOutputType output_type>
@@ -518,6 +520,118 @@ public:
}
};
+template<ExprNodeOutputType output_type>
+void
+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
+{
+ for (auto it : tt)
+ {
+ if (dynamic_cast<AbstractExternalFunctionNode *>(it))
+ it->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
+
+ it->writeOutput(output, output_type, tt, tt_idxs, tef_terms);
+ output << " = ";
+ it->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
+
+ if constexpr(isCOutput(output_type) || isMatlabOutput(output_type))
+ output << ";";
+ output << endl;
+
+ temp_term_union.insert(it);
+ }
+}
+
+template<ExprNodeOutputType output_type>
+void
+ModelTree::writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_union,
+ const temporary_terms_idxs_t &tt_idxs,
+ ostream &output, deriv_node_temp_terms_t &tef_terms) const
+{
+ temporary_terms_t tto;
+ for (const auto &[mlv, value] : temporary_terms_mlv)
+ tto.insert(mlv);
+
+ for (const auto &[mlv, value] : temporary_terms_mlv)
+ {
+ value->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
+
+ if constexpr(isJuliaOutput(output_type))
+ output << " const ";
+
+ mlv->writeOutput(output, output_type, tto, tt_idxs, tef_terms);
+ output << " = ";
+ value->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
+
+ if constexpr(isCOutput(output_type) || isMatlabOutput(output_type))
+ output << ";";
+ output << endl;
+
+ /* We put in temp_term_union the VariableNode corresponding to the MLV,
+ not its definition, so that when equations use the MLV,
+ T(XXX) is printed instead of the MLV name */
+ temp_term_union.insert(mlv);
+ }
+}
+
+template<ExprNodeOutputType output_type>
+void
+ModelTree::writeModelEquations(ostream &output, const temporary_terms_t &temporary_terms) const
+{
+ for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
+ {
+ BinaryOpNode *eq_node { equations[eq] };
+ expr_t lhs { eq_node->arg1 }, rhs { eq_node->arg2 };
+
+ // Test if the right hand side of the equation is empty.
+ double vrhs {1.0};
+ try
+ {
+ vrhs = rhs->eval({});
+ }
+ catch (ExprNode::EvalException &e)
+ {
+ }
+
+ if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
+ if constexpr(isJuliaOutput(output_type))
+ {
+ output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << " = (";
+ lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+ output << ") - (";
+ rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+ output << ")" << endl;
+ }
+ else
+ {
+ output << "lhs = ";
+ lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+ output << ";" << endl
+ << "rhs = ";
+ rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+ output << ";" << endl
+ << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << " = lhs - rhs;" << endl;
+ }
+ else // The right hand side of the equation is empty ==> residual=lhs;
+ {
+ output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << RIGHT_ARRAY_SUBSCRIPT(output_type)
+ << " = ";
+ lhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
+ output << ";" << endl;
+ }
+ }
+}
+
template<ExprNodeOutputType output_type>
pair<vector<ostringstream>, vector<ostringstream>>
ModelTree::writeModelFileHelper() const
@@ -528,23 +642,19 @@ ModelTree::writeModelFileHelper() const
deriv_node_temp_terms_t tef_terms;
temporary_terms_t temp_term_union;
- writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs,
- tt_output[0], output_type, tef_terms);
+ writeModelLocalVariableTemporaryTerms<output_type>(temp_term_union, temporary_terms_idxs,
+ tt_output[0], tef_terms);
- writeTemporaryTerms(temporary_terms_derivatives[0],
- temp_term_union,
- temporary_terms_idxs,
- tt_output[0], output_type, tef_terms);
+ writeTemporaryTerms<output_type>(temporary_terms_derivatives[0], temp_term_union,
+ temporary_terms_idxs, tt_output[0], tef_terms);
- writeModelEquations(d_output[0], output_type, temp_term_union);
+ writeModelEquations<output_type>(d_output[0], temp_term_union);
// Writing Jacobian
if (!derivatives[1].empty())
{
- writeTemporaryTerms(temporary_terms_derivatives[1],
- temp_term_union,
- temporary_terms_idxs,
- tt_output[1], output_type, tef_terms);
+ writeTemporaryTerms<output_type>(temporary_terms_derivatives[1], temp_term_union,
+ temporary_terms_idxs, tt_output[1], tef_terms);
for (const auto &[indices, d1] : derivatives[1])
{
@@ -566,10 +676,8 @@ ModelTree::writeModelFileHelper() const
for (size_t i = 2; i < derivatives.size(); i++)
if (!derivatives[i].empty())
{
- writeTemporaryTerms(temporary_terms_derivatives[i],
- temp_term_union,
- temporary_terms_idxs,
- tt_output[i], output_type, tef_terms);
+ writeTemporaryTerms<output_type>(temporary_terms_derivatives[i], temp_term_union,
+ temporary_terms_idxs, tt_output[i], tef_terms);
/* When creating the sparse matrix (in MATLAB or C mode), since storage
is in column-major order, output the first column, then the second,
@@ -680,9 +788,12 @@ ModelTree::writeParamsDerivativesFileHelper() const
temporary_terms_t temp_term_union;
deriv_node_temp_terms_t tef_terms;
- writeModelLocalVariableTemporaryTerms(temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
+ writeModelLocalVariableTemporaryTerms<output_type>(temp_term_union,
+ params_derivs_temporary_terms_idxs,
+ tt_output, tef_terms);
for (const auto &[order, tts] : params_derivs_temporary_terms)
- writeTemporaryTerms(tts, temp_term_union, params_derivs_temporary_terms_idxs, tt_output, output_type, tef_terms);
+ writeTemporaryTerms<output_type>(tts, temp_term_union, params_derivs_temporary_terms_idxs,
+ tt_output, tef_terms);
for (const auto &[indices, d1] : params_derivatives.find({ 0, 1 })->second)
{
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 0a276ad200a11cdf3cd254c13576b849d6f2e1fd..6dc463b7b1e20b6706d09eb07d6893ba66b605e2 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -114,107 +114,6 @@ StaticModel::writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk
code_file << FLDZ_{};
}
-void
-StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms) const
-{
- BlockSimulationType simulation_type = blocks[blk].simulation_type;
- int block_recursive_size = blocks[blk].getRecursiveSize();
-
- // The equations
- deriv_node_temp_terms_t tef_terms;
-
- 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;
- }
- };
-
- for (int eq = 0; eq < blocks[blk].size; eq++)
- {
- write_eq_tt(eq);
-
- EquationType equ_type = getBlockEquationType(blk, eq);
- BinaryOpNode *e = getBlockEquationExpr(blk, eq);
- expr_t lhs = e->arg1, rhs = e->arg2;
- switch (simulation_type)
- {
- case BlockSimulationType::evaluateBackward:
- case BlockSimulationType::evaluateForward:
- evaluation:
- if (equ_type == EquationType::evaluateRenormalized)
- {
- e = getBlockEquationRenormalizedExpr(blk, eq);
- lhs = e->arg1;
- rhs = e->arg2;
- }
- else if (equ_type != EquationType::evaluate)
- {
- cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
- exit(EXIT_FAILURE);
- }
- 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::solveBackwardSimple:
- case BlockSimulationType::solveForwardSimple:
- case BlockSimulationType::solveBackwardComplete:
- case BlockSimulationType::solveForwardComplete:
- if (eq < block_recursive_size)
- goto evaluation;
- output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
- << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
- << 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;
- default:
- cerr << "Incorrect type for block " << blk+1 << endl;
- exit(EXIT_FAILURE);
- }
- }
- // The Jacobian if we have to solve the block
- if (simulation_type != BlockSimulationType::evaluateBackward
- && simulation_type != BlockSimulationType::evaluateForward)
- {
- // Write temporary terms for derivatives
- write_eq_tt(blocks[blk].size);
-
- ostringstream i_output, j_output, v_output;
- for (int line_counter{ARRAY_SUBSCRIPT_OFFSET(output_type)};
- const auto &[indices, d] : blocks_derivatives[blk])
- {
- auto [eq, var, ignore] = indices;
- i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1-block_recursive_size
- << ';' << endl;
- j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << var+1-block_recursive_size
- << ';' << endl;
- v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
- << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
- d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
- v_output << ';' << endl;
- line_counter++;
- }
- output << i_output.str() << j_output.str() << v_output.str();
- }
-}
-
void
StaticModel::writeStaticPerBlockMFiles(const string &basename) const
{
@@ -258,7 +157,7 @@ StaticModel::writeStaticPerBlockMFiles(const string &basename) const
<< " g1_v=zeros(" << blocks_derivatives[blk].size() << ",1);" << endl
<< endl;
- writeStaticPerBlockHelper(blk, output, ExprNodeOutputType::matlabStaticModel, temporary_terms);
+ writeStaticPerBlockHelper<ExprNodeOutputType::matlabStaticModel>(blk, output, temporary_terms);
if (simulation_type != BlockSimulationType::evaluateBackward
&& simulation_type != BlockSimulationType::evaluateForward)
@@ -306,7 +205,7 @@ StaticModel::writeStaticPerBlockCFiles(const string &basename) const
output << "void static_" << blk+1 << "(double *restrict y, const double *restrict x, const double *restrict params, double *restrict T, double *restrict residual, double *restrict g1_i, double *restrict g1_j, double *restrict g1_v)" << endl;
output << '{' << endl;
- writeStaticPerBlockHelper(blk, output, ExprNodeOutputType::CStaticModel, temporary_terms);
+ writeStaticPerBlockHelper<ExprNodeOutputType::CStaticModel>(blk, output, temporary_terms);
output << '}' << endl
<< endl;
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index d781d0a35d532e2987681ead6cc9edbc751dcdfe..fa412814b5872f6a0c6982ffe00927721ccb3ca3 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -50,7 +50,8 @@ private:
void writeStaticBlockCFile(const string &basename) const;
//! Helper for writing a per-block static file of block decomposed model
- void writeStaticPerBlockHelper(int blk, ostream &output, ExprNodeOutputType output_type, temporary_terms_t &temporary_terms) const;
+ template<ExprNodeOutputType output_type>
+ void writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
//! Writes the per-block static files of block decomposed model (MATLAB version)
void writeStaticPerBlockMFiles(const string &basename) const;
@@ -172,6 +173,108 @@ public:
void addAllParamDerivId(set<int> &deriv_id_set) override;
};
+template<ExprNodeOutputType output_type>
+void
+StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
+{
+ BlockSimulationType simulation_type { blocks[blk].simulation_type };
+ int block_recursive_size { blocks[blk].getRecursiveSize() };
+
+ // The equations
+ deriv_node_temp_terms_t tef_terms;
+
+ 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;
+ }
+ };
+
+ for (int eq {0}; eq < blocks[blk].size; eq++)
+ {
+ write_eq_tt(eq);
+
+ EquationType equ_type { getBlockEquationType(blk, eq) };
+ BinaryOpNode *e { getBlockEquationExpr(blk, eq) };
+ expr_t lhs { e->arg1 }, rhs { e->arg2 };
+ switch (simulation_type)
+ {
+ case BlockSimulationType::evaluateBackward:
+ case BlockSimulationType::evaluateForward:
+ evaluation:
+ if (equ_type == EquationType::evaluateRenormalized)
+ {
+ e = getBlockEquationRenormalizedExpr(blk, eq);
+ lhs = e->arg1;
+ rhs = e->arg2;
+ }
+ else if (equ_type != EquationType::evaluate)
+ {
+ cerr << "Type mismatch for equation " << getBlockEquationID(blk, eq)+1 << endl;
+ exit(EXIT_FAILURE);
+ }
+ 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::solveBackwardSimple:
+ case BlockSimulationType::solveForwardSimple:
+ case BlockSimulationType::solveBackwardComplete:
+ case BlockSimulationType::solveForwardComplete:
+ if (eq < block_recursive_size)
+ goto evaluation;
+ output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+ << eq-block_recursive_size+ARRAY_SUBSCRIPT_OFFSET(output_type)
+ << 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;
+ default:
+ cerr << "Incorrect type for block " << blk+1 << endl;
+ exit(EXIT_FAILURE);
+ }
+ }
+ // The Jacobian if we have to solve the block
+ if (simulation_type != BlockSimulationType::evaluateBackward
+ && simulation_type != BlockSimulationType::evaluateForward)
+ {
+ // Write temporary terms for derivatives
+ write_eq_tt(blocks[blk].size);
+
+ ostringstream i_output, j_output, v_output;
+ for (int line_counter { ARRAY_SUBSCRIPT_OFFSET(output_type) };
+ const auto &[indices, d] : blocks_derivatives[blk])
+ {
+ const auto &[eq, var, ignore] {indices};
+ i_output << " g1_i" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << eq+1-block_recursive_size
+ << ';' << endl;
+ j_output << " g1_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=' << var+1-block_recursive_size
+ << ';' << endl;
+ v_output << " g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << line_counter
+ << RIGHT_ARRAY_SUBSCRIPT(output_type) << '=';
+ d->writeOutput(v_output, output_type, temporary_terms, blocks_temporary_terms_idxs);
+ v_output << ';' << endl;
+ line_counter++;
+ }
+ output << i_output.str() << j_output.str() << v_output.str();
+ }
+}
+
template<bool julia>
void
StaticModel::writeParamsDerivativesFile(const string &basename) const