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

Implement new sparse model representation in MATLAB/Octave 

Ref. dynare#1859
master
Sébastien Villemot 2022-09-30 17:26:43 +02:00
parent 47290087f6
commit 4ab3e937ea
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
5 changed files with 215 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/DynamicModel.cc
View File

@ -3285,6 +3285,13 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
create_directories(plusfolder);
if (block && !use_dll)
create_directories(plusfolder / "+block");
auto sparsefolder {plusfolder / "+sparse"};
create_directories(sparsefolder);
if (!use_dll)
create_directories(sparsefolder / "private");
if (block_decomposed)
create_directories(sparsefolder / "+block");
}
create_directories(model_dir / "bytecode");
@ -3321,8 +3328,14 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool use_dll,
}
// Sparse representation
if (julia)
if (use_dll)
{
// TBD
}
else if (julia)
writeSparseModelJuliaFiles<true>(basename);
else // MATLAB/Octave
writeSparseModelMFiles<true>(basename);
writeSetAuxiliaryVariables(basename, julia);
}

17
src/DynamicModel.hh
View File

@ -118,9 +118,10 @@ private:
//! Used for var_expectation and var_model
map<string, set<int>> var_expectation_functions_to_write;
//! Writes dynamic model file (Matlab version)
// Writes dynamic model file (MATLAB/Octave version, legacy representation)
void writeDynamicMFile(const string &basename) const;
//! Writes the main dynamic function of block decomposed model (MATLAB version)
/* Writes the main dynamic function of block decomposed model (MATLAB/Octave
version, legacy representation) */
void writeDynamicBlockMFile(const string &basename) const;
/* Writes the main dynamic functions of block decomposed model (C version),
then compiles it with the per-block functions into a single MEX */
@ -131,7 +132,8 @@ private:
// 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)
/* Writes the per-block dynamic files of block decomposed model (MATLAB/Octave
version, legacy representation) */
void writeDynamicPerBlockMFiles(const string &basename) const;
/* Writes and compiles the per-block dynamic files of block decomposed model
(C version). Returns the list of paths to the compiled object files. */
@ -206,9 +208,11 @@ private:
//! Write reverse cross references
void writeRevXrefs(ostream &output, const map<pair<int, int>, set<int>> &xrefmap, const string &type) const;
// Writes MATLAB/Octave wrapper function for computing residuals and derivatives at the same time
/* Writes MATLAB/Octave wrapper function for computing residuals and
derivatives at the same time (legacy representation) */
void writeDynamicMWrapperFunction(const string &name, const string &ending) const;
// Helper for writing MATLAB/Octave functions for residuals/derivatives and their temporary terms
/* Helper for writing MATLAB/Octave functions for residuals/derivatives and
their temporary terms (legacy representation) */
void writeDynamicMFileHelper(const string &basename,
const string &name, const string &retvalname,
const string &name_tt, size_t ttlen,
@ -216,7 +220,8 @@ private:
const ostringstream &init_s, const ostringstream &end_s,
const ostringstream &s, const ostringstream &s_tt) const;
//! Create a legacy *_dynamic.m file for MATLAB/Octave not yet using the temporary terms array interface
/* Create the compatibility dynamic.m file for MATLAB/Octave not yet using
the temporary terms array interface (legacy representation) */
void writeDynamicMCompatFile(const string &basename) const;
//! Internal helper for the copy constructor and assignment operator

165
src/ModelTree.hh
View File

