diff --git a/src/DataTree.hh b/src/DataTree.hh
index 42305b15595bb296ad41717454d301f1ac07d55f..5e8ad6d70c73200f896c2d2ce13b241443a821c3 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 93e9627e4ea41bf148f433181385ae562cdbb557..aab43047d794b03bb45c9bf0833f404224e88aaa 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<int>(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()) };
+ int ncols{getJacobianColsNbr(false)};
+ int ntt { static_cast<int>(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<int>(NNZDerivatives.size()); i++)
output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
output << "];" << endl;
+
+ writeDriverSparseIndicesHelper<true>(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<int>::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<int>::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<true>(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 92c9c75cd6f2bfe3da2335c310922a7fa72a2dd6..c075ab210019df8b28b4114e753dc454e81ca2c1 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<pair<int, int>> 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<int, int> 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<int> &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 efc38c0de8b3f8145e55dbb7caab816129aedb36..47a2107f2df1e69765aa246e4056d1d7f8cd8d1f 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 58b1024ed9d011c659bbfd090a493e00aa4a85af..f7cbc96ca361aa8b5e40754e91266b7e70c7d3b8 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -83,12 +83,18 @@ using lag_equivalence_table_t = map<expr_t, map<int, expr_t>>;
//! 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 933f0cf7cec21ace7c64332c6e7d00033240faad..0af60f89f09eddaa95376e54a9d323c06bbbc1d3 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 99163df0f15a4c72617980afc646db323e45f4c0..be51592590455c1e2a8f84c4773764249d839519 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<int> &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<int>
+ModelTree::computeCSCColPtr(const SparseColumnMajorOrderMatrix &matrix, int ncols)
+{
+ vector<int> 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 fcc59547e4358869c4f1f374e4f7c1256514ce67..d4cf1f9d3a67a072620be0f67746a60de0ae8a67 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -119,6 +119,24 @@ protected:
derivatives, ... */
vector<int> NNZDerivatives;
+ // Used to order pairs of indices (row, col) according to column-major order
+ struct columnMajorOrderLess
+ {
+ bool
+ operator()(const pair<int, int> &p1, const pair<int, int> &p2) const
+ {
+ return p1.second < p2.second || (p1.second == p2.second && p1.first < p2.first);
+ }
+ };
+ using SparseColumnMajorOrderMatrix = map<pair<int, int>, expr_t, columnMajorOrderLess>;
+ /* The nonzero values of the sparse Jacobian in column-major order (which is
+ the natural order for Compressed Sparse Column (CSC) storage).
+ The pair of indices is (row, column). */
+ SparseColumnMajorOrderMatrix jacobian_sparse_column_major_order;
+ /* Column indices for the sparse Jacobian in Compressed Sparse Column (CSC)
+ storage (corresponds to the “jc” vector in MATLAB terminology) */
+ vector<int> jacobian_sparse_colptr;
+
//! 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<bool dynamic>
+ void writeDriverSparseIndicesHelper(ostream &output) const;
+
+ // Helper for writing sparse derivatives indices in JSON
+ template<bool dynamic>
+ void writeJsonSparseIndicesHelper(ostream &output) const;
+
/* Helper for writing JSON output for residuals and derivatives.
Returns mlv and derivatives output at each derivation order. */
template<bool dynamic>
@@ -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<int> 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<ExprNodeOutputType output_type>
pair<vector<ostringstream>, vector<ostringstream>>
ModelTree::writeModelFileHelper() const
{
+ constexpr bool sparse {isSparseModelOutput(output_type)};
+
vector<ostringstream> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
vector<ostringstream> tt_output(derivatives.size()); // Temp terms output (at all orders)
@@ -716,19 +748,37 @@ ModelTree::writeModelFileHelper() const
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])
+ 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<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,
- 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++;
}
-
- // Output symetric elements at order 2
- if (i == 2 && vidx[1] != vidx[2])
+ }
+ 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;
+
+ 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<bool dynamic>
+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<bool dynamic>
+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 8db084783107037e7111083cffcd81d591467e5e..adc4e3d9b7cb450304fa618b07b596db7c3f6d2f 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<false>(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<false>(output);
}
void
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index 4551ebd9fdf204b8f493d2912190a018576e573b..071c119aaef430e80e22f5592fdae941615900a6 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();
}