From 0278c8577ca5bce72992fb3d09e416dfa1701596 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 14 Sep 2022 17:07:08 +0200 Subject: [PATCH] Provisions for new sparse representation of dynamic/static files The new representation is only supported for MATLAB/Octave, C and Julia output for the time being. Bytecode and JSON are unsupported. This commit adds new fields in M_. This is a preliminary step for dynare#1859. --- src/DataTree.hh | 12 +- src/DynamicModel.cc | 42 +++--- src/DynamicModel.hh | 51 +++++-- src/ExprNode.cc | 33 ++++- src/ExprNode.hh | 35 ++++- src/ModFile.cc | 2 + src/ModelTree.cc | 26 ++++ src/ModelTree.hh | 347 ++++++++++++++++++++++++++++++++------------ src/StaticModel.cc | 7 +- src/StaticModel.hh | 4 +- 10 files changed, 420 insertions(+), 139 deletions(-) diff --git a/src/DataTree.hh b/src/DataTree.hh index 42305b15..5e8ad6d7 100644 --- a/src/DataTree.hh +++ b/src/DataTree.hh @@ -311,16 +311,20 @@ public: return symbol_table.getName(getSymbIDByDerivID(deriv_id)); } - //! Returns the column of the Jacobian associated to a derivation ID + /* Returns the column of the Jacobian associated to a derivation ID. + The “sparse” argument selects between the legacy representation and the + sparse representation. */ virtual int - getJacobianCol([[maybe_unused]] int deriv_id) const + getJacobianCol([[maybe_unused]] int deriv_id, [[maybe_unused]] bool sparse) const { throw UnknownDerivIDException(); } - //! Returns the number of columns of the Jacobian + /* Returns the number of columns of the Jacobian + The “sparse” argument selects between the legacy representation and the + sparse representation. */ virtual int - getJacobianColsNbr() const + getJacobianColsNbr([[maybe_unused]] bool sparse) const { throw UnknownDerivIDException(); } diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 93e9627e..aab43047 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -683,18 +683,18 @@ DynamicModel::writeDynamicMFile(const string &basename) const "", init_output, end_output, d_output[0], tt_output[0]); init_output.str(""); - init_output << "g1 = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ");"; + init_output << "g1 = zeros(" << equations.size() << ", " << getJacobianColsNbr(false) << ");"; writeDynamicMFileHelper(basename, "dynamic_g1", "g1", "dynamic_g1_tt", temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(), "dynamic_resid_tt", init_output, end_output, d_output[1], tt_output[1]); writeDynamicMWrapperFunction(basename, "g1"); // For order ≥ 2 - int ncols{getJacobianColsNbr()}; - int ntt { static_cast(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()) }; + int ncols{getJacobianColsNbr(false)}; + int ntt { static_cast(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())}; for (size_t i{2}; i < derivatives.size(); i++) { - ncols *= getJacobianColsNbr(); + ncols *= getJacobianColsNbr(false); ntt += temporary_terms_derivatives[i].size(); string gname{"g" + to_string(i)}; string gprevname{"g" + to_string(i-1)}; @@ -804,7 +804,7 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl << " @assert length(T) >= " << temporary_terms_derivatives[0].size() << endl << " @assert length(residual) == " << equations.size() << endl - << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr(false) << endl << " @assert length(params) == " << symbol_table.param_nbr() << endl << " if T_flag" << endl << " dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl @@ -832,8 +832,8 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl << " @assert length(T) >= " << temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl - << " @assert size(g1) == (" << equations.size() << ", " << getJacobianColsNbr() << ")" << endl - << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert size(g1) == (" << equations.size() << ", " << getJacobianColsNbr(false) << ")" << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr(false) << endl << " @assert length(params) == " << symbol_table.param_nbr() << endl << " if T_flag" << endl << " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl @@ -857,13 +857,13 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const << "end" << endl << endl; // dynamicG2! - int hessianColsNbr{getJacobianColsNbr() * getJacobianColsNbr()}; + int hessianColsNbr {getJacobianColsNbr(false) * getJacobianColsNbr(false)}; output << "function dynamicG2!(T::Vector{<: Real}, g2::Matrix{<: Real}," << endl << " y::Vector{<: Real}, x::Matrix{<: Real}, " << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl << " @assert length(T) >= " << temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl << " @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl - << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr(false) << endl << " @assert length(params) == " << symbol_table.param_nbr() << endl << " if T_flag" << endl << " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl @@ -887,14 +887,14 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const << "end" << endl << endl; // dynamicG3! - int ncols{hessianColsNbr * getJacobianColsNbr()}; + int ncols {hessianColsNbr * getJacobianColsNbr(false)}; output << "function dynamicG3!(T::Vector{<: Real}, g3::Matrix{<: Real}," << endl << " y::Vector{<: Real}, x::Matrix{<: Real}, " << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl << " @assert length(T) >= " << temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl << " @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl - << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr(false) << endl << " @assert length(params) == " << symbol_table.param_nbr() << endl << " if T_flag" << endl << " dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl @@ -1883,7 +1883,7 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl try { int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag); - output << " " << getJacobianCol(varID) + 1; + output << " " << getJacobianCol(varID, false) + 1; if (lag == -1) { sstatic = 0; @@ -2035,6 +2035,8 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl for (int i = 1; i < static_cast(NNZDerivatives.size()); i++) output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; "; output << "];" << endl; + + writeDriverSparseIndicesHelper(output); } void @@ -3133,8 +3135,11 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_c indices, check that we will not overflow (see #89). Note that such a check is not needed for parameter derivatives, since tensors for those are not stored as matrices. This check cannot be done before since - getJacobianColsNbr() is not yet set.*/ - if (log2(getJacobianColsNbr())*derivsOrder >= numeric_limits::digits) + getJacobianColsNbr() is not yet set. + We only do the check for the legacy representation, since the sparse + representation is not affected by this problem (TODO: thus the check can be + removed once the legacy representation is dropped). */ + if (log2(getJacobianColsNbr(false))*derivsOrder >= numeric_limits::digits) { cerr << "ERROR: The derivatives matrix of the " << modelClassName() << " is too large. Please decrease the approximation order." << endl; exit(EXIT_FAILURE); @@ -3927,12 +3932,13 @@ DynamicModel::computeDynJacobianCols() symbol_table.getType(symb_id) == SymbolType::endogenous) ordered_dyn_endo[{ lag, symbol_table.getTypeSpecificID(symb_id) }] = deriv_id; - // Fill the dynamic jacobian columns for endogenous + // Fill the dynamic jacobian columns for endogenous (legacy representation) for (int sorted_id{0}; const auto &[ignore, deriv_id] : ordered_dyn_endo) dyn_jacobian_cols_table[deriv_id] = sorted_id++; - // Fill the dynamic columns for exogenous and exogenous deterministic + /* Fill the dynamic columns for exogenous and exogenous deterministic (legacy + representation) */ for (const auto &[symb_lag, deriv_id] : deriv_id_table) { int symb_id{symb_lag.first}; @@ -4651,6 +4657,8 @@ DynamicModel::writeJsonOutput(ostream &output) const writeJsonAST(output); output << ", "; writeJsonVariableMapping(output); + output << ", "; + writeJsonSparseIndicesHelper(output); } void @@ -4776,7 +4784,7 @@ DynamicModel::writeJsonDynamicModelInfo(ostream &output) const if (lag != -max_endo_lag) output << ","; int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag); - output << " " << getJacobianCol(varID) + 1; + output << " " << getJacobianCol(varID, false) + 1; if (lag == -1) { sstatic = 0; diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 92c9c75c..c075ab21 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -63,8 +63,9 @@ private: //! Maps a deriv ID to a pair (symbol_id, lag) vector> inv_deriv_id_table; - //! Maps a deriv_id to the column index of the dynamic Jacobian - /*! Contains only endogenous, exogenous and exogenous deterministic */ + /* Maps a deriv_id to the column index of the dynamic Jacobian, in the legacy + representation. + Contains only endogenous, exogenous and exogenous deterministic */ map dyn_jacobian_cols_table; //! Maximum lag and lead over all types of variables (positive values) @@ -465,18 +466,46 @@ public: int getDerivID(int symb_id, int lag) const noexcept(false) override; int - getJacobianCol(int deriv_id) const override + getJacobianCol(int deriv_id, bool sparse) const override { - if (auto it = dyn_jacobian_cols_table.find(deriv_id); - it == dyn_jacobian_cols_table.end()) - throw UnknownDerivIDException(); + if (sparse) + { + SymbolType type {getTypeByDerivID(deriv_id)}; + int tsid {getTypeSpecificIDByDerivID(deriv_id)}; + int lag {getLagByDerivID(deriv_id)}; + if (type == SymbolType::endogenous) + { + assert(lag >= -1 && lag <= 1); + return tsid+(lag+1)*symbol_table.endo_nbr(); + } + else if (type == SymbolType::exogenous) + { + assert(lag == 0); + return tsid+3*symbol_table.endo_nbr(); + } + else if (type == SymbolType::exogenousDet) + { + assert(lag == 0); + return tsid+3*symbol_table.endo_nbr()+symbol_table.exo_nbr(); + } + else + throw UnknownDerivIDException(); + } else - return it->second; + { + if (auto it = dyn_jacobian_cols_table.find(deriv_id); + it == dyn_jacobian_cols_table.end()) + throw UnknownDerivIDException(); + else + return it->second; + } } int - getJacobianColsNbr() const override + getJacobianColsNbr(bool sparse) const override { - return dyn_jacobian_cols_table.size(); + return sparse ? + 3*symbol_table.endo_nbr() + symbol_table.exo_nbr() + symbol_table.exo_det_nbr() : + dyn_jacobian_cols_table.size(); } void addAllParamDerivId(set &deriv_id_set) override; @@ -904,7 +933,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << rp_output.str() - << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr(false) << ", " << symbol_table.param_nbr() << ");" << endl << gp_output.str() << "if nargout >= 3" << endl << "rpp = zeros(" << params_derivatives.at({ 0, 2 }).size() << ",4);" << endl @@ -939,7 +968,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const << "@inbounds begin" << endl << rp_output.str() << "end" << endl - << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr(false) << ", " << symbol_table.param_nbr() << ");" << endl << "@inbounds begin" << endl << gp_output.str() << "end" << endl diff --git a/src/ExprNode.cc b/src/ExprNode.cc index efc38c0d..47a2107f 100644 --- a/src/ExprNode.cc +++ b/src/ExprNode.cc @@ -1072,14 +1072,20 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, output_type) { case ExprNodeOutputType::juliaDynamicModel: + case ExprNodeOutputType::juliaSparseDynamicModel: case ExprNodeOutputType::matlabDynamicModel: + case ExprNodeOutputType::matlabSparseDynamicModel: case ExprNodeOutputType::CDynamicModel: - i = datatree.getJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type); + case ExprNodeOutputType::CSparseDynamicModel: + i = datatree.getJacobianCol(datatree.getDerivID(symb_id, lag), isSparseModelOutput(output_type)) + ARRAY_SUBSCRIPT_OFFSET(output_type); output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); break; case ExprNodeOutputType::CStaticModel: + case ExprNodeOutputType::CSparseStaticModel: case ExprNodeOutputType::juliaStaticModel: + case ExprNodeOutputType::juliaSparseStaticModel: case ExprNodeOutputType::matlabStaticModel: + case ExprNodeOutputType::matlabSparseStaticModel: i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type); output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); break; @@ -1145,9 +1151,17 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, else output << "x[it_" << lag << "+" << i << "*nb_row_x]"; break; + case ExprNodeOutputType::juliaSparseDynamicModel: + case ExprNodeOutputType::matlabSparseDynamicModel: + case ExprNodeOutputType::CSparseDynamicModel: + assert(lag == 0); + [[fallthrough]]; case ExprNodeOutputType::CStaticModel: + case ExprNodeOutputType::CSparseStaticModel: case ExprNodeOutputType::juliaStaticModel: + case ExprNodeOutputType::juliaSparseStaticModel: case ExprNodeOutputType::matlabStaticModel: + case ExprNodeOutputType::matlabSparseStaticModel: output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); break; case ExprNodeOutputType::matlabOutsideModel: @@ -1206,9 +1220,17 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, else output << "x[it_" << lag << "+" << i << "*nb_row_x]"; break; + case ExprNodeOutputType::juliaSparseDynamicModel: + case ExprNodeOutputType::matlabSparseDynamicModel: + case ExprNodeOutputType::CSparseDynamicModel: + assert(lag == 0); + [[fallthrough]]; case ExprNodeOutputType::CStaticModel: + case ExprNodeOutputType::CSparseStaticModel: case ExprNodeOutputType::juliaStaticModel: + case ExprNodeOutputType::juliaSparseStaticModel: case ExprNodeOutputType::matlabStaticModel: + case ExprNodeOutputType::matlabSparseStaticModel: output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); break; case ExprNodeOutputType::matlabOutsideModel: @@ -2830,7 +2852,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, output << "abs"; break; case UnaryOpcode::sign: - if (output_type == ExprNodeOutputType::CDynamicModel || output_type == ExprNodeOutputType::CStaticModel) + if (isCOutput(output_type)) output << "copysign"; else output << "sign"; @@ -2840,6 +2862,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, switch (output_type) { case ExprNodeOutputType::matlabDynamicModel: + case ExprNodeOutputType::matlabSparseDynamicModel: case ExprNodeOutputType::occbinDifferenceFile: new_output_type = ExprNodeOutputType::matlabDynamicSteadyStateOperator; break; @@ -2847,9 +2870,11 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, new_output_type = ExprNodeOutputType::latexDynamicSteadyStateOperator; break; case ExprNodeOutputType::CDynamicModel: + case ExprNodeOutputType::CSparseDynamicModel: new_output_type = ExprNodeOutputType::CDynamicSteadyStateOperator; break; case ExprNodeOutputType::juliaDynamicModel: + case ExprNodeOutputType::juliaSparseDynamicModel: new_output_type = ExprNodeOutputType::juliaDynamicSteadyStateOperator; break; default: @@ -2931,7 +2956,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, && arg->precedence(output_type, temporary_terms) < precedence(output_type, temporary_terms))) { output << LEFT_PAR(output_type); - if (op_code == UnaryOpcode::sign && (output_type == ExprNodeOutputType::CDynamicModel || output_type == ExprNodeOutputType::CStaticModel)) + if (op_code == UnaryOpcode::sign && isCOutput(output_type)) output << "1.0,"; close_parenthesis = true; } @@ -4568,7 +4593,7 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type, unpackPowerDeriv()->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms); else { - if (output_type == ExprNodeOutputType::juliaStaticModel || output_type == ExprNodeOutputType::juliaDynamicModel) + if (isJuliaOutput(output_type)) output << "get_power_deriv("; else output << "getPowerDeriv("; diff --git a/src/ExprNode.hh b/src/ExprNode.hh index 58b1024e..f7cbc96c 100644 --- a/src/ExprNode.hh +++ b/src/ExprNode.hh @@ -83,12 +83,18 @@ using lag_equivalence_table_t = map>; //! Possible types of output when writing ExprNode(s) (not used for bytecode) enum class ExprNodeOutputType { - matlabStaticModel, //!< Matlab code, static model - matlabDynamicModel, //!< Matlab code, dynamic model - CDynamicModel, //!< C code, dynamic model - CStaticModel, //!< C code, static model - juliaStaticModel, //!< Julia code, static model - juliaDynamicModel, //!< Julia code, dynamic model + matlabStaticModel, //!< Matlab code, static model, legacy representation + matlabDynamicModel, //!< Matlab code, dynamic model, legacy representation + matlabSparseStaticModel, //!< Matlab code, static model, sparse representation + matlabSparseDynamicModel, //!< Matlab code, dynamic model, sparse representation + CDynamicModel, //!< C code, dynamic model, legacy representation + CStaticModel, //!< C code, static model, legacy representation + CSparseDynamicModel, //!< C code, dynamic model, sparse representation + CSparseStaticModel, //!< C code, static model, sparse representation + juliaStaticModel, //!< Julia code, static model, legacy representation + juliaDynamicModel, //!< Julia code, dynamic model, legacy representation + juliaSparseStaticModel, //!< Julia code, static model, sparse representation + juliaSparseDynamicModel, //!< Julia code, dynamic model, sparse representation matlabOutsideModel, //!< Matlab code, outside model block (for example in initval) latexStaticModel, //!< LaTeX code, static model latexDynamicModel, //!< LaTeX code, dynamic model @@ -119,6 +125,8 @@ isMatlabOutput(ExprNodeOutputType output_type) { return output_type == ExprNodeOutputType::matlabStaticModel || output_type == ExprNodeOutputType::matlabDynamicModel + || output_type == ExprNodeOutputType::matlabSparseStaticModel + || output_type == ExprNodeOutputType::matlabSparseDynamicModel || output_type == ExprNodeOutputType::matlabOutsideModel || output_type == ExprNodeOutputType::matlabDynamicSteadyStateOperator || output_type == ExprNodeOutputType::steadyStateFile @@ -132,6 +140,8 @@ isJuliaOutput(ExprNodeOutputType output_type) { return output_type == ExprNodeOutputType::juliaStaticModel || output_type == ExprNodeOutputType::juliaDynamicModel + || output_type == ExprNodeOutputType::juliaSparseStaticModel + || output_type == ExprNodeOutputType::juliaSparseDynamicModel || output_type == ExprNodeOutputType::juliaDynamicSteadyStateOperator || output_type == ExprNodeOutputType::juliaSteadyStateFile || output_type == ExprNodeOutputType::juliaTimeDataFrame; @@ -142,6 +152,8 @@ isCOutput(ExprNodeOutputType output_type) { return output_type == ExprNodeOutputType::CDynamicModel || output_type == ExprNodeOutputType::CStaticModel + || output_type == ExprNodeOutputType::CSparseDynamicModel + || output_type == ExprNodeOutputType::CSparseStaticModel || output_type == ExprNodeOutputType::CDynamicSteadyStateOperator; } @@ -162,6 +174,17 @@ isSteadyStateOperatorOutput(ExprNodeOutputType output_type) || output_type == ExprNodeOutputType::juliaDynamicSteadyStateOperator; } +constexpr bool +isSparseModelOutput(ExprNodeOutputType output_type) +{ + return output_type == ExprNodeOutputType::matlabSparseStaticModel + || output_type == ExprNodeOutputType::matlabSparseDynamicModel + || output_type == ExprNodeOutputType::juliaSparseStaticModel + || output_type == ExprNodeOutputType::juliaSparseDynamicModel + || output_type == ExprNodeOutputType::CSparseDynamicModel + || output_type == ExprNodeOutputType::CSparseStaticModel; +} + constexpr bool isAssignmentLHSBytecodeOutput(ExprNodeBytecodeOutputType output_type) { diff --git a/src/ModFile.cc b/src/ModFile.cc index 933f0cf7..0af60f89 100644 --- a/src/ModFile.cc +++ b/src/ModFile.cc @@ -1167,6 +1167,8 @@ ModFile::writeJsonOutputParsingCheck(const string &basename, JsonFileOutputType symbol_table.writeJsonOutput(output); output << ", "; dynamic_model.writeJsonOutput(output); + output << ", "; + static_model.writeJsonOutput(output); if (!statements.empty() || !var_model_table.empty() diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 99163df0..be515925 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -71,6 +71,8 @@ ModelTree::copyHelper(const ModelTree &m) derivatives.push_back(convert_deriv_map(it)); for (const auto &it : m.params_derivatives) params_derivatives.emplace(it.first, convert_deriv_map(it.second)); + for (const auto &it : m.jacobian_sparse_column_major_order) + jacobian_sparse_column_major_order.emplace(it.first, f(it.second)); auto convert_temporary_terms_t = [f](const temporary_terms_t &tt) { @@ -150,6 +152,7 @@ ModelTree::ModelTree(const ModelTree &m) : equation_tags{m.equation_tags}, computed_derivs_order{m.computed_derivs_order}, NNZDerivatives{m.NNZDerivatives}, + jacobian_sparse_colptr{m.jacobian_sparse_colptr}, eq_idx_block2orig{m.eq_idx_block2orig}, endo_idx_block2orig{m.endo_idx_block2orig}, eq_idx_orig2block{m.eq_idx_orig2block}, @@ -177,6 +180,10 @@ ModelTree::operator=(const ModelTree &m) NNZDerivatives = m.NNZDerivatives; derivatives.clear(); + + jacobian_sparse_column_major_order.clear(); + jacobian_sparse_colptr = m.jacobian_sparse_colptr; + params_derivatives.clear(); temporary_terms_derivatives.clear(); @@ -902,6 +909,11 @@ ModelTree::computeDerivatives(int order, const set &vars) ++NNZDerivatives[1]; } + // Compute the sparse representation of the Jacobian + for (const auto &[indices, d1] : derivatives[1]) + jacobian_sparse_column_major_order.emplace(pair{indices[0], getJacobianCol(indices[1], true)}, d1); + jacobian_sparse_colptr = computeCSCColPtr(jacobian_sparse_column_major_order, getJacobianColsNbr(true)); + // Higher-order derivatives for (int o = 2; o <= order; o++) for (const auto &[lower_indices, lower_d] : derivatives[o-1]) @@ -1972,3 +1984,17 @@ ModelTree::computingPassBlock(const eval_context_t &eval_context, bool no_tmp_te computeBlockTemporaryTerms(); block_decomposed = true; } + +vector +ModelTree::computeCSCColPtr(const SparseColumnMajorOrderMatrix &matrix, int ncols) +{ + vector colptr(ncols+1, matrix.size()); + for (int k {0}, current_col {0}; + const auto &[indices, d1] : matrix) + { + while (indices.second >= current_col) + colptr[current_col++] = k; + k++; + } + return colptr; +} diff --git a/src/ModelTree.hh b/src/ModelTree.hh index fcc59547..d4cf1f9d 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -119,6 +119,24 @@ protected: derivatives, ... */ vector NNZDerivatives; + // Used to order pairs of indices (row, col) according to column-major order + struct columnMajorOrderLess + { + bool + operator()(const pair &p1, const pair &p2) const + { + return p1.second < p2.second || (p1.second == p2.second && p1.first < p2.first); + } + }; + using SparseColumnMajorOrderMatrix = map, expr_t, columnMajorOrderLess>; + /* The nonzero values of the sparse Jacobian in column-major order (which is + the natural order for Compressed Sparse Column (CSC) storage). + The pair of indices is (row, column). */ + SparseColumnMajorOrderMatrix jacobian_sparse_column_major_order; + /* Column indices for the sparse Jacobian in Compressed Sparse Column (CSC) + storage (corresponds to the “jc” vector in MATLAB terminology) */ + vector jacobian_sparse_colptr; + //! Derivatives with respect to parameters /*! The key of the outer map is a pair (derivation order w.r.t. endogenous, derivation order w.r.t. parameters). For e.g., { 1, 2 } corresponds to the jacobian @@ -303,6 +321,14 @@ protected: const temporary_terms_t &temporary_terms_union, const deriv_node_temp_terms_t &tef_terms) const; + // Helper for writing sparse derivatives indices in MATLAB/Octave driver file + template + void writeDriverSparseIndicesHelper(ostream &output) const; + + // Helper for writing sparse derivatives indices in JSON + template + void writeJsonSparseIndicesHelper(ostream &output) const; + /* Helper for writing JSON output for residuals and derivatives. Returns mlv and derivatives output at each derivation order. */ template @@ -490,6 +516,10 @@ protected: // Returns a human-readable string describing the model class (e.g. “dynamic model”…) virtual string modelClassName() const = 0; + /* Given a sparse matrix in column major order, returns the colptr pointer for + the CSC storage */ + static vector computeCSCColPtr(const SparseColumnMajorOrderMatrix &matrix, int ncols); + private: //! Internal helper for the copy constructor and assignment operator /*! Copies all the structures that contain ExprNode*, by the converting the @@ -699,6 +729,8 @@ template pair, vector> ModelTree::writeModelFileHelper() const { + constexpr bool sparse {isSparseModelOutput(output_type)}; + vector d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual) vector tt_output(derivatives.size()); // Temp terms output (at all orders) @@ -716,19 +748,37 @@ ModelTree::writeModelFileHelper() const writeTemporaryTerms(temporary_terms_derivatives[1], temp_term_union, temporary_terms_idxs, tt_output[1], tef_terms); - for (const auto &[indices, d1] : derivatives[1]) + if constexpr(sparse) { - auto [eq, var] = vectorToTuple<2>(indices); + // NB: we iterate over the Jacobian reordered in column-major order + // Indices of rows and columns are output in M_ and the JSON file (since they are constant) + for (int k {0}; + const auto &[row_col, d1] : jacobian_sparse_column_major_order) + { + d_output[1] << "g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; + d1->writeOutput(d_output[1], output_type, temp_term_union, temporary_terms_idxs, tef_terms); + d_output[1] << ";" << endl; + k++; + } + } + else // Legacy representation (dense matrix) + { + for (const auto &[indices, d1] : derivatives[1]) + { + auto [eq, var] = vectorToTuple<2>(indices); - d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type); - if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type)) - d_output[1] << eq + 1 << "," << getJacobianCol(var) + 1; - else - d_output[1] << eq + getJacobianCol(var)*equations.size(); - d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; - d1->writeOutput(d_output[1], output_type, - temp_term_union, temporary_terms_idxs, tef_terms); - d_output[1] << ";" << endl; + d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type); + if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type)) + d_output[1] << eq + 1 << "," << getJacobianCol(var, sparse) + 1; + else + d_output[1] << eq + getJacobianCol(var, sparse)*equations.size(); + d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; + d1->writeOutput(d_output[1], output_type, + temp_term_union, temporary_terms_idxs, tef_terms); + d_output[1] << ";" << endl; + } } } @@ -739,58 +789,49 @@ ModelTree::writeModelFileHelper() const writeTemporaryTerms(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, - then the third. This gives a significant performance boost in use_dll - mode (at both compilation and runtime), because it facilitates memory - accesses and expression reusage. */ - ostringstream i_output, j_output, v_output; - - for (int k{0}; // Current line index in the 3-column matrix - const auto &[vidx, d] : derivatives[i]) + if constexpr(sparse) { - int eq{vidx[0]}; - - int col_idx{0}; - for (size_t j = 1; j < vidx.size(); j++) + /* List non-zero elements of the tensor in row-major order (this is + suitable for the k-order solver according to Normann). */ + // Tensor indices are output in M_ and the JSON file (since they are constant) + for (int k {0}; + const auto &[vidx, d] : derivatives[i]) { - col_idx *= getJacobianColsNbr(); - col_idx += getJacobianCol(vidx[j]); - } - - if constexpr(isJuliaOutput(output_type)) - { - d_output[i] << " g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = "; + d_output[i] << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms); - d_output[i] << endl; - } - else - { - i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << eq + 1 << ";" << endl; - j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << col_idx + 1 << ";" << endl; - v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; - d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms); - v_output << ";" << endl; - + d_output[i] << ";" << endl; k++; } + } + else // Legacy representation + { + /* When creating the sparse matrix (in MATLAB or C mode), since storage + is in column-major order, output the first column, then the second, + then the third. This gives a significant performance boost in use_dll + mode (at both compilation and runtime), because it facilitates memory + accesses and expression reusage. */ + ostringstream i_output, j_output, v_output; - // Output symetric elements at order 2 - if (i == 2 && vidx[1] != vidx[2]) + for (int k{0}; // Current line index in the 3-column matrix + const auto &[vidx, d] : derivatives[i]) { - int col_idx_sym{getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1])}; + int eq{vidx[0]}; + + int col_idx{0}; + for (size_t j = 1; j < vidx.size(); j++) + { + col_idx *= getJacobianColsNbr(sparse); + col_idx += getJacobianCol(vidx[j], sparse); + } if constexpr(isJuliaOutput(output_type)) - d_output[2] << " g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = " - << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl; + { + d_output[i] << " g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = "; + d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms); + d_output[i] << endl; + } else { i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) @@ -800,20 +841,48 @@ ModelTree::writeModelFileHelper() const j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) << k + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << col_idx_sym + 1 << ";" << endl; + << "=" << col_idx + 1 << ";" << endl; v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" - << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl; + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; + d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms); + v_output << ";" << endl; k++; } + + // Output symetric elements at order 2 + if (i == 2 && vidx[1] != vidx[2]) + { + int col_idx_sym{getJacobianCol(vidx[2], sparse) * getJacobianColsNbr(sparse) + getJacobianCol(vidx[1], sparse)}; + + if constexpr(isJuliaOutput(output_type)) + d_output[2] << " g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = " + << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl; + else + { + i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) + << "=" << eq + 1 << ";" << endl; + j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) + << "=" << col_idx_sym + 1 << ";" << endl; + v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" + << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl; + + k++; + } + } } + if constexpr(!isJuliaOutput(output_type)) + d_output[i] << i_output.str() << j_output.str() << v_output.str(); } - if constexpr(!isJuliaOutput(output_type)) - d_output[i] << i_output.str() << j_output.str() << v_output.str(); } if constexpr(isMatlabOutput(output_type)) @@ -985,7 +1054,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext, << endl << " if (nlhs >= 2)" << endl << " {" << endl - << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getJacobianColsNbr() << ", mxREAL);" << endl + << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getJacobianColsNbr(false) << ", mxREAL);" << endl << " double *g1 = mxGetPr(plhs[1]);" << endl << " " << prefix << "g1_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl << " " << prefix << "g1(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, g1);" << endl @@ -999,7 +1068,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext, << " " << prefix << "g2_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl << " " << prefix << "g2(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl - << " mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr(false)*getJacobianColsNbr(false) << ");" << endl << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl << " plhs[2] = plhs_sparse[0];" << endl @@ -1019,7 +1088,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext, << " " << prefix << "g3_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << endl << " " << prefix << "g3(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T, mxGetPr(g3_i), mxGetPr(g3_j), mxGetPr(g3_v));" << endl << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl - << " mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr(false)*getJacobianColsNbr(false)*getJacobianColsNbr(false) << ");" << endl << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g3_i, g3_j, g3_v, m, n };" << endl << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl << " plhs[3] = plhs_sparse[0];" << endl @@ -1130,6 +1199,8 @@ ModelTree::writeParamsDerivativesFileHelper() const { static_assert(!isCOutput(output_type), "C output is not implemented"); + constexpr bool sparse {isSparseModelOutput(output_type)}; + ostringstream tt_output; // Used for storing model temp vars and equations ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters ostringstream gp_output; // 1st deriv. of Jacobian w.r.t. parameters @@ -1161,7 +1232,7 @@ ModelTree::writeParamsDerivativesFileHelper() const { auto [eq, var, param] { vectorToTuple<3>(indices) }; - int var_col { getJacobianCol(var) + 1 }; + int var_col { getJacobianCol(var, sparse) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 }; gp_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << var_col @@ -1213,7 +1284,7 @@ ModelTree::writeParamsDerivativesFileHelper() const { auto [eq, var, param1, param2] { vectorToTuple<4>(indices) }; - int var_col { getJacobianCol(var) + 1 }; + int var_col { getJacobianCol(var, sparse) + 1 }; int param1_col { getTypeSpecificIDByDerivID(param1) + 1 }; int param2_col { getTypeSpecificIDByDerivID(param2) + 1 }; @@ -1256,8 +1327,8 @@ ModelTree::writeParamsDerivativesFileHelper() const { auto [eq, var1, var2, param] { vectorToTuple<4>(indices) }; - int var1_col { getJacobianCol(var1) + 1 }; - int var2_col { getJacobianCol(var2) + 1 }; + int var1_col { getJacobianCol(var1, sparse) + 1 }; + int var2_col { getJacobianCol(var2, sparse) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 }; hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1" @@ -1301,9 +1372,9 @@ ModelTree::writeParamsDerivativesFileHelper() const { auto [eq, var1, var2, var3, param] { vectorToTuple<5>(indices) }; - int var1_col { getJacobianCol(var1) + 1 }; - int var2_col { getJacobianCol(var2) + 1 }; - int var3_col { getJacobianCol(var3) + 1 }; + int var1_col { getJacobianCol(var1, sparse) + 1 }; + int var2_col { getJacobianCol(var2, sparse) + 1 }; + int var3_col { getJacobianCol(var3, sparse) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 }; g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1" @@ -1517,7 +1588,7 @@ ModelTree::writeBytecodeHelper(BytecodeWriter &code_file) const if constexpr(dynamic) { // Bytecode MEX uses a separate matrix for exogenous and exodet Jacobians - int jacob_col { type == SymbolType::endogenous ? getJacobianCol(deriv_id) : tsid }; + int jacob_col { type == SymbolType::endogenous ? getJacobianCol(deriv_id, false) : tsid }; code_file << FSTPG3_{eq, tsid, lag, jacob_col}; } else @@ -1763,7 +1834,7 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const d_output[0] << ", "; writeJsonModelEquations(d_output[0], true); - int ncols { getJacobianColsNbr() }; + int ncols { getJacobianColsNbr(false) }; for (size_t i {1}; i < derivatives.size(); i++) { string matrix_name { i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative"}; @@ -1785,8 +1856,8 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const int col_idx {0}; for (size_t j {1}; j < vidx.size(); j++) { - col_idx *= getJacobianColsNbr(); - col_idx += getJacobianCol(vidx[j]); + col_idx *= getJacobianColsNbr(false); + col_idx += getJacobianCol(vidx[j], false); } if (writeDetails) @@ -1798,7 +1869,7 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian { - int col_idx_sym { getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1])}; + int col_idx_sym { getJacobianCol(vidx[2], false) * getJacobianColsNbr(false) + getJacobianCol(vidx[1], false)}; d_output[i] << ", " << col_idx_sym + 1; } if (i > 1) @@ -1818,7 +1889,7 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const } d_output[i] << "]}"; - ncols *= getJacobianColsNbr(); + ncols *= getJacobianColsNbr(false); } return { move(mlv_output), move(d_output) }; @@ -1877,7 +1948,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const gp_output << R"("deriv_jacobian_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvarcols": )" << getJacobianColsNbr() + << R"(, "nvarcols": )" << getJacobianColsNbr(false) << R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; for (bool printed_something {false}; @@ -1888,7 +1959,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const auto [eq, var, param] { vectorToTuple<3>(vidx) }; - int var_col { getJacobianCol(var) + 1 }; + int var_col { getJacobianCol(var, false) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 }; if (writeDetails) @@ -1948,7 +2019,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const gpp_output << R"("second_deriv_jacobian_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvarcols": )" << getJacobianColsNbr() + << R"(, "nvarcols": )" << getJacobianColsNbr(false) << R"(, "nparam1cols": )" << symbol_table.param_nbr() << R"(, "nparam2cols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; @@ -1960,7 +2031,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const auto [eq, var, param1, param2] { vectorToTuple<4>(vidx) }; - int var_col { getJacobianCol(var) + 1 }; + int var_col { getJacobianCol(var, false) + 1 }; int param1_col { getTypeSpecificIDByDerivID(param1) + 1 }; int param2_col { getTypeSpecificIDByDerivID(param2) + 1 }; @@ -1990,8 +2061,8 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const hp_output << R"("derivative_hessian_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvar1cols": )" << getJacobianColsNbr() - << R"(, "nvar2cols": )" << getJacobianColsNbr() + << R"(, "nvar1cols": )" << getJacobianColsNbr(false) + << R"(, "nvar2cols": )" << getJacobianColsNbr(false) << R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; for (bool printed_something {false}; @@ -2002,8 +2073,8 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const auto [eq, var1, var2, param] { vectorToTuple<4>(vidx) }; - int var1_col { getJacobianCol(var1) + 1 }; - int var2_col { getJacobianCol(var2) + 1 }; + int var1_col { getJacobianCol(var1, false) + 1 }; + int var2_col { getJacobianCol(var2, false) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 }; if (writeDetails) @@ -2036,9 +2107,9 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const { g3p_output << R"("derivative_g3_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvar1cols": )" << getJacobianColsNbr() - << R"(, "nvar2cols": )" << getJacobianColsNbr() - << R"(, "nvar3cols": )" << getJacobianColsNbr() + << R"(, "nvar1cols": )" << getJacobianColsNbr(false) + << R"(, "nvar2cols": )" << getJacobianColsNbr(false) + << R"(, "nvar3cols": )" << getJacobianColsNbr(false) << R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; for (bool printed_something {false}; @@ -2049,9 +2120,9 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const auto [eq, var1, var2, var3, param] { vectorToTuple<5>(vidx) }; - int var1_col { getJacobianCol(var1) + 1 }; - int var2_col { getJacobianCol(var2) + 1 }; - int var3_col { getJacobianCol(var3) + 1 }; + int var1_col { getJacobianCol(var1, false) + 1 }; + int var2_col { getJacobianCol(var2, false) + 1 }; + int var3_col { getJacobianCol(var3, false) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 }; if (writeDetails) @@ -2084,4 +2155,98 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output) }; } +template +void +ModelTree::writeDriverSparseIndicesHelper(ostream &output) const +{ + // TODO: when C++20 support is complete, mark this constexpr + const string model_name {dynamic ? "dynamic" : "static"}; + + // Write indices for the sparse Jacobian (both naive and CSC storage) + output << "M_." << model_name << "_g1_sparse_rowval = int32(["; + for (const auto &[indices, d1] : jacobian_sparse_column_major_order) + output << indices.first+1 << ' '; + output << "]);" << endl + << "M_." << model_name << "_g1_sparse_colval = int32(["; + for (const auto &[indices, d1] : jacobian_sparse_column_major_order) + output << indices.second+1 << ' '; + output << "]);" << endl + << "M_." << model_name << "_g1_sparse_colptr = int32(["; + for (int it : jacobian_sparse_colptr) + output << it+1 << ' '; + output << "]);" << endl; + + // Write indices for the sparse higher-order derivatives + for (int i {2}; i < computed_derivs_order; i++) + { + output << "M_." << model_name << "_g" << i << "_sparse_indices = int32(["; + for (const auto &[vidx, d] : derivatives[i]) + { + for (int it : vidx) + output << it+1 << ' '; + output << ';' << endl; + } + output << "]);" << endl; + } +} + +template +void +ModelTree::writeJsonSparseIndicesHelper(ostream &output) const +{ + // TODO: when C++20 support is complete, mark this constexpr + const string model_name {dynamic ? "dynamic" : "static"}; + + // Write indices for the sparse Jacobian (both naive and CSC storage) + output << '"' << model_name << R"(_g1_sparse_rowval": [)"; + for (bool printed_something {false}; + const auto &[indices, d1] : jacobian_sparse_column_major_order) + { + if (exchange(printed_something, true)) + output << ", "; + output << indices.first+1; + } + output << "], " << endl + << '"' << model_name << R"(_g1_sparse_colval": [)"; + for (bool printed_something {false}; + const auto &[indices, d1] : jacobian_sparse_column_major_order) + { + if (exchange(printed_something, true)) + output << ", "; + output << indices.second+1; + } + output << "], " << endl + << '"' << model_name << R"(_g1_sparse_colptr": [)"; + for (bool printed_something {false}; + int it : jacobian_sparse_colptr) + { + if (exchange(printed_something, true)) + output << ", "; + output << it+1; + } + output << ']' << endl; + + // Write indices for the sparse higher-order derivatives + for (int i {2}; i < computed_derivs_order; i++) + { + output << R"(, ")" << model_name << "_g" << i << R"(_sparse_indices": [)"; + for (bool printed_something {false}; + const auto &[vidx, d] : derivatives[i]) + { + if (exchange(printed_something, true)) + output << ", "; + output << '['; + for (bool printed_something2 {false}; + int it : vidx) + { + if (exchange(printed_something2, true)) + output << ", "; + output << it+1; + } + output << ']' << endl; + } + output << ']' << endl; + } +} + #endif diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 8db08478..adc4e3d9 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -1022,6 +1022,8 @@ StaticModel::writeDriverOutput(ostream &output, bool block) const if (block) writeBlockDriverOutput(output); + + writeDriverSparseIndicesHelper(output); } void @@ -1301,10 +1303,7 @@ StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const void StaticModel::writeJsonOutput(ostream &output) const { - deriv_node_temp_terms_t tef_terms; - writeJsonModelLocalVariables(output, false, tef_terms); - output << ", "; - writeJsonModelEquations(output, false); + writeJsonSparseIndicesHelper(output); } void diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 4551ebd9..071c119a 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -77,12 +77,12 @@ private: int getTypeSpecificIDByDerivID(int deriv_id) const override; int - getJacobianCol(int deriv_id) const override + getJacobianCol(int deriv_id, [[maybe_unused]] bool sparse) const override { return getTypeSpecificIDByDerivID(deriv_id); } int - getJacobianColsNbr() const override + getJacobianColsNbr([[maybe_unused]] bool sparse) const override { return symbol_table.endo_nbr(); }