@ -375,6 +375,10 @@ protected:
template<ExprNodeBytecodeOutputType output_type>
void writeBytecodeModelEquations(BytecodeWriter &code_file, const temporary_terms_t &temporary_terms, const deriv_node_temp_terms_t &tef_terms) const;
// Writes the sparse representation of the model in MATLAB/Octave
template<bool dynamic>
void writeSparseModelMFiles(const string &basename) const;
// Writes the sparse representation of the model in Julia
// Assumes that the directory <MODFILE>/model/julia/ already exists
template<bool dynamic>
@ -2445,3 +2449,164 @@ ModelTree::writeSparseModelJuliaFiles(const string &basename) const
}
}
#endif
template<bool dynamic>
void
ModelTree::writeSparseModelMFiles(const string &basename) const
{
constexpr ExprNodeOutputType output_type {dynamic ? ExprNodeOutputType::matlabSparseDynamicModel : ExprNodeOutputType::matlabSparseStaticModel};
auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper<output_type>();
const filesystem::path m_dir {packageDir(basename) / "+sparse"};
const filesystem::path private_m_dir {m_dir / "private"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix { dynamic ? "dynamic_" : "static_" };
const string full_prefix { basename + ".sparse." + prefix };
const string ss_arg { dynamic ? ", steady_state" : "" };
size_t ttlen {temporary_terms_derivatives[0].size()};
ofstream output;
auto open_file = [&output](const filesystem::path &p)
{
output.open(p, ios::out | ios::binary);
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << p.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
};
// Residuals (non-block)
open_file(private_m_dir / (prefix + "resid_tt.m"));
output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << ss_arg << ", T_order, T)" << endl
<< "if isempty(T_order)" << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "else if T_order >= 0" << endl
<< " return" << endl
<< "end" << endl
<< "T_order = 0;" << endl
<< "if size(T, 1) < " << ttlen << endl
<< " T = resize(T, " << ttlen << ", 1);" << endl
<< "end" << endl
<< tt_sparse_output[0].str()
<< "end" << endl;
output.close();
open_file(m_dir / (prefix + "resid.m"));
output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << ss_arg << ", T_order, T)" << endl
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
<< "residual = NaN(" << equations.size() << ", 1);" << endl
<< d_sparse_output[0].str();
if constexpr(!dynamic)
output << "if ~isreal(residual)" << endl
<< " residual = real(residual)+imag(residual).^2;" << endl
<< "end" << endl;
output << "end" << endl;
output.close();
// Jacobian (non-block)
ttlen += temporary_terms_derivatives[1].size();
open_file(private_m_dir / (prefix + "g1_tt.m"));
output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T)" << endl
<< "if isempty(T_order)" << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "else if T_order >= 1" << endl
<< " return" << endl
<< "end" << endl
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
<< "T_order = 1;" << endl
<< "if size(T, 1) < " << ttlen << endl
<< " T = resize(T, " << ttlen << ", 1);" << endl
<< "end" << endl
<< tt_sparse_output[1].str()
<< "end" << endl;
output.close();
open_file(m_dir / (prefix + "g1.m"));
// NB: At first order, sparse indices are passed as extra arguments
output << "function [g1, T_order, T] = " << prefix << "g1(y, x, params" << ss_arg << ", sparse_rowval, sparse_colval, sparse_colptr, T_order, T)" << endl
<< "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
<< "g1_v = NaN(" << jacobian_sparse_column_major_order.size() << ", 1);" << endl
<< d_sparse_output[1].str();
if constexpr(!dynamic)
output << "if ~isreal(g1_v)" << endl
<< " g1_v = real(g1_v)+2*imag(g1_v);" << endl
<< "end" << endl;
output << "g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << equations.size() << ", " << getJacobianColsNbr(true) << ");" << endl
<< "end" << endl;
output.close();
// Higher-order derivatives (non-block)
for (int i {2}; i <= computed_derivs_order; i++)
{
ttlen += temporary_terms_derivatives[i].size();
open_file(private_m_dir / (prefix + "g" + to_string(i) + "_tt.m"));
output << "function T = " << prefix << "g" << i << "_tt(y, x, params" << ss_arg << ")" << endl
<< "if isempty(T_order)" << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "else if T_order >= " << i << endl
<< " return" << endl
<< "end" << endl
<< "[T_order, T] = " << full_prefix << "g" << i-1 << "_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
<< "T_order = " << i << ";" << endl
<< "if size(T, 1) < " << ttlen << endl
<< " T = resize(T, " << ttlen << ", 1);" << endl
<< "end" << endl
<< tt_sparse_output[i].str()
<< "end" << endl;
output.close();
open_file(m_dir / (prefix + "g" + to_string(i) + ".m"));
output << "function [g" << i << "_v, T_order, T] = " << prefix << "g" << i << "(y, x, params" << ss_arg << ", T_order, T)" << endl
<< "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << ss_arg << ", T_order, T);" << endl
<< "g" << i << "_v = NaN(" << derivatives[i].size() << ", 1);" << endl
<< d_sparse_output[i].str()
<< "end" << endl;
output.close();
}
// Block decomposition
if (block_decomposed)
{
const filesystem::path block_dir {m_dir / "+block"};
temporary_terms_t temporary_terms_written;
for (int blk {0}; blk < static_cast<int>(blocks.size()); blk++)
{
const string funcname {prefix + to_string(blk+1)};
const BlockSimulationType simulation_type {blocks[blk].simulation_type};
const bool evaluate {simulation_type == BlockSimulationType::evaluateForward
|| simulation_type == BlockSimulationType::evaluateBackward};
const string resid_g1_arg {evaluate ? "" : ", residual, g1"};
open_file(block_dir / (funcname + ".m"));
output << "function [y, T" << resid_g1_arg << "] = " << funcname << "(y, x, params" << ss_arg << ", sparse_rowval, sparse_colval, sparse_colptr, T)" << endl;
if (!evaluate)
output << "residual=NaN(" << blocks[blk].mfs_size << ", 1);" << endl;
// Write residuals and temporary terms (incl. those for derivatives)
writePerBlockHelper<output_type>(blk, output, temporary_terms_written);
// Write Jacobian
if (!evaluate)
{
const bool one_boundary {simulation_type == BlockSimulationType::solveBackwardSimple
|| simulation_type == BlockSimulationType::solveForwardSimple
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete};
output << "if nargout > 3" << endl
<< " g1_v = NaN(" << blocks_jacobian_sparse_column_major_order[blk].size() << ", 1);" << endl;
writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
output << " g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << blocks[blk].mfs_size << ", "
<< (one_boundary ? 1 : 3)*blocks[blk].mfs_size << ");" << endl
<< "end" << endl;
}
output << "end" << endl;
output.close();
}
}
}

