/*
* Copyright © 2003-2023 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see .
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "DynamicModel.hh"
#include "ParsingDriver.hh"
void
DynamicModel::copyHelper(const DynamicModel& m)
{
auto f = [this](const ExprNode* e) { return e->clone(*this); };
for (const auto& it : m.static_only_equations)
static_only_equations.push_back(dynamic_cast(f(it)));
}
DynamicModel::DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg) :
ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
trend_component_model_table {trend_component_model_table_arg},
var_model_table {var_model_table_arg}
{
}
DynamicModel::DynamicModel(const DynamicModel& m) :
ModelTree {m},
trend_component_model_table {m.trend_component_model_table},
var_model_table {m.var_model_table},
balanced_growth_test_tol {m.balanced_growth_test_tol},
static_only_equations_lineno {m.static_only_equations_lineno},
static_only_equations_equation_tags {m.static_only_equations_equation_tags},
deriv_id_table {m.deriv_id_table},
inv_deriv_id_table {m.inv_deriv_id_table},
dyn_jacobian_cols_table {m.dyn_jacobian_cols_table},
dyn_jacobian_ncols {m.dyn_jacobian_ncols},
max_lag {m.max_lag},
max_lead {m.max_lead},
max_endo_lag {m.max_endo_lag},
max_endo_lead {m.max_endo_lead},
max_exo_lag {m.max_exo_lag},
max_exo_lead {m.max_exo_lead},
max_exo_det_lag {m.max_exo_det_lag},
max_exo_det_lead {m.max_exo_det_lead},
max_lag_orig {m.max_lag_orig},
max_lead_orig {m.max_lead_orig},
max_lag_with_diffs_expanded_orig {m.max_lag_with_diffs_expanded_orig},
max_endo_lag_orig {m.max_endo_lag_orig},
max_endo_lead_orig {m.max_endo_lead_orig},
max_exo_lag_orig {m.max_exo_lag_orig},
max_exo_lead_orig {m.max_exo_lead_orig},
max_exo_det_lag_orig {m.max_exo_det_lag_orig},
max_exo_det_lead_orig {m.max_exo_det_lead_orig},
xrefs {m.xrefs},
xref_param {m.xref_param},
xref_endo {m.xref_endo},
xref_exo {m.xref_exo},
xref_exo_det {m.xref_exo_det},
nonzero_hessian_eqs {m.nonzero_hessian_eqs},
variableMapping {m.variableMapping},
blocks_jacob_cols_endo {m.blocks_jacob_cols_endo},
var_expectation_functions_to_write {m.var_expectation_functions_to_write},
mfs {m.mfs},
static_mfs {m.static_mfs}
{
copyHelper(m);
}
DynamicModel&
DynamicModel::operator=(const DynamicModel& m)
{
ModelTree::operator=(m);
assert(&trend_component_model_table == &m.trend_component_model_table);
assert(&var_model_table == &m.var_model_table);
balanced_growth_test_tol = m.balanced_growth_test_tol;
static_only_equations_lineno = m.static_only_equations_lineno;
static_only_equations_equation_tags = m.static_only_equations_equation_tags;
deriv_id_table = m.deriv_id_table;
inv_deriv_id_table = m.inv_deriv_id_table;
dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
dyn_jacobian_ncols = m.dyn_jacobian_ncols;
max_lag = m.max_lag;
max_lead = m.max_lead;
max_endo_lag = m.max_endo_lag;
max_endo_lead = m.max_endo_lead;
max_exo_lag = m.max_exo_lag;
max_exo_lead = m.max_exo_lead;
max_exo_det_lag = m.max_exo_det_lag;
max_exo_det_lead = m.max_exo_det_lead;
max_lag_orig = m.max_lag_orig;
max_lead_orig = m.max_lead_orig;
max_lag_with_diffs_expanded_orig = m.max_lag_with_diffs_expanded_orig;
max_endo_lag_orig = m.max_endo_lag_orig;
max_endo_lead_orig = m.max_endo_lead_orig;
max_exo_lag_orig = m.max_exo_lag_orig;
max_exo_lead_orig = m.max_exo_lead_orig;
max_exo_det_lag_orig = m.max_exo_det_lag_orig;
max_exo_det_lead_orig = m.max_exo_det_lead_orig;
xrefs = m.xrefs;
xref_param = m.xref_param;
xref_endo = m.xref_endo;
xref_exo = m.xref_exo;
xref_exo_det = m.xref_exo_det;
nonzero_hessian_eqs = m.nonzero_hessian_eqs;
variableMapping = m.variableMapping;
blocks_jacob_cols_endo = m.blocks_jacob_cols_endo;
var_expectation_functions_to_write = m.var_expectation_functions_to_write;
mfs = m.mfs;
static_mfs = m.static_mfs;
copyHelper(m);
return *this;
}
void
DynamicModel::writeDynamicBytecode(const string& basename) const
{
/* Bytecode only works when there are with as many endogenous as equations.
(e.g. the constructor of FBEGINBLOCK_ makes this assumption) */
assert(static_cast(equations.size()) == symbol_table.endo_nbr());
// Determine the type of model (used for typing the single block)
BlockSimulationType simulation_type;
if (max_endo_lag > 0 && max_endo_lead > 0)
simulation_type = BlockSimulationType::solveTwoBoundariesComplete;
else if (max_endo_lag >= 0 && max_endo_lead == 0)
simulation_type = BlockSimulationType::solveForwardComplete;
else
simulation_type = BlockSimulationType::solveBackwardComplete;
// First write the .bin file
int u_count_int {
writeBytecodeBinFile(basename + "/model/bytecode/dynamic.bin",
simulation_type == BlockSimulationType::solveTwoBoundariesComplete)};
BytecodeWriter code_file {basename + "/model/bytecode/dynamic.cod"};
// Declare temporary terms
code_file << FDIMT_ {static_cast(temporary_terms_derivatives[0].size()
+ temporary_terms_derivatives[1].size())};
// Declare the (single) block
vector exo(symbol_table.exo_nbr()), exo_det(symbol_table.exo_det_nbr());
iota(exo.begin(), exo.end(), 0);
iota(exo_det.begin(), exo_det.end(), 0);
int jacobian_ncols_endo {static_cast(count_if(
dyn_jacobian_cols_table.begin(), dyn_jacobian_cols_table.end(),
[this](const auto& v) { return getTypeByDerivID(v.first) == SymbolType::endogenous; }))};
vector eq_idx(equations.size());
iota(eq_idx.begin(), eq_idx.end(), 0);
vector endo_idx(symbol_table.endo_nbr());
iota(endo_idx.begin(), endo_idx.end(), 0);
code_file << FBEGINBLOCK_ {symbol_table.endo_nbr(),
simulation_type,
0,
symbol_table.endo_nbr(),
endo_idx,
eq_idx,
false,
u_count_int,
jacobian_ncols_endo,
symbol_table.exo_det_nbr(),
symbol_table.exo_nbr(),
exo_det,
exo};
writeBytecodeHelper(code_file);
}
void
DynamicModel::writeDynamicBlockBytecode(const string& basename) const
{
BytecodeWriter code_file {basename + "/model/bytecode/block/dynamic.cod"};
const filesystem::path bin_filename {basename + "/model/bytecode/block/dynamic.bin"};
ofstream bin_file {bin_filename, ios::out | ios::binary};
if (!bin_file.is_open())
{
cerr << R"(Error : Can't open file ")" << bin_filename.string() << R"(" for writing)" << endl;
exit(EXIT_FAILURE);
}
// Temporary variables declaration
code_file << FDIMT_ {static_cast(blocks_temporary_terms_idxs.size())};
temporary_terms_t temporary_terms_written;
for (int block {0}; block < static_cast(blocks.size()); block++)
{
const BlockSimulationType simulation_type {blocks[block].simulation_type};
// Write section of .bin file except for evaluate blocks and solve simple blocks
const int u_count {simulation_type == BlockSimulationType::solveTwoBoundariesSimple
|| simulation_type
== BlockSimulationType::solveTwoBoundariesComplete
|| simulation_type == BlockSimulationType::solveBackwardComplete
|| simulation_type == BlockSimulationType::solveForwardComplete
? writeBlockBytecodeBinFile(bin_file, block)
: 0};
code_file << FBEGINBLOCK_ {blocks[block].mfs_size,
simulation_type,
blocks[block].first_equation,
blocks[block].size,
endo_idx_block2orig,
eq_idx_block2orig,
blocks[block].linear,
u_count,
static_cast(blocks_jacob_cols_endo[block].size())};
writeBlockBytecodeHelper(code_file, block, temporary_terms_written);
}
code_file << FEND_ {};
}
void
DynamicModel::writeDynamicMFile(const string& basename) const
{
auto [d_output, tt_output] = writeModelFileHelper();
ostringstream init_output, end_output;
init_output << "residual = zeros(" << equations.size() << ", 1);";
writeDynamicMFileHelper(basename, "dynamic_resid", "residual", "dynamic_resid_tt",
temporary_terms_derivatives[0].size(), "", init_output, end_output,
d_output[0], tt_output[0]);
init_output.str("");
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(false)};
int ntt {static_cast(temporary_terms_derivatives[0].size()
+ temporary_terms_derivatives[1].size())};
for (size_t i {2}; i < derivatives.size(); i++)
{
ncols *= getJacobianColsNbr(false);
ntt += temporary_terms_derivatives[i].size();
string gname {"g" + to_string(i)};
string gprevname {"g" + to_string(i - 1)};
init_output.str("");
end_output.str("");
if (derivatives[i].size())
{
init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl
<< gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl
<< gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl;
end_output << gname << " = sparse(" << gname << "_i," << gname << "_j," << gname << "_v,"
<< equations.size() << "," << ncols << ");";
}
else
init_output << gname << " = sparse([],[],[]," << equations.size() << "," << ncols << ");";
writeDynamicMFileHelper(basename, "dynamic_" + gname, gname, "dynamic_" + gname + "_tt", ntt,
"dynamic_" + gprevname + "_tt", init_output, end_output, d_output[i],
tt_output[i]);
if (i <= 3)
writeDynamicMWrapperFunction(basename, gname);
}
writeDynamicMCompatFile(basename);
}
string
DynamicModel::reform(const string& name1) const
{
string name = name1;
int pos = name.find('\\', 0);
while (pos >= 0)
{
if (name.substr(pos + 1, 1) != R"(\)")
{
name = name.insert(pos, R"(\)");
pos++;
}
pos++;
pos = name.find('\\', pos);
}
return name;
}
void
DynamicModel::printNonZeroHessianEquations(ostream& output) const
{
if (nonzero_hessian_eqs.size() != 1)
output << "[";
for (bool printed_something {false}; int it : nonzero_hessian_eqs)
{
if (exchange(printed_something, true))
output << " ";
output << it + 1;
}
if (nonzero_hessian_eqs.size() != 1)
output << "]";
}
void
DynamicModel::writeDynamicMWrapperFunction(const string& basename, const string& ending) const
{
string name;
if (ending == "g1")
name = "dynamic_resid_g1";
else if (ending == "g2")
name = "dynamic_resid_g1_g2";
else if (ending == "g3")
name = "dynamic_resid_g1_g2_g3";
filesystem::path filename {packageDir(basename) / (name + ".m")};
ofstream output {filename, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
if (ending == "g1")
output << "function [residual, g1] = " << name << "(T, y, x, params, steady_state, it_, T_flag)"
<< endl
<< "% function [residual, g1] = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl;
else if (ending == "g2")
output << "function [residual, g1, g2] = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function [residual, g1, g2] = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl;
else if (ending == "g3")
output << "function [residual, g1, g2, g3] = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function [residual, g1, g2, g3] = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl;
output << "%" << endl
<< "% Wrapper function automatically created by Dynare" << endl
<< "%" << endl
<< endl
<< " if T_flag" << endl
<< " T = " << basename << ".dynamic_" << ending
<< "_tt(T, y, x, params, steady_state, it_);" << endl
<< " end" << endl;
if (ending == "g1")
output << " residual = " << basename
<< ".dynamic_resid(T, y, x, params, steady_state, it_, false);" << endl
<< " g1 = " << basename
<< ".dynamic_g1(T, y, x, params, steady_state, it_, false);" << endl;
else if (ending == "g2")
output << " [residual, g1] = " << basename
<< ".dynamic_resid_g1(T, y, x, params, steady_state, it_, false);" << endl
<< " g2 = " << basename
<< ".dynamic_g2(T, y, x, params, steady_state, it_, false);" << endl;
else if (ending == "g3")
output << " [residual, g1, g2] = " << basename
<< ".dynamic_resid_g1_g2(T, y, x, params, steady_state, it_, false);" << endl
<< " g3 = " << basename
<< ".dynamic_g3(T, y, x, params, steady_state, it_, false);" << endl;
output << endl << "end" << endl;
output.close();
}
void
DynamicModel::writeDynamicMFileHelper(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
{
filesystem::path filename {packageDir(basename) / (name_tt + ".m")};
ofstream output {filename, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
<< "% function T = " << name_tt << "(T, y, x, params, steady_state, it_)" << endl
<< "%" << endl
<< "% File created by Dynare Preprocessor from .mod file" << endl
<< "%" << endl
<< "% Inputs:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms to be "
"filled by function"
<< endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables "
"in the order stored"
<< endl
<< "% in M_.lead_lag_incidence; see "
"the Manual"
<< endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables "
"(in declaration order)"
<< endl
<< "% for all simulation periods"
<< endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values"
<< endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in "
"declaration order"
<< endl
<< "% it_ scalar double time period for exogenous "
"variables for which"
<< endl
<< "% to evaluate the model" << endl
<< "%" << endl
<< "% Output:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms" << endl
<< "%" << endl
<< endl
<< "assert(length(T) >= " << ttlen << ");" << endl
<< endl;
if (!previous_tt_name.empty())
output << "T = " << basename << "." << previous_tt_name
<< "(T, y, x, params, steady_state, it_);" << endl
<< endl;
output << s_tt.str() << endl << "end" << endl;
output.close();
filename = packageDir(basename) / (name + ".m");
output.open(filename, ios::out | ios::binary);
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function " << retvalname << " = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "% function " << retvalname << " = " << name
<< "(T, y, x, params, steady_state, it_, T_flag)" << endl
<< "%" << endl
<< "% File created by Dynare Preprocessor from .mod file" << endl
<< "%" << endl
<< "% Inputs:" << endl
<< "% T [#temp variables by 1] double vector of temporary terms to be "
"filled by function"
<< endl
<< "% y [#dynamic variables by 1] double vector of endogenous variables "
"in the order stored"
<< endl
<< "% in M_.lead_lag_incidence; see "
"the Manual"
<< endl
<< "% x [nperiods by M_.exo_nbr] double matrix of exogenous variables "
"(in declaration order)"
<< endl
<< "% for all simulation periods"
<< endl
<< "% steady_state [M_.endo_nbr by 1] double vector of steady state values"
<< endl
<< "% params [M_.param_nbr by 1] double vector of parameter values in "
"declaration order"
<< endl
<< "% it_ scalar double time period for exogenous "
"variables for which"
<< endl
<< "% to evaluate the model" << endl
<< "% T_flag boolean boolean flag saying whether or not to "
"calculate temporary terms"
<< endl
<< "%" << endl
<< "% Output:" << endl
<< "% " << retvalname << endl
<< "%" << endl
<< endl;
if (!name_tt.empty())
output << "if T_flag" << endl
<< " T = " << basename << "." << name_tt << "(T, y, x, params, steady_state, it_);"
<< endl
<< "end" << endl;
output << init_s.str() << endl << s.str() << end_s.str() << endl << "end" << endl;
output.close();
}
void
DynamicModel::writeDynamicMCompatFile(const string& basename) const
{
filesystem::path filename {packageDir(basename) / "dynamic.m"};
ofstream output {filename, ios::out | ios::binary};
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
exit(EXIT_FAILURE);
}
int ntt {static_cast(
temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size()
+ temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size())};
output << "function [residual, g1, g2, g3] = dynamic(y, x, params, steady_state, it_)" << endl
<< " T = NaN(" << ntt << ", 1);" << endl
<< " if nargout <= 1" << endl
<< " residual = " << basename
<< ".dynamic_resid(T, y, x, params, steady_state, it_, true);" << endl
<< " elseif nargout == 2" << endl
<< " [residual, g1] = " << basename
<< ".dynamic_resid_g1(T, y, x, params, steady_state, it_, true);" << endl
<< " elseif nargout == 3" << endl
<< " [residual, g1, g2] = " << basename
<< ".dynamic_resid_g1_g2(T, y, x, params, steady_state, it_, true);" << endl
<< " else" << endl
<< " [residual, g1, g2, g3] = " << basename
<< ".dynamic_resid_g1_g2_g3(T, y, x, params, steady_state, it_, true);" << endl
<< " end" << endl
<< "end" << endl;
output.close();
}
vector