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

Provisions for new sparse representation of block-decomposed dynamic/static files 

The stochastic mode in currently unsupported.

This commit adds new fields in M_.

This is a preliminary step for dynare#1859.
master
Sébastien Villemot 2022-09-28 19:18:15 +02:00
parent 1ed72f6da2
commit 47290087f6
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 136 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
src/DynamicModel.cc
View File

@ -1461,6 +1461,8 @@ DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
<< "M_.block_structure.block(" << blk+1 << ").NNZDerivatives = " << blocks_derivatives[blk].size() << ';' << endl;
}
writeBlockDriverSparseIndicesHelper<true>(output);
output << "M_.block_structure.variable_reordered = [";
for (int i = 0; i < symbol_table.endo_nbr(); i++)
output << " " << endo_idx_block2orig[i]+1;
@ -3101,10 +3103,14 @@ DynamicModel::computeChainRuleJacobian()
size_t nb_blocks { blocks.size() };
blocks_derivatives.resize(nb_blocks);
blocks_jacobian_sparse_column_major_order.resize(nb_blocks);
blocks_jacobian_sparse_colptr.resize(nb_blocks);
for (int blk {0}; blk < static_cast<int>(nb_blocks); blk++)
{
int nb_recursives = blocks[blk].getRecursiveSize();
int mfs_size { blocks[blk].mfs_size };
BlockSimulationType simulation_type { blocks[blk].simulation_type };
// Create a map from recursive vars to their defining (normalized) equation
map<int, BinaryOpNode *> recursive_vars;
@ -3144,6 +3150,25 @@ DynamicModel::computeChainRuleJacobian()
if (d != Zero)
blocks_derivatives[blk][{ eq, var, lag }] = d;
}
// Compute the sparse representation of the Jacobian
if (simulation_type != BlockSimulationType::evaluateForward
&& simulation_type != BlockSimulationType::evaluateBackward)
{
const bool one_boundary {simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete};
for (const auto &[indices, d1] : blocks_derivatives[blk])
{
auto &[eq, var, lag] { indices };
assert(lag >= -1 && lag <= 1);
if (eq >= nb_recursives && var >= nb_recursives
&& !(one_boundary && lag != 0))
blocks_jacobian_sparse_column_major_order[blk].emplace(pair{eq-nb_recursives, var-nb_recursives+static_cast<int>(!one_boundary)*(lag+1)*mfs_size}, d1);
}
blocks_jacobian_sparse_colptr[blk] = computeCSCColPtr(blocks_jacobian_sparse_column_major_order[blk], (one_boundary ? 1 : 3)*mfs_size);
}
}
/* Also store information and derivatives w.r.t. other types of variables

4
src/DynamicModel.hh
View File

@ -128,7 +128,7 @@ private:
/* Computes the number of nonzero elements in deterministic Jacobian of
block-decomposed model */
int nzeDeterministicJacobianForBlock(int blk) const;
//! Helper for writing the per-block dynamic files of block decomposed models
// Helper for writing the per-block dynamic files of block decomposed models (legacy representation)
template<ExprNodeOutputType output_type>
void writeDynamicPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms, int nze_stochastic, int nze_deterministic, int nze_exo, int nze_exo_det, int nze_other_endo) const;
//! Writes the per-block dynamic files of block decomposed model (MATLAB version)
@ -695,6 +695,8 @@ DynamicModel::writeDynamicPerBlockHelper(int blk, ostream &output, temporary_ter
int nze_stochastic, int nze_deterministic, int nze_exo,
int nze_exo_det, int nze_other_endo) const
{
static_assert(!isSparseModelOutput(output_type));
BlockSimulationType simulation_type { blocks[blk].simulation_type };
int block_mfs_size { blocks[blk].mfs_size };
int block_recursive_size { blocks[blk].getRecursiveSize() };

11
src/ModelTree.cc
View File

@ -125,6 +125,14 @@ ModelTree::copyHelper(const ModelTree &m)
blocks_temporary_terms.push_back(convert_vector_tt(it));
for (const auto &it : m.blocks_temporary_terms_idxs)
blocks_temporary_terms_idxs.emplace(f(it.first), it.second);
for (const auto &it : m.blocks_jacobian_sparse_column_major_order)
{
map<pair<int, int>, expr_t, columnMajorOrderLess> v;
for (const auto &it2 : it)
v.emplace(it2.first, f(it2.second));
blocks_jacobian_sparse_column_major_order.push_back(v);
}
}
ModelTree::ModelTree(SymbolTable &symbol_table_arg,
@ -160,6 +168,7 @@ ModelTree::ModelTree(const ModelTree &m) :
blocks{m.blocks},
endo2block{m.endo2block},
eq2block{m.eq2block},
blocks_jacobian_sparse_colptr{m.blocks_jacobian_sparse_colptr},
endo2eq{m.endo2eq},
cutoff{m.cutoff},
mfs{m.mfs}
@ -204,6 +213,8 @@ ModelTree::operator=(const ModelTree &m)
eq2block = m.eq2block;
blocks_temporary_terms.clear();
blocks_temporary_terms_idxs.clear();
blocks_jacobian_sparse_column_major_order.clear();
blocks_jacobian_sparse_colptr = m.blocks_jacobian_sparse_colptr;
endo2eq = m.endo2eq;
cutoff = m.cutoff;
mfs = m.mfs;

72
src/ModelTree.hh
View File

@ -238,6 +238,18 @@ protected:
the vector of all temporary terms */
temporary_terms_idxs_t blocks_temporary_terms_idxs;
/* The nonzero values of the sparse block Jacobians in column-major order
(which is the natural order for Compressed Sparse Column (CSC) storage).
The pair of indices is (row, column). Nothing is stored for blocks of
type “evaluate”, since their Jacobian is not used in the sparse
representation (since the latter does not support the stochastic mode,
which only exists in the legacy representation). */
vector<SparseColumnMajorOrderMatrix> blocks_jacobian_sparse_column_major_order;
/* Column indices for the sparse block Jacobian in Compressed Sparse Column
(CSC) storage (corresponds to the “jc” vector in MATLAB terminology).
Same remark as above regarding blocks of type “evaluate”. */
vector<vector<int>> blocks_jacobian_sparse_colptr;
//! Computes derivatives
/*! \param order the derivation order
\param vars the derivation IDs w.r.t. which compute the derivatives */
@ -296,9 +308,15 @@ protected:
void writeModelCFile(const string &basename, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
// Writes per-block residuals and temporary terms (incl. for derivatives)
// This part is common to the legacy and sparse representations
template<ExprNodeOutputType output_type>
void writePerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
// Writes per-block Jacobian (sparse representation)
// Assumes temporary terms for derivatives are already set.
template<ExprNodeOutputType output_type>
void writeSparsePerBlockJacobianHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
/* Helper for writing derivatives w.r.t. parameters.
Returns { tt, rp, gp, rpp, gpp, hp, g3p }.
g3p is empty if requesting a static output type. */
@ -329,6 +347,10 @@ protected:
template<bool dynamic>
void writeJsonSparseIndicesHelper(ostream &output) const;
// Helper for writing sparse block derivatives indices in MATLAB/Octave driver file
template<bool dynamic>
void writeBlockDriverSparseIndicesHelper(ostream &output) const;
/* Helper for writing JSON output for residuals and derivatives.
Returns mlv and derivatives output at each derivation order. */
template<bool dynamic>
@ -2254,6 +2276,56 @@ ModelTree::writeJsonSparseIndicesHelper(ostream &output) const
}
}
template<bool dynamic>
void
ModelTree::writeBlockDriverSparseIndicesHelper(ostream &output) const
{
for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
{
const string struct_name { "M_.block_structure"s + (dynamic ? "" : "_stat") + ".block(" + to_string(blk+1) + ")." };
// Write indices for the sparse Jacobian (both naive and CSC storage)
output << struct_name << "g1_sparse_rowval = int32([";
for (const auto &[indices, d1] : blocks_jacobian_sparse_column_major_order.at(blk))
output << indices.first+1 << ' ';
output << "]);" << endl
<< struct_name << "g1_sparse_colval = int32([";
for (const auto &[indices, d1] : blocks_jacobian_sparse_column_major_order.at(blk))
output << indices.second+1 << ' ';
output << "]);" << endl
<< struct_name << "g1_sparse_colptr = int32([";
for (int it : blocks_jacobian_sparse_colptr.at(blk))
output << it+1 << ' ';
output << "]);" << endl;
}
}
template<ExprNodeOutputType output_type>
void
ModelTree::writeSparsePerBlockJacobianHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
{
static_assert(isSparseModelOutput(output_type));
// NB: stochastic mode is currently unsupported by sparse representation
/* See also the comment above the definition of
blocks_jacobian_sparse_column_major_order and
blocks_jacobian_sparse_column_major_colptr */
assert(blocks[blk].simulation_type != BlockSimulationType::evaluateForward
&& blocks[blk].simulation_type != BlockSimulationType::evaluateBackward);
int k {0};
for (const auto &[row_col, d1] : blocks_jacobian_sparse_column_major_order[blk])
{
output << "g1_v" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< k + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
d1->writeOutput(output, output_type, temporary_terms, blocks_temporary_terms_idxs);
output << ";" << endl;
k++;
}
}
template<bool dynamic>
void
ModelTree::writeSparseModelJuliaFiles(const string &basename) const