15
src/StaticModel.cc
View File

@ -628,6 +628,13 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
create_directories(plusfolder);
if (block && !use_dll)
create_directories(plusfolder / "+block");
auto sparsefolder {plusfolder / "+sparse"};
create_directories(sparsefolder);
if (!use_dll)
create_directories(sparsefolder / "private");
if (block_decomposed)
create_directories(sparsefolder / "+block");
}
create_directories(model_dir / "bytecode");
@ -664,8 +671,14 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool use_dll, c
}
// Sparse representation
if (julia)
if (use_dll)
{
// TBD
}
else if (julia)
writeSparseModelJuliaFiles<false>(basename);
else // MATLAB/Octave
writeSparseModelMFiles<false>(basename);
writeSetAuxiliaryVariables(basename, julia);
}

17
src/StaticModel.hh
View File

@ -34,10 +34,11 @@ class DynamicModel;
class StaticModel : public ModelTree
{
private:
//! Writes static model file (standard Matlab version)
// Writes static model file (MATLAB/Octave version, legacy representation)
void writeStaticMFile(const string &basename) const;
//! Writes the main static function of block decomposed model (MATLAB version)
/* Writes the main static function of block decomposed model (MATLAB/Octave
version, legacy representation) */
void writeStaticBlockMFile(const string &basename) const;
/* Writes the main static functions of block decomposed model (C version),
@ -48,7 +49,8 @@ private:
template<ExprNodeOutputType output_type>
void writeStaticPerBlockHelper(int blk, ostream &output, temporary_terms_t &temporary_terms) const;
//! Writes the per-block static files of block decomposed model (MATLAB version)
/* Writes the per-block static files of block decomposed model (MATLAB/Octave
version, legacy representation) */
void writeStaticPerBlockMFiles(const string &basename) const;
/* Writes and compiles the per-block static files of block decomposed model
@ -86,17 +88,20 @@ private:
void computeChainRuleJacobian() override;
// Helper for writing MATLAB/Octave functions for residuals/derivatives and their temporary terms
/* Helper for writing MATLAB/Octave functions for residuals/derivatives and
their temporary terms (legacy representation) */
void writeStaticMFileHelper(const string &basename,
const string &name, const string &retvalname,
const string &name_tt, size_t ttlen,
const string &previous_tt_name,
const ostringstream &init_s, const ostringstream &end_s,
const ostringstream &s, const ostringstream &s_tt) const;
// Writes MATLAB/Octave wrapper function for computing residuals and derivatives at the same time
/* Writes MATLAB/Octave wrapper function for computing residuals and
derivatives at the same time (legacy representation) */
void writeStaticMWrapperFunction(const string &basename, const string &ending) const;
//! Create a legacy *_static.m file for MATLAB/Octave not yet using the temporary terms array interface
/* Create the compatibility static.m file for MATLAB/Octave not yet using the
temporary terms array interface (legacy representation) */
void writeStaticMCompatFile(const string &name) const;
int