Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 4.6
  • 5.x
  • 6.x
  • aux_vars_fix
  • julia
  • llvm-15
  • master
  • python-codegen
  • rework_pac
  • uop
  • julia-6.2.0
  • julia-6.3.0
  • julia-6.4.0
  • julia-7.0.0
14 results

Target

Select target project
  • normann/preprocessor
  • Dynare/preprocessor
  • FerhatMihoubi/preprocessor
  • MichelJuillard/preprocessor
  • sebastien/preprocessor
  • lnsongxf/preprocessor
  • albop/preprocessor
  • DoraK/preprocessor
  • amg/preprocessor
  • wmutschl/preprocessor
  • JohannesPfeifer/preprocessor
11 results
Select Git revision
  • julia-6.3.0
  • julia-6.4.0
  • master
  • created_preprocessor_repo
4 results
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 2726 additions and 795 deletions
src/ExprNode.cc
View file @ 6618d7d1
/*
* Copyright © 2007-2023 Dynare Team
* Copyright © 2007-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -24,6 +24,7 @@
#include <limits>
#include <numbers>
#include <numeric>
#include <ranges>
#include <utility>
#include "DataTree.hh"
......@@ -194,8 +195,7 @@ ExprNode::collectVariables(SymbolType type, set<int>& result) const
{
set<pair<int, int>> symbs_lags;
collectDynamicVariables(type, symbs_lags);
transform(symbs_lags.begin(), symbs_lags.end(), inserter(result, result.begin()),
[](const auto& x) { return x.first; });
ranges::copy(views::keys(symbs_lags), inserter(result, result.begin()));
}
void
......@@ -395,9 +395,8 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
{
auto [orig_var_id, orig_lag] = datatree.symbol_table.unrollDiffLeadLagChain(var_id, lag);
not_ec = not_ec
|| (find(target_lhs.begin(), target_lhs.end(), orig_var_id) == target_lhs.end()
&& find(nontarget_lhs.begin(), nontarget_lhs.end(), orig_var_id)
== nontarget_lhs.end());
|| (ranges::find(target_lhs, orig_var_id) == target_lhs.end()
&& ranges::find(nontarget_lhs, orig_var_id) == nontarget_lhs.end());
}
if (not_ec)
continue;
......@@ -405,7 +404,7 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
// Now fill the matrices
for (const auto& [var_id, lag, param_id, constant] : error_linear_combination)
if (auto [orig_vid, orig_lag] = datatree.symbol_table.unrollDiffLeadLagChain(var_id, lag);
find(target_lhs.begin(), target_lhs.end(), orig_vid) == target_lhs.end())
ranges::find(target_lhs, orig_vid) == target_lhs.end())
{
if (orig_lag != -1)
{
......@@ -422,15 +421,15 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
<< endl;
exit(EXIT_FAILURE);
}
if (*param_id)
if (param_id)
{
cerr
<< "ERROR in trend component model: spurious parameter in error correction term"
<< endl;
exit(EXIT_FAILURE);
}
int colidx = static_cast<int>(distance(
nontarget_lhs.begin(), find(nontarget_lhs.begin(), nontarget_lhs.end(), orig_vid)));
int colidx = static_cast<int>(
ranges::distance(nontarget_lhs.begin(), ranges::find(nontarget_lhs, orig_vid)));
if (A0.contains({eqn, colidx}))
{
cerr << "ExprNode::fillErrorCorrection: Error filling A0 matrix: "
......@@ -443,7 +442,7 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
{
// This is a target, so fill A0star
int colidx = static_cast<int>(
distance(target_lhs.begin(), find(target_lhs.begin(), target_lhs.end(), orig_vid)));
ranges::distance(target_lhs.begin(), ranges::find(target_lhs, orig_vid)));
expr_t e = datatree.AddTimes(datatree.AddVariable(speed_of_adjustment_param),
datatree.AddPossiblyNegativeConstant(-constant));
if (param_id)
......@@ -804,10 +803,7 @@ NumConstNode::differentiateForwardVars([[maybe_unused]] const vector<string>& su
bool
NumConstNode::isNumConstNodeEqualTo(double value) const
{
if (datatree.num_constants.getDouble(id) == value)
return true;
else
return false;
return datatree.num_constants.getDouble(id) == value;
}
bool
......@@ -874,15 +870,20 @@ NumConstNode::substituteLogTransform([[maybe_unused]] int orig_symb_id,
return const_cast<NumConstNode*>(this);
}
expr_t
NumConstNode::substituteAggregationOperators([[maybe_unused]] subst_table_t& subst_table,
[[maybe_unused]] vector<BinaryOpNode*>& neweqs) const
{
return const_cast<NumConstNode*>(this);
}
VariableNode::VariableNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg, int lag_arg) :
ExprNode {datatree_arg, idx_arg}, symb_id {symb_id_arg}, lag {lag_arg}
{
// It makes sense to allow a lead/lag on parameters: during steady state calibration, endogenous
// and parameters can be swapped
assert(get_type() != SymbolType::externalFunction
&& (lag == 0
|| (get_type() != SymbolType::modelLocalVariable
&& get_type() != SymbolType::modFileLocalVariable)));
&& (lag == 0 || get_type() != SymbolType::modFileLocalVariable));
}
void
......@@ -901,18 +902,21 @@ VariableNode::prepareForDerivation()
case SymbolType::trend:
case SymbolType::logTrend:
// In static models, exogenous and trends do not have deriv IDs
if (!datatree.isDynamic())
if (!datatree.is_dynamic)
break;
[[fallthrough]];
case SymbolType::endogenous:
case SymbolType::parameter:
non_null_derivatives.insert(datatree.getDerivID(symb_id, lag));
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
non_null_derivatives.insert(getDerivID());
break;
case SymbolType::modelLocalVariable:
datatree.getLocalVariable(symb_id)->prepareForDerivation();
datatree.getLocalVariable(symb_id, lag)->prepareForDerivation();
// Non null derivatives are those of the value of the local parameter
non_null_derivatives = datatree.getLocalVariable(symb_id)->non_null_derivatives;
non_null_derivatives = datatree.getLocalVariable(symb_id, lag)->non_null_derivatives;
break;
case SymbolType::heterogeneousParameter:
case SymbolType::modFileLocalVariable:
case SymbolType::statementDeclaredVariable:
case SymbolType::unusedEndogenous:
......@@ -926,7 +930,7 @@ VariableNode::prepareForDerivation()
cerr << "VariableNode::prepareForDerivation: impossible case: "
<< "You are trying to derive a variable that has been excluded via "
"model_remove/var_remove/include_eqs/exclude_eqs: "
<< datatree.symbol_table.getName(symb_id) << endl;
<< getName() << endl;
exit(EXIT_FAILURE);
}
}
......@@ -944,14 +948,13 @@ VariableNode::prepareForChainRuleDerivation(
case SymbolType::endogenous:
{
set<int>& nnd {non_null_chain_rule_derivatives[const_cast<VariableNode*>(this)]};
int my_deriv_id {datatree.getDerivID(symb_id, lag)};
if (auto it = recursive_variables.find(my_deriv_id); it != recursive_variables.end())
if (auto it = recursive_variables.find(getDerivID()); it != recursive_variables.end())
{
it->second->arg2->prepareForChainRuleDerivation(recursive_variables,
non_null_chain_rule_derivatives);
nnd = non_null_chain_rule_derivatives.at(it->second->arg2);
}
nnd.insert(my_deriv_id);
nnd.insert(getDerivID());
}
break;
case SymbolType::exogenous:
......@@ -962,12 +965,15 @@ VariableNode::prepareForChainRuleDerivation(
case SymbolType::modFileLocalVariable:
case SymbolType::statementDeclaredVariable:
case SymbolType::unusedEndogenous:
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
case SymbolType::heterogeneousParameter:
// Those variables are never derived using chain rule
non_null_chain_rule_derivatives.try_emplace(const_cast<VariableNode*>(this));
break;
case SymbolType::modelLocalVariable:
{
expr_t def {datatree.getLocalVariable(symb_id)};
expr_t def {datatree.getLocalVariable(symb_id, lag)};
// Non null derivatives are those of the value of the model local variable
def->prepareForChainRuleDerivation(recursive_variables, non_null_chain_rule_derivatives);
non_null_chain_rule_derivatives.emplace(const_cast<VariableNode*>(this),
......@@ -992,17 +998,21 @@ VariableNode::computeDerivative(int deriv_id)
case SymbolType::trend:
case SymbolType::logTrend:
// In static models, exogenous and trends do not have deriv IDs
if (!datatree.isDynamic())
if (!datatree.is_dynamic)
return datatree.Zero;
[[fallthrough]];
case SymbolType::endogenous:
case SymbolType::parameter:
if (deriv_id == datatree.getDerivID(symb_id, lag))
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
if (deriv_id == getDerivID())
return datatree.One;
else
return datatree.Zero;
case SymbolType::heterogeneousParameter:
return datatree.Zero;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->getDerivative(deriv_id);
return datatree.getLocalVariable(symb_id, lag)->getDerivative(deriv_id);
case SymbolType::modFileLocalVariable:
cerr << "modFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
......@@ -1025,7 +1035,7 @@ bool
VariableNode::containsExternalFunction() const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->containsExternalFunction();
return datatree.getLocalVariable(symb_id, lag)->containsExternalFunction();
return false;
}
......@@ -1034,7 +1044,7 @@ void
VariableNode::writeJsonAST(ostream& output) const
{
output << R"({"node_type" : "VariableNode", )"
<< R"("name" : ")" << datatree.symbol_table.getName(symb_id) << R"(", "type" : ")";
<< R"("name" : ")" << getName() << R"(", "type" : ")";
switch (get_type())
{
case SymbolType::endogenous:
......@@ -1076,8 +1086,23 @@ VariableNode::writeJsonAST(ostream& output) const
case SymbolType::excludedVariable:
cerr << "VariableNode::computeDerivative: Impossible case!" << endl;
exit(EXIT_FAILURE);
case SymbolType::heterogeneousEndogenous:
output << "heterogeneousEndogenous";
break;
case SymbolType::heterogeneousExogenous:
output << "heterogeneousExogenous";
break;
case SymbolType::heterogeneousParameter:
output << "heterogeneousParameter";
break;
}
output << R"(", "lag" : )" << lag << "}";
output << '"';
if (isHeterogeneous(get_type()))
output << R"(, "heterogeneity_dimension" : ")"
<< datatree.heterogeneity_table.getName(
datatree.symbol_table.getHeterogeneityDimension(symb_id))
<< '"';
output << R"(, "lag" : )" << lag << "}";
}
void
......@@ -1091,7 +1116,7 @@ VariableNode::writeJsonOutput(ostream& output, const temporary_terms_t& temporar
return;
}
output << datatree.symbol_table.getName(symb_id);
output << getName();
if (isdynamic && lag != 0)
output << "(" << lag << ")";
}
......@@ -1131,7 +1156,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
auto juliaTimeDataFrameHelper = [&] {
if (lag != 0)
output << "lag(";
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
if (lag != 0)
{
if (lag != -1)
......@@ -1144,13 +1169,13 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
switch (type)
{
case SymbolType::parameter:
if (int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
output_type == ExprNodeOutputType::matlabOutsideModel)
if (output_type == ExprNodeOutputType::matlabOutsideModel)
output << "M_.params"
<< "(" << tsid + 1 << ")";
<< "(" << getTypeSpecificID() + 1 << ")";
else
output << "params" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< tsid + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type);
<< getTypeSpecificID() + ARRAY_SUBSCRIPT_OFFSET(output_type)
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case SymbolType::modelLocalVariable:
......@@ -1158,31 +1183,36 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
|| output_type == ExprNodeOutputType::CDynamicSteadyStateOperator)
{
output << "(";
datatree.getLocalVariable(symb_id)->writeOutput(output, output_type, temporary_terms,
temporary_terms_idxs, tef_terms);
datatree.getLocalVariable(symb_id, lag)
->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
output << ")";
}
else
/* We append underscores to avoid name clashes with "g1" or "oo_".
But we probably never arrive here because MLV are temporary terms… */
output << datatree.symbol_table.getName(symb_id) << "__";
output << getName() << "__";
break;
case SymbolType::modFileLocalVariable:
output << datatree.symbol_table.getName(symb_id);
output << getName();
break;
case SymbolType::endogenous:
switch (int tsid = datatree.symbol_table.getTypeSpecificID(symb_id); output_type)
switch (int tsid {getTypeSpecificID()}; output_type)
{
case ExprNodeOutputType::juliaDynamicModel:
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
i = datatree.getJacobianCol(datatree.getDerivID(symb_id, lag),
isSparseModelOutput(output_type))
assert(lag >= -1 && lag <= 1);
i = tsid + (lag + 1) * datatree.symbol_table.endo_nbr()
+ ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case ExprNodeOutputType::juliaDynamicModel:
case ExprNodeOutputType::matlabDynamicModel:
case ExprNodeOutputType::CDynamicModel:
i = datatree.getJacobianCol(getDerivID(), isSparseModelOutput(output_type))
+ ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
......@@ -1214,7 +1244,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case ExprNodeOutputType::matlabDseries:
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
if (lag != 0)
output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
......@@ -1223,7 +1253,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
juliaTimeDataFrameHelper();
break;
case ExprNodeOutputType::epilogueFile:
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
if (lag != 0)
output << lag;
......@@ -1239,7 +1269,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
break;
case SymbolType::exogenous:
i = datatree.symbol_table.getTypeSpecificID(symb_id) + ARRAY_SUBSCRIPT_OFFSET(output_type);
i = getTypeSpecificID() + ARRAY_SUBSCRIPT_OFFSET(output_type);
switch (output_type)
{
case ExprNodeOutputType::juliaDynamicModel:
......@@ -1289,7 +1319,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case ExprNodeOutputType::matlabDseries:
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
if (lag != 0)
output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
......@@ -1298,7 +1328,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
juliaTimeDataFrameHelper();
break;
case ExprNodeOutputType::epilogueFile:
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
if (lag != 0)
output << lag;
......@@ -1311,7 +1341,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
break;
case SymbolType::exogenousDet:
i = datatree.symbol_table.getTypeSpecificID(symb_id) + datatree.symbol_table.exo_nbr()
i = getTypeSpecificID() + datatree.symbol_table.exo_nbr()
+ ARRAY_SUBSCRIPT_OFFSET(output_type);
switch (output_type)
{
......@@ -1351,12 +1381,10 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
break;
case ExprNodeOutputType::matlabOutsideModel:
assert(lag == 0);
output << "oo_.exo_det_steady_state("
<< datatree.symbol_table.getTypeSpecificID(symb_id) + 1 << ")";
output << "oo_.exo_det_steady_state(" << getTypeSpecificID() + 1 << ")";
break;
case ExprNodeOutputType::matlabDynamicSteadyStateOperator:
output << "oo_.exo_det_steady_state("
<< datatree.symbol_table.getTypeSpecificID(symb_id) + 1 << ")";
output << "oo_.exo_det_steady_state(" << getTypeSpecificID() + 1 << ")";
break;
case ExprNodeOutputType::juliaSteadyStateFile:
case ExprNodeOutputType::steadyStateFile:
......@@ -1364,7 +1392,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case ExprNodeOutputType::matlabDseries:
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
if (lag != 0)
output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
......@@ -1373,7 +1401,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
juliaTimeDataFrameHelper();
break;
case ExprNodeOutputType::epilogueFile:
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
if (lag != 0)
output << lag;
......@@ -1387,7 +1415,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
case SymbolType::epilogue:
if (output_type == ExprNodeOutputType::epilogueFile)
{
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
if (lag != 0)
output << lag;
......@@ -1396,7 +1424,7 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
else if (output_type == ExprNodeOutputType::matlabDseries
|| output_type == ExprNodeOutputType::juliaTimeDataFrame)
// Only writing dseries for epilogue_static, hence no need to check lag
output << "ds." << datatree.symbol_table.getName(symb_id);
output << "ds." << getName();
else
{
cerr << "VariableNode::writeOutput: Impossible case" << endl;
......@@ -1415,6 +1443,64 @@ VariableNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
case SymbolType::excludedVariable:
cerr << "VariableNode::writeOutput: Impossible case" << endl;
exit(EXIT_FAILURE);
case SymbolType::heterogeneousEndogenous:
switch (int tsid = datatree.symbol_table.getTypeSpecificID(symb_id); output_type)
{
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
assert(lag >= -1 && lag <= 1);
i = tsid
+ (lag + 1)
* datatree.symbol_table.het_endo_nbr(
datatree.symbol_table.getHeterogeneityDimension(symb_id))
+ ARRAY_SUBSCRIPT_OFFSET(output_type);
output << "yh" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
default:
cerr << "VariableNode::writeOutput: unsupported output type for heterogeneousEndogenous "
"symbol"
<< endl;
exit(EXIT_FAILURE);
}
break;
case SymbolType::heterogeneousExogenous:
i = datatree.symbol_table.getTypeSpecificID(symb_id) + ARRAY_SUBSCRIPT_OFFSET(output_type);
switch (output_type)
{
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
assert(lag == 0);
output << "xh" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
default:
cerr << "VariableNode::writeOutput: unsupported output type for heterogeneousExogenous "
"symbol"
<< endl;
exit(EXIT_FAILURE);
}
break;
case SymbolType::heterogeneousParameter:
i = datatree.symbol_table.getTypeSpecificID(symb_id) + ARRAY_SUBSCRIPT_OFFSET(output_type);
switch (output_type)
{
case ExprNodeOutputType::juliaSparseDynamicModel:
case ExprNodeOutputType::matlabSparseDynamicModel:
case ExprNodeOutputType::CSparseDynamicModel:
output << "paramsh" << LEFT_ARRAY_SUBSCRIPT(output_type) << i
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
default:
cerr << "VariableNode::writeOutput: unsupported output type for heterogeneousParameter "
"symbol"
<< endl;
exit(EXIT_FAILURE);
}
break;
}
}
......@@ -1422,7 +1508,7 @@ double
VariableNode::eval(const eval_context_t& eval_context) const noexcept(false)
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->eval(eval_context);
return datatree.getLocalVariable(symb_id, lag)->eval(eval_context);
auto it = eval_context.find(symb_id);
if (it == eval_context.end())
......@@ -1442,33 +1528,30 @@ VariableNode::writeBytecodeOutput(Bytecode::Writer& code_file,
temporary_terms_idxs))
return;
auto type = get_type();
if (type == SymbolType::modelLocalVariable || type == SymbolType::modFileLocalVariable)
datatree.getLocalVariable(symb_id)->writeBytecodeOutput(code_file, output_type, temporary_terms,
temporary_terms_idxs, tef_terms);
if (auto type = get_type(); type == SymbolType::modelLocalVariable)
datatree.getLocalVariable(symb_id, lag)
->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs,
tef_terms);
else
{
int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
switch (output_type)
{
case ExprNodeBytecodeOutputType::dynamicModel:
code_file << Bytecode::FLDV {type, tsid, lag};
code_file << Bytecode::FLDV {type, getTypeSpecificID(), lag};
break;
case ExprNodeBytecodeOutputType::staticModel:
code_file << Bytecode::FLDSV {type, tsid};
code_file << Bytecode::FLDSV {type, getTypeSpecificID()};
break;
case ExprNodeBytecodeOutputType::dynamicSteadyStateOperator:
code_file << Bytecode::FLDVS {type, tsid};
code_file << Bytecode::FLDVS {type, getTypeSpecificID()};
break;
case ExprNodeBytecodeOutputType::dynamicAssignmentLHS:
code_file << Bytecode::FSTPV {type, tsid, lag};
code_file << Bytecode::FSTPV {type, getTypeSpecificID(), lag};
break;
case ExprNodeBytecodeOutputType::staticAssignmentLHS:
code_file << Bytecode::FSTPSV {type, tsid};
code_file << Bytecode::FSTPSV {type, getTypeSpecificID()};
break;
}
}
}
void
VariableNode::collectVARLHSVariable(set<expr_t>& result) const
......@@ -1488,7 +1571,7 @@ VariableNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int>>&
if (get_type() == type_arg)
result.emplace(symb_id, lag);
if (get_type() == SymbolType::modelLocalVariable)
datatree.getLocalVariable(symb_id)->collectDynamicVariables(type_arg, result);
datatree.getLocalVariable(symb_id, lag)->collectDynamicVariables(type_arg, result);
}
void
......@@ -1498,8 +1581,8 @@ VariableNode::computeSubExprContainingVariable(int symb_id_arg, int lag_arg,
if (symb_id == symb_id_arg && lag == lag_arg)
contain_var.insert(const_cast<VariableNode*>(this));
if (get_type() == SymbolType::modelLocalVariable)
datatree.getLocalVariable(symb_id)->computeSubExprContainingVariable(symb_id_arg, lag_arg,
contain_var);
datatree.getLocalVariable(symb_id, lag)
->computeSubExprContainingVariable(symb_id_arg, lag_arg, contain_var);
}
BinaryOpNode*
......@@ -1508,7 +1591,7 @@ VariableNode::normalizeEquationHelper(const set<expr_t>& contain_var, expr_t rhs
assert(contain_var.contains(const_cast<VariableNode*>(this)));
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->normalizeEquationHelper(contain_var, rhs);
return datatree.getLocalVariable(symb_id, lag)->normalizeEquationHelper(contain_var, rhs);
// This the LHS variable: we have finished the normalization
return datatree.AddEqual(const_cast<VariableNode*>(this), rhs);
......@@ -1527,23 +1610,28 @@ VariableNode::computeChainRuleDerivative(
case SymbolType::trend:
case SymbolType::logTrend:
// In static models, exogenous and trends do not have deriv IDs
if (!datatree.isDynamic())
if (!datatree.is_dynamic)
return datatree.Zero;
[[fallthrough]];
case SymbolType::endogenous:
case SymbolType::parameter:
if (int my_deriv_id {datatree.getDerivID(symb_id, lag)}; deriv_id == my_deriv_id)
case SymbolType::heterogeneousEndogenous:
if (deriv_id == getDerivID())
return datatree.One;
// If there is in the equation a recursive variable we could use a chaine rule derivation
else if (auto it = recursive_variables.find(my_deriv_id); it != recursive_variables.end())
else if (auto it = recursive_variables.find(getDerivID()); it != recursive_variables.end())
return it->second->arg2->getChainRuleDerivative(deriv_id, recursive_variables,
non_null_chain_rule_derivatives, cache);
else
return datatree.Zero;
case SymbolType::heterogeneousExogenous:
case SymbolType::heterogeneousParameter:
return datatree.Zero;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->getChainRuleDerivative(
deriv_id, recursive_variables, non_null_chain_rule_derivatives, cache);
return datatree.getLocalVariable(symb_id, lag)
->getChainRuleDerivative(deriv_id, recursive_variables, non_null_chain_rule_derivatives,
cache);
case SymbolType::modFileLocalVariable:
cerr << "modFileLocalVariable is not derivable" << endl;
exit(EXIT_FAILURE);
......@@ -1585,17 +1673,20 @@ VariableNode::computeXrefs(EquationInfo& ei) const
case SymbolType::parameter:
ei.param.emplace(symb_id, 0);
break;
case SymbolType::modFileLocalVariable:
datatree.getLocalVariable(symb_id)->computeXrefs(ei);
case SymbolType::modelLocalVariable:
datatree.getLocalVariable(symb_id, lag)->computeXrefs(ei);
break;
case SymbolType::trend:
case SymbolType::logTrend:
case SymbolType::modelLocalVariable:
case SymbolType::modFileLocalVariable:
case SymbolType::statementDeclaredVariable:
case SymbolType::unusedEndogenous:
case SymbolType::externalFunction:
case SymbolType::epilogue:
case SymbolType::excludedVariable:
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
case SymbolType::heterogeneousParameter:
break;
}
}
......@@ -1606,6 +1697,24 @@ VariableNode::get_type() const
return datatree.symbol_table.getType(symb_id);
}
string
VariableNode::getName() const
{
return datatree.symbol_table.getName(symb_id);
}
int
VariableNode::getDerivID() const
{
return datatree.getDerivID(symb_id, lag);
}
int
VariableNode::getTypeSpecificID() const
{
return datatree.symbol_table.getTypeSpecificID(symb_id);
}
expr_t
VariableNode::clone(DataTree& alt_datatree) const
{
......@@ -1620,7 +1729,7 @@ VariableNode::maxEndoLead() const
case SymbolType::endogenous:
return max(lag, 0);
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxEndoLead();
return datatree.getLocalVariable(symb_id, lag)->maxEndoLead();
default:
return 0;
}
......@@ -1634,7 +1743,7 @@ VariableNode::maxExoLead() const
case SymbolType::exogenous:
return max(lag, 0);
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxExoLead();
return datatree.getLocalVariable(symb_id, lag)->maxExoLead();
default:
return 0;
}
......@@ -1648,7 +1757,7 @@ VariableNode::maxEndoLag() const
case SymbolType::endogenous:
return max(-lag, 0);
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxEndoLag();
return datatree.getLocalVariable(symb_id, lag)->maxEndoLag();
default:
return 0;
}
......@@ -1662,7 +1771,7 @@ VariableNode::maxExoLag() const
case SymbolType::exogenous:
return max(-lag, 0);
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxExoLag();
return datatree.getLocalVariable(symb_id, lag)->maxExoLag();
default:
return 0;
}
......@@ -1676,9 +1785,12 @@ VariableNode::maxLead() const
case SymbolType::endogenous:
case SymbolType::exogenous:
case SymbolType::exogenousDet:
case SymbolType::epilogue:
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
return lag;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxLead();
return datatree.getLocalVariable(symb_id, lag)->maxLead();
default:
return 0;
}
......@@ -1692,9 +1804,12 @@ VariableNode::maxLag() const
case SymbolType::endogenous:
case SymbolType::exogenous:
case SymbolType::exogenousDet:
case SymbolType::epilogue:
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
return -lag;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxLag();
return datatree.getLocalVariable(symb_id, lag)->maxLag();
default:
return 0;
}
......@@ -1709,9 +1824,11 @@ VariableNode::maxLagWithDiffsExpanded() const
case SymbolType::exogenous:
case SymbolType::exogenousDet:
case SymbolType::epilogue:
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
return -lag;
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->maxLagWithDiffsExpanded();
return datatree.getLocalVariable(symb_id, lag)->maxLagWithDiffsExpanded();
default:
return 0;
}
......@@ -1745,7 +1862,7 @@ expr_t
VariableNode::substituteModelLocalVariables() const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id);
return datatree.getLocalVariable(symb_id, lag);
return const_cast<VariableNode*>(this);
}
......@@ -1754,7 +1871,7 @@ expr_t
VariableNode::substituteVarExpectation(const map<string, expr_t>& subst_table) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substituteVarExpectation(subst_table);
return datatree.getLocalVariable(symb_id, lag)->substituteVarExpectation(subst_table);
return const_cast<VariableNode*>(this);
}
......@@ -1763,21 +1880,21 @@ void
VariableNode::findDiffNodes(lag_equivalence_table_t& nodes) const
{
if (get_type() == SymbolType::modelLocalVariable)
datatree.getLocalVariable(symb_id)->findDiffNodes(nodes);
datatree.getLocalVariable(symb_id, lag)->findDiffNodes(nodes);
}
void
VariableNode::findUnaryOpNodesForAuxVarCreation(lag_equivalence_table_t& nodes) const
{
if (get_type() == SymbolType::modelLocalVariable)
datatree.getLocalVariable(symb_id)->findUnaryOpNodesForAuxVarCreation(nodes);
datatree.getLocalVariable(symb_id, lag)->findUnaryOpNodesForAuxVarCreation(nodes);
}
optional<int>
VariableNode::findTargetVariable(int lhs_symb_id) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->findTargetVariable(lhs_symb_id);
return datatree.getLocalVariable(symb_id, lag)->findTargetVariable(lhs_symb_id);
return nullopt;
}
......@@ -1787,7 +1904,7 @@ VariableNode::substituteDiff(const lag_equivalence_table_t& nodes, subst_table_t
vector<BinaryOpNode*>& neweqs) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substituteDiff(nodes, subst_table, neweqs);
return datatree.getLocalVariable(symb_id, lag)->substituteDiff(nodes, subst_table, neweqs);
return const_cast<VariableNode*>(this);
}
......@@ -1798,7 +1915,8 @@ VariableNode::substituteUnaryOpNodes(const lag_equivalence_table_t& nodes,
vector<BinaryOpNode*>& neweqs) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substituteUnaryOpNodes(nodes, subst_table, neweqs);
return datatree.getLocalVariable(symb_id, lag)
->substituteUnaryOpNodes(nodes, subst_table, neweqs);
return const_cast<VariableNode*>(this);
}
......@@ -1807,7 +1925,7 @@ expr_t
VariableNode::substitutePacExpectation(const string& name, expr_t subexpr)
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substitutePacExpectation(name, subexpr);
return datatree.getLocalVariable(symb_id, lag)->substitutePacExpectation(name, subexpr);
return const_cast<VariableNode*>(this);
}
......@@ -1816,7 +1934,7 @@ expr_t
VariableNode::substitutePacTargetNonstationary(const string& name, expr_t subexpr)
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substitutePacTargetNonstationary(name, subexpr);
return datatree.getLocalVariable(symb_id, lag)->substitutePacTargetNonstationary(name, subexpr);
return const_cast<VariableNode*>(this);
}
......@@ -1829,11 +1947,14 @@ VariableNode::decreaseLeadsLags(int n) const
case SymbolType::endogenous:
case SymbolType::exogenous:
case SymbolType::exogenousDet:
case SymbolType::epilogue:
case SymbolType::trend:
case SymbolType::logTrend:
case SymbolType::heterogeneousEndogenous:
case SymbolType::heterogeneousExogenous:
return datatree.AddVariable(symb_id, lag - n);
case SymbolType::modelLocalVariable:
return datatree.getLocalVariable(symb_id)->decreaseLeadsLags(n);
return datatree.getLocalVariable(symb_id, lag)->decreaseLeadsLags(n);
default:
return const_cast<VariableNode*>(this);
}
......@@ -1864,7 +1985,7 @@ VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t& subst_table,
else
return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
case SymbolType::modelLocalVariable:
if (expr_t value = datatree.getLocalVariable(symb_id); value->maxEndoLead() <= 1)
if (expr_t value = datatree.getLocalVariable(symb_id, lag); value->maxEndoLead() <= 1)
return const_cast<VariableNode*>(this);
else
return value->substituteEndoLeadGreaterThanTwo(subst_table, neweqs, deterministic_model);
......@@ -1913,7 +2034,7 @@ VariableNode::substituteEndoLagGreaterThanTwo(subst_table_t& subst_table,
return substexpr;
case SymbolType::modelLocalVariable:
if (expr_t value = datatree.getLocalVariable(symb_id); value->maxEndoLag() <= 1)
if (expr_t value = datatree.getLocalVariable(symb_id, lag); value->maxEndoLag() <= 1)
return const_cast<VariableNode*>(this);
else
return value->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
......@@ -1934,7 +2055,7 @@ VariableNode::substituteExoLead(subst_table_t& subst_table, vector<BinaryOpNode*
else
return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
case SymbolType::modelLocalVariable:
if (expr_t value = datatree.getLocalVariable(symb_id); value->maxExoLead() == 0)
if (expr_t value = datatree.getLocalVariable(symb_id, lag); value->maxExoLead() == 0)
return const_cast<VariableNode*>(this);
else
return value->substituteExoLead(subst_table, neweqs, deterministic_model);
......@@ -1982,7 +2103,7 @@ VariableNode::substituteExoLag(subst_table_t& subst_table, vector<BinaryOpNode*>
return substexpr;
case SymbolType::modelLocalVariable:
if (expr_t value = datatree.getLocalVariable(symb_id); value->maxExoLag() == 0)
if (expr_t value = datatree.getLocalVariable(symb_id, lag); value->maxExoLag() == 0)
return const_cast<VariableNode*>(this);
else
return value->substituteExoLag(subst_table, neweqs);
......@@ -1996,8 +2117,8 @@ VariableNode::substituteExpectation(subst_table_t& subst_table, vector<BinaryOpN
bool partial_information_model) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substituteExpectation(subst_table, neweqs,
partial_information_model);
return datatree.getLocalVariable(symb_id, lag)
->substituteExpectation(subst_table, neweqs, partial_information_model);
return const_cast<VariableNode*>(this);
}
......@@ -2010,10 +2131,7 @@ VariableNode::differentiateForwardVars(const vector<string>& subset, subst_table
{
case SymbolType::endogenous:
assert(lag <= 1);
if (lag <= 0
|| (subset.size() > 0
&& find(subset.begin(), subset.end(), datatree.symbol_table.getName(symb_id))
== subset.end()))
if (lag <= 0 || (subset.size() > 0 && ranges::find(subset, getName()) == subset.end()))
return const_cast<VariableNode*>(this);
else
{
......@@ -2033,7 +2151,7 @@ VariableNode::differentiateForwardVars(const vector<string>& subset, subst_table
return datatree.AddPlus(datatree.AddVariable(symb_id, 0), diffvar);
}
case SymbolType::modelLocalVariable:
if (expr_t value = datatree.getLocalVariable(symb_id); value->maxEndoLead() <= 0)
if (expr_t value = datatree.getLocalVariable(symb_id, lag); value->maxEndoLead() <= 0)
return const_cast<VariableNode*>(this);
else
return value->differentiateForwardVars(subset, subst_table, neweqs);
......@@ -2051,18 +2169,14 @@ VariableNode::isNumConstNodeEqualTo([[maybe_unused]] double value) const
bool
VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
{
if (get_type() == type_arg && datatree.symbol_table.getTypeSpecificID(symb_id) == variable_id
&& lag == lag_arg)
return true;
else
return false;
return get_type() == type_arg && getTypeSpecificID() == variable_id && lag == lag_arg;
}
bool
VariableNode::containsPacExpectation(const string& pac_model_name) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->containsPacExpectation(pac_model_name);
return datatree.getLocalVariable(symb_id, lag)->containsPacExpectation(pac_model_name);
return false;
}
......@@ -2071,7 +2185,7 @@ bool
VariableNode::containsPacTargetNonstationary(const string& pac_model_name) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->containsPacTargetNonstationary(pac_model_name);
return datatree.getLocalVariable(symb_id, lag)->containsPacTargetNonstationary(pac_model_name);
return false;
}
......@@ -2080,7 +2194,7 @@ expr_t
VariableNode::replaceTrendVar() const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->replaceTrendVar();
return datatree.getLocalVariable(symb_id, lag)->replaceTrendVar();
if (get_type() == SymbolType::trend)
return datatree.One;
......@@ -2094,7 +2208,7 @@ expr_t
VariableNode::detrend(int symb_id, bool log_trend, expr_t trend) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->detrend(symb_id, log_trend, trend);
return datatree.getLocalVariable(symb_id, lag)->detrend(symb_id, log_trend, trend);
if (this->symb_id != symb_id)
return const_cast<VariableNode*>(this);
......@@ -2119,7 +2233,7 @@ int
VariableNode::countDiffs() const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->countDiffs();
return datatree.getLocalVariable(symb_id, lag)->countDiffs();
return 0;
}
......@@ -2128,7 +2242,7 @@ expr_t
VariableNode::removeTrendLeadLag(const map<int, expr_t>& trend_symbols_map) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->removeTrendLeadLag(trend_symbols_map);
return datatree.getLocalVariable(symb_id, lag)->removeTrendLeadLag(trend_symbols_map);
if ((get_type() != SymbolType::trend && get_type() != SymbolType::logTrend) || lag == 0)
return const_cast<VariableNode*>(this);
......@@ -2180,7 +2294,7 @@ bool
VariableNode::isInStaticForm() const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->isInStaticForm();
return datatree.getLocalVariable(symb_id, lag)->isInStaticForm();
return lag == 0;
}
......@@ -2189,7 +2303,7 @@ bool
VariableNode::isParamTimesEndogExpr() const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->isParamTimesEndogExpr();
return datatree.getLocalVariable(symb_id, lag)->isParamTimesEndogExpr();
return false;
}
......@@ -2214,8 +2328,7 @@ VariableNode::matchMatchedMoment(vector<int>& symb_ids, vector<int>& lags,
model local variables */
if (get_type() != SymbolType::endogenous)
throw MatchFailureException {"Variable " + datatree.symbol_table.getName(symb_id)
+ " is not an endogenous"};
throw MatchFailureException {"Variable " + getName() + " is not an endogenous"};
symb_ids.push_back(symb_id);
lags.push_back(lag);
......@@ -2226,7 +2339,8 @@ expr_t
VariableNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id)->substituteLogTransform(orig_symb_id, aux_symb_id);
return datatree.getLocalVariable(symb_id, lag)
->substituteLogTransform(orig_symb_id, aux_symb_id);
if (symb_id == orig_symb_id)
return datatree.AddExp(datatree.AddVariable(aux_symb_id, lag));
......@@ -2234,6 +2348,17 @@ VariableNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return const_cast<VariableNode*>(this);
}
expr_t
VariableNode::substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const
{
if (get_type() == SymbolType::modelLocalVariable)
return datatree.getLocalVariable(symb_id, lag)
->substituteAggregationOperators(subst_table, neweqs);
return const_cast<VariableNode*>(this);
}
UnaryOpNode::UnaryOpNode(DataTree& datatree_arg, int idx_arg, UnaryOpcode op_code_arg,
const expr_t arg_arg, int expectation_information_set_arg,
int param1_symb_id_arg, int param2_symb_id_arg, string adl_param_name_arg,
......@@ -2259,12 +2384,12 @@ UnaryOpNode::prepareForDerivation()
/* Non-null derivatives are those of the argument (except for STEADY_STATE in
a dynamic context, in which case the potentially non-null derivatives are
all the parameters) */
all the parameters; and for SUM, which will be substituted out before deriving) */
if ((op_code == UnaryOpcode::steadyState || op_code == UnaryOpcode::steadyStateParamDeriv
|| op_code == UnaryOpcode::steadyStateParam2ndDeriv)
&& datatree.isDynamic())
&& datatree.is_dynamic)
datatree.addAllParamDerivId(non_null_derivatives);
else
else if (op_code != UnaryOpcode::sum)
{
arg->prepareForDerivation();
non_null_derivatives = arg->non_null_derivatives;
......@@ -2285,7 +2410,7 @@ UnaryOpNode::prepareForChainRuleDerivation(
set<int>& nnd {non_null_chain_rule_derivatives[const_cast<UnaryOpNode*>(this)]};
if ((op_code == UnaryOpcode::steadyState || op_code == UnaryOpcode::steadyStateParamDeriv
|| op_code == UnaryOpcode::steadyStateParam2ndDeriv)
&& datatree.isDynamic())
&& datatree.is_dynamic)
datatree.addAllParamDerivId(nnd);
else
{
......@@ -2367,7 +2492,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
case UnaryOpcode::sign:
return datatree.Zero;
case UnaryOpcode::steadyState:
if (datatree.isDynamic())
if (datatree.is_dynamic)
{
if (datatree.getTypeByDerivID(deriv_id) == SymbolType::parameter)
{
......@@ -2392,7 +2517,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
else
return darg;
case UnaryOpcode::steadyStateParamDeriv:
assert(datatree.isDynamic());
assert(datatree.is_dynamic);
if (datatree.getTypeByDerivID(deriv_id) == SymbolType::parameter)
{
auto varg = dynamic_cast<VariableNode*>(arg);
......@@ -2404,7 +2529,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
else
return datatree.Zero;
case UnaryOpcode::steadyStateParam2ndDeriv:
assert(datatree.isDynamic());
assert(datatree.is_dynamic);
if (datatree.getTypeByDerivID(deriv_id) == SymbolType::parameter)
{
cerr << "3rd derivative of STEADY_STATE node w.r.t. three parameters not implemented"
......@@ -2441,6 +2566,9 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
case UnaryOpcode::adl:
cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
case UnaryOpcode::sum:
cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::sum" << endl;
exit(EXIT_FAILURE);
}
__builtin_unreachable(); // Silence GCC warning
}
......@@ -2469,9 +2597,8 @@ UnaryOpNode::cost(const vector<vector<unordered_set<expr_t>>>& blocks_temporary_
bool is_matlab) const
{
// For a temporary term, the cost is null
for (const auto& blk_tt : blocks_temporary_terms)
for (const auto& eq_tt : blk_tt)
if (eq_tt.contains(const_cast<UnaryOpNode*>(this)))
if (ranges::any_of(views::join(blocks_temporary_terms),
[this](auto& tt) { return tt.contains(const_cast<UnaryOpNode*>(this)); }))
return 0;
return cost(arg->cost(blocks_temporary_terms, is_matlab), is_matlab);
......@@ -2535,6 +2662,8 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
case UnaryOpcode::adl:
cerr << "UnaryOpNode::cost: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
case UnaryOpcode::sum:
return 0; // In the generated files, the SUM() operator behaves like a variable
}
else
// Cost for C files
......@@ -2585,6 +2714,8 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
case UnaryOpcode::adl:
cerr << "UnaryOpNode::cost: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
case UnaryOpcode::sum:
return 0; // In the generated files, the SUM() operator behaves like a variable
}
__builtin_unreachable(); // Silence GCC warning
}
......@@ -2729,6 +2860,9 @@ UnaryOpNode::writeJsonAST(ostream& output) const
case UnaryOpcode::erfc:
output << "erfc";
break;
case UnaryOpcode::sum:
output << "sum";
break;
}
output << R"(", "arg" : )";
arg->writeJsonAST(output);
......@@ -2843,10 +2977,8 @@ UnaryOpNode::writeJsonOutput(ostream& output, const temporary_terms_t& temporary
output << "])";
return;
case UnaryOpcode::steadyState:
output << "(";
arg->writeJsonOutput(output, temporary_terms, tef_terms, isdynamic);
output << ")";
return;
output << "STEADY_STATE";
break;
case UnaryOpcode::steadyStateParamDeriv:
{
auto varg = dynamic_cast<VariableNode*>(arg);
......@@ -2881,6 +3013,9 @@ UnaryOpNode::writeJsonOutput(ostream& output, const temporary_terms_t& temporary
case UnaryOpcode::erfc:
output << "erfc";
break;
case UnaryOpcode::sum:
output << "sum";
break;
}
bool close_parenthesis = false;
......@@ -3111,6 +3246,20 @@ UnaryOpNode::writeOutput(ostream& output, ExprNodeOutputType output_type,
case UnaryOpcode::adl:
output << "adl";
break;
case UnaryOpcode::sum:
if (isLatexOutput(output_type))
output << "sum";
else
{
auto varg {dynamic_cast<VariableNode*>(arg)};
assert(varg && varg->lag == 0);
int idx {
datatree.heterogeneity_table.getSummedHeterogenousEndogenousIndex(varg->symb_id)};
output << "yagg" << LEFT_ARRAY_SUBSCRIPT(output_type) << idx + 1
<< RIGHT_ARRAY_SUBSCRIPT(output_type);
return;
}
break;
}
if (output_type == ExprNodeOutputType::juliaTimeDataFrame && op_code != UnaryOpcode::uminus)
......@@ -3240,6 +3389,9 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) noexcept(false)
case UnaryOpcode::adl:
cerr << "UnaryOpNode::eval_opcode: not implemented on UnaryOpcode::adl" << endl;
exit(EXIT_FAILURE);
case UnaryOpcode::sum:
cerr << "UnaryOpNode::eval_opcode: not implemented on UnaryOpcode::sum" << endl;
exit(EXIT_FAILURE);
}
__builtin_unreachable(); // Silence GCC warning
}
......@@ -3497,6 +3649,8 @@ UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree& alt_datatree) con
return alt_datatree.AddDiff(alt_arg);
case UnaryOpcode::adl:
return alt_datatree.AddAdl(alt_arg, adl_param_name, adl_lags);
case UnaryOpcode::sum:
return alt_datatree.AddSum(alt_arg);
}
__builtin_unreachable(); // Silence GCC warning
}
......@@ -4089,6 +4243,34 @@ UnaryOpNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
expr_t
UnaryOpNode::substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const
{
if (op_code == UnaryOpcode::sum)
{
if (auto it = subst_table.find(const_cast<UnaryOpNode*>(this)); it != subst_table.end())
return const_cast<VariableNode*>(it->second);
// Arriving here, we need to create an auxiliary variable for this operator
const VariableNode* varg {dynamic_cast<const VariableNode*>(arg)};
assert(varg && varg->lag == 0);
const string varname {"SUM_" + varg->getName()};
int symb_id {datatree.symbol_table.addAggregationOpAuxiliaryVar(
varname, const_cast<UnaryOpNode*>(this))};
expr_t newAuxE = datatree.AddVariable(symb_id, 0);
subst_table[this] = dynamic_cast<VariableNode*>(newAuxE);
neweqs.push_back(datatree.AddEqual(newAuxE, const_cast<UnaryOpNode*>(this)));
datatree.heterogeneity_table.addSummedHeterogeneousEndogenous(varg->symb_id);
return newAuxE;
}
else
return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
}
BinaryOpNode::BinaryOpNode(DataTree& datatree_arg, int idx_arg, const expr_t arg1_arg,
BinaryOpcode op_code_arg, const expr_t arg2_arg,
int powerDerivOrder_arg) :
......@@ -4114,8 +4296,7 @@ BinaryOpNode::prepareForDerivation()
// Non-null derivatives are the union of those of the arguments
// Compute set union of arg1->non_null_derivatives and arg2->non_null_derivatives
set_union(arg1->non_null_derivatives.begin(), arg1->non_null_derivatives.end(),
arg2->non_null_derivatives.begin(), arg2->non_null_derivatives.end(),
ranges::set_union(arg1->non_null_derivatives, arg2->non_null_derivatives,
inserter(non_null_derivatives, non_null_derivatives.begin()));
}
......@@ -4131,10 +4312,8 @@ BinaryOpNode::prepareForChainRuleDerivation(
arg2->prepareForChainRuleDerivation(recursive_variables, non_null_chain_rule_derivatives);
set<int>& nnd {non_null_chain_rule_derivatives[const_cast<BinaryOpNode*>(this)]};
set_union(non_null_chain_rule_derivatives.at(arg1).begin(),
non_null_chain_rule_derivatives.at(arg1).end(),
non_null_chain_rule_derivatives.at(arg2).begin(),
non_null_chain_rule_derivatives.at(arg2).end(), inserter(nnd, nnd.begin()));
ranges::set_union(non_null_chain_rule_derivatives.at(arg1),
non_null_chain_rule_derivatives.at(arg2), inserter(nnd, nnd.begin()));
}
expr_t
......@@ -4368,9 +4547,8 @@ BinaryOpNode::cost(const vector<vector<unordered_set<expr_t>>>& blocks_temporary
bool is_matlab) const
{
// For a temporary term, the cost is null
for (const auto& blk_tt : blocks_temporary_terms)
for (const auto& eq_tt : blk_tt)
if (eq_tt.contains(const_cast<BinaryOpNode*>(this)))
if (ranges::any_of(views::join(blocks_temporary_terms),
[this](auto& tt) { return tt.contains(const_cast<BinaryOpNode*>(this)); }))
return 0;
int arg_cost = arg1->cost(blocks_temporary_terms, is_matlab)
......@@ -5938,7 +6116,7 @@ BinaryOpNode::fillAutoregressiveRow(int eqn, const vector<int>& lhs,
tie(vid, lag) = datatree.symbol_table.unrollDiffLeadLagChain(*vid, lag);
if (find(lhs.begin(), lhs.end(), *vid) == lhs.end())
if (ranges::find(lhs, *vid) == lhs.end())
continue;
if (AR.contains({eqn, -lag, *vid}))
......@@ -6017,6 +6195,13 @@ BinaryOpNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
expr_t
BinaryOpNode::substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const
{
return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
}
TrinaryOpNode::TrinaryOpNode(DataTree& datatree_arg, int idx_arg, const expr_t arg1_arg,
TrinaryOpcode op_code_arg, const expr_t arg2_arg,
const expr_t arg3_arg) :
......@@ -6043,11 +6228,9 @@ TrinaryOpNode::prepareForDerivation()
// Non-null derivatives are the union of those of the arguments
// Compute set union of arg{1,2,3}->non_null_derivatives
set<int> non_null_derivatives_tmp;
set_union(arg1->non_null_derivatives.begin(), arg1->non_null_derivatives.end(),
arg2->non_null_derivatives.begin(), arg2->non_null_derivatives.end(),
ranges::set_union(arg1->non_null_derivatives, arg2->non_null_derivatives,
inserter(non_null_derivatives_tmp, non_null_derivatives_tmp.begin()));
set_union(non_null_derivatives_tmp.begin(), non_null_derivatives_tmp.end(),
arg3->non_null_derivatives.begin(), arg3->non_null_derivatives.end(),
ranges::set_union(non_null_derivatives_tmp, arg3->non_null_derivatives,
inserter(non_null_derivatives, non_null_derivatives.begin()));
}
......@@ -6065,12 +6248,9 @@ TrinaryOpNode::prepareForChainRuleDerivation(
set<int>& nnd {non_null_chain_rule_derivatives[const_cast<TrinaryOpNode*>(this)]};
set<int> nnd_tmp;
set_union(non_null_chain_rule_derivatives.at(arg1).begin(),
non_null_chain_rule_derivatives.at(arg1).end(),
non_null_chain_rule_derivatives.at(arg2).begin(),
non_null_chain_rule_derivatives.at(arg2).end(), inserter(nnd_tmp, nnd_tmp.begin()));
set_union(nnd_tmp.begin(), nnd_tmp.end(), non_null_chain_rule_derivatives.at(arg3).begin(),
non_null_chain_rule_derivatives.at(arg3).end(), inserter(nnd, nnd.begin()));
ranges::set_union(non_null_chain_rule_derivatives.at(arg1),
non_null_chain_rule_derivatives.at(arg2), inserter(nnd_tmp, nnd_tmp.begin()));
ranges::set_union(nnd_tmp, non_null_chain_rule_derivatives.at(arg3), inserter(nnd, nnd.begin()));
}
expr_t
......@@ -6186,9 +6366,8 @@ TrinaryOpNode::cost(const vector<vector<unordered_set<expr_t>>>& blocks_temporar
bool is_matlab) const
{
// For a temporary term, the cost is null
for (const auto& blk_tt : blocks_temporary_terms)
for (const auto& eq_tt : blk_tt)
if (eq_tt.contains(const_cast<TrinaryOpNode*>(this)))
if (ranges::any_of(views::join(blocks_temporary_terms),
[this](auto& tt) { return tt.contains(const_cast<TrinaryOpNode*>(this)); }))
return 0;
int arg_cost = arg1->cost(blocks_temporary_terms, is_matlab)
......@@ -6274,9 +6453,9 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
switch (op_code)
{
case TrinaryOpcode::normcdf:
return (0.5 * (1 + erf((v1 - v2) / v3 / numbers::sqrt2)));
return 0.5 * (1 + erf((v1 - v2) / v3 / numbers::sqrt2));
case TrinaryOpcode::normpdf:
return (1 / (v3 * sqrt(2 * numbers::pi) * exp(pow((v1 - v2) / v3, 2) / 2)));
return 1 / (v3 * sqrt(2 * numbers::pi) * exp(pow((v1 - v2) / v3, 2) / 2));
}
__builtin_unreachable(); // Silence GCC warning
}
......@@ -6805,9 +6984,9 @@ TrinaryOpNode::isVariableNodeEqualTo([[maybe_unused]] SymbolType type_arg,
bool
TrinaryOpNode::containsPacExpectation(const string& pac_model_name) const
{
return (arg1->containsPacExpectation(pac_model_name)
return arg1->containsPacExpectation(pac_model_name)
|| arg2->containsPacExpectation(pac_model_name)
|| arg3->containsPacExpectation(pac_model_name));
|| arg3->containsPacExpectation(pac_model_name);
}
bool
......@@ -6861,6 +7040,13 @@ TrinaryOpNode::substituteLogTransform(int orig_symb_id, int aux_symb_id) const
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
expr_t
TrinaryOpNode::substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const
{
return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
}
AbstractExternalFunctionNode::AbstractExternalFunctionNode(DataTree& datatree_arg, int idx_arg,
int symb_id_arg,
vector<expr_t> arguments_arg) :
......@@ -6868,6 +7054,12 @@ AbstractExternalFunctionNode::AbstractExternalFunctionNode(DataTree& datatree_ar
{
}
string
AbstractExternalFunctionNode::getName() const
{
return datatree.symbol_table.getName(symb_id);
}
void
AbstractExternalFunctionNode::prepareForDerivation()
{
......@@ -6881,9 +7073,7 @@ AbstractExternalFunctionNode::prepareForDerivation()
for (int i = 1; i < static_cast<int>(arguments.size()); i++)
{
set<int> non_null_derivatives_tmp;
set_union(non_null_derivatives.begin(), non_null_derivatives.end(),
arguments.at(i)->non_null_derivatives.begin(),
arguments.at(i)->non_null_derivatives.end(),
ranges::set_union(non_null_derivatives, arguments.at(i)->non_null_derivatives,
inserter(non_null_derivatives_tmp, non_null_derivatives_tmp.begin()));
non_null_derivatives = move(non_null_derivatives_tmp);
}
......@@ -6909,8 +7099,7 @@ AbstractExternalFunctionNode::prepareForChainRuleDerivation(
for (int i {1}; i < static_cast<int>(arguments.size()); i++)
{
set<int> nnd_tmp;
set_union(nnd.begin(), nnd.end(), non_null_chain_rule_derivatives.at(arguments.at(i)).begin(),
non_null_chain_rule_derivatives.at(arguments.at(i)).end(),
ranges::set_union(nnd, non_null_chain_rule_derivatives.at(arguments.at(i)),
inserter(nnd_tmp, nnd_tmp.begin()));
nnd = move(nnd_tmp);
}
......@@ -6985,43 +7174,43 @@ AbstractExternalFunctionNode::maxHelper(const function<int(expr_t)>& f) const
int
AbstractExternalFunctionNode::maxEndoLead() const
{
return maxHelper([](expr_t e) { return e->maxEndoLead(); });
return maxHelper(&ExprNode::maxEndoLead);
}
int
AbstractExternalFunctionNode::maxExoLead() const
{
return maxHelper([](expr_t e) { return e->maxExoLead(); });
return maxHelper(&ExprNode::maxExoLead);
}
int
AbstractExternalFunctionNode::maxEndoLag() const
{
return maxHelper([](expr_t e) { return e->maxEndoLag(); });
return maxHelper(&ExprNode::maxEndoLag);
}
int
AbstractExternalFunctionNode::maxExoLag() const
{
return maxHelper([](expr_t e) { return e->maxExoLag(); });
return maxHelper(&ExprNode::maxExoLag);
}
int
AbstractExternalFunctionNode::maxLead() const
{
return maxHelper([](expr_t e) { return e->maxLead(); });
return maxHelper(&ExprNode::maxLead);
}
int
AbstractExternalFunctionNode::maxLag() const
{
return maxHelper([](expr_t e) { return e->maxLag(); });
return maxHelper(&ExprNode::maxLag);
}
int
AbstractExternalFunctionNode::maxLagWithDiffsExpanded() const
{
return maxHelper([](expr_t e) { return e->maxLagWithDiffsExpanded(); });
return maxHelper(&ExprNode::maxLagWithDiffsExpanded);
}
expr_t
......@@ -7159,7 +7348,7 @@ AbstractExternalFunctionNode::substituteUnaryOpNodes(const lag_equivalence_table
int
AbstractExternalFunctionNode::countDiffs() const
{
return maxHelper([](expr_t e) { return e->countDiffs(); });
return maxHelper(&ExprNode::countDiffs);
}
expr_t
......@@ -7218,7 +7407,7 @@ AbstractExternalFunctionNode::computeTemporaryTerms(
expr_t this2 = const_cast<AbstractExternalFunctionNode*>(this);
for (auto& tt : temp_terms_map)
if (find_if(tt.second.cbegin(), tt.second.cend(), sameTefTermPredicate()) != tt.second.cend())
if (ranges::find_if(tt.second, sameTefTermPredicate()) != tt.second.cend())
{
tt.second.insert(this2);
return;
......@@ -7236,7 +7425,7 @@ AbstractExternalFunctionNode::computeBlockTemporaryTerms(
expr_t this2 = const_cast<AbstractExternalFunctionNode*>(this);
for (auto& btt : blocks_temporary_terms)
for (auto& tt : btt)
if (find_if(tt.cbegin(), tt.cend(), sameTefTermPredicate()) != tt.cend())
if (ranges::find_if(tt, sameTefTermPredicate()) != tt.cend())
{
tt.insert(this2);
return;
......@@ -7262,15 +7451,15 @@ AbstractExternalFunctionNode::isVariableNodeEqualTo([[maybe_unused]] SymbolType
bool
AbstractExternalFunctionNode::containsPacExpectation(const string& pac_model_name) const
{
return any_of(arguments.begin(), arguments.end(),
return ranges::any_of(arguments,
[&](expr_t e) { return e->containsPacExpectation(pac_model_name); });
}
bool
AbstractExternalFunctionNode::containsPacTargetNonstationary(const string& pac_model_name) const
{
return any_of(arguments.begin(), arguments.end(),
[&](expr_t e) { return e->containsPacTargetNonstationary(pac_model_name); });
return ranges::any_of(
arguments, [&](expr_t e) { return e->containsPacTargetNonstationary(pac_model_name); });
}
expr_t
......@@ -7294,7 +7483,7 @@ AbstractExternalFunctionNode::removeTrendLeadLag(const map<int, expr_t>& trend_s
bool
AbstractExternalFunctionNode::isInStaticForm() const
{
return all_of(arguments.begin(), arguments.end(), [](expr_t e) { return e->isInStaticForm(); });
return ranges::all_of(arguments, &ExprNode::isInStaticForm);
}
bool
......@@ -7407,6 +7596,13 @@ AbstractExternalFunctionNode::substituteLogTransform(int orig_symb_id, int aux_s
return recurseTransform(&ExprNode::substituteLogTransform, orig_symb_id, aux_symb_id);
}
expr_t
AbstractExternalFunctionNode::substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const
{
return recurseTransform(&ExprNode::substituteAggregationOperators, subst_table, neweqs);
}
expr_t
AbstractExternalFunctionNode::toStatic(DataTree& static_datatree) const
{
......@@ -7507,7 +7703,7 @@ ExternalFunctionNode::writeBytecodeExternalFunctionOutput(
}
code_file << Bytecode::FCALL {nb_output_arguments, static_cast<int>(arguments.size()),
datatree.symbol_table.getName(symb_id), indx, call_type}
getName(), indx, call_type}
<< Bytecode::FSTPTEF {indx};
}
}
......@@ -7516,7 +7712,7 @@ void
ExternalFunctionNode::writeJsonAST(ostream& output) const
{
output << R"({"node_type" : "ExternalFunctionNode", )"
<< R"("name" : ")" << datatree.symbol_table.getName(symb_id) << R"(", "args" : [)";
<< R"("name" : ")" << getName() << R"(", "args" : [)";
writeJsonASTExternalFunctionArguments(output);
output << "]}";
}
......@@ -7540,7 +7736,7 @@ ExternalFunctionNode::writeJsonOutput(ostream& output, const temporary_terms_t&
catch (UnknownFunctionNameAndArgs&)
{
// When writing the JSON output at parsing pass, we don’t use TEF terms
output << datatree.symbol_table.getName(symb_id) << "(";
output << getName() << "(";
writeJsonExternalFunctionArguments(output, temporary_terms, tef_terms, isdynamic);
output << ")";
}
......@@ -7558,8 +7754,8 @@ ExternalFunctionNode::writeOutput(ostream& output, ExprNodeOutputType output_typ
|| output_type == ExprNodeOutputType::epilogueFile
|| output_type == ExprNodeOutputType::occbinDifferenceFile || isLatexOutput(output_type))
{
string name = isLatexOutput(output_type) ? datatree.symbol_table.getTeXName(symb_id)
: datatree.symbol_table.getName(symb_id);
string name
= isLatexOutput(output_type) ? datatree.symbol_table.getTeXName(symb_id) : getName();
output << name << "(";
writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs,
tef_terms);
......@@ -7622,7 +7818,7 @@ ExternalFunctionNode::writeExternalFunctionOutput(
writePrhs(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
output << " mexCallMATLAB(" << nlhs << ", plhs, " << arguments.size() << ", prhs, "
<< R"(")" << datatree.symbol_table.getName(symb_id) << R"(");)" << endl;
<< R"(")" << getName() << R"(");)" << endl;
output << " TEF_" << indx << " = mxGetPr(plhs[0]);" << endl;
if (symb_id == first_deriv_symb_id)
......@@ -7643,7 +7839,7 @@ ExternalFunctionNode::writeExternalFunctionOutput(
else
output << "TEF_" << indx << " = ";
output << datatree.symbol_table.getName(symb_id) << "(";
output << getName() << "(";
writeExternalFunctionArguments(output, output_type, temporary_terms, temporary_terms_idxs,
tef_terms);
output << ");" << endl;
......@@ -7680,7 +7876,7 @@ ExternalFunctionNode::writeJsonExternalFunctionOutput(vector<string>& efout,
if (symb_id == second_deriv_symb_id)
ef << R"(, "external_function_term_dd": "TEFDD_)" << indx << R"(")";
ef << R"(, "value": ")" << datatree.symbol_table.getName(symb_id) << "(";
ef << R"(, "value": ")" << getName() << "(";
writeJsonExternalFunctionArguments(ef, temporary_terms, tef_terms, isdynamic);
ef << R"lit()"}})lit";
efout.push_back(ef.str());
......@@ -7690,7 +7886,6 @@ ExternalFunctionNode::writeJsonExternalFunctionOutput(vector<string>& efout,
void
ExternalFunctionNode::computeXrefs(EquationInfo& ei) const
{
vector<expr_t> dynamic_arguments;
for (auto argument : arguments)
argument->computeXrefs(ei);
}
......@@ -7707,7 +7902,7 @@ ExternalFunctionNode::sameTefTermPredicate() const
{
return [this](expr_t e) {
auto e2 = dynamic_cast<ExternalFunctionNode*>(e);
return (e2 != nullptr && e2->symb_id == symb_id && e2->arguments == arguments);
return e2 != nullptr && e2->symb_id == symb_id && e2->arguments == arguments;
};
}
......@@ -7736,7 +7931,7 @@ void
FirstDerivExternalFunctionNode::writeJsonAST(ostream& output) const
{
output << R"({"node_type" : "FirstDerivExternalFunctionNode", )"
<< R"("name" : ")" << datatree.symbol_table.getName(symb_id) << R"(", "args" : [)";
<< R"("name" : ")" << getName() << R"(", "args" : [)";
writeJsonASTExternalFunctionArguments(output);
output << "]}";
}
......@@ -7874,8 +8069,7 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(
<< "{" << endl
<< " const mwSize dims[2] = {1, " << arguments.size() << "};" << endl
<< " mxArray *plhs[1], *prhs[3];" << endl
<< R"( prhs[0] = mxCreateString(")" << datatree.symbol_table.getName(symb_id)
<< R"(");)" << endl
<< R"( prhs[0] = mxCreateString(")" << getName() << R"(");)" << endl
<< " prhs[1] = mxCreateDoubleScalar(" << inputIndex << ");" << endl
<< " prhs[2] = mxCreateCellArray(2, dims);" << endl;
......@@ -7913,8 +8107,7 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(
{
if (first_deriv_symb_id == ExternalFunctionsTable::IDNotSet)
output << "TEFD_fdd_" << getIndxInTefTerms(symb_id, tef_terms) << "_" << inputIndex
<< " = jacob_element('" << datatree.symbol_table.getName(symb_id) << "',"
<< inputIndex << ",{";
<< " = jacob_element('" << getName() << "'," << inputIndex << ",{";
else
{
tef_terms[{first_deriv_symb_id, arguments}] = static_cast<int>(tef_terms.size());
......@@ -7957,8 +8150,7 @@ FirstDerivExternalFunctionNode::writeJsonExternalFunctionOutput(
<< R"("external_function_term": "TEFD_fdd_)" << getIndxInTefTerms(symb_id, tef_terms) << "_"
<< inputIndex << R"(")"
<< R"(, "analytic_derivative": false)"
<< R"(, "wrt": )" << inputIndex << R"(, "value": ")"
<< datatree.symbol_table.getName(symb_id) << "(";
<< R"(, "wrt": )" << inputIndex << R"(, "value": ")" << getName() << "(";
else
{
tef_terms[{first_deriv_symb_id, arguments}] = static_cast<int>(tef_terms.size());
......@@ -8006,7 +8198,7 @@ FirstDerivExternalFunctionNode::writeBytecodeExternalFunctionOutput(
int nb_output_arguments {1};
Bytecode::FCALL fcall {nb_output_arguments, nb_input_arguments, "jacob_element", indx,
Bytecode::ExternalFunctionCallType::numericalFirstDerivative};
fcall.set_arg_func_name(datatree.symbol_table.getName(symb_id));
fcall.set_arg_func_name(getName());
fcall.set_row(inputIndex);
fcall.set_nb_add_input_arguments(static_cast<int>(arguments.size()));
code_file << fcall << Bytecode::FSTPTEFD {indx, inputIndex};
......@@ -8037,7 +8229,6 @@ FirstDerivExternalFunctionNode::buildSimilarExternalFunctionNode(vector<expr_t>&
void
FirstDerivExternalFunctionNode::computeXrefs(EquationInfo& ei) const
{
vector<expr_t> dynamic_arguments;
for (auto argument : arguments)
argument->computeXrefs(ei);
}
......@@ -8049,12 +8240,12 @@ FirstDerivExternalFunctionNode::sameTefTermPredicate() const
if (first_deriv_symb_id == symb_id)
return [this](expr_t e) {
auto e2 = dynamic_cast<ExternalFunctionNode*>(e);
return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
};
else
return [this](expr_t e) {
auto e2 = dynamic_cast<FirstDerivExternalFunctionNode*>(e);
return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
};
}
......@@ -8078,7 +8269,7 @@ void
SecondDerivExternalFunctionNode::writeJsonAST(ostream& output) const
{
output << R"({"node_type" : "SecondDerivExternalFunctionNode", )"
<< R"("name" : ")" << datatree.symbol_table.getName(symb_id) << R"(", "args" : [)";
<< R"("name" : ")" << getName() << R"(", "args" : [)";
writeJsonASTExternalFunctionArguments(output);
output << "]}";
}
......@@ -8206,8 +8397,7 @@ SecondDerivExternalFunctionNode::writeExternalFunctionOutput(
<< "{" << endl
<< " const mwSize dims[2]= {1, " << arguments.size() << "};" << endl
<< " mxArray *plhs[1], *prhs[4];" << endl
<< R"( prhs[0] = mxCreateString(")" << datatree.symbol_table.getName(symb_id)
<< R"(");)" << endl
<< R"( prhs[0] = mxCreateString(")" << getName() << R"(");)" << endl
<< " prhs[1] = mxCreateDoubleScalar(" << inputIndex1 << ");" << endl
<< " prhs[2] = mxCreateDoubleScalar(" << inputIndex2 << ");" << endl
<< " prhs[3] = mxCreateCellArray(2, dims);" << endl;
......@@ -8246,8 +8436,8 @@ SecondDerivExternalFunctionNode::writeExternalFunctionOutput(
{
if (second_deriv_symb_id == ExternalFunctionsTable::IDNotSet)
output << "TEFDD_fdd_" << getIndxInTefTerms(symb_id, tef_terms) << "_" << inputIndex1 << "_"
<< inputIndex2 << " = hess_element('" << datatree.symbol_table.getName(symb_id)
<< "'," << inputIndex1 << "," << inputIndex2 << ",{";
<< inputIndex2 << " = hess_element('" << getName() << "'," << inputIndex1 << ","
<< inputIndex2 << ",{";
else
{
tef_terms[{second_deriv_symb_id, arguments}] = static_cast<int>(tef_terms.size());
......@@ -8291,7 +8481,7 @@ SecondDerivExternalFunctionNode::writeJsonExternalFunctionOutput(
<< inputIndex1 << "_" << inputIndex2 << R"(")"
<< R"(, "analytic_derivative": false)"
<< R"(, "wrt1": )" << inputIndex1 << R"(, "wrt2": )" << inputIndex2 << R"(, "value": ")"
<< datatree.symbol_table.getName(symb_id) << "(";
<< getName() << "(";
else
{
tef_terms[{second_deriv_symb_id, arguments}] = static_cast<int>(tef_terms.size());
......@@ -8317,7 +8507,6 @@ SecondDerivExternalFunctionNode::buildSimilarExternalFunctionNode(vector<expr_t>
void
SecondDerivExternalFunctionNode::computeXrefs(EquationInfo& ei) const
{
vector<expr_t> dynamic_arguments;
for (auto argument : arguments)
argument->computeXrefs(ei);
}
......@@ -8381,7 +8570,7 @@ SecondDerivExternalFunctionNode::writeBytecodeExternalFunctionOutput(
{
Bytecode::FCALL fcall {1, 0, "hess_element", indx,
Bytecode::ExternalFunctionCallType::numericalSecondDerivative};
fcall.set_arg_func_name(datatree.symbol_table.getName(symb_id));
fcall.set_arg_func_name(getName());
fcall.set_row(inputIndex1);
fcall.set_col(inputIndex2);
fcall.set_nb_add_input_arguments(static_cast<int>(arguments.size()));
......@@ -8406,12 +8595,12 @@ SecondDerivExternalFunctionNode::sameTefTermPredicate() const
if (second_deriv_symb_id == symb_id)
return [this](expr_t e) {
auto e2 = dynamic_cast<ExternalFunctionNode*>(e);
return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
};
else
return [this](expr_t e) {
auto e2 = dynamic_cast<SecondDerivExternalFunctionNode*>(e);
return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
};
}
......@@ -8772,6 +8961,13 @@ SubModelNode::substituteLogTransform([[maybe_unused]] int orig_symb_id,
return const_cast<SubModelNode*>(this);
}
expr_t
SubModelNode::substituteAggregationOperators([[maybe_unused]] subst_table_t& subst_table,
[[maybe_unused]] vector<BinaryOpNode*>& neweqs) const
{
return const_cast<SubModelNode*>(this);
}
VarExpectationNode::VarExpectationNode(DataTree& datatree_arg, int idx_arg, string model_name_arg) :
SubModelNode {datatree_arg, idx_arg, move(model_name_arg)}
{
......@@ -9162,8 +9358,7 @@ VariableNode::matchVTCTPHelper(optional<int>& var_id, int& lag, optional<int>& p
param_id = symb_id;
}
else
throw MatchFailureException {"Symbol " + datatree.symbol_table.getName(symb_id)
+ " not allowed here"};
throw MatchFailureException {"Symbol " + getName() + " not allowed here"};
}
void
......@@ -9291,9 +9486,9 @@ ExprNode::matchParamTimesTargetMinusVariable(int symb_id) const
auto& avi = datatree.symbol_table.getAuxVarInfo(target->symb_id);
if (avi.type == AuxVarType::pacTargetNonstationary && target->lag == -1)
return true;
return (avi.type == AuxVarType::unaryOp && avi.unary_op == "log" && avi.orig_symb_id
return avi.type == AuxVarType::unaryOp && avi.unary_op == "log" && avi.orig_symb_id
&& !datatree.symbol_table.isAuxiliaryVariable(*avi.orig_symb_id)
&& target->lag + avi.orig_lead_lag.value() == -1);
&& target->lag + avi.orig_lead_lag.value() == -1;
}
else
return target->lag == -1;
......@@ -9370,3 +9565,83 @@ ExprNode::toString() const
writeJsonOutput(ss, {}, {});
return ss.str();
}
tuple<int, expr_t, expr_t>
ExprNode::matchComplementarityCondition(
[[maybe_unused]] const optional<int>& heterogeneity_dimension) const
{
throw MatchFailureException {"This expression is not an inequality"};
}
tuple<int, expr_t, expr_t>
BinaryOpNode::matchComplementarityCondition(const optional<int>& heterogeneity_dimension) const
{
bool is_greater {[&] {
switch (op_code)
{
case BinaryOpcode::less:
case BinaryOpcode::lessEqual:
return false;
case BinaryOpcode::greater:
case BinaryOpcode::greaterEqual:
return true;
default:
throw MatchFailureException {"This expression is not an inequality"};
}
}()};
auto match_contemporaneous_endogenous = [&](expr_t e) -> optional<int> {
auto* ve = dynamic_cast<VariableNode*>(e);
if (ve && ve->lag == 0
&& ((!heterogeneity_dimension
&& datatree.symbol_table.getType(ve->symb_id) == SymbolType::endogenous)
|| (heterogeneity_dimension
&& datatree.symbol_table.getType(ve->symb_id) == SymbolType::heterogeneousEndogenous
&& datatree.symbol_table.getHeterogeneityDimension(ve->symb_id)
== *heterogeneity_dimension)))
return ve->symb_id;
else
return nullopt;
};
auto check_bound_constant = [](expr_t e) {
if (!e->isConstant())
throw MatchFailureException {"Bounds must not contain any endogenous or exogenous variable"};
};
// Match “endogenous ≶ bound”
if (auto id = match_contemporaneous_endogenous(arg1); id)
{
check_bound_constant(arg2);
return {*id, is_greater ? arg2 : nullptr, is_greater ? nullptr : arg2};
}
// Match “bound ≶ endogenous”
if (auto id = match_contemporaneous_endogenous(arg2); id)
{
check_bound_constant(arg1);
return {*id, is_greater ? nullptr : arg1, is_greater ? arg1 : nullptr};
}
// Match: “bound < endogenous < bound” or “bound > endogenous > bound”
/* NB: we exploit the fact that inequality signs are left-associative, hence the constraint is
necessarily “(bound < endogenous) < bound” or “(bound > endogenous) > bound” (assuming the user
did not add a spurious parenthesis). */
auto barg1 = dynamic_cast<BinaryOpNode*>(arg1);
if (!(barg1
&& ((is_greater
&& (barg1->op_code == BinaryOpcode::greater
|| barg1->op_code == BinaryOpcode::greaterEqual))
|| (!is_greater
&& (barg1->op_code == BinaryOpcode::less
|| barg1->op_code == BinaryOpcode::lessEqual)))))
throw MatchFailureException {};
auto id = match_contemporaneous_endogenous(barg1->arg2);
if (!id)
throw MatchFailureException {};
check_bound_constant(barg1->arg1);
check_bound_constant(arg2);
return {*id, is_greater ? arg2 : barg1->arg1, is_greater ? barg1->arg1 : arg2};
}
src/ExprNode.hh
View file @ 6618d7d1
/*
* Copyright © 2007-2024 Dynare Team
* Copyright © 2007-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -219,13 +219,13 @@ RIGHT_ARRAY_SUBSCRIPT(ExprNodeOutputType output_type)
}
// Left and right parentheses
inline string
constexpr string
LEFT_PAR(ExprNodeOutputType output_type)
{
return isLatexOutput(output_type) ? "\\left(" : "(";
}
inline string
constexpr string
RIGHT_PAR(ExprNodeOutputType output_type)
{
return isLatexOutput(output_type) ? "\\right)" : ")";
......@@ -283,7 +283,7 @@ protected:
min_cost(bool is_matlab)
{
return is_matlab ? min_cost_matlab : min_cost_c;
};
}
//! Initializes data member non_null_derivatives
virtual void prepareForDerivation() = 0;
......@@ -572,38 +572,35 @@ public:
{
};
//! Returns the maximum lead of endogenous in this expression
//! Returns the maximum lead of endogenous in this expression (not incl. heterogeneous endo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxEndoLead() const = 0;
//! Returns the maximum lead of exogenous in this expression
//! Returns the maximum lead of exogenous in this expression (not incl. heterogeneous exo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxExoLead() const = 0;
//! Returns the maximum lag of endogenous in this expression
//! Returns the maximum lag of endogenous in this expression (not incl. heterogeneous endo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxEndoLag() const = 0;
//! Returns the maximum lag of exogenous in this expression
//! Returns the maximum lag of exogenous in this expression (not incl. heterogeneous exo)
/*! Always returns a non-negative value */
[[nodiscard]] virtual int maxExoLag() const = 0;
//! Returns the maximum lead of endo/exo/exodet in this expression
/*! A negative value means that the expression contains only lagged
variables. A value of numeric_limits<int>::min() means that there is
no variable. */
/* Returns the maximum lead of endo/exo/exodet in this expression (including heterogeneous
endo/exo). A negative value means that the expression contains only lagged variables. A value
of numeric_limits<int>::min() means that there is no variable. */
[[nodiscard]] virtual int maxLead() const = 0;
//! Returns the maximum lag of endo/exo/exodet in this expression
/*! A negative value means that the expression contains only leaded
variables. A value of numeric_limits<int>::min() means that there is
no variable. */
/* Returns the maximum lag of endo/exo/exodet in this expression (including heterogeneous
endo/exo). A negative value means that the expression contains only leaded variables. A value
of numeric_limits<int>::min() means that there is no variable. */
[[nodiscard]] virtual int maxLag() const = 0;
//! Returns the maximum lag of endo/exo/exodet, as if diffs were expanded
/*! This function behaves as maxLag(), except that it treats diff()
differently. For e.g., on diff(diff(x(-1))), maxLag() returns 1 while
maxLagWithDiffsExpanded() returns 3. */
/* Returns the maximum lag of endo/exo/exodet (including heterogeneous endo/exo), as if diffs were
expanded. This function behaves as maxLag(), except that it treats diff() differently. For
e.g., on diff(diff(x(-1))), maxLag() returns 1 while maxLagWithDiffsExpanded() returns 3. */
[[nodiscard]] virtual int maxLagWithDiffsExpanded() const = 0;
[[nodiscard]] virtual expr_t undiff() const = 0;
......@@ -942,6 +939,19 @@ public:
// Substitutes orig_symb_id(±l) with exp(aux_symb_id(±l)) (used for “var(log)”)
[[nodiscard]] virtual expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const = 0;
/* Matches an expression that constitutes a complementarity condition.
If successful, returns a triplet (endo_symb_id, lower_bound, upper_bound).
Otherwise, throws a MatchFailureException. */
[[nodiscard]] virtual tuple<int, expr_t, expr_t>
matchComplementarityCondition(const optional<int>& heterogeneity_dimension = nullopt) const;
/* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
Also declares those aggregation operators in the HeterogeneityTable, so as to
compute their index in the dedicated vector in argument of the dynamic/static files. */
[[nodiscard]] virtual expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const
= 0;
};
//! Object used to compare two nodes (using their indexes)
......@@ -1053,6 +1063,8 @@ public:
= "") const override;
[[nodiscard]] bool isParamTimesEndogExpr() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
};
//! Symbol or variable node
......@@ -1085,6 +1097,10 @@ protected:
public:
VariableNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg, int lag_arg);
[[nodiscard]] SymbolType get_type() const;
[[nodiscard]] string getName() const;
[[nodiscard]] int getDerivID() const;
[[nodiscard]] int getTypeSpecificID() const;
void writeOutput(ostream& output, ExprNodeOutputType output_type,
const temporary_terms_t& temporary_terms,
const temporary_terms_idxs_t& temporary_terms_idxs,
......@@ -1102,7 +1118,6 @@ public:
const deriv_node_temp_terms_t& tef_terms) const override;
expr_t toStatic(DataTree& static_datatree) const override;
void computeXrefs(EquationInfo& ei) const override;
[[nodiscard]] SymbolType get_type() const;
BinaryOpNode* normalizeEquationHelper(const set<expr_t>& contain_var, expr_t rhs) const override;
[[nodiscard]] int maxEndoLead() const override;
[[nodiscard]] int maxExoLead() const override;
......@@ -1157,6 +1172,8 @@ public:
vector<int>& powers) const override;
[[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
};
//! Unary operator node
......@@ -1304,6 +1321,8 @@ public:
[[nodiscard]] bool isParamTimesEndogExpr() const override;
void decomposeAdditiveTerms(vector<pair<expr_t, int>>& terms, int current_sign) const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
};
//! Binary operator node
......@@ -1454,7 +1473,7 @@ public:
[[nodiscard]] expr_t unpackPowerDeriv() const;
//! Returns MULT_i*(lhs-rhs) = 0, creating multiplier MULT_i
expr_t addMultipliersToConstraints(int i);
//! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
//! Returns the non-zero hand side of an equation (that must have a hand side equal to zero)
[[nodiscard]] expr_t getNonZeroPartofEquation() const;
[[nodiscard]] bool isInStaticForm() const override;
void fillAutoregressiveRow(int eqn, const vector<int>& lhs,
......@@ -1499,6 +1518,11 @@ public:
vector<int>& powers) const override;
[[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
[[nodiscard]] tuple<int, expr_t, expr_t>
matchComplementarityCondition(const optional<int>& heterogeneity_dimension
= nullopt) const override;
};
//! Trinary operator node
......@@ -1642,6 +1666,8 @@ public:
[[nodiscard]] bool containsPacTargetNonstationary(const string& pac_model_name
= "") const override;
[[nodiscard]] bool isParamTimesEndogExpr() const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
};
......@@ -1714,6 +1740,7 @@ protected:
public:
AbstractExternalFunctionNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg,
vector<expr_t> arguments_arg);
[[nodiscard]] string getName() const;
void computeTemporaryTerms(const pair<int, int>& derivOrder,
map<pair<int, int>, unordered_set<expr_t>>& temp_terms_map,
unordered_map<expr_t, pair<int, pair<int, int>>>& reference_count,
......@@ -1816,6 +1843,8 @@ public:
= "") const override;
[[nodiscard]] bool isParamTimesEndogExpr() const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
};
class ExternalFunctionNode : public AbstractExternalFunctionNode
......@@ -2018,6 +2047,8 @@ public:
expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
[[nodiscard]] expr_t removeTrendLeadLag(const map<int, expr_t>& trend_symbols_map) const override;
[[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
vector<BinaryOpNode*>& neweqs) const override;
protected:
void prepareForDerivation() override;
......
src/ExtendedPreprocessorTypes.hh
View file @ 6618d7d1
......@@ -24,6 +24,7 @@
enum class OutputType
{
standard, // Default value, infer the derivation order from .mod file only
first, // Output only 1st dynamic derivatives with no other computations
second, // Output at least 2nd dynamic derivatives
third, // Output at least 3rd dynamic derivatives
};
......
src/HeterogeneityTable.cc 0 → 100644
View file @ 6618d7d1
/*
* Copyright © 2024 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 <https://www.gnu.org/licenses/>.
*/
#include <cassert>
#include <utility>
#include "HeterogeneityTable.hh"
#include "SymbolTable.hh"
void
HeterogeneityTable::setSymbolTable(SymbolTable* symbol_table_arg)
{
symbol_table = symbol_table_arg;
}
int
HeterogeneityTable::addDimension(string name)
{
if (name_to_id.contains(name))
throw AlreadyDeclaredDimensionException {move(name)};
int id {static_cast<int>(id_to_name.size())};
name_to_id.emplace(name, id);
id_to_name.push_back(move(name));
return id;
}
bool
HeterogeneityTable::exists(const string& name) const
{
return name_to_id.contains(name);
}
int
HeterogeneityTable::getID(const string& name) const
{
if (auto it = name_to_id.find(name); it != name_to_id.end())
return it->second;
else
throw UnknownDimensionNameException {name};
}
string
HeterogeneityTable::getName(int id) const
{
if (id < 0 || id >= static_cast<int>(id_to_name.size()))
throw UnknownDimensionIDException {id};
else
return id_to_name[id];
}
bool
HeterogeneityTable::empty() const
{
return name_to_id.empty();
}
vector<string>
HeterogeneityTable::getDimensions() const
{
return id_to_name;
}
int
HeterogeneityTable::size() const
{
return static_cast<int>(name_to_id.size());
}
void
HeterogeneityTable::addSummedHeterogeneousEndogenous(int symb_id)
{
assert(symbol_table->getType(symb_id) == SymbolType::heterogeneousEndogenous);
if (summed_het_endo_to_index.contains(symb_id))
throw AlreadyDeclaredSummedHeterogeneousEndogenousException {symb_id};
int index {static_cast<int>(index_to_summed_het_endo.size())};
summed_het_endo_to_index.emplace(symb_id, index);
index_to_summed_het_endo.push_back(symb_id);
}
int
HeterogeneityTable::getSummedHeterogenousEndogenousIndex(int symb_id) const
{
if (auto it = summed_het_endo_to_index.find(symb_id); it != summed_het_endo_to_index.end())
return it->second;
else
throw UnknownSummedHeterogeneousEndogenousException {symb_id};
}
int
HeterogeneityTable::aggregateEndoSize() const
{
return index_to_summed_het_endo.size();
}
void
HeterogeneityTable::writeOutput(ostream& output) const
{
for (size_t id {0}; id < id_to_name.size(); id++)
output << "M_.heterogeneity(" << id + 1 << ").dimension_name = '" << id_to_name[id] << "';"
<< endl;
output << "M_.heterogeneity_aggregates = {" << endl;
for (int symb_id : index_to_summed_het_endo)
output << "'sum', " << symbol_table->getHeterogeneityDimension(symb_id) + 1 << ", "
<< symbol_table->getTypeSpecificID(symb_id) + 1 << ";" << endl;
output << "};" << endl;
}
void
HeterogeneityTable::writeJsonOutput(ostream& output) const
{
assert(!empty());
output << R"("heterogeneity_dimension": [)";
for (bool first_written {false}; const auto& dim : id_to_name)
{
if (exchange(first_written, true))
output << ", ";
output << '"' << dim << '"';
}
output << "]" << endl;
}
src/HeterogeneityTable.hh 0 → 100644
View file @ 6618d7d1
/*
* Copyright © 2024 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 <https://www.gnu.org/licenses/>.
*/
#ifndef HETEROGENEITY_TABLE_HH
#define HETEROGENEITY_TABLE_HH
#include <map>
#include <ostream>
#include <string>
#include <vector>
using namespace std;
class SymbolTable; // Forward declaration, to avoid circularity
/*
There is a guarantee that heterogeneity IDs are increasing, i.e. if dimension A is added after
dimension B, then the ID of A is greater than the ID of B.
Moreover, the IDs form a contiguous interval starting at 0.
*/
class HeterogeneityTable
{
private:
// Maps dimension names to IDs
map<string, int> name_to_id;
// Maps dimension IDs to names
vector<string> id_to_name;
SymbolTable* symbol_table {nullptr}; // Cannot be a reference, because of circularity
/* Keeps track of the SUM() operator instances.
Maps a symbol ID that appears inside a SUM() operator into an index in
M_.heterogeneity_aggregates */
map<int, int> summed_het_endo_to_index;
// Maps an index in M_.heterogeneity_aggregates into a symbol ID
vector<int> index_to_summed_het_endo;
public:
void setSymbolTable(SymbolTable* symbol_table_arg);
struct AlreadyDeclaredDimensionException
{
// Dimension name
const string name;
};
struct UnknownDimensionNameException
{
// Dimension name
const string name;
};
struct UnknownDimensionIDException
{
// Dimension ID
const int id;
};
// Returns the dimension ID
int addDimension(string name);
[[nodiscard]] bool exists(const string& name) const;
[[nodiscard]] int getID(const string& name) const;
[[nodiscard]] string getName(int id) const;
[[nodiscard]] bool empty() const;
[[nodiscard]] vector<string> getDimensions() const;
[[nodiscard]] int size() const;
struct AlreadyDeclaredSummedHeterogeneousEndogenousException
{
const int symb_id;
};
struct UnknownSummedHeterogeneousEndogenousException
{
const int symb_id;
};
void addSummedHeterogeneousEndogenous(int symb_id);
[[nodiscard]] int getSummedHeterogenousEndogenousIndex(int symb_id) const;
[[nodiscard]] int aggregateEndoSize() const;
void writeOutput(ostream& output) const;
void writeJsonOutput(ostream& output) const;
};
#endif
src/HeterogeneousModel.cc 0 → 100644
View file @ 6618d7d1
/*
* Copyright © 2024-2025 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 <https://www.gnu.org/licenses/>.
*/
#include <cstdlib>
#include <iostream>
#include "HeterogeneousModel.hh"
HeterogeneousModel::HeterogeneousModel(SymbolTable& symbol_table_arg,
NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
int heterogeneity_dimension_arg) :
ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
heterogeneity_table_arg, true},
heterogeneity_dimension {heterogeneity_dimension_arg}
{
}
HeterogeneousModel&
HeterogeneousModel::operator=(const HeterogeneousModel& m)
{
ModelTree::operator=(m);
assert(heterogeneity_dimension == m.heterogeneity_dimension);
deriv_id_table = m.deriv_id_table;
inv_deriv_id_table = m.inv_deriv_id_table;
return *this;
}
void
HeterogeneousModel::computeChainRuleJacobian()
{
cerr << "Heterogeneous::computeChainRuleJacobian(): unimplemented" << endl;
exit(EXIT_FAILURE);
}
int
HeterogeneousModel::getBlockJacobianEndoCol([[maybe_unused]] int blk, [[maybe_unused]] int var,
[[maybe_unused]] int lead_lag) const
{
cerr << "Heterogeneous::getBlockJacobianEndoCol(): unimplemented" << endl;
exit(EXIT_FAILURE);
}
int
HeterogeneousModel::getMFS() const
{
cerr << "Heterogeneous::getMFS(): unimplemented" << endl;
exit(EXIT_FAILURE);
}
void
HeterogeneousModel::computeDerivIDs()
{
set<pair<int, int>> dynvars;
for (auto& equation : equations)
{
equation->collectDynamicVariables(SymbolType::heterogeneousEndogenous, dynvars);
equation->collectDynamicVariables(SymbolType::heterogeneousExogenous, dynvars);
equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
equation->collectDynamicVariables(SymbolType::parameter, dynvars);
}
for (const auto& [symb_id, lead_lag] : dynvars)
{
auto type {symbol_table.getType(symb_id)};
if (isHeterogeneous(type))
assert(symbol_table.getHeterogeneityDimension(symb_id) == heterogeneity_dimension);
if (type == SymbolType::heterogeneousEndogenous || type == SymbolType::endogenous)
assert(abs(lead_lag) <= 1);
if (type == SymbolType::heterogeneousExogenous || type == SymbolType::exogenous)
assert(lead_lag == 0);
int deriv_id {static_cast<int>(deriv_id_table.size())};
deriv_id_table.emplace(pair {symb_id, lead_lag}, deriv_id);
inv_deriv_id_table.emplace_back(symb_id, lead_lag);
}
}
/*
* Unfold complementarity conditions: (i) declare the multipliers associated
* with each bound constraint μ_l and μ_u ; (ii) add or substract the
* multiplier into the associated condition; (iii) add the the complementarity
* slackness conditions into the set of equations. For example,
* households choose {cₜ, aₜ₊₁} to maximize expected lifetime utility:
* max 𝐸ₜ [∑ₛ₌₀^∞ βˢ · u(cₜ₊ₛ)]
*
* Subject to:
* 1. Budget constraint: cₜ + aₜ₊₁ = yₜ + (1 + rₜ) · aₜ
* 2. Borrowing constraint: aₜ₊₁ ≥ aₘᵢₙ
*
* Let u'(cₜ) denote the marginal utility of consumption.
* Let μₜ ≥ 0 be the Lagrange multiplier on the borrowing constraint.
*
* Then, the Euler equation becomes:
* u′(cₜ) = β · (1 + rₜ₊₁) · u′(cₜ₊₁) − μₜ
*
* Together with:
* aₜ₊₁ ≥ aₘᵢₙ [primal feasibility]
* μₜ ≥ 0 [dual feasibility]
* μₜ · (aₜ₊₁ − aₘᵢₙ) = 0 [complementarity slackness]
* Note that the primal feasibility and dual feasibility constraints are not
* introduced here, but Bhandari et al. (2023) show in Appendix B.1 that they
* are redundant.
*/
void
HeterogeneousModel::transformPass()
{
for (int i = 0; i < static_cast<int>(equations.size()); ++i)
{
if (!complementarity_conditions[i])
continue;
/*
* `const auto& [symb_id, lb, ub] = *complementarity_conditions[i];` was not used here because
* the call to `addEquation` may eventually lead to a resize of the
* `complementarity_conditions` vector, which may invalidate the reference to its element. We
* take a copy instead for safety.
*/
auto [symb_id, lb, ub] = *complementarity_conditions[i];
VariableNode* var = getVariable(symb_id);
if (lb)
{
int mu_id = symbol_table.addHeterogeneousMultiplierAuxiliaryVar(
heterogeneity_dimension, i, "MULT_L_" + symbol_table.getName(symb_id));
expr_t mu_L = AddVariable(mu_id);
auto substeq = AddEqual(AddPlus(equations[i]->arg1, mu_L), equations[i]->arg2);
assert(substeq);
equations[i] = substeq;
addEquation(AddEqual(AddTimes(mu_L, AddMinus(var, lb)), Zero), nullopt);
}
if (ub)
{
int mu_id = symbol_table.addHeterogeneousMultiplierAuxiliaryVar(
heterogeneity_dimension, i, "MULT_U_" + symbol_table.getName(symb_id));
auto mu_U = AddVariable(mu_id);
auto substeq = AddEqual(AddMinus(equations[i]->arg1, mu_U), equations[i]->arg2);
assert(substeq);
equations[i] = substeq;
addEquation(AddEqual(AddTimes(mu_U, AddMinus(ub, var)), Zero), nullopt);
}
}
}
void
HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll)
{
assert(!use_dll); // Not yet implemented
computeDerivIDs();
set<int> vars;
for (auto& [symb_lag, deriv_id] : deriv_id_table)
if (symbol_table.getType(symb_lag.first) != SymbolType::parameter)
vars.insert(deriv_id);
cout << "Computing " << modelClassName() << " derivatives (order " << derivsOrder << ")." << endl;
computeDerivatives(derivsOrder, vars);
computeTemporaryTerms(!use_dll, no_tmp_terms);
computeMCPEquationsReordering(heterogeneity_dimension);
}
void
HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
{
assert(!julia); // Not yet implemented
writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
writeComplementarityConditionsFile<true>(basename, heterogeneity_dimension);
}
int
HeterogeneousModel::getJacobianCol(int deriv_id, bool sparse) const
{
assert(sparse);
SymbolType type {getTypeByDerivID(deriv_id)};
int tsid {getTypeSpecificIDByDerivID(deriv_id)};
int lag {getLagByDerivID(deriv_id)};
if (type == SymbolType::heterogeneousEndogenous)
return tsid + (lag + 1) * symbol_table.het_endo_nbr(heterogeneity_dimension);
int shift {3 * symbol_table.het_endo_nbr(heterogeneity_dimension)};
if (type == SymbolType::heterogeneousExogenous)
return shift + tsid;
shift += symbol_table.het_exo_nbr(heterogeneity_dimension);
if (type == SymbolType::endogenous)
return shift + tsid + (lag + 1) * symbol_table.endo_nbr();
shift += symbol_table.endo_nbr();
if (type == SymbolType::exogenous)
return shift + tsid;
throw UnknownDerivIDException();
}
int
HeterogeneousModel::getJacobianColsNbr(bool sparse) const
{
assert(sparse);
return 3 * (symbol_table.het_endo_nbr(heterogeneity_dimension) + symbol_table.endo_nbr())
+ symbol_table.het_exo_nbr(heterogeneity_dimension) + symbol_table.exo_nbr();
}
SymbolType
HeterogeneousModel::getTypeByDerivID(int deriv_id) const noexcept(false)
{
return symbol_table.getType(getSymbIDByDerivID(deriv_id));
}
int
HeterogeneousModel::getLagByDerivID(int deriv_id) const noexcept(false)
{
if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
throw UnknownDerivIDException();
return inv_deriv_id_table[deriv_id].second;
}
int
HeterogeneousModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
{
if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
throw UnknownDerivIDException();
return inv_deriv_id_table[deriv_id].first;
}
int
HeterogeneousModel::getTypeSpecificIDByDerivID(int deriv_id) const
{
return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id));
}
int
HeterogeneousModel::getDerivID(int symb_id, int lead_lag) const noexcept(false)
{
if (auto it = deriv_id_table.find({symb_id, lead_lag}); it == deriv_id_table.end())
throw UnknownDerivIDException();
else
return it->second;
}
void
HeterogeneousModel::writeDriverOutput(ostream& output) const
{
std::vector<int> state_var;
for (int endoID = 0; endoID < symbol_table.het_endo_nbr(heterogeneity_dimension); endoID++)
try
{
getDerivID(symbol_table.getID(SymbolType::heterogeneousEndogenous, endoID,
heterogeneity_dimension),
-1);
if (ranges::find(state_var, endoID) == state_var.end())
state_var.push_back(endoID);
}
catch (UnknownDerivIDException& e)
{
}
output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").state_var = [";
for (int it : state_var)
output << it + 1 << " ";
output << "];" << endl;
output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").dynamic_tmp_nbr = [";
for (const auto& it : temporary_terms_derivatives)
output << it.size() << "; ";
output << "];" << endl;
writeDriverSparseIndicesHelper(
"heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").dynamic", output);
output << "M_.heterogeneity(" << heterogeneity_dimension + 1
<< ").dynamic_mcp_equations_reordering = [";
for (auto i : mcp_equations_reordering)
output << i + 1 << "; ";
output << "];" << endl;
}
src/HeterogeneousModel.hh 0 → 100644
View file @ 6618d7d1
/*
* Copyright © 2024 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 <https://www.gnu.org/licenses/>.
*/
#ifndef HETEROGENEOUS_MODEL_HH
#define HETEROGENEOUS_MODEL_HH
#include <string>
#include "ModelTree.hh"
using namespace std;
class HeterogeneousModel : public ModelTree
{
public:
const int heterogeneity_dimension;
HeterogeneousModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg, int heterogeneity_dimension_arg);
HeterogeneousModel(const HeterogeneousModel& m) = default;
HeterogeneousModel& operator=(const HeterogeneousModel& m);
void transformPass();
void computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll);
void writeModelFiles(const string& basename, bool julia) const;
void writeDriverOutput(ostream& output) const;
[[nodiscard]] int getJacobianCol(int deriv_id, bool sparse) const override;
[[nodiscard]] int getJacobianColsNbr(bool sparse) const override;
#if 0
void substituteEndoLeadGreaterThanTwo();
//! Transforms the model by removing all lags greater or equal than 2 on endos
void substituteEndoLagGreaterThanTwo();
//! Transforms the model by removing all leads on exos
/*! Note that this can create new lags on endos and exos */
void substituteExoLead();
//! Transforms the model by removing all lags on exos
void substituteExoLag();
//! Transforms the model by removing all UnaryOpcode::expectation
void substituteExpectation(bool partial_information_model);
//! Transforms the model by decreasing the lead/lag of predetermined variables in model equations
//! by one
void transformPredeterminedVariables();
//! Substitutes out all model-local variables
void substituteModelLocalVariables();
#endif
// FIXME: the following 5 functions are identical to those in DynamicModel. Factorization?
[[nodiscard]] int getDerivID(int symb_id, int lead_lag) const noexcept(false) override;
[[nodiscard]] SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
[[nodiscard]] int getLagByDerivID(int deriv_id) const noexcept(false) override;
[[nodiscard]] int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
[[nodiscard]] int getTypeSpecificIDByDerivID(int deriv_id) const override;
protected:
void computeChainRuleJacobian() override;
int getBlockJacobianEndoCol(int blk, int var, int lead_lag) const override;
string
modelClassName() const override
{
return "dynamic model for heterogeneity dimension '"
+ heterogeneity_table.getName(heterogeneity_dimension) + "'";
}
int getMFS() const override;
private:
// Maps a pair (symbol ID, lead/lag) to a deriv ID
map<pair<int, int>, int> deriv_id_table;
// Maps a deriv ID to a pair (symbol ID, lead/lag)
vector<pair<int, int>> inv_deriv_id_table;
// Allocates the derivation IDs for all endogenous variables for this heterogeneity dimension
void computeDerivIDs();
};
#endif
src/MacroExpandModFile.cc
View file @ 6618d7d1
/*
* Copyright © 2015-2023 Dynare Team
* Copyright © 2015-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -41,7 +41,7 @@ macroExpandModFile(const filesystem::path& filename, const istream& modfile, boo
if (save_macro)
{
if (save_macro_file.empty())
save_macro_file = filename.stem().string() + "-macroexp.mod";
save_macro_file = filename.stem().string() + "_macroexp.mod";
ofstream macro_output_file {save_macro_file};
if (macro_output_file.fail())
{
......@@ -52,15 +52,8 @@ macroExpandModFile(const filesystem::path& filename, const istream& modfile, boo
string str(macro_output.str());
if (!line_macro)
{
/* Remove the @#line directives.
Unfortunately GCC 11 does not yet support std::regex::multiline
(despite it being in the C++17 standard), so we are forced to use
a trick to emulate the “usual” behaviour of the caret ^;
here, the latter only matches the beginning of file.
This also means that we are forced to remove the EOL before the
@#line, and not the one after it (matching the EOL before and the
EOL after in the same regexp does not work). */
str = regex_replace(str, regex(R"((^|\r?\n)@#line.*)"), "");
// Remove the @#line directives.
str = regex_replace(str, regex(R"(^@#line.*$)", std::regex::multiline), "");
/* Remove the EOLs at the beginning of the output, the first one
being a remnant of the first @#line directive. */
str = regex_replace(str, regex(R"(^(\r?\n)+)"), "");
......
src/ModFile.cc
View file @ 6618d7d1
/*
* Copyright © 2006-2023 Dynare Team
* Copyright © 2006-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -31,25 +31,48 @@
#include "Shocks.hh"
ModFile::ModFile(WarningConsolidation& warnings_arg) :
symbol_table {heterogeneity_table},
var_model_table {symbol_table},
trend_component_model_table {symbol_table},
var_expectation_model_table {symbol_table},
pac_model_table {symbol_table},
expressions_tree {symbol_table, num_constants, external_functions_table},
original_model {symbol_table, num_constants, external_functions_table,
trend_component_model_table, var_model_table},
dynamic_model {symbol_table, num_constants, external_functions_table,
trend_component_model_table, var_model_table},
trend_dynamic_model {symbol_table, num_constants, external_functions_table,
trend_component_model_table, var_model_table},
orig_ramsey_dynamic_model {symbol_table, num_constants, external_functions_table,
trend_component_model_table, var_model_table},
epilogue {symbol_table, num_constants, external_functions_table, trend_component_model_table,
expressions_tree {symbol_table, num_constants, external_functions_table, heterogeneity_table},
original_model {symbol_table,
num_constants,
external_functions_table,
heterogeneity_table,
trend_component_model_table,
var_model_table},
static_model {symbol_table, num_constants, external_functions_table},
steady_state_model {symbol_table, num_constants, external_functions_table, static_model},
dynamic_model {symbol_table,
num_constants,
external_functions_table,
heterogeneity_table,
trend_component_model_table,
var_model_table},
trend_dynamic_model {symbol_table,
num_constants,
external_functions_table,
heterogeneity_table,
trend_component_model_table,
var_model_table},
orig_ramsey_dynamic_model {symbol_table,
num_constants,
external_functions_table,
heterogeneity_table,
trend_component_model_table,
var_model_table},
epilogue {symbol_table,
num_constants,
external_functions_table,
heterogeneity_table,
trend_component_model_table,
var_model_table},
static_model {symbol_table, num_constants, external_functions_table, heterogeneity_table},
steady_state_model {symbol_table, num_constants, external_functions_table, heterogeneity_table,
static_model},
warnings {warnings_arg}
{
heterogeneity_table.setSymbolTable(&symbol_table);
}
void
......@@ -167,7 +190,7 @@ ModFile::checkPass(bool nostrict, bool stochastic)
exit(EXIT_FAILURE);
}
if (mod_file_struct.ramsey_constraints_present && !mod_file_struct.ramsey_model_present)
if (!ramsey_constraints.empty() && !mod_file_struct.ramsey_model_present)
{
cerr << "ERROR: A ramsey_constraints block requires the presence of a ramsey_model or "
"ramsey_policy statement"
......@@ -355,10 +378,8 @@ ModFile::checkPass(bool nostrict, bool stochastic)
// Test if some estimated parameters are used within the values of shocks
// statements (see issue #469)
set<int> parameters_intersect;
set_intersection(mod_file_struct.parameters_within_shocks_values.begin(),
mod_file_struct.parameters_within_shocks_values.end(),
mod_file_struct.estimated_parameters.begin(),
mod_file_struct.estimated_parameters.end(),
ranges::set_intersection(mod_file_struct.parameters_within_shocks_values,
mod_file_struct.estimated_parameters,
inserter(parameters_intersect, parameters_intersect.begin()));
if (parameters_intersect.size() > 0)
{
......@@ -378,8 +399,8 @@ ModFile::checkPass(bool nostrict, bool stochastic)
// Check if some exogenous is not used in the model block, Issue #841
set<int> unusedExo0 = dynamic_model.findUnusedExogenous();
set<int> unusedExo;
set_difference(unusedExo0.begin(), unusedExo0.end(), mod_file_struct.pac_params.begin(),
mod_file_struct.pac_params.end(), inserter(unusedExo, unusedExo.begin()));
ranges::set_difference(unusedExo0, mod_file_struct.pac_params,
inserter(unusedExo, unusedExo.begin()));
if (unusedExo.size() > 0)
{
ostringstream unused_exos;
......@@ -398,6 +419,120 @@ ModFile::checkPass(bool nostrict, bool stochastic)
exit(EXIT_FAILURE);
}
}
if (!heterogeneity_table.empty())
{
if (block)
{
cerr << "ERROR: the 'block' option of the 'model' block is not supported for "
"heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.check_present)
{
cerr << "ERROR: The 'check' command is not supported for heterogeneous models" << endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.steady_present)
{
cerr << "ERROR: The 'steady' command is not supported for heterogeneous models" << endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.perfect_foresight_solver_present)
{
cerr << "ERROR: The 'perfect_foresight_solver' command is not supported for "
"heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.perfect_foresight_with_expectation_errors_solver_present)
{
cerr << "ERROR: The 'perfect_foresight_with_expectation_errors_solver' command is not "
"supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.stoch_simul_present)
{
cerr << "ERROR: The 'stoch_simul' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.estimation_present)
{
cerr << "ERROR: The 'estimation' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.osr_present)
{
cerr << "ERROR: The 'osr' command is not supported for heterogeneous models" << endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.osr_params_present)
{
cerr << "ERROR: The 'osr_params' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.optim_weights_present)
{
cerr << "ERROR: The 'optim_weights' block is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.ramsey_model_present)
{
cerr << "ERROR: The 'ramsey_model' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.discretionary_policy_present)
{
cerr << "ERROR: The 'discretionary_policy' command is not supported for heterogeneous "
"models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.planner_objective_present)
{
cerr << "ERROR: The 'planner_objective' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.extended_path_present)
{
cerr << "ERROR: The 'extended_path' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.identification_present)
{
cerr << "ERROR: The 'identification' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.sensitivity_present)
{
cerr << "ERROR: The 'sensitivity' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.mom_estimation_present)
{
cerr
<< "ERROR: The 'methods_of_moments' command is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
if (mod_file_struct.occbin_constraints_present)
{
cerr << "ERROR: The 'occbin_constraints' block is not supported for heterogeneous models"
<< endl;
exit(EXIT_FAILURE);
}
}
}
void
......@@ -531,12 +666,22 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
*/
if (linear)
orig_ramsey_dynamic_model = dynamic_model;
DynamicModel ramsey_FOC_equations_dynamic_model {
symbol_table, num_constants, external_functions_table, trend_component_model_table,
DynamicModel ramsey_FOC_equations_dynamic_model {symbol_table,
num_constants,
external_functions_table,
heterogeneity_table,
trend_component_model_table,
var_model_table};
ramsey_FOC_equations_dynamic_model = dynamic_model;
auto clone_if_not_null
= [&](expr_t e) { return e ? e->clone(ramsey_FOC_equations_dynamic_model) : nullptr; };
map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints;
for (const auto& [symb_id, bounds] : ramsey_constraints)
cloned_ramsey_constraints.try_emplace(symb_id, clone_if_not_null(bounds.first),
clone_if_not_null(bounds.second));
mod_file_struct.ramsey_orig_endo_nbr
= ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective);
= ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective,
cloned_ramsey_constraints);
ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
}
......@@ -545,6 +690,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
// Must come after detrending of variables and Ramsey policy transformation
dynamic_model.substituteLogTransform();
if (!heterogeneity_table.empty())
dynamic_model.substituteAggregationOperators();
/* Create auxiliary vars for leads and lags greater than 2, on both endos and
exos. The transformation is not exactly the same on stochastic and
deterministic models, because there is no need to take into account the
......@@ -583,6 +731,10 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
dynamic_model.reorderAuxiliaryEquations();
symbol_table.resizeHetAuxVars();
for (auto& hm : heterogeneous_models)
hm.transformPass();
// Freeze the symbol table
symbol_table.freeze();
......@@ -697,6 +849,16 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
exit(EXIT_FAILURE);
}
}
for (int dim {0}; dim < heterogeneity_table.size(); dim++)
if (heterogeneous_models.at(dim).equation_number() != symbol_table.het_endo_nbr(dim))
{
cerr << "ERROR: There are " << heterogeneous_models.at(dim).equation_number()
<< " equations but " << symbol_table.het_endo_nbr(dim)
<< " endogenous variables in the model for heterogeneity dimension '"
<< heterogeneity_table.getName(dim) << "'!" << endl;
exit(EXIT_FAILURE);
}
}
void
......@@ -802,8 +964,22 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
no_tmp_terms, block, use_dll);
}
}
else // No computing task requested, compute derivatives up to 2nd order by default
else // No computing task requested, compute derivatives up to 2nd order by default unless
// output=first (preprocessor#100) or third (preprocessor#121) is requested
{
switch (output)
{
case OutputType::first:
dynamic_model.computingPass(1, 0, global_eval_context, no_tmp_terms, block, use_dll);
break;
case OutputType::third:
dynamic_model.computingPass(3, 0, global_eval_context, no_tmp_terms, block, use_dll);
break;
default:
dynamic_model.computingPass(2, 0, global_eval_context, no_tmp_terms, block, use_dll);
break;
}
}
if (linear)
{
......@@ -817,6 +993,9 @@ ModFile::computingPass(bool no_tmp_terms, OutputType output, int params_derivs_o
// Those matrices can only be filled here, because we use derivatives
dynamic_model.fillVarModelTableMatrices();
for (auto& hm : heterogeneous_models)
hm.computingPass(mod_file_struct.order_option, no_tmp_terms, use_dll);
for (auto& statement : statements)
statement->computingPass(mod_file_struct);
......@@ -859,6 +1038,11 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
auto plusfolder {DataTree::packageDir(basename)};
if (check_model_changes && !heterogeneity_table.empty())
{
cerr << "ERROR: the 'fast' option is not supported for heterogeneous models" << endl;
exit(EXIT_FAILURE);
}
bool hasModelChanged = !dynamic_model.isChecksumMatching(basename) || !check_model_changes;
if (hasModelChanged)
{
......@@ -899,8 +1083,8 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
mOutputFile << "clearvars -global" << endl
<< "clear_persistent_variables(fileparts(which('dynare')), false)" << endl;
else if (clear_global)
mOutputFile
<< "clear M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info estimation_info;"
mOutputFile << "clearvars -global M_ options_ oo_ estim_params_ bayestopt_ dataset_ "
"dataset_info estimation_info;"
<< endl;
if (!notime)
......@@ -934,6 +1118,7 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
mOutputFile << "M_.parameter_used_with_lead_lag = true;" << endl;
symbol_table.writeOutput(mOutputFile);
heterogeneity_table.writeOutput(mOutputFile);
var_model_table.writeOutput(basename, mOutputFile);
trend_component_model_table.writeOutput(basename, mOutputFile);
......@@ -945,6 +1130,10 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
<< ");" << endl
<< "M_.Correlation_matrix = eye(" << symbol_table.exo_nbr() << ", "
<< symbol_table.exo_nbr() << ");" << endl;
for (int hd {0}; hd < heterogeneity_table.size(); hd++)
mOutputFile << "M_.heterogeneity(" << hd + 1 << ").Sigma_e = zeros("
<< symbol_table.het_exo_nbr(hd) << ", " << symbol_table.het_exo_nbr(hd) << ");"
<< endl;
if (mod_file_struct.calibrated_measurement_errors)
mOutputFile << "M_.H = zeros(" << symbol_table.observedVariablesNbr() << ", "
......@@ -966,7 +1155,8 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
<< "M_.heteroskedastic_shocks.Qvalue_orig = [];" << endl
<< "M_.heteroskedastic_shocks.Qscale_orig = [];" << endl
<< "M_.matched_irfs = {};" << endl
<< "M_.matched_irfs_weights = {};" << endl;
<< "M_.matched_irfs_weights = {};" << endl
<< "M_.perfect_foresight_controlled_paths = [];" << endl;
// NB: options_.{ramsey,discretionary}_policy should rather be fields of M_
mOutputFile << boolalpha << "options_.linear = " << linear << ";" << endl
......@@ -1033,6 +1223,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
}
}
for (const auto& hm : heterogeneous_models)
hm.writeDriverOutput(mOutputFile);
if (onlymodel || gui)
for (const auto& statement : statements)
{
......@@ -1137,8 +1330,8 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
if (!no_warn)
{
if (warnings.countWarnings() > 0)
mOutputFile << "disp('Note: " << warnings.countWarnings()
if (int num_warnings {warnings.numWarnings()}; num_warnings > 0)
mOutputFile << "disp('Note: " << num_warnings
<< " warning(s) encountered in the preprocessor')" << endl;
mOutputFile << "if ~isempty(lastwarn)" << endl
......@@ -1181,6 +1374,9 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
epilogue.writeEpilogueFile(basename);
pac_model_table.writeTargetCoefficientsFile(basename);
for (const auto& hm : heterogeneous_models)
hm.writeModelFiles(basename, false);
}
}
......@@ -1258,6 +1454,11 @@ ModFile::writeJsonOutputParsingCheck(const string& basename, JsonFileOutputType
symbol_table.writeJsonOutput(output);
output << ", ";
if (!heterogeneity_table.empty())
{
heterogeneity_table.writeJsonOutput(output);
output << ", ";
}
dynamic_model.writeJsonOutput(output);
output << ", ";
static_model.writeJsonOutput(output);
......
src/ModFile.hh
View file @ 6618d7d1
/*
* Copyright © 2006-2023 Dynare Team
* Copyright © 2006-2024 Dynare Team
*
* This file is part of Dynare.
*
......@@ -30,6 +30,8 @@
#include "DynamicModel.hh"
#include "ExtendedPreprocessorTypes.hh"
#include "ExternalFunctionsTable.hh"
#include "HeterogeneityTable.hh"
#include "HeterogeneousModel.hh"
#include "ModelEquationBlock.hh"
#include "NumericalConstants.hh"
#include "NumericalInitialization.hh"
......@@ -46,6 +48,8 @@ class ModFile
{
public:
explicit ModFile(WarningConsolidation& warnings_arg);
// For heterogeneity dimensions
HeterogeneityTable heterogeneity_table;
//! Symbol table
SymbolTable symbol_table;
//! External Functions table
......@@ -76,6 +80,8 @@ public:
StaticModel static_model;
//! Static model, as declared in the "steady_state_model" block if present
SteadyStateModel steady_state_model;
// Heterogeneous model blocks, ordered per heterogeneity dimension ID
vector<HeterogeneousModel> heterogeneous_models;
//! Option linear
bool linear {false};
......@@ -117,6 +123,11 @@ public:
/*! (i.e. option parallel_local_files of model block) */
vector<string> parallel_local_files;
/* Contents of the ramsey_constraints block.
Maps symb_id → (lower_bound, upper_bound).
NB: The two expr_t live in dynamic_model */
map<int, pair<expr_t, expr_t>> ramsey_constraints;
private:
//! List of statements
vector<unique_ptr<Statement>> statements;
......
src/ModelEquationBlock.cc
View file @ 6618d7d1
......@@ -26,8 +26,10 @@
PlannerObjective::PlannerObjective(SymbolTable& symbol_table_arg,
NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg) :
StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg}
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg) :
StaticModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
heterogeneity_table_arg}
{
}
......@@ -38,7 +40,7 @@ PlannerObjective::writeDriverOutput(ostream& output) const
for (const auto& it : temporary_terms_derivatives)
output << it.size() << "; ";
output << "];" << endl;
writeDriverSparseIndicesHelper<false, true>(output);
writeDriverSparseIndicesHelper("objective", output);
}
void
......@@ -51,9 +53,14 @@ PlannerObjective::computingPassBlock([[maybe_unused]] const eval_context_t& eval
OrigRamseyDynamicModel::OrigRamseyDynamicModel(
SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg, VarModelTable& var_model_table_arg) :
DynamicModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
trend_component_model_table_arg, var_model_table_arg}
DynamicModel {symbol_table_arg,
num_constants_arg,
external_functions_table_arg,
heterogeneity_table_arg,
trend_component_model_table_arg,
var_model_table_arg}
{
}
......@@ -67,8 +74,10 @@ OrigRamseyDynamicModel::operator=(const DynamicModel& m)
SteadyStateModel::SteadyStateModel(SymbolTable& symbol_table_arg,
NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
const StaticModel& static_model_arg) :
DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg},
DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
heterogeneity_table_arg},
static_model {static_model_arg}
{
}
......@@ -327,10 +336,15 @@ SteadyStateModel::writeJsonSteadyStateFile(ostream& output, bool transformComput
Epilogue::Epilogue(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg) :
DynamicModel {symbol_table_arg, num_constants_arg, external_functions_table_arg,
trend_component_model_table_arg, var_model_table_arg}
DynamicModel {symbol_table_arg,
num_constants_arg,
external_functions_table_arg,
heterogeneity_table_arg,
trend_component_model_table_arg,
var_model_table_arg}
{
}
......@@ -377,7 +391,7 @@ Epilogue::checkPass(ModFileStructure& mod_file_struct) const
for (const auto& [symb_id, expr] : dynamic_def_table)
if (so_far_defined.contains(symb_id))
{
cerr << "WARNING: in the 'epilogue' block, variable '" << symbol_table.getName(symb_id)
cerr << "ERROR: in the 'epilogue' block, variable '" << symbol_table.getName(symb_id)
<< "' is declared twice" << endl;
exit(EXIT_FAILURE);
}
......@@ -396,10 +410,10 @@ void
Epilogue::detrend(const map<int, expr_t>& trend_symbols_map,
const nonstationary_symbols_map_t& nonstationary_symbols_map)
{
for (const auto& it : ranges::reverse_view(nonstationary_symbols_map))
for (const auto& [symb_id, deflator] : ranges::reverse_view(nonstationary_symbols_map))
for (auto& [symb_id, expr] : dynamic_def_table)
{
expr = expr->detrend(it.first, it.second.first, it.second.second);
expr = expr->detrend(symb_id, deflator.first, deflator.second);
assert(expr);
}
......
src/ModelEquationBlock.hh
View file @ 6618d7d1
/*
* Copyright © 2010-2024 Dynare Team
* Copyright © 2010-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -30,7 +30,8 @@ class PlannerObjective : public StaticModel
{
public:
PlannerObjective(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg);
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg);
// NB: masks the method with the same name in StaticModel (not in a virtual fashion)
void writeDriverOutput(ostream& output) const;
......@@ -50,6 +51,7 @@ class OrigRamseyDynamicModel : public DynamicModel
public:
OrigRamseyDynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg);
OrigRamseyDynamicModel& operator=(const DynamicModel& m);
......@@ -60,6 +62,17 @@ protected:
{
return "original Ramsey model";
}
int
getJacobianCol(int deriv_id, [[maybe_unused]] bool sparse) const override
{
/* Override the DynamicModel method by returning a dummy Jacobian column number.
The override is necessary because the method from DynamicModel fails with
endos with lag/lead greater than 1 or exos with a lag/lead, while substitutions
are by definition not done for an original model.
In particular, this fixes dynare#1960 (equation derivatives are computed for models declared
as linear, to check whether they are truly linear). */
return deriv_id;
}
};
class SteadyStateModel : public DataTree
......@@ -75,6 +88,7 @@ private:
public:
SteadyStateModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
const StaticModel& static_model_arg);
SteadyStateModel(const SteadyStateModel& m);
......@@ -108,6 +122,7 @@ private:
public:
Epilogue(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg,
TrendComponentModelTable& trend_component_model_table_arg,
VarModelTable& var_model_table_arg);
......
src/ModelTree.cc
View file @ 6618d7d1
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -34,6 +34,7 @@
#endif
#include <algorithm>
#include <numeric>
#include <regex>
#include <utility>
......@@ -51,13 +52,21 @@ vector<jthread> ModelTree::mex_compilation_workers;
void
ModelTree::copyHelper(const ModelTree& m)
{
auto f = [this](expr_t e) { return e->clone(*this); };
auto f = [this](expr_t e) { return e ? e->clone(*this) : nullptr; };
// Equations
for (const auto& it : m.equations)
equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
for (const auto& it : m.aux_equations)
aux_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
for (const auto& it : m.complementarity_conditions)
if (it)
{
const auto& [symb_id, lb, ub] = *it;
complementarity_conditions.emplace_back(in_place, symb_id, f(lb), f(ub));
}
else
complementarity_conditions.emplace_back(nullopt);
auto convert_deriv_map = [f](const map<vector<int>, expr_t>& dm) {
map<vector<int>, expr_t> dm2;
......@@ -135,8 +144,10 @@ ModelTree::copyHelper(const ModelTree& m)
}
ModelTree::ModelTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg},
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
heterogeneity_table_arg, is_dynamic_arg},
derivatives(4),
NNZDerivatives(4, 0),
temporary_terms_derivatives(4)
......@@ -184,6 +195,8 @@ ModelTree::operator=(const ModelTree& m)
equations_lineno = m.equations_lineno;
aux_equations.clear();
equation_tags = m.equation_tags;
complementarity_conditions.clear();
computed_derivs_order = m.computed_derivs_order;
NNZDerivatives = m.NNZDerivatives;
......@@ -264,7 +277,7 @@ ModelTree::computeNormalization(const jacob_map_t& contemporaneous_jacobian)
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to
// them
endo2eq.resize(equations.size());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(),
ranges::transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(),
[=](vertex_descriptor_t i) { return i - n; });
}
......@@ -882,7 +895,7 @@ ModelTree::determineLinearBlocks()
int
ModelTree::equation_number() const
{
return (equations.size());
return equations.size();
}
void
......@@ -972,7 +985,7 @@ ModelTree::computeTemporaryTerms(bool is_matlab, bool no_tmp_terms)
temporary_terms_derivatives.clear();
temporary_terms_derivatives.resize(derivatives.size());
for (int order = 0; order < static_cast<int>(derivatives.size()); order++)
copy(temp_terms_map[{order, 0}].begin(), temp_terms_map[{order, 0}].end(),
ranges::copy(temp_terms_map[{order, 0}],
inserter(temporary_terms_derivatives.at(order),
temporary_terms_derivatives.at(order).begin()));
......@@ -1027,7 +1040,7 @@ ModelTree::computeBlockTemporaryTerms(bool no_tmp_terms)
{
blocks_temporary_terms.at(blk).resize(temp_terms.at(blk).size());
for (size_t i {0}; i < temp_terms.at(blk).size(); i++)
copy(temp_terms.at(blk).at(i).begin(), temp_terms.at(blk).at(i).end(),
ranges::copy(temp_terms.at(blk).at(i),
inserter(blocks_temporary_terms.at(blk).at(i),
blocks_temporary_terms.at(blk).at(i).begin()));
}
......@@ -1385,8 +1398,7 @@ ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_b
content_output << endl << R"(\end{dmath})" << endl;
}
output << R"(\include{)" << latex_basename + "_content"
<< "}" << endl
output << R"(\include{)" << latex_basename + "_content" << "}" << endl
<< R"(\end{document})" << endl;
output.close();
......@@ -1394,28 +1406,33 @@ ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_b
}
void
ModelTree::addEquation(expr_t eq, const optional<int>& lineno)
ModelTree::addEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition)
{
auto beq = dynamic_cast<BinaryOpNode*>(eq);
assert(beq && beq->op_code == BinaryOpcode::equal);
equations.push_back(beq);
equations_lineno.push_back(lineno);
complementarity_conditions.push_back(move(complementarity_condition));
}
void
ModelTree::findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const
ModelTree::findConstantEquationsWithoutComplementarityCondition(
map<VariableNode*, NumConstNode*>& subst_table) const
{
for (size_t i = 0; i < equations.size(); i++)
if (!equation_tags.exists(i, "mcp"))
if (!complementarity_conditions[i])
equations[i]->findConstantEquations(subst_table);
}
void
ModelTree::addEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags)
ModelTree::addEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition,
map<string, string> eq_tags)
{
equation_tags.add(equations.size(), move(eq_tags));
addEquation(eq, lineno);
addEquation(eq, lineno, move(complementarity_condition));
}
void
......@@ -1528,8 +1545,7 @@ ModelTree::computeParamsDerivativesTemporaryTerms()
d->computeTemporaryTerms(order, temp_terms_map, reference_count, true);
for (const auto& [order, tts] : temp_terms_map)
copy(temp_terms_map[order].begin(), temp_terms_map[order].end(),
inserter(params_derivs_temporary_terms[order],
ranges::copy(temp_terms_map[order], inserter(params_derivs_temporary_terms[order],
params_derivs_temporary_terms[order].begin()));
for (int idx {0}; const auto& [order, tts] : params_derivs_temporary_terms)
......@@ -1593,6 +1609,27 @@ ModelTree::writeJsonModelEquations(ostream& output, bool residuals) const
eqtags.clear();
}
}
if (complementarity_conditions[eq])
{
auto& [symb_id, lower_bound, upper_bound] = *complementarity_conditions[eq];
output << R"(, "complementarity_condition": {"variable": ")"
<< symbol_table.getName(symb_id) << '"';
if (lower_bound)
{
output << R"(, "lower_bound": ")";
lower_bound->writeJsonOutput(output, {}, {});
output << '"';
}
if (upper_bound)
{
output << R"(, "upper_bound": ")";
upper_bound->writeJsonOutput(output, {}, {});
output << '"';
}
output << "}";
}
output << "}" << endl;
}
output << endl << "]" << endl;
......@@ -1637,7 +1674,7 @@ ModelTree::findCompilerOnMacos(const string& mexext)
Apple’s clang is located both in /usr/bin/gcc and /usr/bin/clang, it
automatically selects x86_64 or arm64 depending on the compile-time
environment. */
const string macos_gcc_version {"13"};
const string macos_gcc_version {"15"};
if (filesystem::path global_gcc_path {"/usr/local/bin/gcc-" + macos_gcc_version};
exists(global_gcc_path) && mexext == "mexmaci64")
......@@ -1649,14 +1686,14 @@ ModelTree::findCompilerOnMacos(const string& mexext)
return {global_clang_path, true};
else
{
cerr << "ERROR: You must install gcc-" << macos_gcc_version
cerr << "ERROR: You must install gcc@" << macos_gcc_version
<< " on your system before using the `use_dll` option of Dynare. "
<< "You should install Homebrew";
if (mexext == "mexmaca64")
cerr << " for arm64";
else if (mexext == "mexmaci64")
cerr << " for x86_64";
cerr << " and run `brew install gcc-" << macos_gcc_version << "` in a terminal." << endl;
cerr << " and run `brew install gcc@" << macos_gcc_version << "` in a terminal." << endl;
exit(EXIT_FAILURE);
}
}
......@@ -1783,7 +1820,7 @@ ModelTree::compileMEX(const filesystem::path& output_dir, const string& output_b
// The prerequisites are the object files among the input files
set<filesystem::path> prerequisites;
copy_if(input_files.begin(), input_files.end(), inserter(prerequisites, prerequisites.end()),
ranges::copy_if(input_files, inserter(prerequisites, prerequisites.end()),
[](const auto& p) { return p.extension() == ".o"; });
unique_lock<mutex> lk {mex_compilation_mut};
......@@ -1944,8 +1981,7 @@ ModelTree::initializeMEXCompilationWorkers(int numworkers, const filesystem::pat
auto pick_job = [&cmd, &output] {
for (auto it {mex_compilation_queue.begin()}; it != mex_compilation_queue.end(); ++it)
if (const auto& prerequisites {get<1>(*it)}; // Will become dangling after erase
includes(mex_compilation_done.begin(), mex_compilation_done.end(),
prerequisites.begin(), prerequisites.end()))
ranges::includes(mex_compilation_done, prerequisites))
{
output = get<0>(*it);
cmd = get<2>(*it);
......@@ -2037,6 +2073,8 @@ ModelTree::waitForMEXCompilationWorkers()
void
ModelTree::computingPassBlock(const eval_context_t& eval_context, bool no_tmp_terms)
{
if (!heterogeneity_table.empty())
return;
if (!computeNonSingularNormalization(eval_context))
return;
auto [prologue, epilogue] = computePrologueAndEpilogue();
......@@ -2078,3 +2116,70 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o
output << ";" << endl;
}
}
void
ModelTree::computeMCPEquationsReordering(const optional<int>& heterogeneous_dimension)
{
/* Optimal policy models (discretionary, or Ramsey before computing FOCs) do not have as many
equations as variables. Do not even try to compute the reordering. */
if (static_cast<int>(equations.size())
!= (heterogeneous_dimension ? symbol_table.het_endo_nbr(*heterogeneous_dimension)
: symbol_table.endo_nbr()))
return;
assert(equations.size() == complementarity_conditions.size());
mcp_equations_reordering.resize(equations.size());
iota(mcp_equations_reordering.begin(), mcp_equations_reordering.end(), 0);
set<int> endos;
for (int eq {0}; eq < static_cast<int>(equations.size()); eq++)
if (complementarity_conditions.at(eq))
{
int symb_id {get<0>(*complementarity_conditions[eq])};
auto [ignore, inserted] = endos.insert(symb_id);
if (!inserted)
{
cerr << "ERROR: variable " << symbol_table.getName(symb_id)
<< " appears in two complementarity conditions" << endl;
exit(EXIT_FAILURE);
}
int endo_id {symbol_table.getTypeSpecificID(symb_id)};
auto it = ranges::find(mcp_equations_reordering, eq);
assert(it != mcp_equations_reordering.end());
swap(mcp_equations_reordering[endo_id], *it);
}
}
void
ModelTree::writeDriverSparseIndicesHelper(const string& prefix, ostream& output) const
{
// Write indices for the sparse Jacobian (both naive and CSC storage)
output << "M_." << prefix << "_g1_sparse_rowval = int32([";
for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
output << indices.first + 1 << ' ';
output << "]);" << endl << "M_." << prefix << "_g1_sparse_colval = int32([";
for (const auto& [indices, d1] : jacobian_sparse_column_major_order)
output << indices.second + 1 << ' ';
output << "]);" << endl << "M_." << prefix << "_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_." << prefix << "_g" << i << "_sparse_indices = int32([";
for (const auto& [vidx, d] : derivatives[i])
{
for (bool row_number {true}; // First element of vidx is row number
int it : vidx)
output << (exchange(row_number, false) ? it : getJacobianCol(it, true)) + 1 << ' ';
output << ';' << endl;
}
output << "]);" << endl;
}
}
src/ModelTree.hh
View file @ 6618d7d1
/*
* Copyright © 2003-2024 Dynare Team
* Copyright © 2003-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -91,6 +91,9 @@ protected:
vector<optional<int>> equations_lineno;
//! Stores equation tags
EquationTags equation_tags;
/* The tuple is: endogenous symbol ID, lower bound (possibly nullptr), upper bound (possibly
nullptr). */
vector<optional<tuple<int, expr_t, expr_t>>> complementarity_conditions;
/*
* ************** END **************
*/
......@@ -208,7 +211,7 @@ protected:
getRecursiveSize() const
{
return size - mfs_size;
};
}
};
// Whether block decomposition has been successfully computed
......@@ -271,6 +274,12 @@ protected:
Same remark as above regarding blocks of type “evaluate”. */
vector<vector<int>> blocks_jacobian_sparse_colptr;
/* Indices of reordered equations for use with an MCP solver. Contains a permutation of the
equation indices, so that reordered equations appear at the (type-specific) index of the
endogenous to which they are associated through the complementarity condition.
Also see computeMCPEquationsReordering() method. */
vector<int> mcp_equations_reordering;
//! Computes derivatives
/*! \param order the derivation order
\param vars the derivation IDs w.r.t. which compute the derivatives */
......@@ -284,6 +293,11 @@ protected:
//! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms();
/* Computes the mcp_equations_reordering vector.
Also checks that a variable does not appear as constrained in two different equations. */
void computeMCPEquationsReordering(const optional<int>& heterogeneous_dimension = nullopt);
//! Writes temporary terms
template<ExprNodeOutputType output_type>
void writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union,
......@@ -352,9 +366,8 @@ protected:
temporary_terms_t& temporary_terms_union) const;
/* Helper for writing sparse derivatives indices in MATLAB/Octave driver file.
Also supports the planner objective through the corresponding boolean. */
template<bool dynamic, bool objective>
void writeDriverSparseIndicesHelper(ostream& output) const;
The “prefix” will be prepended to variable names to construct objects under M_. */
void writeDriverSparseIndicesHelper(const string& prefix, ostream& output) const;
// Helper for writing sparse derivatives indices in JSON
template<bool dynamic>
......@@ -395,7 +408,8 @@ protected:
// Writes the sparse representation of the model in MATLAB/Octave
template<bool dynamic>
void writeSparseModelMFiles(const string& basename) const;
void writeSparseModelMFiles(const string& basename,
const optional<int>& heterogeneous_dimension = nullopt) const;
// Writes and compiles the sparse representation of the model in C
template<bool dynamic>
......@@ -424,6 +438,11 @@ protected:
template<bool dynamic>
void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
template<bool dynamic>
void writeComplementarityConditionsFile(const string& basename,
const optional<int>& heterogeneous_dimension
= nullopt) const;
private:
//! Sparse matrix of double to store the values of the static Jacobian
/*! First index is equation number, second index is endogenous type specific ID */
......@@ -542,7 +561,7 @@ protected:
{
return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
.first;
};
}
//! Return true if the equation has been normalized
bool
isBlockEquationRenormalized(int blk, int eq) const
......@@ -550,44 +569,44 @@ protected:
return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
.first
== EquationType::evaluateRenormalized;
};
}
//! Return the expr_t of equation belonging to the block
BinaryOpNode*
getBlockEquationExpr(int blk, int eq) const
{
return equations[eq_idx_block2orig[blocks[blk].first_equation + eq]];
};
}
//! Return the expr_t of renormalized equation belonging to the block
BinaryOpNode*
getBlockEquationRenormalizedExpr(int blk, int eq) const
{
return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
.second;
};
}
//! Return the original number of equation belonging to the block
int
getBlockEquationID(int blk, int eq) const
{
return eq_idx_block2orig[blocks[blk].first_equation + eq];
};
}
//! Return the original number of variable belonging to the block
int
getBlockVariableID(int blk, int var) const
{
return endo_idx_block2orig[blocks[blk].first_equation + var];
};
}
//! Return the position of an equation (given by its original index) inside its block
int
getBlockInitialEquationID(int blk, int eq) const
{
return eq_idx_orig2block[eq] - blocks[blk].first_equation;
};
}
//! Return the position of a variable (given by its original index) inside its block
int
getBlockInitialVariableID(int blk, int var) const
{
return endo_idx_orig2block[var] - blocks[blk].first_equation;
};
}
//! Initialize equation_reordered & variable_reordered
void initializeVariablesAndEquations();
......@@ -642,7 +661,8 @@ private:
public:
ModelTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg = false);
ExternalFunctionsTable& external_functions_table_arg,
HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);
protected:
ModelTree(const ModelTree& m);
......@@ -651,10 +671,14 @@ protected:
public:
//! Absolute value under which a number is considered to be zero
double cutoff {1e-15};
//! Declare a node as an equation of the model; also give its line number
void addEquation(expr_t eq, const optional<int>& lineno);
/* Declare a node as an equation of the model; also give its line number and complementarity
condition */
void addEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition = nullopt);
//! Declare a node as an equation of the model, also giving its tags
void addEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags);
void addEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition,
map<string, string> eq_tags);
//! Declare a node as an auxiliary equation of the model, adding it at the end of the list of
//! auxiliary equations
void addAuxEquation(expr_t eq);
......@@ -671,9 +695,10 @@ public:
/*! Reorder auxiliary variables so that they appear in recursive order in
set_auxiliary_variables.m and dynamic_set_auxiliary_series.m */
void reorderAuxiliaryEquations();
//! Find equations of the form “variable=constant”, excluding equations with “mcp” tag (see
//! dynare#1697)
void findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const;
/* Find equations of the form “variable=constant”, excluding equations with a complementarity
condition (see dynare#1697) */
void findConstantEquationsWithoutComplementarityCondition(
map<VariableNode*, NumConstNode*>& subst_table) const;
/* Given an expression, searches for the first equation that has exactly this
expression on the LHS, and returns the RHS of that equation.
If no such equation can be found, throws an ExprNode::MatchFailureExpression */
......@@ -754,7 +779,7 @@ ModelTree::writeModelEquations(ostream& output, const temporary_terms_t& tempora
BinaryOpNode* eq_node {equations[eq]};
expr_t lhs {eq_node->arg1}, rhs {eq_node->arg2};
// Test if the right hand side of the equation is empty.
// Test if the right-hand side of the equation is empty.
double vrhs {1.0};
try
{
......@@ -764,7 +789,7 @@ ModelTree::writeModelEquations(ostream& output, const temporary_terms_t& tempora
{
}
if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
if (vrhs != 0) // The right-hand side of the equation is not empty ==> residual=lhs-rhs;
{
output << " residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
......@@ -774,7 +799,7 @@ ModelTree::writeModelEquations(ostream& output, const temporary_terms_t& tempora
rhs->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs);
output << ");" << endl;
}
else // The right hand side of the equation is empty ==> residual=lhs;
else // The right-hand side of the equation is empty ==> residual=lhs;
{
output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
<< eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type)
......@@ -981,12 +1006,10 @@ ModelTree::writeModelCFile(const string& basename, const string& mexext,
const filesystem::path model_src_dir {filesystem::path {basename} / "model" / "src"};
auto [d_output, tt_output] = writeModelFileHelper < dynamic
? ExprNodeOutputType::CDynamicModel
auto [d_output, tt_output] = writeModelFileHelper<dynamic ? ExprNodeOutputType::CDynamicModel
: ExprNodeOutputType::CStaticModel>();
vector<filesystem::path> header_files, object_files;
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "dynamic_" : "static_"};
const string ss_it_argin {dynamic ? ", const double *restrict steady_state, int it_" : ""};
const string ss_it_argout {dynamic ? ", steady_state, it_" : ""};
......@@ -1535,7 +1558,7 @@ ModelTree::writeBytecodeModelEquations(Bytecode::Writer& code_file,
BinaryOpNode* eq_node {equations[eq]};
expr_t lhs {eq_node->arg1}, rhs {eq_node->arg2};
code_file << Bytecode::FNUMEXPR {Bytecode::ExpressionType::ModelEquation, eq};
// Test if the right hand side of the equation is empty.
// Test if the right-hand side of the equation is empty.
double vrhs {1.0};
try
{
......@@ -1545,7 +1568,7 @@ ModelTree::writeBytecodeModelEquations(Bytecode::Writer& code_file,
{
}
if (vrhs != 0) // The right hand side of the equation is not empty ⇒ residual=lhs-rhs
if (vrhs != 0) // The right-hand side of the equation is not empty ⇒ residual=lhs-rhs
{
lhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs,
tef_terms);
......@@ -1554,7 +1577,7 @@ ModelTree::writeBytecodeModelEquations(Bytecode::Writer& code_file,
code_file << Bytecode::FBINARY {BinaryOpcode::minus} << Bytecode::FSTPR {eq};
}
else // The right hand side of the equation is empty ⇒ residual=lhs
else // The right-hand side of the equation is empty ⇒ residual=lhs
{
lhs->writeBytecodeOutput(code_file, output_type, temporary_terms, temporary_terms_idxs,
tef_terms);
......@@ -2242,46 +2265,10 @@ ModelTree::writeJsonParamsDerivativesHelper(bool writeDetails) const
move(rpp_output), move(gpp_output), move(hp_output), move(g3p_output)};
}
template<bool dynamic, bool objective>
void
ModelTree::writeDriverSparseIndicesHelper(ostream& output) const
{
static_assert(!(objective && dynamic), "There is no such thing as a dynamic planner objective");
// TODO: when C++20 support is complete, mark this constexpr
const string model_name {objective ? "objective" : (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 (bool row_number {true}; // First element of vidx is row number
int it : vidx)
output << (exchange(row_number, false) ? it : getJacobianCol(it, true)) + 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)
......@@ -2386,12 +2373,13 @@ template<bool dynamic>
void
ModelTree::writeSparseModelJuliaFiles(const string& basename) const
{
auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper < dynamic
? ExprNodeOutputType::juliaSparseDynamicModel
assert(heterogeneity_table.empty());
auto [d_sparse_output, tt_sparse_output]
= writeModelFileHelper<dynamic ? ExprNodeOutputType::juliaSparseDynamicModel
: ExprNodeOutputType::juliaSparseStaticModel>();
filesystem::path julia_dir {filesystem::path {basename} / "model" / "julia"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "SparseDynamic" : "SparseStatic"};
const string ss_argin {dynamic ? ", steady_state::Vector{<: Real}" : ""};
const string ss_argout {dynamic ? ", steady_state" : ""};
......@@ -2496,17 +2484,25 @@ ModelTree::writeSparseModelJuliaFiles(const string& basename) const
template<bool dynamic>
void
ModelTree::writeSparseModelMFiles(const string& basename) const
ModelTree::writeSparseModelMFiles(const string& basename,
const optional<int>& heterogeneous_dimension) 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"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "dynamic_" : "static_"};
const string prefix {
(dynamic ? "dynamic_"s : "static_"s)
+ (heterogeneous_dimension ? "het"s + to_string(*heterogeneous_dimension + 1) + "_"s : ""s)};
const string full_prefix {basename + ".sparse." + prefix};
const string ss_arg {dynamic ? ", steady_state" : ""};
const string extra_args {(dynamic ? ", steady_state"s : ""s)
+ (heterogeneous_dimension
? ", yh, xh, paramsh"s
: (heterogeneity_table.empty() ? ""s : ", yagg"s))};
const int nextra_args {
static_cast<int>(dynamic)
+ (heterogeneous_dimension ? 3 : static_cast<int>(!heterogeneity_table.empty()))};
size_t ttlen {temporary_terms_derivatives[0].size()};
......@@ -2522,7 +2518,7 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
// Residuals (non-block)
open_file(m_dir / (prefix + "resid_tt.m"));
output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << ss_arg
output << "function [T_order, T] = " << prefix << "resid_tt(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "if T_order >= 0" << endl
<< " return" << endl
......@@ -2535,13 +2531,13 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
output.close();
open_file(m_dir / (prefix + "resid.m"));
output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << ss_arg
output << "function [residual, T_order, T] = " << prefix << "resid(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "if nargin < " << 5 + static_cast<int>(dynamic) << endl
<< "if nargin < " << 5 + nextra_args << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "end" << endl
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "residual = NaN(" << equations.size() << ", 1);" << endl
<< d_sparse_output[0].str();
......@@ -2552,12 +2548,12 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
ttlen += temporary_terms_derivatives[1].size();
open_file(m_dir / (prefix + "g1_tt.m"));
output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << ss_arg
output << "function [T_order, T] = " << prefix << "g1_tt(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "if T_order >= 1" << endl
<< " return" << endl
<< "end" << endl
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << ss_arg
<< "[T_order, T] = " << full_prefix << "resid_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "T_order = 1;" << endl
<< "if size(T, 1) < " << ttlen << endl
......@@ -2568,22 +2564,17 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
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
output << "function [g1, T_order, T] = " << prefix << "g1(y, x, params" << extra_args
<< ", sparse_rowval, sparse_colval, sparse_colptr, T_order, T)" << endl
<< "if nargin < " << 8 + static_cast<int>(dynamic) << endl
<< "if nargin < " << 8 + nextra_args << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "end" << endl
<< "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << ss_arg << ", T_order, T);"
<< endl
<< "[T_order, T] = " << full_prefix << "g1_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "g1_v = NaN(" << jacobian_sparse_column_major_order.size() << ", 1);" << endl
<< d_sparse_output[1].str();
// On MATLAB < R2020a, sparse() does not accept int32 indices
output << "if ~isoctave && matlab_ver_less_than('9.8')" << endl
<< " sparse_rowval = double(sparse_rowval);" << endl
<< " sparse_colval = double(sparse_colval);" << endl
<< "end" << endl
<< "g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << equations.size() << ", "
output << "g1 = sparse(sparse_rowval, sparse_colval, g1_v, " << equations.size() << ", "
<< getJacobianColsNbr(true) << ");" << endl
<< "end" << endl;
output.close();
......@@ -2594,12 +2585,12 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
ttlen += temporary_terms_derivatives[i].size();
open_file(m_dir / (prefix + "g" + to_string(i) + "_tt.m"));
output << "function [T_order, T] = " << prefix << "g" << i << "_tt(y, x, params" << ss_arg
output << "function [T_order, T] = " << prefix << "g" << i << "_tt(y, x, params" << extra_args
<< ", T_order, T)" << endl
<< "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] = " << full_prefix << "g" << i - 1 << "_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "T_order = " << i << ";" << endl
<< "if size(T, 1) < " << ttlen << endl
......@@ -2610,12 +2601,12 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
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
<< "if nargin < " << 5 + static_cast<int>(dynamic) << endl
<< extra_args << ", T_order, T)" << endl
<< "if nargin < " << 5 + nextra_args << endl
<< " T_order = -1;" << endl
<< " T = NaN(" << ttlen << ", 1);" << endl
<< "end" << endl
<< "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << ss_arg
<< "[T_order, T] = " << full_prefix << "g" << i << "_tt(y, x, params" << extra_args
<< ", T_order, T);" << endl
<< "g" << i << "_v = NaN(" << derivatives[i].size() << ", 1);" << endl
<< d_sparse_output[i].str() << "end" << endl;
......@@ -2636,7 +2627,7 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
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;
<< extra_args << ", sparse_rowval, sparse_colval, sparse_colptr, T)" << endl;
if (!evaluate)
output << "residual=NaN(" << blocks[blk].mfs_size << ", 1);" << endl;
......@@ -2655,12 +2646,7 @@ ModelTree::writeSparseModelMFiles(const string& basename) const
<< " g1_v = NaN(" << blocks_jacobian_sparse_column_major_order[blk].size()
<< ", 1);" << endl;
writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
// On MATLAB < R2020a, sparse() does not accept int32 indices
output << " if ~isoctave && matlab_ver_less_than('9.8')" << endl
<< " sparse_rowval = double(sparse_rowval);" << endl
<< " sparse_colval = double(sparse_colval);" << endl
<< " end" << endl
<< " g1 = sparse(sparse_rowval, sparse_colval, g1_v, "
output << " g1 = sparse(sparse_rowval, sparse_colval, g1_v, "
<< blocks[blk].mfs_size << ", "
<< (one_boundary ? 1 : 3) * blocks[blk].mfs_size << ");" << endl
<< "end" << endl;
......@@ -2682,10 +2668,14 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
const filesystem::path mex_dir {packageDir(basename) / "+sparse"};
const filesystem::path model_src_dir {filesystem::path {basename} / "model" / "src" / "sparse"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "dynamic_" : "static_"};
const string ss_argin {dynamic ? ", const double *restrict steady_state" : ""};
const string ss_argout {dynamic ? ", steady_state" : ""};
const string extra_argin {
(dynamic ? ", const double *restrict steady_state"s : ""s)
+ (heterogeneity_table.empty() ? ""s : ", const double *restrict yagg"s)};
const string extra_argout {(dynamic ? ", steady_state"s : ""s)
+ (heterogeneity_table.empty() ? ""s : ", yagg"s)};
const int nextra_args {static_cast<int>(dynamic)
+ static_cast<int>(!heterogeneity_table.empty())};
const int ylen {(dynamic ? 3 : 1) * symbol_table.endo_nbr()};
const int xlen {symbol_table.exo_nbr() + symbol_table.exo_det_nbr()};
......@@ -2703,8 +2693,8 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
size_t ttlen {0};
// Helper for dealing with y, x, params and steady_state inputs (shared with block case)
auto y_x_params_ss_inputs = [&](bool assign_y) {
// Helper for dealing with y, x, params, steady_state and yagg inputs (shared with block case)
auto y_x_params_ss_yagg_inputs = [&](bool assign_y) {
output << " if (!(mxIsDouble(prhs[0]) && !mxIsComplex(prhs[0]) && !mxIsSparse(prhs[0]) && "
"mxGetNumberOfElements(prhs[0]) == "
<< ylen << "))" << endl
......@@ -2731,12 +2721,22 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
<< R"( mexErrMsgTxt("steady_state must be a real dense numeric array with )"
<< symbol_table.endo_nbr() << R"( elements");)" << endl
<< " const double *restrict steady_state = mxGetPr(prhs[3]);" << endl;
if (!heterogeneity_table.empty())
{
const int idx {3 + static_cast<int>(dynamic)};
output << " if (!(mxIsDouble(prhs[" << idx << "]) && !mxIsComplex(prhs[" << idx
<< "]) && !mxIsSparse(prhs[" << idx << "]) && mxGetNumberOfElements(prhs[" << idx
<< "]) == " << heterogeneity_table.aggregateEndoSize() << "))" << endl
<< R"( mexErrMsgTxt("yagg must be a real dense numeric array with )"
<< heterogeneity_table.aggregateEndoSize() << R"( elements");)" << endl
<< " const double *restrict yagg = mxGetPr(prhs[" << idx << "]);" << endl;
}
};
// Helper for dealing with sparse_rowval and sparse_colptr inputs (shared with block case)
auto sparse_indices_inputs = [&](int ncols, int nzval) {
// We use sparse_rowval and sparse_colptr (sparse_colval is unused)
const int row_idx {3 + static_cast<int>(dynamic)}, col_idx {row_idx + 2};
const int row_idx {3 + nextra_args}, col_idx {row_idx + 2};
output << " if (!(mxIsInt32(prhs[" << row_idx << "]) && mxGetNumberOfElements(prhs[" << row_idx
<< "]) == " << nzval << "))" << endl
<< R"( mexErrMsgTxt("sparse_rowval must be an int32 array with )" << nzval
......@@ -2778,7 +2778,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
const string prototype_tt {
"void " + funcname
+ "_tt(const double *restrict y, const double *restrict x, const double *restrict params"
+ ss_argin + ", double *restrict T)"};
+ extra_argin + ", double *restrict T)"};
const filesystem::path header_tt {model_src_dir / (funcname + "_tt.h")};
open_file(header_tt);
......@@ -2803,7 +2803,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
const string prototype_main {
"void " + funcname
+ "(const double *restrict y, const double *restrict x, const double *restrict params"
+ ss_argin + ", const double *restrict T, double *restrict "
+ extra_argin + ", const double *restrict T, double *restrict "
+ (i == 0 ? "residual" : "g" + to_string(i) + "_v") + ")"};
const filesystem::path header_main {model_src_dir / (funcname + ".h")};
......@@ -2828,7 +2828,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
compileMEX(model_src_dir, funcname, mexext, {source_main}, matlabroot, false)};
const filesystem::path source_mex {model_src_dir / (funcname + "_mex.c")};
int nargin {5 + static_cast<int>(dynamic) + 3 * static_cast<int>(i == 1)};
int nargin {5 + nextra_args + 3 * static_cast<int>(i == 1)};
open_file(source_mex);
output << "#include <string.h>" << endl // For memcpy()
<< R"(#include "mex.h")" << endl
......@@ -2848,7 +2848,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
<< " if (nlhs != 1 && nlhs != 3)" << endl
<< R"( mexErrMsgTxt("Accepts exactly 1 or 3 output arguments");)" << endl;
y_x_params_ss_inputs(true);
y_x_params_ss_yagg_inputs(true);
if (i == 1)
sparse_indices_inputs(getJacobianColsNbr(true), jacobian_sparse_column_major_order.size());
......@@ -2896,7 +2896,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
output << " default:" << endl << " " << prefix << "resid";
else
output << " case " << j - 1 << ":" << endl << " " << prefix << "g" << j;
output << "_tt(y, x, params" << ss_argout << ", T);" << endl;
output << "_tt(y, x, params" << extra_argout << ", T);" << endl;
}
output << " }" << endl;
if (i == 1)
......@@ -2906,7 +2906,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
output << " plhs[0] = mxCreateDoubleMatrix("
<< (i == 0 ? equations.size() : derivatives[i].size()) << ", 1, mxREAL);" << endl;
output << " " << prefix << (i == 0 ? "resid" : "g" + to_string(i)) << "(y, x, params"
<< ss_argout << ", T, mxGetPr(plhs[0]));" << endl
<< extra_argout << ", T, mxGetPr(plhs[0]));" << endl
<< " if (nlhs == 3)" << endl
<< " {" << endl
<< " plhs[1] = T_order_mx;" << endl
......@@ -2944,6 +2944,32 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
output << "#include <math.h>" << endl << R"(#include "mex.h")" << endl << endl;
writeCHelpersDefinition(output);
writeCHelpersDeclaration(output); // Provide external definition of helpers
// Function computing residuals and/or endogenous, and temporary terms (incl. those for
// derivatives)
output << endl
<< "void " << funcname
<< "_resid(double *restrict y, const double *restrict x, const double *restrict "
"params"
<< extra_argin << ", double *restrict T"
<< (evaluate ? "" : ", double *restrict residual") << ")" << endl
<< "{" << endl;
writePerBlockHelper<output_type>(blk, output, temporary_terms_written);
output << "}" << endl;
// Function computing the Jacobian
if (!evaluate)
{
output << endl
<< "void " << funcname
<< "_g1(const double *restrict y, const double *restrict x, const double "
"*restrict params"
<< extra_argin << ", double *restrict T, double *restrict g1_v)" << endl
<< "{" << endl;
writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
output << "}" << endl;
}
output << endl
<< "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])"
<< endl
......@@ -2959,7 +2985,7 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
evaluate the recursive variables. */
output << " if (nlhs < 2 || nlhs > 4)" << endl
<< R"( mexErrMsgTxt("Accepts 2 to 4 output arguments");)" << endl;
y_x_params_ss_inputs(false);
y_x_params_ss_yagg_inputs(false);
/* We’d like to avoid copying y if this is a “solve” block without
recursive variables. Unfortunately “plhs[0]=prhs[0]” leads to
......@@ -2998,22 +3024,26 @@ ModelTree::writeSparseModelCFiles(const string& basename, const string& mexext,
if (!evaluate)
output << " mxArray *residual_mx = mxCreateDoubleMatrix(" << blocks[blk].mfs_size
<< ", 1, mxREAL);" << endl
<< " double *restrict residual = mxGetPr(residual_mx);" << endl
<< " if (nlhs > 2)" << endl
<< " plhs[2] = residual_mx;" << endl;
<< " double *restrict residual = mxGetPr(residual_mx);" << endl;
// Write residuals and temporary terms (incl. those for derivatives)
writePerBlockHelper<output_type>(blk, output, temporary_terms_written);
output << " " << funcname << "_resid(y, x, params" << extra_argout << ", T"
<< (evaluate ? "" : ", residual") << ");" << endl;
if (!evaluate)
{
output << " if (nlhs > 2)" << endl
<< " plhs[2] = residual_mx;" << endl
<< " else" << endl
<< " mxDestroyArray(residual_mx);" << endl;
// Write Jacobian
output << " if (nlhs > 3)" << endl << " {" << endl;
sparse_jacobian_create(3, blocks[blk].mfs_size, g1_ncols,
blocks_jacobian_sparse_column_major_order[blk].size());
output << " double *restrict g1_v = mxGetPr(plhs[3]);" << endl;
writeSparsePerBlockJacobianHelper<output_type>(blk, output, temporary_terms_written);
output << " }" << endl;
output << " double *restrict g1_v = mxGetPr(plhs[3]);" << endl
<< " " << funcname << "_g1(y, x, params" << extra_argout << ", T, g1_v);"
<< endl
<< " }" << endl;
}
output << "}" << endl;
output.close();
......@@ -3030,9 +3060,7 @@ ModelTree::writeDebugModelMFiles(const string& basename) const
: ExprNodeOutputType::matlabSparseStaticModel};
const filesystem::path m_dir {packageDir(basename) / "+debug"};
// TODO: when C++20 support is complete, mark the following strings constexpr
const string prefix {dynamic ? "dynamic_" : "static_"};
const string ss_arg {dynamic ? ", steady_state" : ""};
const filesystem::path resid_filename {m_dir / (prefix + "resid.m")};
ofstream output {resid_filename, ios::out | ios::binary};
......@@ -3042,7 +3070,13 @@ ModelTree::writeDebugModelMFiles(const string& basename) const
exit(EXIT_FAILURE);
}
output << "function [lhs, rhs] = " << prefix << "resid(y, x, params" << ss_arg << ")" << endl
output << "function [lhs, rhs] = " << prefix << "resid(y, x, params";
if (dynamic)
output << ", steady_state";
if (!heterogeneity_table.empty())
output << ", yagg";
output << ")" << endl
<< "T = NaN(" << temporary_terms_derivatives[0].size() << ", 1);" << endl
<< "lhs = NaN(" << equations.size() << ", 1);" << endl
<< "rhs = NaN(" << equations.size() << ", 1);" << endl;
......@@ -3133,4 +3167,53 @@ ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) co
}
}
template<bool dynamic>
void
ModelTree::writeComplementarityConditionsFile(const string& basename,
const optional<int>& heterogeneous_dimension) const
{
const string funcname {
(dynamic ? "dynamic"s : "static"s)
+ (heterogeneous_dimension ? "_het"s + to_string(*heterogeneous_dimension + 1) : ""s)
+ "_complementarity_conditions"};
const filesystem::path filename {packageDir(basename) / (funcname + ".m")};
/* Can’t use matlabOutsideModel for output type, since it uses M_.
Static is ok even for the dynamic model, since there are no leads/lags. */
constexpr ExprNodeOutputType output_type {ExprNodeOutputType::matlabStaticModel};
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 [lb, ub] = " << funcname << "(params)" << endl
<< "ub = inf(" << equations.size() << ",1);" << endl
<< "lb = -ub;" << endl;
for (const auto& it : complementarity_conditions)
if (it)
{
const auto& [symb_id, lb, ub] = *it;
int endo_id {symbol_table.getTypeSpecificID(symb_id)};
if (lb)
{
output << "lb(" << endo_id + 1 << ")=";
lb->writeOutput(output, output_type);
output << ";" << endl;
}
if (ub)
{
output << "ub(" << endo_id + 1 << ")=";
ub->writeOutput(output, output_type);
output << ";" << endl;
}
}
output << "end" << endl;
output.close();
}
#endif
src/NumericalInitialization.cc
View file @ 6618d7d1
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2024 Dynare Team
*
* This file is part of Dynare.
*
......@@ -303,10 +303,10 @@ EndValStatement::writeJsonOutput(ostream& output) const
output << "]}";
}
EndValLearntInStatement::EndValLearntInStatement(int learnt_in_period_arg,
EndValLearntInStatement::EndValLearntInStatement(variant<int, string> learnt_in_period_arg,
learnt_end_values_t learnt_end_values_arg,
const SymbolTable& symbol_table_arg) :
learnt_in_period {learnt_in_period_arg},
learnt_in_period {move(learnt_in_period_arg)},
learnt_end_values {move(learnt_end_values_arg)},
symbol_table {symbol_table_arg}
{
......@@ -343,9 +343,10 @@ EndValLearntInStatement::writeOutput(ostream& output, [[maybe_unused]] const str
{
if (symbol_table.getType(symb_id) == SymbolType::unusedEndogenous) // See #82
continue;
output << "struct('learnt_in'," << learnt_in_period << ",'exo_id',"
<< symbol_table.getTypeSpecificID(symb_id) + 1 << ",'type','" << typeToString(type)
<< "'"
output << "struct('learnt_in',";
visit([&](const auto& p) { output << p; }, learnt_in_period);
output << ",'exo_id'," << symbol_table.getTypeSpecificID(symb_id) + 1 << ",'type','"
<< typeToString(type) << "'"
<< ",'value',";
value->writeOutput(output);
output << ");" << endl;
......@@ -356,7 +357,18 @@ EndValLearntInStatement::writeOutput(ostream& output, [[maybe_unused]] const str
void
EndValLearntInStatement::writeJsonOutput(ostream& output) const
{
output << R"({"statementName": "endval", "learnt_in": )" << learnt_in_period << R"(, "vals": [)";
output << R"({"statementName": "endval", "learnt_in": )";
visit(
[&]<class T>(const T& p) {
if constexpr (is_same_v<T, int>)
output << p;
else if constexpr (is_same_v<T, string>)
output << '"' << p << '"';
else
static_assert(always_false_v<T>, "Non-exhaustive visitor!");
},
learnt_in_period);
output << R"(, "vals": [)";
for (bool printed_something {false}; auto& [type, symb_id, value] : learnt_end_values)
{
if (symbol_table.getType(symb_id) == SymbolType::unusedEndogenous) // See #82
......
src/NumericalInitialization.hh
View file @ 6618d7d1
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2024 Dynare Team
*
* This file is part of Dynare.
*
......@@ -23,6 +23,7 @@
#include <filesystem>
#include <map>
#include <string>
#include <variant>
#include <vector>
#include "ExprNode.hh"
......@@ -101,7 +102,7 @@ public:
class EndValLearntInStatement : public Statement
{
public:
const int learnt_in_period;
const variant<int, string> learnt_in_period;
enum class LearntEndValType
{
level,
......@@ -117,7 +118,8 @@ private:
static string typeToString(LearntEndValType type);
public:
EndValLearntInStatement(int learnt_in_period_arg, learnt_end_values_t learnt_end_values_arg,
EndValLearntInStatement(variant<int, string> learnt_in_period_arg,
learnt_end_values_t learnt_end_values_arg,
const SymbolTable& symbol_table_arg);
void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
......
src/ParsingDriver.cc
View file @ 6618d7d1
/*
* Copyright © 2003-2024 Dynare Team
* Copyright © 2003-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -22,11 +22,14 @@
#include <fstream>
#include <iostream>
#include <numeric>
#include <ranges>
#include <sstream>
#include "ExprNode.hh"
#include "ParsingDriver.hh"
#include "Statement.hh"
/* NB: the following also imports our specialization of operator<< for location class,
used below in error() and undeclared_model_variable_error() */
#include "WarningConsolidation.hh"
bool
......@@ -71,6 +74,7 @@ ParsingDriver::set_current_data_tree(DataTree* data_tree_arg)
data_tree = data_tree_arg;
model_tree = dynamic_cast<ModelTree*>(data_tree_arg);
dynamic_model = dynamic_cast<DynamicModel*>(data_tree_arg);
heterogeneous_model = dynamic_cast<HeterogeneousModel*>(data_tree_arg);
}
void
......@@ -111,7 +115,7 @@ ParsingDriver::parse(istream& in, bool debug)
void
ParsingDriver::error(const Dynare::parser::location_type& l, const string& m)
{
create_error_string(l, m, cerr);
cerr << "ERROR: " << l << ": " << m << endl;
exit(EXIT_FAILURE);
}
......@@ -121,54 +125,12 @@ ParsingDriver::error(const string& m)
error(location, m);
}
void
ParsingDriver::create_error_string(const Dynare::parser::location_type& l, const string& m,
ostream& stream)
{
stream << "ERROR: " << *l.begin.filename << ": line " << l.begin.line;
if (l.begin.line == l.end.line)
if (l.begin.column == l.end.column - 1)
stream << ", col " << l.begin.column;
else
stream << ", cols " << l.begin.column << "-" << l.end.column - 1;
else
stream << ", col " << l.begin.column << " -"
<< " line " << l.end.line << ", col " << l.end.column - 1;
stream << ": " << m << endl;
}
void
ParsingDriver::create_error_string(const Dynare::parser::location_type& l, const string& m,
const string& var)
{
ostringstream stream;
create_error_string(l, m, stream);
model_errors.emplace_back(var, stream.str());
}
void
ParsingDriver::model_error(const string& m, const string& var)
{
create_error_string(location, m, var);
}
void
ParsingDriver::undeclared_model_variable_error(const string& m, const string& var)
{
ostringstream stream;
if (!nostrict)
{
stream << "ERROR: " << *location.begin.filename << ": line " << location.begin.line;
if (location.begin.line == location.end.line)
if (location.begin.column == location.end.column - 1)
stream << ", col " << location.begin.column;
else
stream << ", cols " << location.begin.column << "-" << location.end.column - 1;
else
stream << ", col " << location.begin.column << " -"
<< " line " << location.end.line << ", col " << location.end.column - 1;
stream << ": ";
}
stream << "ERROR: " << location << ": ";
stream << m;
if (nostrict)
stream << " automatically declared exogenous.";
......@@ -183,12 +145,14 @@ ParsingDriver::warning(const string& m)
int
ParsingDriver::declare_symbol(const string& name, SymbolType type, const string& tex_name,
const vector<pair<string, string>>& partition_value)
const vector<pair<string, string>>& partition_value,
const optional<int>& heterogeneity_dimension)
{
int symb_id;
try
{
symb_id = mod_file->symbol_table.addSymbol(name, type, tex_name, partition_value);
symb_id = mod_file->symbol_table.addSymbol(name, type, tex_name, partition_value,
heterogeneity_dimension);
}
catch (SymbolTable::AlreadyDeclaredException& e)
{
......@@ -207,16 +171,30 @@ int
ParsingDriver::declare_endogenous(const string& name, const string& tex_name,
const vector<pair<string, string>>& partition_value)
{
return declare_symbol(name, SymbolType::endogenous, tex_name, partition_value);
return declare_symbol(name, SymbolType::endogenous, tex_name, partition_value, {});
}
void
ParsingDriver::var(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
bool log_option)
const optional<string>& heterogeneity_dimension, bool log_option)
{
for (auto& [name, tex_name, partition] : symbol_list)
{
int symb_id = declare_endogenous(name, tex_name, partition);
int symb_id {[&] {
if (heterogeneity_dimension)
try
{
return declare_symbol(name, SymbolType::heterogeneousEndogenous, tex_name, partition,
mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
}
catch (HeterogeneityTable::UnknownDimensionNameException&)
{
error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
}
else
return declare_endogenous(name, tex_name, partition);
}()};
if (log_option)
mod_file->symbol_table.markWithLogTransform(symb_id);
}
......@@ -226,14 +204,26 @@ int
ParsingDriver::declare_exogenous(const string& name, const string& tex_name,
const vector<pair<string, string>>& partition_value)
{
return declare_symbol(name, SymbolType::exogenous, tex_name, partition_value);
return declare_symbol(name, SymbolType::exogenous, tex_name, partition_value, {});
}
void
ParsingDriver::varexo(
const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
const optional<string>& heterogeneity_dimension)
{
for (auto& [name, tex_name, partition] : symbol_list)
if (heterogeneity_dimension)
try
{
declare_symbol(name, SymbolType::heterogeneousExogenous, tex_name, partition,
mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
}
catch (HeterogeneityTable::UnknownDimensionNameException&)
{
error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
}
else
declare_exogenous(name, tex_name, partition);
}
......@@ -242,21 +232,33 @@ ParsingDriver::varexo_det(
const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
{
for (auto& [name, tex_name, partition] : symbol_list)
declare_symbol(name, SymbolType::exogenousDet, tex_name, partition);
declare_symbol(name, SymbolType::exogenousDet, tex_name, partition, {});
}
int
ParsingDriver::declare_parameter(const string& name, const string& tex_name,
const vector<pair<string, string>>& partition_value)
{
return declare_symbol(name, SymbolType::parameter, tex_name, partition_value);
return declare_symbol(name, SymbolType::parameter, tex_name, partition_value, {});
}
void
ParsingDriver::parameters(
const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list)
const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
const optional<string>& heterogeneity_dimension)
{
for (auto& [name, tex_name, partition] : symbol_list)
if (heterogeneity_dimension)
try
{
declare_symbol(name, SymbolType::heterogeneousParameter, tex_name, partition,
mod_file->heterogeneity_table.getID(*heterogeneity_dimension));
}
catch (HeterogeneityTable::UnknownDimensionNameException&)
{
error("Unknown heterogeneity dimension: " + *heterogeneity_dimension);
}
else
declare_parameter(name, tex_name, partition);
}
......@@ -266,7 +268,7 @@ ParsingDriver::declare_statement_local_variable(const string& name)
if (mod_file->symbol_table.exists(name))
error("Symbol " + name + " cannot be assigned within a statement "
+ "while being assigned elsewhere in the modfile");
declare_symbol(name, SymbolType::statementDeclaredVariable, "", {});
declare_symbol(name, SymbolType::statementDeclaredVariable, "", {}, {});
}
void
......@@ -293,7 +295,7 @@ ParsingDriver::end_trend_var(bool log_trend, expr_t growth_factor,
for (auto& [name, tex_name] : symbol_list)
{
int symb_id = declare_symbol(name, log_trend ? SymbolType::logTrend : SymbolType::trend,
tex_name, {});
tex_name, {}, {});
declared_trend_vars.push_back(symb_id);
}
......@@ -418,10 +420,6 @@ ParsingDriver::add_model_variable(int symb_id, int lag)
error("Exogenous deterministic variable " + mod_file->symbol_table.getName(symb_id)
+ " cannot be given a lead or a lag.");
if (type == SymbolType::modelLocalVariable && lag != 0)
error("Model local variable " + mod_file->symbol_table.getName(symb_id)
+ " cannot be given a lead or a lag.");
if (data_tree == planner_objective.get())
{
if (lag != 0)
......@@ -788,22 +786,26 @@ ParsingDriver::end_endval(bool all_values_required)
}
void
ParsingDriver::end_endval_learnt_in(const string& learnt_in_period)
ParsingDriver::end_endval_learnt_in(variant<int, string> learnt_in_period)
{
int learnt_in_period_int = stoi(learnt_in_period);
if (holds_alternative<int>(learnt_in_period))
{
int learnt_in_period_int = get<int>(learnt_in_period);
if (learnt_in_period_int < 1)
error("endval: value '" + learnt_in_period + "' is not allowed for 'learnt_in' option");
error("endval: value '" + to_string(learnt_in_period_int)
+ "' is not allowed for 'learnt_in' option");
if (learnt_in_period_int == 1)
{
end_endval(false);
return;
}
}
for (auto [type, symb_id, value] : end_values)
if (mod_file->symbol_table.getType(symb_id) != SymbolType::exogenous)
error("endval(learnt_in=...): " + mod_file->symbol_table.getName(symb_id)
+ " is not an exogenous variable");
mod_file->addStatement(make_unique<EndValLearntInStatement>(
learnt_in_period_int, move(end_values), mod_file->symbol_table));
move(learnt_in_period), move(end_values), mod_file->symbol_table));
end_values.clear();
}
......@@ -840,7 +842,7 @@ ParsingDriver::end_epilogue()
void
ParsingDriver::add_epilogue_variable(const string& name)
{
declare_symbol(name, SymbolType::epilogue, "", {});
declare_symbol(name, SymbolType::epilogue, "", {}, {});
}
void
......@@ -856,17 +858,26 @@ ParsingDriver::begin_model()
}
void
ParsingDriver::end_model()
ParsingDriver::begin_heterogeneous_model(const string& heterogeneity_dimension)
{
bool exit_after_write = false;
if (model_errors.size() > 0)
for (auto& it : model_errors)
int het_dim_id {[&] {
try
{
if (it.first.empty())
exit_after_write = true;
cerr << it.second;
return mod_file->heterogeneity_table.getID(heterogeneity_dimension);
}
catch (HeterogeneityTable::UnknownDimensionNameException&)
{
error("Unknown heterogeneity dimension: " + heterogeneity_dimension);
}
}()};
set_current_data_tree(&mod_file->heterogeneous_models.at(het_dim_id));
}
void
ParsingDriver::end_model()
{
bool exit_after_write = false;
if (undeclared_model_variable_errors.size() > 0)
for (auto& it : undeclared_model_variable_errors)
{
......@@ -904,6 +915,28 @@ ParsingDriver::end_shocks(bool overwrite)
corr_shocks.clear();
}
void
ParsingDriver::end_heterogeneous_shocks(const string& heterogeneity_dimension, bool overwrite)
{
int het_dim_id {[&] {
try
{
return mod_file->heterogeneity_table.getID(heterogeneity_dimension);
}
catch (HeterogeneityTable::UnknownDimensionNameException&)
{
error("Unknown heterogeneity dimension: " + heterogeneity_dimension);
}
}()};
mod_file->addStatement(make_unique<HeterogeneousShocksStatement>(
het_dim_id, overwrite, move(var_shocks), move(std_shocks), move(covar_shocks),
move(corr_shocks), mod_file->symbol_table, mod_file->heterogeneity_table));
var_shocks.clear();
std_shocks.clear();
covar_shocks.clear();
corr_shocks.clear();
}
void
ParsingDriver::end_mshocks(bool overwrite, bool relative_to_initval)
{
......@@ -919,6 +952,11 @@ ParsingDriver::end_mshocks(bool overwrite, bool relative_to_initval)
void
ParsingDriver::end_shocks_surprise(bool overwrite)
{
if (ranges::any_of(views::values(det_shocks), [](auto& v) {
return ranges::any_of(views::keys(v),
[](auto& p) { return !holds_alternative<pair<int, int>>(p); });
}))
error("shocks(surprise): dates are not allowed in the 'periods' keyword");
mod_file->addStatement(
make_unique<ShocksSurpriseStatement>(overwrite, move(det_shocks), mod_file->symbol_table));
det_shocks.clear();
......@@ -929,11 +967,14 @@ ParsingDriver::end_shocks_surprise(bool overwrite)
}
void
ParsingDriver::end_shocks_learnt_in(const string& learnt_in_period, bool overwrite)
ParsingDriver::end_shocks_learnt_in(variant<int, string> learnt_in_period, bool overwrite)
{
int learnt_in_period_int = stoi(learnt_in_period);
if (holds_alternative<int>(learnt_in_period))
{
int learnt_in_period_int = get<int>(learnt_in_period);
if (learnt_in_period_int < 1)
error("shocks: value '" + learnt_in_period + "' is not allowed for 'learnt_in' option");
error("shocks: value '" + to_string(learnt_in_period_int)
+ "' is not allowed for 'learnt_in' option");
if (learnt_in_period_int == 1)
{
end_shocks(overwrite);
......@@ -941,65 +982,78 @@ ParsingDriver::end_shocks_learnt_in(const string& learnt_in_period, bool overwri
}
for (auto& storage : {det_shocks, learnt_shocks_add, learnt_shocks_multiply})
for (auto& [symb_id, vals] : storage)
for (auto [period1, period2, expr] : vals)
if (period1 < learnt_in_period_int)
for (const auto& [period_range, expr] : vals)
if (holds_alternative<pair<int, int>>(period_range))
if (int period1 = get<pair<int, int>>(period_range).first;
period1 < learnt_in_period_int)
error("shocks: for variable " + mod_file->symbol_table.getName(symb_id)
+ ", shock period (" + to_string(period1)
+ ") is earlier than the period in which the shock is learnt (" + learnt_in_period
+ ")");
+ ") is earlier than the period in which the shock is learnt ("
+ to_string(learnt_in_period_int) + ")");
}
// Aggregate the three types of shocks
ShocksLearntInStatement::learnt_shocks_t learnt_shocks;
for (const auto& [id, v] : det_shocks)
{
vector<tuple<ShocksLearntInStatement::LearntShockType, int, int, expr_t>> v2;
for (auto [period1, period2, value] : v)
v2.emplace_back(ShocksLearntInStatement::LearntShockType::level, period1, period2, value);
vector<tuple<ShocksLearntInStatement::LearntShockType,
AbstractShocksStatement::period_range_t, expr_t>>
v2;
for (const auto& [period_range, value] : v)
v2.emplace_back(ShocksLearntInStatement::LearntShockType::level, period_range, value);
learnt_shocks[id] = v2;
}
for (const auto& [id, v] : learnt_shocks_add)
{
vector<tuple<ShocksLearntInStatement::LearntShockType, int, int, expr_t>> v2;
for (auto [period1, period2, value] : v)
v2.emplace_back(ShocksLearntInStatement::LearntShockType::add, period1, period2, value);
vector<tuple<ShocksLearntInStatement::LearntShockType,
AbstractShocksStatement::period_range_t, expr_t>>
v2;
for (const auto& [period_range, value] : v)
v2.emplace_back(ShocksLearntInStatement::LearntShockType::add, period_range, value);
learnt_shocks[id] = v2;
}
for (const auto& [id, v] : learnt_shocks_multiply)
{
vector<tuple<ShocksLearntInStatement::LearntShockType, int, int, expr_t>> v2;
for (auto [period1, period2, value] : v)
v2.emplace_back(ShocksLearntInStatement::LearntShockType::multiply, period1, period2,
value);
vector<tuple<ShocksLearntInStatement::LearntShockType,
AbstractShocksStatement::period_range_t, expr_t>>
v2;
for (const auto& [period_range, value] : v)
v2.emplace_back(ShocksLearntInStatement::LearntShockType::multiply, period_range, value);
learnt_shocks[id] = v2;
}
mod_file->addStatement(make_unique<ShocksLearntInStatement>(
learnt_in_period_int, overwrite, move(learnt_shocks), mod_file->symbol_table));
move(learnt_in_period), overwrite, move(learnt_shocks), mod_file->symbol_table));
det_shocks.clear();
learnt_shocks_add.clear();
learnt_shocks_multiply.clear();
}
void
ParsingDriver::end_mshocks_learnt_in(const string& learnt_in_period, bool overwrite,
ParsingDriver::end_mshocks_learnt_in(variant<int, string> learnt_in_period, bool overwrite,
bool relative_to_initval)
{
int learnt_in_period_int = stoi(learnt_in_period);
if (holds_alternative<int>(learnt_in_period))
{
int learnt_in_period_int = get<int>(learnt_in_period);
if (learnt_in_period_int < 1)
error("mshocks: value '" + learnt_in_period + "' is not allowed for 'learnt_in' option");
error("mshocks: value '" + to_string(learnt_in_period_int)
+ "' is not allowed for 'learnt_in' option");
if (learnt_in_period_int == 1)
{
end_mshocks(overwrite, relative_to_initval);
return;
}
for (auto& [symb_id, vals] : det_shocks)
for (auto [period1, period2, expr] : vals)
if (period1 < learnt_in_period_int)
for (const auto& [period_range, expr] : vals)
if (holds_alternative<pair<int, int>>(period_range))
if (int period1 = get<pair<int, int>>(period_range).first;
period1 < learnt_in_period_int)
error("mshocks: for variable " + mod_file->symbol_table.getName(symb_id)
+ ", shock period (" + to_string(period1)
+ ") is earlier than the period in which the shock is learnt (" + learnt_in_period
+ ")");
+ ") is earlier than the period in which the shock is learnt ("
+ to_string(learnt_in_period_int) + ")");
}
ShocksLearntInStatement::learnt_shocks_t learnt_shocks;
const auto type {relative_to_initval
......@@ -1007,14 +1061,16 @@ ParsingDriver::end_mshocks_learnt_in(const string& learnt_in_period, bool overwr
: ShocksLearntInStatement::LearntShockType::multiplySteadyState};
for (const auto& [id, v] : det_shocks)
{
vector<tuple<ShocksLearntInStatement::LearntShockType, int, int, expr_t>> v2;
for (auto [period1, period2, value] : v)
v2.emplace_back(type, period1, period2, value);
vector<tuple<ShocksLearntInStatement::LearntShockType,
AbstractShocksStatement::period_range_t, expr_t>>
v2;
for (const auto& [period_range, value] : v)
v2.emplace_back(type, period_range, value);
learnt_shocks[id] = v2;
}
mod_file->addStatement(make_unique<ShocksLearntInStatement>(
learnt_in_period_int, overwrite, move(learnt_shocks), mod_file->symbol_table));
move(learnt_in_period), overwrite, move(learnt_shocks), mod_file->symbol_table));
det_shocks.clear();
if (!learnt_shocks_add.empty())
error("mshocks: 'add' keyword not allowed");
......@@ -1033,7 +1089,8 @@ ParsingDriver::end_heteroskedastic_shocks(bool overwrite)
}
void
ParsingDriver::add_det_shock(const string& var, const vector<pair<int, int>>& periods,
ParsingDriver::add_det_shock(const string& var,
const vector<AbstractShocksStatement::period_range_t>& periods,
const vector<expr_t>& values, DetShockType type)
{
switch (type)
......@@ -1061,10 +1118,10 @@ ParsingDriver::add_det_shock(const string& var, const vector<pair<int, int>>& pe
error("shocks/conditional_forecast_paths: variable " + var
+ ": number of periods is different from number of shock values");
vector<tuple<int, int, expr_t>> v;
vector<pair<AbstractShocksStatement::period_range_t, expr_t>> v;
for (size_t i = 0; i < periods.size(); i++)
v.emplace_back(periods[i].first, periods[i].second, values[i]);
v.emplace_back(periods[i], values[i]);
switch (type)
{
......@@ -1082,11 +1139,15 @@ ParsingDriver::add_det_shock(const string& var, const vector<pair<int, int>>& pe
}
void
ParsingDriver::add_heteroskedastic_shock(const string& var, const vector<pair<int, int>>& periods,
ParsingDriver::add_heteroskedastic_shock(
const string& var, const vector<AbstractShocksStatement::period_range_t>& periods,
const vector<expr_t>& values, bool scales)
{
check_symbol_is_exogenous(var, false);
if (ranges::any_of(periods, [](auto& p) { return !holds_alternative<pair<int, int>>(p); }))
error("heteroskedastic_shocks: dates are not allowed in the 'periods' keyword");
int symb_id = mod_file->symbol_table.getID(var);
if ((!scales && heteroskedastic_shocks_values.contains(symb_id))
......@@ -1099,7 +1160,10 @@ ParsingDriver::add_heteroskedastic_shock(const string& var, const vector<pair<in
vector<tuple<int, int, expr_t>> v;
for (size_t i = 0; i < periods.size(); i++)
v.emplace_back(periods[i].first, periods[i].second, values[i]);
{
auto [period1, period2] = get<pair<int, int>>(periods[i]);
v.emplace_back(period1, period2, values[i]);
}
if (scales)
heteroskedastic_shocks_scales[symb_id] = v;
......@@ -1292,7 +1356,7 @@ void
ParsingDriver::add_positive_restriction_element(expr_t value, const string& variable,
const string& lag)
{
// if the expression is not on the left handside, change its sign
// if the expression is not on the left-hand side, change its sign
if (!svar_left_handside)
value = add_uminus(value);
......@@ -1304,7 +1368,7 @@ ParsingDriver::add_positive_restriction_element(const string& variable, const st
{
expr_t value(data_tree->One);
// if the expression is not on the left handside, change its sign
// if the expression is not on the left-hand side, change its sign
if (!svar_left_handside)
value = add_uminus(value);
......@@ -1315,7 +1379,7 @@ void
ParsingDriver::add_negative_restriction_element(expr_t value, const string& variable,
const string& lag)
{
// if the expression is on the left handside, change its sign
// if the expression is on the left-hand side, change its sign
if (svar_left_handside)
value = add_uminus(value);
......@@ -1327,7 +1391,7 @@ ParsingDriver::add_negative_restriction_element(const string& variable, const st
{
expr_t value(data_tree->One);
// if the expression is on the left handside, change its sign
// if the expression is on the left-hand side, change its sign
if (svar_left_handside)
value = add_uminus(value);
......@@ -1710,11 +1774,9 @@ ParsingDriver::set_unit_root_vars()
}
void
ParsingDriver::set_time(const string& arg)
ParsingDriver::set_time(string period)
{
option_date("initial_period", arg);
mod_file->addStatement(make_unique<SetTimeStatement>(move(options_list)));
options_list.clear();
mod_file->addStatement(make_unique<SetTimeStatement>(move(period)));
}
void
......@@ -2224,7 +2286,8 @@ void
ParsingDriver::begin_planner_objective()
{
planner_objective = make_unique<PlannerObjective>(mod_file->symbol_table, mod_file->num_constants,
mod_file->external_functions_table);
mod_file->external_functions_table,
mod_file->heterogeneity_table);
set_current_data_tree(planner_objective.get());
}
......@@ -2601,6 +2664,11 @@ ParsingDriver::plot_conditional_forecast(const optional<string>& periods,
void
ParsingDriver::conditional_forecast_paths()
{
if (ranges::any_of(views::values(det_shocks), [](auto& v) {
return ranges::any_of(views::keys(v),
[](auto& p) { return !holds_alternative<pair<int, int>>(p); });
}))
error("conditional_forecast_paths: dates are not allowed in the 'periods' keyword");
mod_file->addStatement(
make_unique<ConditionalForecastPathsStatement>(move(det_shocks), mod_file->symbol_table));
det_shocks.clear();
......@@ -2635,6 +2703,9 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
if (key == "endogenous")
declare_or_change_type(SymbolType::endogenous, value);
optional<tuple<int, expr_t, expr_t>> matched_complementarity_condition;
// Handle legacy “mcp” tags, for backward compatibility
if (eq_tags.contains("mcp"))
{
if (complementarity_condition)
......@@ -2643,18 +2714,64 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
else
warning("Specifying complementarity conditions with the 'mcp' tag is obsolete. Please "
"consider switching to the new syntax using the perpendicular symbol.");
auto [var_name, is_lower_bound, constant] {[&]() -> tuple<string, bool, string> {
auto tagval = eq_tags.at("mcp");
if (auto less_pos = tagval.find('<'); less_pos != string::npos)
return {tagval.substr(0, less_pos), false, tagval.substr(less_pos + 1)};
else if (auto greater_pos = tagval.find('>'); greater_pos != string::npos)
return {tagval.substr(0, greater_pos), true, tagval.substr(greater_pos + 1)};
else
error("'mcp' tag does not contain an inequality");
}()};
// Trim whitespace
var_name.erase(var_name.find_last_not_of(" \n\t") + 1);
var_name.erase(0, var_name.find_first_not_of(" \n\t"));
constant.erase(constant.find_last_not_of(" \n\t") + 1);
constant.erase(0, constant.find_first_not_of(" \n\t"));
int symb_id {[&] {
try
{
return mod_file->symbol_table.getID(var_name);
}
catch (SymbolTable::UnknownSymbolNameException&)
{
error("Left-hand side of expression in 'mcp' tag is not a variable");
}
}()};
if (heterogeneous_model)
error("'mcp' tags are not allowed in heterogeneous model blocks");
if (mod_file->symbol_table.getType(symb_id) != SymbolType::endogenous)
error("Left-hand side of expression in 'mcp' tag is not an endogenous variable");
expr_t matched_constant {[&] {
char* str_end;
double d = strtod(constant.c_str(), &str_end);
if (str_end == constant.c_str())
error("Right-hand side of expression in 'mcp' tag should be a constant");
return data_tree->AddPossiblyNegativeConstant(d);
}()};
matched_complementarity_condition = {symb_id, is_lower_bound ? matched_constant : nullptr,
is_lower_bound ? nullptr : matched_constant};
}
if (complementarity_condition)
try
{
if (auto bcomp = dynamic_cast<BinaryOpNode*>(complementarity_condition);
!(bcomp
&& (bcomp->op_code == BinaryOpcode::less || bcomp->op_code == BinaryOpcode::lessEqual
|| bcomp->op_code == BinaryOpcode::greater
|| bcomp->op_code == BinaryOpcode::greaterEqual)))
error("The complementarity constraint must be an inequality.");
eq_tags.emplace("mcp", complementarity_condition->toString());
matched_complementarity_condition
= heterogeneous_model ? complementarity_condition->matchComplementarityCondition(
heterogeneous_model->heterogeneity_dimension)
: complementarity_condition->matchComplementarityCondition();
}
catch (ExprNode::MatchFailureException& e)
{
error("Complementarity condition has an incorrect form"s
+ (e.message.empty() ? ""s : ": "s + e.message));
}
if (eq_tags.contains("static"))
......@@ -2664,7 +2781,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
error("An equation tagged [static] cannot contain leads, lags, expectations or "
"STEADY_STATE operators");
dynamic_model->addStaticOnlyEquation(id, location.begin.line, eq_tags);
dynamic_model->addStaticOnlyEquation(id, location.begin.line,
move(matched_complementarity_condition), eq_tags);
}
else if (eq_tags.contains("bind") || eq_tags.contains("relax"))
{
......@@ -2740,7 +2858,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
}
}
else // General case
model_tree->addEquation(id, location.begin.line, move(eq_tags));
model_tree->addEquation(id, location.begin.line, move(matched_complementarity_condition),
move(eq_tags));
return id;
}
......@@ -2756,7 +2875,7 @@ void
ParsingDriver::model_local_variable(const vector<pair<string, string>>& symbol_list)
{
for (auto& [name, tex_name] : symbol_list)
declare_symbol(name, SymbolType::modelLocalVariable, tex_name, {});
declare_symbol(name, SymbolType::modelLocalVariable, tex_name, {}, {});
}
void
......@@ -2993,9 +3112,8 @@ ParsingDriver::add_diff(expr_t arg1)
expr_t
ParsingDriver::add_adl(expr_t arg1, const string& name, const string& lag)
{
vector<int> lags(stoi(lag));
iota(lags.begin(), lags.end(), 1);
return add_adl(arg1, name, lags);
auto lags = views::iota(1, stoi(lag) + 1);
return add_adl(arg1, name, vector<int> {lags.begin(), lags.end()});
}
expr_t
......@@ -3176,6 +3294,23 @@ ParsingDriver::add_steady_state(expr_t arg1)
return data_tree->AddSteadyState(arg1);
}
expr_t
ParsingDriver::add_sum(expr_t arg)
{
if (heterogeneous_model)
error("The SUM() operator cannot be used inside a model(heterogeneity=...) block");
VariableNode* varg {dynamic_cast<VariableNode*>(arg)};
if (!varg)
error("The argument to the SUM() operator must be a single variable");
if (varg->lag != 0)
error("The argument to the SUM() operator must not have a lead or lag");
if (mod_file->symbol_table.getType(varg->symb_id) != SymbolType::heterogeneousEndogenous)
error("The argument to the SUM() operator must be a heterogeneous endogenous variable");
return data_tree->AddSum(arg);
}
void
ParsingDriver::external_function_option(const string& name_option, const string& opt)
{
......@@ -3184,7 +3319,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
if (opt.empty())
error("An argument must be passed to the 'name' option of the external_function() "
"statement.");
declare_symbol(opt, SymbolType::externalFunction, "", {});
declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
current_external_function_id = mod_file->symbol_table.getID(opt);
}
else if (name_option == "first_deriv_provided")
......@@ -3194,7 +3329,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
= ExternalFunctionsTable::IDSetButNoNameProvided;
else
{
int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {});
int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
current_external_function_options.firstDerivSymbID = symb_id;
}
}
......@@ -3205,7 +3340,7 @@ ParsingDriver::external_function_option(const string& name_option, const string&
= ExternalFunctionsTable::IDSetButNoNameProvided;
else
{
int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {});
int symb_id = declare_symbol(opt, SymbolType::externalFunction, "", {}, {});
current_external_function_options.secondDerivSymbID = symb_id;
}
}
......@@ -3315,10 +3450,9 @@ ParsingDriver::add_model_var_or_external_function(const string& function_name, b
optional<int> rv {is_there_one_integer_argument()};
if (!rv)
model_error("Symbol " + function_name
+ " is being treated as if it were a function (i.e., takes an argument "
"that is not an integer).",
"");
error("Symbol " + function_name
+ " is being treated as if it were a function (i.e., takes an argument that is "
"not an integer).");
nid = add_model_variable(mod_file->symbol_table.getID(function_name), *rv);
stack_external_function_args.pop();
......@@ -3371,7 +3505,7 @@ ParsingDriver::add_model_var_or_external_function(const string& function_name, b
+ ") within the model block, you must first declare it via the "
"external_function() statement.");
}
int symb_id = declare_symbol(function_name, SymbolType::externalFunction, "", {});
int symb_id = declare_symbol(function_name, SymbolType::externalFunction, "", {}, {});
current_external_function_options.nargs = stack_external_function_args.top().size();
mod_file->external_functions_table.addExternalFunction(
symb_id, current_external_function_options, in_model_block);
......@@ -3613,6 +3747,29 @@ ParsingDriver::perfect_foresight_with_expectation_errors_solver()
options_list.clear();
}
void
ParsingDriver::perfect_foresight_controlled_paths(
const vector<tuple<string, vector<AbstractShocksStatement::period_range_t>, vector<expr_t>,
string>>& paths,
variant<int, string> learnt_in_period)
{
PerfectForesightControlledPathsStatement::paths_t paths_transformed;
for (const auto& [exogenize, periods, values, endogenize] : paths)
{
int exogenize_id = mod_file->symbol_table.getID(exogenize);
int endogenize_id = mod_file->symbol_table.getID(endogenize);
vector<pair<AbstractShocksStatement::period_range_t, expr_t>> v;
for (size_t i = 0; i < periods.size(); i++)
v.emplace_back(periods[i], values[i]);
paths_transformed.emplace_back(exogenize_id, move(v), endogenize_id);
}
mod_file->addStatement(make_unique<PerfectForesightControlledPathsStatement>(
move(paths_transformed), move(learnt_in_period), mod_file->symbol_table));
}
void
ParsingDriver::method_of_moments()
{
......@@ -3629,48 +3786,31 @@ ParsingDriver::prior_posterior_function(bool prior_func)
}
void
ParsingDriver::add_ramsey_constraints_statement()
ParsingDriver::begin_ramsey_constraints()
{
mod_file->addStatement(
make_unique<RamseyConstraintsStatement>(mod_file->symbol_table, move(ramsey_constraints)));
ramsey_constraints.clear();
set_current_data_tree(&mod_file->dynamic_model);
}
void
ParsingDriver::ramsey_constraint_add_less(const string& name, const expr_t rhs)
ParsingDriver::end_ramsey_constraints(const vector<expr_t>& constraints)
{
add_ramsey_constraint(name, BinaryOpcode::less, rhs);
}
void
ParsingDriver::ramsey_constraint_add_greater(const string& name, const expr_t rhs)
for (expr_t c : constraints)
try
{
add_ramsey_constraint(name, BinaryOpcode::greater, rhs);
}
auto [symb_id, lower_bound, upper_bound] = c->matchComplementarityCondition();
void
ParsingDriver::ramsey_constraint_add_less_equal(const string& name, const expr_t rhs)
{
add_ramsey_constraint(name, BinaryOpcode::lessEqual, rhs);
auto [it, success]
= mod_file->ramsey_constraints.try_emplace(symb_id, lower_bound, upper_bound);
if (!success)
error("The ramsey_constraints block contains two constraints for variable "
+ mod_file->symbol_table.getName(symb_id));
}
void
ParsingDriver::ramsey_constraint_add_greater_equal(const string& name, const expr_t rhs)
catch (ExprNode::MatchFailureException& e)
{
add_ramsey_constraint(name, BinaryOpcode::greaterEqual, rhs);
error("Ramsey constraint has an incorrect form: " + e.message);
}
void
ParsingDriver::add_ramsey_constraint(const string& name, BinaryOpcode op_code, const expr_t rhs)
{
check_symbol_is_endogenous(name);
int symb_id = mod_file->symbol_table.getID(name);
RamseyConstraintsStatement::Constraint C;
C.endo = symb_id;
C.code = op_code;
C.expression = rhs;
ramsey_constraints.push_back(C);
reset_data_tree();
}
void
......@@ -3814,7 +3954,8 @@ ParsingDriver::begin_occbin_constraints()
if added to the main DynamicModel tree. It also simplifies the
enforcement of various constraints at parsing time. */
occbin_constraints_tree = make_unique<DataTree>(mod_file->symbol_table, mod_file->num_constants,
mod_file->external_functions_table, false);
mod_file->external_functions_table,
mod_file->heterogeneity_table, false);
set_current_data_tree(occbin_constraints_tree.get());
}
......@@ -3949,3 +4090,26 @@ ParsingDriver::matched_irfs_weights(MatchedIrfsWeightsStatement::matched_irfs_we
{
mod_file->addStatement(make_unique<MatchedIrfsWeightsStatement>(move(weights), overwrite));
}
void
ParsingDriver::heterogeneity_dimension(const vector<string>& dims)
{
for (const auto& dim : dims)
{
int het_dim_id {[&] {
try
{
return mod_file->heterogeneity_table.addDimension(dim);
}
catch (HeterogeneityTable::AlreadyDeclaredDimensionException&)
{
error("Heterogeneity dimension '" + dim + "' already declared");
}
}()};
assert(static_cast<int>(mod_file->heterogeneous_models.size()) == het_dim_id);
mod_file->heterogeneous_models.emplace_back(mod_file->symbol_table, mod_file->num_constants,
mod_file->external_functions_table,
mod_file->heterogeneity_table, het_dim_id);
}
}
src/ParsingDriver.hh
View file @ 6618d7d1
/*
* Copyright © 2003-2024 Dynare Team
* Copyright © 2003-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -29,12 +29,12 @@
#include <stack>
#include <string>
#include <string_view>
#include <variant>
#include <vector>
#include "ModFile.hh"
#include "SymbolList.hh"
class ParsingDriver;
#include "DynareBison.hh"
#include "ExprNode.hh"
......@@ -74,9 +74,6 @@ public:
//! Increment the location counter given a token
static void location_increment(Dynare::parser::location_type* yylloc, const char* yytext);
//! Count parens in dates statement
int dates_parens_nb;
};
//! Drives the scanning and parsing of the .mod file, and constructs its abstract representation
......@@ -114,7 +111,8 @@ private:
//! Helper to add a symbol declaration (returns its symbol ID)
int declare_symbol(const string& name, SymbolType type, const string& tex_name,
const vector<pair<string, string>>& partition_value);
const vector<pair<string, string>>& partition_value,
const optional<int>& heterogeneity_dimension);
//! Temporary store for the planner objective
unique_ptr<PlannerObjective> planner_objective;
......@@ -137,6 +135,11 @@ private:
* DynamicModel instance */
DynamicModel* dynamic_model;
//! The heterogeneous model tree in which to add expressions currently parsed
/*! It is only a dynamic cast of data_tree pointer, and is therefore null if data_tree is not a
* HeterogeneousModel instance */
HeterogeneousModel* heterogeneous_model;
//! Sets data_tree and model_tree pointers
void set_current_data_tree(DataTree* data_tree_arg);
......@@ -186,8 +189,6 @@ private:
MomentCalibration::constraints_t moment_calibration_constraints;
//! Temporary storage for irf_calibration
IrfCalibration::constraints_t irf_calibration_constraints;
//! Temporary storage for ramsey_constraints
RamseyConstraintsStatement::constraints_t ramsey_constraints;
//! Temporary storage for svar_identification blocks
SvarIdentificationStatement::svar_identification_restrictions_t svar_ident_restrictions;
//! Temporary storage for mapping the equation number to the restrictions within an
......@@ -269,7 +270,6 @@ private:
bool nostrict;
vector<pair<string, string>> model_errors;
vector<pair<string, string>> undeclared_model_variable_errors;
//! True when parsing the epilogue block
......@@ -316,23 +316,16 @@ public:
expr_t var_expectation_model_discount {nullptr};
//! Error handler with explicit location
void error(const Dynare::parser::location_type& l, const string& m) __attribute__((noreturn));
[[noreturn]] void error(const Dynare::parser::location_type& l, const string& m);
//! Error handler using saved location
void error(const string& m) __attribute__((noreturn));
[[noreturn]] void error(const string& m);
//! Warning handler using saved location
void warning(const string& m);
//! Error handler with explicit location (used in model block, accumulating error messages to be
//! printed later)
void model_error(const string& m, const string& var);
void undeclared_model_variable_error(const string& m, const string& var);
//! Code shared between model_error() and error()
void create_error_string(const Dynare::parser::location_type& l, const string& m,
const string& var);
void create_error_string(const Dynare::parser::location_type& l, const string& m,
ostream& stream);
//! Check if a given symbol exists in the parsing context, and is not a mod file local variable
bool symbol_exists_and_is_not_modfile_local_or_external_function(const string& s);
//! Sets mode of ModelTree class to use C output
......@@ -373,19 +366,21 @@ public:
const vector<pair<string, string>>& partition_value = {});
// Handles a “var” or “var(log)” statement (without “deflator” or “log_deflator” options)
void var(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
bool log_option);
const optional<string>& heterogeneity_dimension, bool log_option);
//! Declares an exogenous variable (and returns its symbol ID)
int declare_exogenous(const string& name, const string& tex_name = "",
const vector<pair<string, string>>& partition_value = {});
// Handles a “varexo” statement
void varexo(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
void varexo(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
const optional<string>& heterogeneity_dimension);
// Handles a “varexo_det” statement
void varexo_det(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
//! Declares a parameter (and returns its symbol ID)
int declare_parameter(const string& name, const string& tex_name = "",
const vector<pair<string, string>>& partition_value = {});
// Handles a “parameters” statement
void parameters(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list);
void parameters(const vector<tuple<string, string, vector<pair<string, string>>>>& symbol_list,
const optional<string>& heterogeneity_dimension);
// Handles a “model_local_variable” statement
void model_local_variable(const vector<pair<string, string>>& symbol_list);
//! Declares a statement local variable
......@@ -443,7 +438,7 @@ public:
//! Writes end of an endval block
void end_endval(bool all_values_required);
//! Writes end of an endval(learnt_in=…) block
void end_endval_learnt_in(const string& learnt_in_period);
void end_endval_learnt_in(variant<int, string> learnt_in_period);
//! Writes end of an histval block
void end_histval(bool all_values_required);
//! Writes end of an homotopy_setup block
......@@ -456,10 +451,11 @@ public:
void add_epilogue_variable(const string& varname);
//! Add equation in epilogue block
void add_epilogue_equal(const string& varname, expr_t expr);
/* Begin a model or model_replace block, or an expression as an option value
of some statement.
Must be followed by a call to reset_data_tree(). */
/* Begin a model block (without heterogeneity option), or an expression as an option value of some
statement. Must be followed by a call to reset_data_tree(). */
void begin_model();
// Begin a model(heterogeneity=…) block
void begin_heterogeneous_model(const string& heterogeneity_dimension);
//! End a model or model_replace block, printing errors that were encountered in parsing
void end_model();
//! Writes a shocks statement
......@@ -469,9 +465,11 @@ public:
//! Writes a shocks(surprise) statement
void end_shocks_surprise(bool overwrite);
//! Writes a shocks(learnt_in=…) block
void end_shocks_learnt_in(const string& learnt_in_period, bool overwrite);
void end_shocks_learnt_in(variant<int, string> learnt_in_period, bool overwrite);
// For a shocks(heterogeneity=…) block
void end_heterogeneous_shocks(const string& heterogeneity_dimension, bool overwrite);
//! Writes a mshocks(learnt_in=…) block
void end_mshocks_learnt_in(const string& learnt_in_period, bool overwrite,
void end_mshocks_learnt_in(variant<int, string> learnt_in_period, bool overwrite,
bool relative_to_initval);
//! Writes a heteroskedastic_shocks statement
void end_heteroskedastic_shocks(bool overwrite);
......@@ -484,10 +482,12 @@ public:
multiply, // for “multiply” in “shocks(learnt_in)”
conditional_forecast
};
void add_det_shock(const string& var, const vector<pair<int, int>>& periods,
void add_det_shock(const string& var,
const vector<AbstractShocksStatement::period_range_t>& periods,
const vector<expr_t>& values, DetShockType type);
//! Adds a heteroskedastic shock (either values or scales)
void add_heteroskedastic_shock(const string& var, const vector<pair<int, int>>& periods,
void add_heteroskedastic_shock(const string& var,
const vector<AbstractShocksStatement::period_range_t>& periods,
const vector<expr_t>& values, bool scales);
//! Adds a std error shock
void add_stderr_shock(const string& var, expr_t value);
......@@ -562,8 +562,8 @@ public:
void osr_params_bounds();
//! Add a line in an osr params block
void add_osr_params_element();
//! Sets the frequency of the data
void set_time(const string& arg);
// Sets the initial period for estimation
void set_time(string period);
//! Estimation Data
void estimation_data();
//! Sets the prior for a parameter
......@@ -667,18 +667,9 @@ public:
void end_planner_objective(expr_t expr);
//! Ramsey model statement
void ramsey_model();
//! Ramsey constraints statement
void add_ramsey_constraints_statement();
//! Ramsey less constraint
void ramsey_constraint_add_less(const string& name, const expr_t rhs);
//! Ramsey greater constraint
void ramsey_constraint_add_greater(const string& name, const expr_t rhs);
//! Ramsey less or equal constraint
void ramsey_constraint_add_less_equal(const string& name, const expr_t rhs);
//! Ramsey greater or equal constraint
void ramsey_constraint_add_greater_equal(const string& name, const expr_t rhs);
//! Ramsey constraint helper function
void add_ramsey_constraint(const string& name, BinaryOpcode op_code, const expr_t rhs);
// Ramsey constraints statement
void begin_ramsey_constraints();
void end_ramsey_constraints(const vector<expr_t>& constraints);
//! Ramsey policy statement
void ramsey_policy(vector<string> symbol_list);
//! Evaluate Planner Objective
......@@ -855,6 +846,8 @@ public:
expr_t add_erfc(expr_t arg);
//! Writes token "steadyState(arg1)" to model tree
expr_t add_steady_state(expr_t arg1);
// Add a “sum(arg)” node to model tree
expr_t add_sum(expr_t arg);
//! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)
void push_external_function_arg_vector_onto_stack();
//! Adds an external function argument
......@@ -928,6 +921,10 @@ public:
void perfect_foresight_solver();
void perfect_foresight_with_expectation_errors_setup();
void perfect_foresight_with_expectation_errors_solver();
void perfect_foresight_controlled_paths(
const vector<tuple<string, vector<AbstractShocksStatement::period_range_t>, vector<expr_t>,
string>>& paths,
variant<int, string> learnt_in_period);
void prior_posterior_function(bool prior_func);
//! Method of Moments estimation statement
void method_of_moments();
......@@ -963,6 +960,8 @@ public:
// Add a matched_irfs_weights block
void matched_irfs_weights(MatchedIrfsWeightsStatement::matched_irfs_weights_t weights,
bool overwrite);
void heterogeneity_dimension(const vector<string>& dims);
// Returns true iff the string is a legal symbol identifier (see NAME token in lexer)
static bool isSymbolIdentifier(const string& str);
// Given an Occbin regime name, returns the corresponding auxiliary parameter
......
src/Shocks.cc
View file @ 6618d7d1
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -19,11 +19,43 @@
#include <cassert>
#include <cstdlib>
#include <functional>
#include <iostream>
#include <utility>
#include "Shocks.hh"
static auto print_matlab_period_range = []<class T>(ostream& output, const T& arg) {
if constexpr (is_same_v<T, pair<int, int>>)
output << arg.first << ":" << arg.second;
else if constexpr (is_same_v<T, pair<string, string>>)
output << "(" << arg.first << "):(" << arg.second << ")";
else
static_assert(always_false_v<T>, "Non-exhaustive visitor!");
};
static auto print_json_period_range = []<class T>(ostream& output, const T& arg) {
if constexpr (is_same_v<T, pair<int, int>>)
output << R"("period1": )" << arg.first << ", "
<< R"("period2": )" << arg.second;
else if constexpr (is_same_v<T, pair<string, string>>)
output << R"("period1": ")" << arg.first << R"(", )"
<< R"("period2": ")" << arg.second << '"';
else
static_assert(always_false_v<T>, "Non-exhaustive visitor!");
};
static auto print_matlab_learnt_in = [](ostream& output, const auto& p) { output << p; };
static auto print_json_learnt_in = []<class T>(ostream& output, const T& p) {
if constexpr (is_same_v<T, int>)
output << p;
else if constexpr (is_same_v<T, string>)
output << '"' << p << '"';
else
static_assert(always_false_v<T>, "Non-exhaustive visitor!");
};
AbstractShocksStatement::AbstractShocksStatement(bool overwrite_arg, ShockType type_arg,
det_shocks_t det_shocks_arg,
const SymbolTable& symbol_table_arg) :
......@@ -37,23 +69,19 @@ AbstractShocksStatement::AbstractShocksStatement(bool overwrite_arg, ShockType t
void
AbstractShocksStatement::writeDetShocks(ostream& output) const
{
int exo_det_length = 0;
for (const auto& [id, shock_vec] : det_shocks)
for (bool exo_det = (symbol_table.getType(id) == SymbolType::exogenousDet);
const auto& [period1, period2, value] : shock_vec)
const auto& [period_range, value] : shock_vec)
{
output << "M_.det_shocks = [ M_.det_shocks;" << endl
<< boolalpha << "struct('exo_det'," << exo_det << ",'exo_id',"
<< symbol_table.getTypeSpecificID(id) + 1 << ",'type','" << typeToString(type) << "'"
<< ",'periods'," << period1 << ":" << period2 << ",'value',";
<< ",'periods',";
visit(bind(print_matlab_period_range, ref(output), placeholders::_1), period_range);
output << ",'value',";
value->writeOutput(output);
output << ") ];" << endl;
if (exo_det && period2 > exo_det_length)
exo_det_length = period2;
}
output << "M_.exo_det_length = " << exo_det_length << ";\n";
}
void
......@@ -66,13 +94,13 @@ AbstractShocksStatement::writeJsonDetShocks(ostream& output) const
output << ", ";
output << R"({"var": ")" << symbol_table.getName(id) << R"(", )"
<< R"("values": [)";
for (bool printed_something2 {false}; const auto& [period1, period2, value] : shock_vec)
for (bool printed_something2 {false}; const auto& [period_range, value] : shock_vec)
{
if (exchange(printed_something2, true))
output << ", ";
output << R"({"period1": )" << period1 << ", "
<< R"("period2": )" << period2 << ", "
<< R"("value": ")";
output << "{";
visit(bind(print_json_period_range, ref(output), placeholders::_1), period_range);
output << R"(, "value": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
......@@ -467,8 +495,9 @@ ShocksSurpriseStatement::writeOutput(ostream& output, [[maybe_unused]] const str
else
output << "M_.surprise_shocks = [ M_.surprise_shocks;" << endl;
for (const auto& [id, shock_vec] : surprise_shocks)
for (const auto& [period1, period2, value] : shock_vec)
for (const auto& [period_range, value] : shock_vec)
{
auto [period1, period2] = get<pair<int, int>>(period_range);
output << "struct('exo_id'," << symbol_table.getTypeSpecificID(id) + 1 << ",'periods',"
<< period1 << ":" << period2 << ",'value',";
value->writeOutput(output);
......@@ -489,8 +518,9 @@ ShocksSurpriseStatement::writeJsonOutput(ostream& output) const
output << ", ";
output << R"({"var": ")" << symbol_table.getName(id) << R"(", )"
<< R"("values": [)";
for (bool printed_something2 {false}; const auto& [period1, period2, value] : shock_vec)
for (bool printed_something2 {false}; const auto& [period_range, value] : shock_vec)
{
auto [period1, period2] = get<pair<int, int>>(period_range);
if (exchange(printed_something2, true))
output << ", ";
output << R"({"period1": )" << period1 << ", "
......@@ -504,10 +534,11 @@ ShocksSurpriseStatement::writeJsonOutput(ostream& output) const
output << "]}";
}
ShocksLearntInStatement::ShocksLearntInStatement(int learnt_in_period_arg, bool overwrite_arg,
ShocksLearntInStatement::ShocksLearntInStatement(variant<int, string> learnt_in_period_arg,
bool overwrite_arg,
learnt_shocks_t learnt_shocks_arg,
const SymbolTable& symbol_table_arg) :
learnt_in_period {learnt_in_period_arg},
learnt_in_period {move(learnt_in_period_arg)},
overwrite {overwrite_arg},
learnt_shocks {move(learnt_shocks_arg)},
symbol_table {symbol_table_arg}
......@@ -545,18 +576,29 @@ ShocksLearntInStatement::writeOutput(ostream& output, [[maybe_unused]] const str
[[maybe_unused]] bool minimal_workspace) const
{
if (overwrite)
{
output << "if ~isempty(M_.learnt_shocks)" << endl
<< " M_.learnt_shocks = M_.learnt_shocks([M_.learnt_shocks.learnt_in] ~= "
<< learnt_in_period << ");" << endl
<< "end" << endl;
<< " M_.learnt_shocks = M_.learnt_shocks(cellfun(@(x) ~isa(x, '";
if (holds_alternative<int>(learnt_in_period))
output << "numeric";
else
output << "dates";
output << "') || x ~= ";
/* NB: date expression not parenthesized since it can only contain a + operator, which has
higher precedence than ~= and || */
visit(bind(print_matlab_learnt_in, ref(output), placeholders::_1), learnt_in_period);
output << ", {M_.learnt_shocks.learnt_in}));" << endl << "end" << endl;
}
output << "M_.learnt_shocks = [ M_.learnt_shocks;" << endl;
for (const auto& [id, shock_vec] : learnt_shocks)
for (const auto& [type, period1, period2, value] : shock_vec)
for (const auto& [type, period_range, value] : shock_vec)
{
output << "struct('learnt_in'," << learnt_in_period << ",'exo_id',"
<< symbol_table.getTypeSpecificID(id) + 1 << ",'periods'," << period1 << ":"
<< period2 << ",'type','" << typeToString(type) << "'"
output << "struct('learnt_in',";
visit(bind(print_matlab_learnt_in, ref(output), placeholders::_1), learnt_in_period);
output << ",'exo_id'," << symbol_table.getTypeSpecificID(id) + 1 << ",'periods',";
visit(bind(print_matlab_period_range, ref(output), placeholders::_1), period_range);
output << ",'type','" << typeToString(type) << "'"
<< ",'value',";
value->writeOutput(output);
output << ");" << endl;
......@@ -568,21 +610,22 @@ void
ShocksLearntInStatement::writeJsonOutput(ostream& output) const
{
output << R"({"statementName": "shocks")"
<< R"(, "learnt_in": )" << learnt_in_period << R"(, "overwrite": )" << boolalpha
<< overwrite << R"(, "learnt_shocks": [)";
<< R"(, "learnt_in": )";
visit(bind(print_json_learnt_in, ref(output), placeholders::_1), learnt_in_period);
output << R"(, "overwrite": )" << boolalpha << overwrite << R"(, "learnt_shocks": [)";
for (bool printed_something {false}; const auto& [id, shock_vec] : learnt_shocks)
{
if (exchange(printed_something, true))
output << ", ";
output << R"({"var": ")" << symbol_table.getName(id) << R"(", )"
<< R"("values": [)";
for (bool printed_something2 {false}; const auto& [type, period1, period2, value] : shock_vec)
for (bool printed_something2 {false}; const auto& [type, period_range, value] : shock_vec)
{
if (exchange(printed_something2, true))
output << ", ";
output << R"({"period1": )" << period1 << ", "
<< R"("period2": )" << period2 << ", "
<< R"("type": ")" << typeToString(type) << R"(", )"
output << "{";
visit(bind(print_json_period_range, ref(output), placeholders::_1), period_range);
output << R"(, "type": ")" << typeToString(type) << R"(", )"
<< R"("value": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
......@@ -592,6 +635,213 @@ ShocksLearntInStatement::writeJsonOutput(ostream& output) const
output << "]}";
}
HeterogeneousShocksStatement::HeterogeneousShocksStatement(
int heterogeneity_dimension_arg, bool overwrite_arg, var_and_std_shocks_t var_shocks_arg,
var_and_std_shocks_t std_shocks_arg, covar_and_corr_shocks_t covar_shocks_arg,
covar_and_corr_shocks_t corr_shocks_arg, const SymbolTable& symbol_table_arg,
const HeterogeneityTable& heterogeneity_table_arg) :
heterogeneity_dimension {heterogeneity_dimension_arg},
overwrite {overwrite_arg},
var_shocks {move(var_shocks_arg)},
std_shocks {move(std_shocks_arg)},
covar_shocks {move(covar_shocks_arg)},
corr_shocks {move(corr_shocks_arg)},
symbol_table {symbol_table_arg},
heterogeneity_table {heterogeneity_table_arg}
{
}
void
HeterogeneousShocksStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
[[maybe_unused]] bool minimal_workspace) const
{
if (overwrite)
output << sigmaeName() << " = zeros(" << symbol_table.het_exo_nbr(heterogeneity_dimension)
<< ", " << symbol_table.het_exo_nbr(heterogeneity_dimension) << ");" << endl;
writeVarAndStdShocks(output);
writeCovarAndCorrShocks(output);
}
void
HeterogeneousShocksStatement::writeJsonOutput(ostream& output) const
{
output << R"({"statementName": "shocks")"
<< R"(, "heterogeneity": ")" << heterogeneity_table.getName(heterogeneity_dimension)
<< R"(", "overwrite": )" << boolalpha << overwrite;
output << R"(, "variance": [)";
for (bool printed_something {false}; auto& [id, value] : var_shocks)
{
if (exchange(printed_something, true))
output << ", ";
output << R"({"name": ")" << symbol_table.getName(id) << R"(", )"
<< R"("variance": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
output << "]"
<< R"(, "stderr": [)";
for (bool printed_something {false}; auto& [id, value] : std_shocks)
{
if (exchange(printed_something, true))
output << ", ";
output << R"({"name": ")" << symbol_table.getName(id) << R"(", )"
<< R"("stderr": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
output << "]"
<< R"(, "covariance": [)";
for (bool printed_something {false}; auto& [ids, value] : covar_shocks)
{
if (exchange(printed_something, true))
output << ", ";
output << "{"
<< R"("name": ")" << symbol_table.getName(ids.first) << R"(", )"
<< R"("name2": ")" << symbol_table.getName(ids.second) << R"(", )"
<< R"("covariance": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
output << "]"
<< R"(, "correlation": [)";
for (bool printed_something {false}; auto& [ids, value] : corr_shocks)
{
if (exchange(printed_something, true))
output << ", ";
output << "{"
<< R"("name": ")" << symbol_table.getName(ids.first) << R"(", )"
<< R"("name2": ")" << symbol_table.getName(ids.second) << R"(", )"
<< R"("correlation": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
output << "]"
<< "}";
}
void
HeterogeneousShocksStatement::writeVarOrStdShock(ostream& output, const pair<int, expr_t>& it,
bool stddev) const
{
SymbolType type = symbol_table.getType(it.first);
assert(type == SymbolType::heterogeneousExogenous);
int id {symbol_table.getTypeSpecificID(it.first) + 1};
output << sigmaeName() << "(" << id << ", " << id << ") = ";
if (stddev)
output << "(";
it.second->writeOutput(output);
if (stddev)
output << ")^2";
output << ";" << endl;
}
void
HeterogeneousShocksStatement::writeVarAndStdShocks(ostream& output) const
{
for (const auto& it : var_shocks)
writeVarOrStdShock(output, it, false);
for (const auto& it : std_shocks)
writeVarOrStdShock(output, it, true);
}
void
HeterogeneousShocksStatement::writeCovarOrCorrShock(ostream& output,
const pair<pair<int, int>, expr_t>& it,
bool corr) const
{
assert(symbol_table.getType(it.first.first) == SymbolType::heterogeneousExogenous
&& symbol_table.getType(it.first.second) == SymbolType::heterogeneousExogenous);
int id1 {symbol_table.getTypeSpecificID(it.first.first) + 1},
id2 {symbol_table.getTypeSpecificID(it.first.second) + 1};
output << sigmaeName() << "(" << id1 << ", " << id2 << ") = ";
it.second->writeOutput(output);
if (corr)
output << "*sqrt(" << sigmaeName() << "(" << id1 << ", " << id1 << ")*" << sigmaeName() << "("
<< id2 << ", " << id2 << "))";
output << ";" << endl
<< sigmaeName() << "(" << id2 << ", " << id1 << ") = " << sigmaeName() << "(" << id1
<< ", " << id2 << ");" << endl;
}
void
HeterogeneousShocksStatement::writeCovarAndCorrShocks(ostream& output) const
{
for (const auto& it : covar_shocks)
writeCovarOrCorrShock(output, it, false);
for (const auto& it : corr_shocks)
writeCovarOrCorrShock(output, it, true);
}
void
HeterogeneousShocksStatement::checkPass(ModFileStructure& mod_file_struct,
[[maybe_unused]] WarningConsolidation& warnings)
{
/* Error out if variables are not of the right type. This must be done here
and not at parsing time (see #448). */
for (auto [id, val] : var_shocks)
if (symbol_table.getType(id) != SymbolType::heterogeneousExogenous)
{
cerr << "shocks: setting a variance on '" << symbol_table.getName(id)
<< "' is not allowed, because it is not a heterogeneous exogenous variable" << endl;
exit(EXIT_FAILURE);
}
for (auto [id, val] : std_shocks)
if (symbol_table.getType(id) != SymbolType::heterogeneousExogenous)
{
cerr << "shocks: setting a standard error on '" << symbol_table.getName(id)
<< "' is not allowed, because it is not a heterogeneous exogenous variable" << endl;
exit(EXIT_FAILURE);
}
for (const auto& [ids, val] : covar_shocks)
{
auto& [symb_id1, symb_id2] = ids;
if (!(symbol_table.getType(symb_id1) == SymbolType::heterogeneousExogenous
&& symbol_table.getType(symb_id2) == SymbolType::heterogeneousExogenous))
{
cerr << "shocks: setting a covariance between '" << symbol_table.getName(symb_id1)
<< "' and '" << symbol_table.getName(symb_id2)
<< "'is not allowed; covariances can only be specified for heterogeneous exogenous "
"variables"
<< endl;
exit(EXIT_FAILURE);
}
}
for (const auto& [ids, val] : corr_shocks)
{
auto& [symb_id1, symb_id2] = ids;
if (!(symbol_table.getType(symb_id1) == SymbolType::heterogeneousExogenous
&& symbol_table.getType(symb_id2) == SymbolType::heterogeneousExogenous))
{
cerr << "shocks: setting a correlation between '" << symbol_table.getName(symb_id1)
<< "' and '" << symbol_table.getName(symb_id2)
<< "'is not allowed; covariances can only be specified for heterogeneous exogenous "
"variables"
<< endl;
exit(EXIT_FAILURE);
}
}
// Fill in mod_file_struct.parameters_with_shocks_values (related to #469)
for (auto [id, val] : var_shocks)
val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
for (auto [id, val] : std_shocks)
val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
for (const auto& [ids, val] : covar_shocks)
val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
for (const auto& [ids, val] : corr_shocks)
val->collectVariables(SymbolType::parameter, mod_file_struct.parameters_within_shocks_values);
}
ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(
AbstractShocksStatement::det_shocks_t paths_arg, const SymbolTable& symbol_table_arg) :
paths {move(paths_arg)}, symbol_table {symbol_table_arg}, path_length {computePathLength(paths)}
......@@ -604,9 +854,12 @@ ConditionalForecastPathsStatement::computePathLength(
{
int length {0};
for (const auto& [ignore, elems] : paths)
for (auto& [period1, period2, value] : elems)
for (auto& [period_range, value] : elems)
{
auto [period1, period2] = get<pair<int, int>>(period_range);
// Period1 < Period2, as enforced in ParsingDriver::add_period()
length = max(length, period2);
}
return length;
}
......@@ -626,10 +879,10 @@ ConditionalForecastPathsStatement::writeOutput(ostream& output,
else
output << "constrained_vars_ = [constrained_vars_; "
<< symbol_table.getTypeSpecificID(id) + 1 << "];" << endl;
for (const auto& [period1, period2, value] : elems)
for (int j = period1; j <= period2; j++)
for (const auto& [period_range, value] : elems)
{
output << "constrained_paths_(" << k << "," << j << ")=";
auto [period1, period2] = get<pair<int, int>>(period_range);
output << "constrained_paths_(" << k << "," << period1 << ":" << period2 << ")=";
value->writeOutput(output);
output << ";" << endl;
}
......@@ -648,13 +901,13 @@ ConditionalForecastPathsStatement::writeJsonOutput(ostream& output) const
output << ", ";
output << R"({"var": ")" << symbol_table.getName(id) << R"(", )"
<< R"("values": [)";
for (bool printed_something2 {false}; const auto& [period1, period2, value] : elems)
for (bool printed_something2 {false}; const auto& [period_range, value] : elems)
{
if (exchange(printed_something2, true))
output << ", ";
output << R"({"period1": )" << period1 << ", "
<< R"("period2": )" << period2 << ", "
<< R"("value": ")";
output << "{";
visit(bind(print_json_period_range, ref(output), placeholders::_1), period_range);
output << R"(, "value": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
......@@ -663,6 +916,65 @@ ConditionalForecastPathsStatement::writeJsonOutput(ostream& output) const
output << "]}";
}
PerfectForesightControlledPathsStatement::PerfectForesightControlledPathsStatement(
paths_t paths_arg, variant<int, string> learnt_in_period_arg,
const SymbolTable& symbol_table_arg) :
paths {move(paths_arg)},
learnt_in_period {move(learnt_in_period_arg)},
symbol_table {symbol_table_arg}
{
}
void
PerfectForesightControlledPathsStatement::writeOutput(ostream& output,
[[maybe_unused]] const string& basename,
[[maybe_unused]] bool minimal_workspace) const
{
for (const auto& [exogenize_id, constraints, endogenize_id] : paths)
for (const auto& [period_range, value] : constraints)
{
output << "M_.perfect_foresight_controlled_paths = [ "
"M_.perfect_foresight_controlled_paths;"
<< endl
<< "struct('exogenize_id'," << symbol_table.getTypeSpecificID(exogenize_id) + 1
<< ",'periods',";
visit(bind(print_matlab_period_range, ref(output), placeholders::_1), period_range);
output << ",'value',";
value->writeOutput(output);
output << ",'endogenize_id'," << symbol_table.getTypeSpecificID(endogenize_id) + 1
<< ",'learnt_in',";
visit(bind(print_matlab_learnt_in, ref(output), placeholders::_1), learnt_in_period);
output << ") ];" << endl;
}
}
void
PerfectForesightControlledPathsStatement::writeJsonOutput(ostream& output) const
{
output << R"({"statementName": "perfect_foresight_controlled_paths")"
<< R"(, "paths": [)";
for (bool printed_something {false};
const auto& [exogenize_id, constraints, endogenize_id] : paths)
{
if (exchange(printed_something, true))
output << ", ";
output << R"({"exogenize": ")" << symbol_table.getName(exogenize_id) << R"(", )"
<< R"("values": [)";
for (bool printed_something2 {false}; const auto& [period_range, value] : constraints)
{
if (exchange(printed_something2, true))
output << ", ";
output << "{";
visit(bind(print_json_period_range, ref(output), placeholders::_1), period_range);
output << R"(, "value": ")";
value->writeJsonOutput(output, {}, {});
output << R"("})";
}
output << R"(], "endogenize": ")" << symbol_table.getName(endogenize_id) << R"("})";
}
output << "]}";
}
MomentCalibration::MomentCalibration(constraints_t constraints_arg,
const SymbolTable& symbol_table_arg) :
constraints {move(constraints_arg)}, symbol_table {symbol_table_arg}
......
src/Shocks.hh
View file @ 6618d7d1
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2025 Dynare Team
*
* This file is part of Dynare.
*
......@@ -22,9 +22,11 @@
#include <map>
#include <string>
#include <variant>
#include <vector>
#include "ExprNode.hh"
#include "HeterogeneityTable.hh"
#include "Statement.hh"
#include "SymbolTable.hh"
......@@ -33,8 +35,10 @@ using namespace std;
class AbstractShocksStatement : public Statement
{
public:
// The tuple is (period1, period2, value)
using det_shocks_t = map<int, vector<tuple<int, int, expr_t>>>;
// A period range is either two indices (1-based), or two dates (from dseries)
using period_range_t = variant<pair<int, int>, pair<string, string>>;
// The pair is (period range, value)
using det_shocks_t = map<int, vector<pair<period_range_t, expr_t>>>;
enum class ShockType
{
level, // The value is the level of the exogenous (“values” statement in “shocks”)
......@@ -122,7 +126,7 @@ public:
class ShocksLearntInStatement : public Statement
{
public:
const int learnt_in_period;
const variant<int, string> learnt_in_period;
//! Does this “shocks(learnt_in=…)” or “mshocks(learnt_in=…)” block replace the previous ones?
const bool overwrite;
enum class LearntShockType
......@@ -139,8 +143,9 @@ public:
// state as anticipated in the same informational period (“values”
// statement in “mshocks(learnt_in=…, relative_to_initval)”)
};
// The tuple is (type, period1, period2, value)
using learnt_shocks_t = map<int, vector<tuple<LearntShockType, int, int, expr_t>>>;
// The tuple is (type, period range, value)
using learnt_shocks_t
= map<int, vector<tuple<LearntShockType, AbstractShocksStatement::period_range_t, expr_t>>>;
const learnt_shocks_t learnt_shocks;
private:
......@@ -148,13 +153,54 @@ private:
static string typeToString(LearntShockType type);
public:
ShocksLearntInStatement(int learnt_in_period_arg, bool overwrite_arg,
ShocksLearntInStatement(variant<int, string> learnt_in_period_arg, bool overwrite_arg,
learnt_shocks_t learnt_shocks_arg, const SymbolTable& symbol_table_arg);
void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
void writeJsonOutput(ostream& output) const override;
};
class HeterogeneousShocksStatement : public Statement
{
public:
const int heterogeneity_dimension;
const bool overwrite;
using var_and_std_shocks_t = map<int, expr_t>;
using covar_and_corr_shocks_t = map<pair<int, int>, expr_t>;
const var_and_std_shocks_t var_shocks, std_shocks;
const covar_and_corr_shocks_t covar_shocks, corr_shocks;
private:
const SymbolTable& symbol_table;
const HeterogeneityTable& heterogeneity_table;
void writeVarOrStdShock(ostream& output, const pair<int, expr_t>& it, bool stddev) const;
void writeVarAndStdShocks(ostream& output) const;
void writeCovarOrCorrShock(ostream& output, const pair<pair<int, int>, expr_t>& it,
bool corr) const;
void writeCovarAndCorrShocks(ostream& output) const;
string
sigmaeName() const
{
return "M_.heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").Sigma_e"s;
}
public:
HeterogeneousShocksStatement(int heterogeneity_dimension_arg, bool overwrite_arg,
var_and_std_shocks_t var_shocks_arg,
var_and_std_shocks_t std_shocks_arg,
covar_and_corr_shocks_t covar_shocks_arg,
covar_and_corr_shocks_t corr_shocks_arg,
const SymbolTable& symbol_table_arg,
const HeterogeneityTable& heterogeneity_table_arg);
void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
void writeJsonOutput(ostream& output) const override;
};
class ConditionalForecastPathsStatement : public Statement
{
private:
......@@ -170,6 +216,26 @@ public:
static int computePathLength(const AbstractShocksStatement::det_shocks_t& paths);
};
class PerfectForesightControlledPathsStatement : public Statement
{
public:
// (exogenize_id, vector of (period range, value), endogenize_id)
using paths_t
= vector<tuple<int, vector<pair<AbstractShocksStatement::period_range_t, expr_t>>, int>>;
private:
const paths_t paths;
const variant<int, string> learnt_in_period;
const SymbolTable& symbol_table;
public:
PerfectForesightControlledPathsStatement(paths_t paths_arg,
variant<int, string> learnt_in_period_arg,
const SymbolTable& symbol_table_arg);
void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
void writeJsonOutput(ostream& output) const override;
};
class MomentCalibration : public Statement
{
public:
......