stepan
/
preprocessor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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.
master
Sébastien Villemot 2022-09-14 17:07:08 +02:00
parent 30853e7360
commit 0278c8577c
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
10 changed files with 420 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/DataTree.hh
View File

@ -311,16 +311,20 @@ public:
return symbol_table.getName(getSymbIDByDerivID(deriv_id)); 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 virtual int
getJacobianCol([[maybe_unused]] int deriv_id) const getJacobianCol([[maybe_unused]] int deriv_id, [[maybe_unused]] bool sparse) const
{ {
throw UnknownDerivIDException(); 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 virtual int
getJacobianColsNbr() const getJacobianColsNbr([[maybe_unused]] bool sparse) const
{ {
throw UnknownDerivIDException(); throw UnknownDerivIDException();
} }

42
src/DynamicModel.cc
View File

@ -683,18 +683,18 @@ DynamicModel::writeDynamicMFile(const string &basename) const
"", init_output, end_output, d_output[0], tt_output[0]); "", init_output, end_output, d_output[0], tt_output[0]);
init_output.str(""); 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", writeDynamicMFileHelper(basename, "dynamic_g1", "g1", "dynamic_g1_tt",
temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(), temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(),
"dynamic_resid_tt", init_output, end_output, d_output[1], tt_output[1]); "dynamic_resid_tt", init_output, end_output, d_output[1], tt_output[1]);
writeDynamicMWrapperFunction(basename, "g1"); writeDynamicMWrapperFunction(basename, "g1");
// For order ≥ 2 // For order ≥ 2
int ncols{getJacobianColsNbr()}; int ncols{getJacobianColsNbr(false)};
int ntt { static_cast<int>(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()) }; int ntt { static_cast<int>(temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())};
for (size_t i{2}; i < derivatives.size(); i++) for (size_t i{2}; i < derivatives.size(); i++)
{ {
ncols *= getJacobianColsNbr(); ncols *= getJacobianColsNbr(false);
ntt += temporary_terms_derivatives[i].size(); ntt += temporary_terms_derivatives[i].size();
string gname{"g" + to_string(i)}; string gname{"g" + to_string(i)};
string gprevname{"g" + to_string(i-1)}; 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 << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
<< " @assert length(T) >= " << temporary_terms_derivatives[0].size() << endl << " @assert length(T) >= " << temporary_terms_derivatives[0].size() << endl
<< " @assert length(residual) == " << equations.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 << " @assert length(params) == " << symbol_table.param_nbr() << endl
<< " if T_flag" << endl << " if T_flag" << endl
<< " dynamicResidTT!(T, y, x, params, steady_state, it_)" << 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 << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
<< " @assert length(T) >= " << " @assert length(T) >= "
<< temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl << temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
<< " @assert size(g1) == (" << equations.size() << ", " << getJacobianColsNbr() << ")" << endl << " @assert size(g1) == (" << equations.size() << ", " << getJacobianColsNbr(false) << ")" << 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 << " @assert length(params) == " << symbol_table.param_nbr() << endl
<< " if T_flag" << endl << " if T_flag" << endl
<< " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl << " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl
@ -857,13 +857,13 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
<< "end" << endl << endl; << "end" << endl << endl;
// dynamicG2! // dynamicG2!
int hessianColsNbr{getJacobianColsNbr() * getJacobianColsNbr()}; int hessianColsNbr {getJacobianColsNbr(false) * getJacobianColsNbr(false)};
output << "function dynamicG2!(T::Vector{<: Real}, g2::Matrix{<: Real}," << endl output << "function dynamicG2!(T::Vector{<: Real}, g2::Matrix{<: Real}," << endl
<< " y::Vector{<: Real}, x::Matrix{<: Real}, " << " y::Vector{<: Real}, x::Matrix{<: Real}, "
<< "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl << "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 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 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 << " @assert length(params) == " << symbol_table.param_nbr() << endl
<< " if T_flag" << endl << " if T_flag" << endl
<< " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl << " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl
@ -887,14 +887,14 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
<< "end" << endl << endl; << "end" << endl << endl;
// dynamicG3! // dynamicG3!
int ncols{hessianColsNbr * getJacobianColsNbr()}; int ncols {hessianColsNbr * getJacobianColsNbr(false)};
output << "function dynamicG3!(T::Vector{<: Real}, g3::Matrix{<: Real}," << endl output << "function dynamicG3!(T::Vector{<: Real}, g3::Matrix{<: Real}," << endl
<< " y::Vector{<: Real}, x::Matrix{<: Real}, " << " y::Vector{<: Real}, x::Matrix{<: Real}, "
<< "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
<< " @assert length(T) >= " << " @assert length(T) >= "
<< temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl << 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 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 << " @assert length(params) == " << symbol_table.param_nbr() << endl
<< " if T_flag" << endl << " if T_flag" << endl
<< " dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl << " dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl
@ -1883,7 +1883,7 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl
try try
{ {
int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag); int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag);
output << " " << getJacobianCol(varID) + 1; output << " " << getJacobianCol(varID, false) + 1;
if (lag == -1) if (lag == -1)
{ {
sstatic = 0; 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++) for (int i = 1; i < static_cast<int>(NNZDerivatives.size()); i++)
output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; "; output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
output << "];" << endl; output << "];" << endl;
writeDriverSparseIndicesHelper<true>(output);
} }
void 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 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 is not needed for parameter derivatives, since tensors for those are not
stored as matrices. This check cannot be done before since stored as matrices. This check cannot be done before since
getJacobianColsNbr() is not yet set.*/ getJacobianColsNbr() is not yet set.
if (log2(getJacobianColsNbr())*derivsOrder >= numeric_limits<int>::digits) 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; cerr << "ERROR: The derivatives matrix of the " << modelClassName() << " is too large. Please decrease the approximation order." << endl;
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
@ -3927,12 +3932,13 @@ DynamicModel::computeDynJacobianCols()
symbol_table.getType(symb_id) == SymbolType::endogenous) symbol_table.getType(symb_id) == SymbolType::endogenous)
ordered_dyn_endo[{ lag, symbol_table.getTypeSpecificID(symb_id) }] = deriv_id; 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}; for (int sorted_id{0};
const auto &[ignore, deriv_id] : ordered_dyn_endo) const auto &[ignore, deriv_id] : ordered_dyn_endo)
dyn_jacobian_cols_table[deriv_id] = sorted_id++; 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) for (const auto &[symb_lag, deriv_id] : deriv_id_table)
{ {
int symb_id{symb_lag.first}; int symb_id{symb_lag.first};
@ -4651,6 +4657,8 @@ DynamicModel::writeJsonOutput(ostream &output) const
writeJsonAST(output); writeJsonAST(output);
output << ", "; output << ", ";
writeJsonVariableMapping(output); writeJsonVariableMapping(output);
output << ", ";
writeJsonSparseIndicesHelper<true>(output);
} }
void void
@ -4776,7 +4784,7 @@ DynamicModel::writeJsonDynamicModelInfo(ostream &output) const
if (lag != -max_endo_lag) if (lag != -max_endo_lag)
output << ","; output << ",";
int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag); int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag);
output << " " << getJacobianCol(varID) + 1; output << " " << getJacobianCol(varID, false) + 1;
if (lag == -1) if (lag == -1)
{ {
sstatic = 0; sstatic = 0;

51
src/DynamicModel.hh
View File

@ -63,8 +63,9 @@ private:
//! Maps a deriv ID to a pair (symbol_id, lag) //! Maps a deriv ID to a pair (symbol_id, lag)
vector<pair<int, int>> inv_deriv_id_table; vector<pair<int, int>> inv_deriv_id_table;
//! Maps a deriv_id to the column index of the dynamic Jacobian /* Maps a deriv_id to the column index of the dynamic Jacobian, in the legacy
/*! Contains only endogenous, exogenous and exogenous deterministic */ representation.
Contains only endogenous, exogenous and exogenous deterministic */
map<int, int> dyn_jacobian_cols_table; map<int, int> dyn_jacobian_cols_table;
//! Maximum lag and lead over all types of variables (positive values) //! 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 getDerivID(int symb_id, int lag) const noexcept(false) override;
int 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); if (sparse)
it == dyn_jacobian_cols_table.end()) {
throw UnknownDerivIDException(); 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 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 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; void addAllParamDerivId(set<int> &deriv_id_set) override;
@ -904,7 +933,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
<< "rp = zeros(" << equations.size() << ", " << "rp = zeros(" << equations.size() << ", "
<< symbol_table.param_nbr() << ");" << endl << symbol_table.param_nbr() << ");" << endl
<< rp_output.str() << 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() << gp_output.str()
<< "if nargout >= 3" << endl << "if nargout >= 3" << endl
<< "rpp = zeros(" << params_derivatives.at({ 0, 2 }).size() << ",4);" << endl << "rpp = zeros(" << params_derivatives.at({ 0, 2 }).size() << ",4);" << endl
@ -939,7 +968,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
<< "@inbounds begin" << endl << "@inbounds begin" << endl
<< rp_output.str() << rp_output.str()
<< "end" << endl << "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 << "@inbounds begin" << endl
<< gp_output.str() << gp_output.str()
<< "end" << endl << "end" << endl

33
src/ExprNode.cc
View File

@ -1072,14 +1072,20 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output_type) output_type)
{ {
case ExprNodeOutputType::juliaDynamicModel: case ExprNodeOutputType::juliaDynamicModel:
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabDynamicModel: case ExprNodeOutputType::matlabDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CDynamicModel: 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); output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break; break;
case ExprNodeOutputType::CStaticModel: case ExprNodeOutputType::CStaticModel:
case ExprNodeOutputType::CSparseStaticModel:
case ExprNodeOutputType::juliaStaticModel: case ExprNodeOutputType::juliaStaticModel:
case ExprNodeOutputType::juliaSparseStaticModel:
case ExprNodeOutputType::matlabStaticModel: case ExprNodeOutputType::matlabStaticModel:
case ExprNodeOutputType::matlabSparseStaticModel:
i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type); i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break; break;
@ -1145,9 +1151,17 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
else else
output << "x[it_" << lag << "+" << i << "*nb_row_x]"; output << "x[it_" << lag << "+" << i << "*nb_row_x]";
break; break;
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
assert(lag == 0);
[[fallthrough]];
case ExprNodeOutputType::CStaticModel: case ExprNodeOutputType::CStaticModel:
case ExprNodeOutputType::CSparseStaticModel:
case ExprNodeOutputType::juliaStaticModel: case ExprNodeOutputType::juliaStaticModel:
case ExprNodeOutputType::juliaSparseStaticModel:
case ExprNodeOutputType::matlabStaticModel: case ExprNodeOutputType::matlabStaticModel:
case ExprNodeOutputType::matlabSparseStaticModel:
output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break; break;
case ExprNodeOutputType::matlabOutsideModel: case ExprNodeOutputType::matlabOutsideModel:
@ -1206,9 +1220,17 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
else else
output << "x[it_" << lag << "+" << i << "*nb_row_x]"; output << "x[it_" << lag << "+" << i << "*nb_row_x]";
break; break;
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
assert(lag == 0);
[[fallthrough]];
case ExprNodeOutputType::CStaticModel: case ExprNodeOutputType::CStaticModel:
case ExprNodeOutputType::CSparseStaticModel:
case ExprNodeOutputType::juliaStaticModel: case ExprNodeOutputType::juliaStaticModel:
case ExprNodeOutputType::juliaSparseStaticModel:
case ExprNodeOutputType::matlabStaticModel: case ExprNodeOutputType::matlabStaticModel:
case ExprNodeOutputType::matlabSparseStaticModel:
output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
break; break;
case ExprNodeOutputType::matlabOutsideModel: case ExprNodeOutputType::matlabOutsideModel:
@ -2830,7 +2852,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << "abs"; output << "abs";
break; break;
case UnaryOpcode::sign: case UnaryOpcode::sign:
if (output_type == ExprNodeOutputType::CDynamicModel || output_type == ExprNodeOutputType::CStaticModel) if (isCOutput(output_type))
output << "copysign"; output << "copysign";
else else
output << "sign"; output << "sign";
@ -2840,6 +2862,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
switch (output_type) switch (output_type)
{ {
case ExprNodeOutputType::matlabDynamicModel: case ExprNodeOutputType::matlabDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::occbinDifferenceFile: case ExprNodeOutputType::occbinDifferenceFile:
new_output_type = ExprNodeOutputType::matlabDynamicSteadyStateOperator; new_output_type = ExprNodeOutputType::matlabDynamicSteadyStateOperator;
break; break;
@ -2847,9 +2870,11 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
new_output_type = ExprNodeOutputType::latexDynamicSteadyStateOperator; new_output_type = ExprNodeOutputType::latexDynamicSteadyStateOperator;
break; break;
case ExprNodeOutputType::CDynamicModel: case ExprNodeOutputType::CDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
new_output_type = ExprNodeOutputType::CDynamicSteadyStateOperator; new_output_type = ExprNodeOutputType::CDynamicSteadyStateOperator;
break; break;
case ExprNodeOutputType::juliaDynamicModel: case ExprNodeOutputType::juliaDynamicModel:
case ExprNodeOutputType::juliaSparseDynamicModel:
new_output_type = ExprNodeOutputType::juliaDynamicSteadyStateOperator; new_output_type = ExprNodeOutputType::juliaDynamicSteadyStateOperator;
break; break;
default: default:
@ -2931,7 +2956,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
&& arg->precedence(output_type, temporary_terms) < precedence(output_type, temporary_terms))) && arg->precedence(output_type, temporary_terms) < precedence(output_type, temporary_terms)))
{ {
output << LEFT_PAR(output_type); 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,"; output << "1.0,";
close_parenthesis = true; 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); unpackPowerDeriv()->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
else else
{ {
if (output_type == ExprNodeOutputType::juliaStaticModel || output_type == ExprNodeOutputType::juliaDynamicModel) if (isJuliaOutput(output_type))
output << "get_power_deriv("; output << "get_power_deriv(";
else else
output << "getPowerDeriv("; output << "getPowerDeriv(";

35
src/ExprNode.hh
View File

@ -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) //! Possible types of output when writing ExprNode(s) (not used for bytecode)
enum class ExprNodeOutputType enum class ExprNodeOutputType
{ {
matlabStaticModel, //!< Matlab code, static model matlabStaticModel, //!< Matlab code, static model, legacy representation
matlabDynamicModel, //!< Matlab code, dynamic model matlabDynamicModel, //!< Matlab code, dynamic model, legacy representation
CDynamicModel, //!< C code, dynamic model matlabSparseStaticModel, //!< Matlab code, static model, sparse representation
CStaticModel, //!< C code, static model matlabSparseDynamicModel, //!< Matlab code, dynamic model, sparse representation
juliaStaticModel, //!< Julia code, static model CDynamicModel, //!< C code, dynamic model, legacy representation
juliaDynamicModel, //!< Julia code, dynamic model 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) matlabOutsideModel, //!< Matlab code, outside model block (for example in initval)
latexStaticModel, //!< LaTeX code, static model latexStaticModel, //!< LaTeX code, static model
latexDynamicModel, //!< LaTeX code, dynamic model latexDynamicModel, //!< LaTeX code, dynamic model
@ -119,6 +125,8 @@ isMatlabOutput(ExprNodeOutputType output_type)
{ {
return output_type == ExprNodeOutputType::matlabStaticModel return output_type == ExprNodeOutputType::matlabStaticModel
|| output_type == ExprNodeOutputType::matlabDynamicModel || output_type == ExprNodeOutputType::matlabDynamicModel
|| output_type == ExprNodeOutputType::matlabSparseStaticModel
|| output_type == ExprNodeOutputType::matlabSparseDynamicModel
|| output_type == ExprNodeOutputType::matlabOutsideModel || output_type == ExprNodeOutputType::matlabOutsideModel
|| output_type == ExprNodeOutputType::matlabDynamicSteadyStateOperator || output_type == ExprNodeOutputType::matlabDynamicSteadyStateOperator
|| output_type == ExprNodeOutputType::steadyStateFile || output_type == ExprNodeOutputType::steadyStateFile
@ -132,6 +140,8 @@ isJuliaOutput(ExprNodeOutputType output_type)
{ {
return output_type == ExprNodeOutputType::juliaStaticModel return output_type == ExprNodeOutputType::juliaStaticModel
|| output_type == ExprNodeOutputType::juliaDynamicModel || output_type == ExprNodeOutputType::juliaDynamicModel
|| output_type == ExprNodeOutputType::juliaSparseStaticModel
|| output_type == ExprNodeOutputType::juliaSparseDynamicModel
|| output_type == ExprNodeOutputType::juliaDynamicSteadyStateOperator || output_type == ExprNodeOutputType::juliaDynamicSteadyStateOperator
|| output_type == ExprNodeOutputType::juliaSteadyStateFile || output_type == ExprNodeOutputType::juliaSteadyStateFile
|| output_type == ExprNodeOutputType::juliaTimeDataFrame; || output_type == ExprNodeOutputType::juliaTimeDataFrame;
@ -142,6 +152,8 @@ isCOutput(ExprNodeOutputType output_type)
{ {
return output_type == ExprNodeOutputType::CDynamicModel return output_type == ExprNodeOutputType::CDynamicModel
|| output_type == ExprNodeOutputType::CStaticModel || output_type == ExprNodeOutputType::CStaticModel
|| output_type == ExprNodeOutputType::CSparseDynamicModel
|| output_type == ExprNodeOutputType::CSparseStaticModel
|| output_type == ExprNodeOutputType::CDynamicSteadyStateOperator; || output_type == ExprNodeOutputType::CDynamicSteadyStateOperator;
} }
@ -162,6 +174,17 @@ isSteadyStateOperatorOutput(ExprNodeOutputType output_type)
|| output_type == ExprNodeOutputType::juliaDynamicSteadyStateOperator; || 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 constexpr bool
isAssignmentLHSBytecodeOutput(ExprNodeBytecodeOutputType output_type) isAssignmentLHSBytecodeOutput(ExprNodeBytecodeOutputType output_type)
{ {

2
src/ModFile.cc
View File

@ -1167,6 +1167,8 @@ ModFile::writeJsonOutputParsingCheck(const string &basename, JsonFileOutputType
symbol_table.writeJsonOutput(output); symbol_table.writeJsonOutput(output);
output << ", "; output << ", ";
dynamic_model.writeJsonOutput(output); dynamic_model.writeJsonOutput(output);
output << ", ";
static_model.writeJsonOutput(output);
if (!statements.empty() if (!statements.empty()
|| !var_model_table.empty() || !var_model_table.empty()

26
src/ModelTree.cc
View File

@ -71,6 +71,8 @@ ModelTree::copyHelper(const ModelTree &m)
derivatives.push_back(convert_deriv_map(it)); derivatives.push_back(convert_deriv_map(it));
for (const auto &it : m.params_derivatives) for (const auto &it : m.params_derivatives)
params_derivatives.emplace(it.first, convert_deriv_map(it.second)); 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) 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}, equation_tags{m.equation_tags},
computed_derivs_order{m.computed_derivs_order}, computed_derivs_order{m.computed_derivs_order},
NNZDerivatives{m.NNZDerivatives}, NNZDerivatives{m.NNZDerivatives},
jacobian_sparse_colptr{m.jacobian_sparse_colptr},
eq_idx_block2orig{m.eq_idx_block2orig}, eq_idx_block2orig{m.eq_idx_block2orig},
endo_idx_block2orig{m.endo_idx_block2orig}, endo_idx_block2orig{m.endo_idx_block2orig},
eq_idx_orig2block{m.eq_idx_orig2block}, eq_idx_orig2block{m.eq_idx_orig2block},
@ -177,6 +180,10 @@ ModelTree::operator=(const ModelTree &m)
NNZDerivatives = m.NNZDerivatives; NNZDerivatives = m.NNZDerivatives;
derivatives.clear(); derivatives.clear();
jacobian_sparse_column_major_order.clear();
jacobian_sparse_colptr = m.jacobian_sparse_colptr;
params_derivatives.clear(); params_derivatives.clear();
temporary_terms_derivatives.clear(); temporary_terms_derivatives.clear();
@ -902,6 +909,11 @@ ModelTree::computeDerivatives(int order, const set<int> &vars)
++NNZDerivatives[1]; ++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 // Higher-order derivatives
for (int o = 2; o <= order; o++) for (int o = 2; o <= order; o++)
for (const auto &[lower_indices, lower_d] : derivatives[o-1]) 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(); computeBlockTemporaryTerms();
block_decomposed = true; 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;
}

347
src/ModelTree.hh
View File

@ -119,6 +119,24 @@ protected:
derivatives, ... */ derivatives, ... */
vector<int> NNZDerivatives; 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 //! Derivatives with respect to parameters
/*! The key of the outer map is a pair (derivation order w.r.t. endogenous, /*! 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 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 temporary_terms_t &temporary_terms_union,
const deriv_node_temp_terms_t &tef_terms) const; 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. /* Helper for writing JSON output for residuals and derivatives.
Returns mlv and derivatives output at each derivation order. */ Returns mlv and derivatives output at each derivation order. */
template<bool dynamic> template<bool dynamic>
@ -490,6 +516,10 @@ protected:
// Returns a human-readable string describing the model class (e.g. “dynamic model”…) // Returns a human-readable string describing the model class (e.g. “dynamic model”…)
virtual string modelClassName() const = 0; 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: private:
//! Internal helper for the copy constructor and assignment operator //! Internal helper for the copy constructor and assignment operator
/*! Copies all the structures that contain ExprNode*, by the converting the /*! Copies all the structures that contain ExprNode*, by the converting the
@ -699,6 +729,8 @@ template<ExprNodeOutputType output_type>
pair<vector<ostringstream>, vector<ostringstream>> pair<vector<ostringstream>, vector<ostringstream>>
ModelTree::writeModelFileHelper() const 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> d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual)
vector<ostringstream> tt_output(derivatives.size()); // Temp terms output (at all orders) 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, writeTemporaryTerms<output_type>(temporary_terms_derivatives[1], temp_term_union,
temporary_terms_idxs, tt_output[1], tef_terms); 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); d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type)) if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type))
d_output[1] << eq + 1 << "," << getJacobianCol(var) + 1; d_output[1] << eq + 1 << "," << getJacobianCol(var, sparse) + 1;
else else
d_output[1] << eq + getJacobianCol(var)*equations.size(); d_output[1] << eq + getJacobianCol(var, sparse)*equations.size();
d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
d1->writeOutput(d_output[1], output_type, d1->writeOutput(d_output[1], output_type,
temp_term_union, temporary_terms_idxs, tef_terms); temp_term_union, temporary_terms_idxs, tef_terms);
d_output[1] << ";" << endl; d_output[1] << ";" << endl;
}
} }
} }
@ -739,58 +789,49 @@ ModelTree::writeModelFileHelper() const
writeTemporaryTerms<output_type>(temporary_terms_derivatives[i], temp_term_union, writeTemporaryTerms<output_type>(temporary_terms_derivatives[i], temp_term_union,
temporary_terms_idxs, tt_output[i], tef_terms); temporary_terms_idxs, tt_output[i], tef_terms);
/* When creating the sparse matrix (in MATLAB or C mode), since storage if constexpr(sparse)
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 eq{vidx[0]}; /* List non-zero elements of the tensor in row-major order (this is
suitable for the k-order solver according to Normann). */
int col_idx{0}; // Tensor indices are output in M_ and the JSON file (since they are constant)
for (size_t j = 1; j < vidx.size(); j++) for (int k {0};
const auto &[vidx, d] : derivatives[i])
{ {
col_idx *= getJacobianColsNbr(); d_output[i] << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
col_idx += getJacobianCol(vidx[j]); << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
} << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
if constexpr(isJuliaOutput(output_type))
{
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->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms);
d_output[i] << endl; 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;
k++; 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 for (int k{0}; // Current line index in the 3-column matrix
if (i == 2 && vidx[1] != vidx[2]) 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)) 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 else
{ {
i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) 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) j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type) << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(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) v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type) << k + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
<< "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms);
<< k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type) v_output << ";" << endl;
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl;
k++; 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)) if constexpr(isMatlabOutput(output_type))
@ -985,7 +1054,7 @@ ModelTree::writeModelCFile(const string &basename, const string &mexext,
<< endl << endl
<< " if (nlhs >= 2)" << endl << " if (nlhs >= 2)" << endl
<< " {" << 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 << " double *g1 = mxGetPr(plhs[1]);" << endl
<< " " << prefix << "g1_tt(y, x" << nb_row_x_argout << ", params" << ss_it_argout << ", T);" << 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 << " " << 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_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 << " " << 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 *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 << " 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 << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
<< " plhs[2] = plhs_sparse[0];" << 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_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 << " " << 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 *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 << " 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 << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl
<< " plhs[3] = plhs_sparse[0];" << endl << " plhs[3] = plhs_sparse[0];" << endl
@ -1130,6 +1199,8 @@ ModelTree::writeParamsDerivativesFileHelper() const
{ {
static_assert(!isCOutput(output_type), "C output is not implemented"); 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 tt_output; // Used for storing model temp vars and equations
ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters ostringstream rp_output; // 1st deriv. of residuals w.r.t. parameters
ostringstream gp_output; // 1st deriv. of Jacobian 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) }; 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 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 };
gp_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << var_col 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) }; 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 param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
int param2_col { getTypeSpecificIDByDerivID(param2) + 1 }; int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
@ -1256,8 +1327,8 @@ ModelTree::writeParamsDerivativesFileHelper() const
{ {
auto [eq, var1, var2, param] { vectorToTuple<4>(indices) }; auto [eq, var1, var2, param] { vectorToTuple<4>(indices) };
int var1_col { getJacobianCol(var1) + 1 }; int var1_col { getJacobianCol(var1, sparse) + 1 };
int var2_col { getJacobianCol(var2) + 1 }; int var2_col { getJacobianCol(var2, sparse) + 1 };
int param_col { getTypeSpecificIDByDerivID(param) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 };
hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",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) }; auto [eq, var1, var2, var3, param] { vectorToTuple<5>(indices) };
int var1_col { getJacobianCol(var1) + 1 }; int var1_col { getJacobianCol(var1, sparse) + 1 };
int var2_col { getJacobianCol(var2) + 1 }; int var2_col { getJacobianCol(var2, sparse) + 1 };
int var3_col { getJacobianCol(var3) + 1 }; int var3_col { getJacobianCol(var3, sparse) + 1 };
int param_col { getTypeSpecificIDByDerivID(param) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 };
g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1" g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
@ -1517,7 +1588,7 @@ ModelTree::writeBytecodeHelper(BytecodeWriter &code_file) const
if constexpr(dynamic) if constexpr(dynamic)
{ {
// Bytecode MEX uses a separate matrix for exogenous and exodet Jacobians // 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}; code_file << FSTPG3_{eq, tsid, lag, jacob_col};
} }
else else
@ -1763,7 +1834,7 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
d_output[0] << ", "; d_output[0] << ", ";
writeJsonModelEquations(d_output[0], true); writeJsonModelEquations(d_output[0], true);
int ncols { getJacobianColsNbr() }; int ncols { getJacobianColsNbr(false) };
for (size_t i {1}; i < derivatives.size(); i++) 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"}; 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}; int col_idx {0};
for (size_t j {1}; j < vidx.size(); j++) for (size_t j {1}; j < vidx.size(); j++)
{ {
col_idx *= getJacobianColsNbr(); col_idx *= getJacobianColsNbr(false);
col_idx += getJacobianCol(vidx[j]); col_idx += getJacobianCol(vidx[j], false);
} }
if (writeDetails) if (writeDetails)
@ -1798,7 +1869,7 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian 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; d_output[i] << ", " << col_idx_sym + 1;
} }
if (i > 1) if (i > 1)
@ -1818,7 +1889,7 @@ ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const
} }
d_output[i] << "]}"; d_output[i] << "]}";
ncols *= getJacobianColsNbr(); ncols *= getJacobianColsNbr(false);
} }
return { move(mlv_output), move(d_output) }; return { move(mlv_output), move(d_output) };
@ -1877,7 +1948,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
gp_output << R"("deriv_jacobian_wrt_params": {)" gp_output << R"("deriv_jacobian_wrt_params": {)"
<< R"( "neqs": )" << equations.size() << R"( "neqs": )" << equations.size()
<< R"(, "nvarcols": )" << getJacobianColsNbr() << R"(, "nvarcols": )" << getJacobianColsNbr(false)
<< R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "nparamcols": )" << symbol_table.param_nbr()
<< R"(, "entries": [)"; << R"(, "entries": [)";
for (bool printed_something {false}; for (bool printed_something {false};
@ -1888,7 +1959,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
auto [eq, var, param] { vectorToTuple<3>(vidx) }; 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 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 };
if (writeDetails) if (writeDetails)
@ -1948,7 +2019,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
gpp_output << R"("second_deriv_jacobian_wrt_params": {)" gpp_output << R"("second_deriv_jacobian_wrt_params": {)"
<< R"( "neqs": )" << equations.size() << R"( "neqs": )" << equations.size()
<< R"(, "nvarcols": )" << getJacobianColsNbr() << R"(, "nvarcols": )" << getJacobianColsNbr(false)
<< R"(, "nparam1cols": )" << symbol_table.param_nbr() << R"(, "nparam1cols": )" << symbol_table.param_nbr()
<< R"(, "nparam2cols": )" << symbol_table.param_nbr() << R"(, "nparam2cols": )" << symbol_table.param_nbr()
<< R"(, "entries": [)"; << R"(, "entries": [)";
@ -1960,7 +2031,7 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
auto [eq, var, param1, param2] { vectorToTuple<4>(vidx) }; 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 param1_col { getTypeSpecificIDByDerivID(param1) + 1 };
int param2_col { getTypeSpecificIDByDerivID(param2) + 1 }; int param2_col { getTypeSpecificIDByDerivID(param2) + 1 };
@ -1990,8 +2061,8 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
hp_output << R"("derivative_hessian_wrt_params": {)" hp_output << R"("derivative_hessian_wrt_params": {)"
<< R"( "neqs": )" << equations.size() << R"( "neqs": )" << equations.size()
<< R"(, "nvar1cols": )" << getJacobianColsNbr() << R"(, "nvar1cols": )" << getJacobianColsNbr(false)
<< R"(, "nvar2cols": )" << getJacobianColsNbr() << R"(, "nvar2cols": )" << getJacobianColsNbr(false)
<< R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "nparamcols": )" << symbol_table.param_nbr()
<< R"(, "entries": [)"; << R"(, "entries": [)";
for (bool printed_something {false}; for (bool printed_something {false};
@ -2002,8 +2073,8 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
auto [eq, var1, var2, param] { vectorToTuple<4>(vidx) }; auto [eq, var1, var2, param] { vectorToTuple<4>(vidx) };
int var1_col { getJacobianCol(var1) + 1 }; int var1_col { getJacobianCol(var1, false) + 1 };
int var2_col { getJacobianCol(var2) + 1 }; int var2_col { getJacobianCol(var2, false) + 1 };
int param_col { getTypeSpecificIDByDerivID(param) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 };
if (writeDetails) if (writeDetails)
@ -2036,9 +2107,9 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
{ {
g3p_output << R"("derivative_g3_wrt_params": {)" g3p_output << R"("derivative_g3_wrt_params": {)"
<< R"( "neqs": )" << equations.size() << R"( "neqs": )" << equations.size()
<< R"(, "nvar1cols": )" << getJacobianColsNbr() << R"(, "nvar1cols": )" << getJacobianColsNbr(false)
<< R"(, "nvar2cols": )" << getJacobianColsNbr() << R"(, "nvar2cols": )" << getJacobianColsNbr(false)
<< R"(, "nvar3cols": )" << getJacobianColsNbr() << R"(, "nvar3cols": )" << getJacobianColsNbr(false)
<< R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "nparamcols": )" << symbol_table.param_nbr()
<< R"(, "entries": [)"; << R"(, "entries": [)";
for (bool printed_something {false}; for (bool printed_something {false};
@ -2049,9 +2120,9 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
auto [eq, var1, var2, var3, param] { vectorToTuple<5>(vidx) }; auto [eq, var1, var2, var3, param] { vectorToTuple<5>(vidx) };
int var1_col { getJacobianCol(var1) + 1 }; int var1_col { getJacobianCol(var1, false) + 1 };
int var2_col { getJacobianCol(var2) + 1 }; int var2_col { getJacobianCol(var2, false) + 1 };
int var3_col { getJacobianCol(var3) + 1 }; int var3_col { getJacobianCol(var3, false) + 1 };
int param_col { getTypeSpecificIDByDerivID(param) + 1 }; int param_col { getTypeSpecificIDByDerivID(param) + 1 };
if (writeDetails) if (writeDetails)
@ -2084,4 +2155,98 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output) }; 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 #endif

7
src/StaticModel.cc
View File

@ -1022,6 +1022,8 @@ StaticModel::writeDriverOutput(ostream &output, bool block) const
if (block) if (block)
writeBlockDriverOutput(output); writeBlockDriverOutput(output);
writeDriverSparseIndicesHelper<false>(output);
} }
void void
@ -1301,10 +1303,7 @@ StaticModel::writeJsonAuxVarRecursiveDefinitions(ostream &output) const
void void
StaticModel::writeJsonOutput(ostream &output) const StaticModel::writeJsonOutput(ostream &output) const
{ {
deriv_node_temp_terms_t tef_terms; writeJsonSparseIndicesHelper<false>(output);
writeJsonModelLocalVariables(output, false, tef_terms);
output << ", ";
writeJsonModelEquations(output, false);
} }
void void

4
src/StaticModel.hh
View File

@ -77,12 +77,12 @@ private:
int getTypeSpecificIDByDerivID(int deriv_id) const override; int getTypeSpecificIDByDerivID(int deriv_id) const override;
int int
getJacobianCol(int deriv_id) const override getJacobianCol(int deriv_id, [[maybe_unused]] bool sparse) const override
{ {
return getTypeSpecificIDByDerivID(deriv_id); return getTypeSpecificIDByDerivID(deriv_id);
} }
int int
getJacobianColsNbr() const override getJacobianColsNbr([[maybe_unused]] bool sparse) const override
{ {
return symbol_table.endo_nbr(); return symbol_table.endo_nbr();
} }