24
src/StaticModel.cc
View File

@ -846,6 +846,8 @@ StaticModel::writeBlockDriverOutput(ostream &output) const
output << "];" << endl
<< "M_.block_structure_stat.tmp_nbr = " << blocks_temporary_terms_idxs.size()
<< ";" << endl;
writeBlockDriverSparseIndicesHelper<false>(output);
}
SymbolType
@ -911,10 +913,15 @@ void
StaticModel::computeChainRuleJacobian()
{
int nb_blocks = blocks.size();
blocks_derivatives.resize(nb_blocks);
blocks_jacobian_sparse_column_major_order.resize(nb_blocks);
blocks_jacobian_sparse_colptr.resize(nb_blocks);
for (int blk = 0; blk < nb_blocks; blk++)
{
int nb_recursives = blocks[blk].getRecursiveSize();
BlockSimulationType simulation_type { blocks[blk].simulation_type };
map<int, BinaryOpNode *> recursive_vars;
for (int i = 0; i < nb_recursives; i++)
@ -926,8 +933,8 @@ StaticModel::computeChainRuleJacobian()
recursive_vars[deriv_id] = getBlockEquationExpr(blk, i);
}
assert(blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesSimple
&& blocks[blk].simulation_type != BlockSimulationType::solveTwoBoundariesComplete);
assert(simulation_type != BlockSimulationType::solveTwoBoundariesSimple
&& simulation_type != BlockSimulationType::solveTwoBoundariesComplete);
int size = blocks[blk].size;
@ -942,6 +949,19 @@ StaticModel::computeChainRuleJacobian()
blocks_derivatives[blk][{ eq, var, 0 }] = d1;
}
}
// Compute the sparse representation of the Jacobian
if (simulation_type != BlockSimulationType::evaluateForward
&& simulation_type != BlockSimulationType::evaluateBackward)
{
for (const auto &[indices, d1] : blocks_derivatives[blk])
{
auto &[eq, var, lag] { indices };
assert(lag == 0);
blocks_jacobian_sparse_column_major_order[blk].emplace(pair{eq, var}, d1);
}
blocks_jacobian_sparse_colptr[blk] = computeCSCColPtr(blocks_jacobian_sparse_column_major_order[blk], size);
}
}
}

4
src/StaticModel.hh
View File

@ -44,7 +44,7 @@ private:
then compiles it with the per-block functions into a single MEX */
void writeStaticBlockCFile(const string &basename, vector<filesystem::path> per_block_object_files, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot) const;
//! Helper for writing a per-block static file of block decomposed model
// Helper for writing a per-block static file of block decomposed model (legacy representation)
template<ExprNodeOutputType output_type>
void writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
@ -178,6 +178,8 @@ template<ExprNodeOutputType output_type>
void
StaticModel::writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const
{
static_assert(!isSparseModelOutput(output_type));
BlockSimulationType simulation_type { blocks[blk].simulation_type };
int block_recursive_size { blocks[blk].getRecursiveSize() };