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

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
Show changes
Commits on Source (10)
Hide whitespace changes
Inline Side-by-side
Showing
with 431 additions and 237 deletions
src/ComputingTasks.cc
View file @ 253c1866
...@@ -508,90 +508,6 @@ RamseyModelStatement::writeJsonOutput(ostream& output) const ...@@ -508,90 +508,6 @@ RamseyModelStatement::writeJsonOutput(ostream& output) const
output << "}"; output << "}";
} }
RamseyConstraintsStatement::RamseyConstraintsStatement(const SymbolTable& symbol_table_arg,
constraints_t constraints_arg) :
symbol_table {symbol_table_arg}, constraints {move(constraints_arg)}
{
}
void
RamseyConstraintsStatement::checkPass(ModFileStructure& mod_file_struct,
[[maybe_unused]] WarningConsolidation& warnings)
{
mod_file_struct.ramsey_constraints_present = true;
}
void
RamseyConstraintsStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
[[maybe_unused]] bool minimal_workspace) const
{
output << "M_.ramsey_model_constraints = {" << endl;
for (bool printed_something {false}; const auto& it : constraints)
{
if (exchange(printed_something, true))
output << ", ";
output << "{" << it.endo + 1 << ", '";
switch (it.code)
{
case BinaryOpcode::less:
output << '<';
break;
case BinaryOpcode::greater:
output << '>';
break;
case BinaryOpcode::lessEqual:
output << "<=";
break;
case BinaryOpcode::greaterEqual:
output << ">=";
break;
default:
cerr << "Ramsey constraints: this shouldn't happen." << endl;
exit(EXIT_FAILURE);
}
output << "', '";
it.expression->writeOutput(output);
output << "'}" << endl;
}
output << "};" << endl;
}
void
RamseyConstraintsStatement::writeJsonOutput(ostream& output) const
{
output << R"({"statementName": "ramsey_constraints")"
<< R"(, "ramsey_model_constraints": [)" << endl;
for (bool printed_something {false}; const auto& it : constraints)
{
if (exchange(printed_something, true))
output << ", ";
output << R"({"constraint": ")" << symbol_table.getName(it.endo) << " ";
switch (it.code)
{
case BinaryOpcode::less:
output << '<';
break;
case BinaryOpcode::greater:
output << '>';
break;
case BinaryOpcode::lessEqual:
output << "<=";
break;
case BinaryOpcode::greaterEqual:
output << ">=";
break;
default:
cerr << "Ramsey constraints: this shouldn't happen." << endl;
exit(EXIT_FAILURE);
}
output << " ";
it.expression->writeJsonOutput(output, {}, {});
output << R"("})" << endl;
}
output << "]" << endl;
output << "}";
}
RamseyPolicyStatement::RamseyPolicyStatement(SymbolList symbol_list_arg, RamseyPolicyStatement::RamseyPolicyStatement(SymbolList symbol_list_arg,
OptionsList options_list_arg, OptionsList options_list_arg,
const SymbolTable& symbol_table_arg) : const SymbolTable& symbol_table_arg) :
...@@ -2014,7 +1930,7 @@ void ...@@ -2014,7 +1930,7 @@ void
OsrParamsStatement::checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) OsrParamsStatement::checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings)
{ {
if (mod_file_struct.osr_params_present) if (mod_file_struct.osr_params_present)
cerr << "WARNING: You have more than one osr_params statement in the .mod file." << endl; warnings << "WARNING: You have more than one osr_params statement in the .mod file." << endl;
mod_file_struct.osr_params_present = true; mod_file_struct.osr_params_present = true;
try try
......
src/ComputingTasks.hh
View file @ 253c1866
/* /*
* Copyright © 2003-2023 Dynare Team * Copyright © 2003-2024 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -179,28 +179,6 @@ public: ...@@ -179,28 +179,6 @@ public:
void writeJsonOutput(ostream& output) const override; void writeJsonOutput(ostream& output) const override;
}; };
class RamseyConstraintsStatement : public Statement
{
public:
struct Constraint
{
int endo;
BinaryOpcode code;
expr_t expression;
};
using constraints_t = vector<Constraint>;
private:
const SymbolTable& symbol_table;
const constraints_t constraints;
public:
RamseyConstraintsStatement(const SymbolTable& symbol_table_arg, constraints_t constraints_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 RamseyPolicyStatement : public Statement class RamseyPolicyStatement : public Statement
{ {
private: private:
......
src/DynamicModel.cc
View file @ 253c1866
...@@ -34,10 +34,18 @@ ...@@ -34,10 +34,18 @@
void void
DynamicModel::copyHelper(const DynamicModel& m) DynamicModel::copyHelper(const DynamicModel& m)
{ {
auto f = [this](const ExprNode* e) { return e->clone(*this); }; auto f = [this](const ExprNode* e) { return e ? e->clone(*this) : nullptr; };
for (const auto& it : m.static_only_equations) for (const auto& it : m.static_only_equations)
static_only_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it))); static_only_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
for (const auto& it : m.static_only_complementarity_conditions)
if (it)
{
const auto& [symb_id, lb, ub] = *it;
static_only_complementarity_conditions.emplace_back(in_place, symb_id, f(lb), f(ub));
}
else
static_only_complementarity_conditions.emplace_back(nullopt);
} }
DynamicModel::DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg, DynamicModel::DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
...@@ -104,6 +112,7 @@ DynamicModel::operator=(const DynamicModel& m) ...@@ -104,6 +112,7 @@ DynamicModel::operator=(const DynamicModel& m)
static_only_equations_lineno = m.static_only_equations_lineno; static_only_equations_lineno = m.static_only_equations_lineno;
static_only_equations_equation_tags = m.static_only_equations_equation_tags; static_only_equations_equation_tags = m.static_only_equations_equation_tags;
static_only_complementarity_conditions.clear();
deriv_id_table = m.deriv_id_table; deriv_id_table = m.deriv_id_table;
inv_deriv_id_table = m.inv_deriv_id_table; inv_deriv_id_table = m.inv_deriv_id_table;
dyn_jacobian_cols_table = m.dyn_jacobian_cols_table; dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
...@@ -651,11 +660,11 @@ DynamicModel::parseIncludeExcludeEquations(const string& inc_exc_option_value, b ...@@ -651,11 +660,11 @@ DynamicModel::parseIncludeExcludeEquations(const string& inc_exc_option_value, b
} }
vector<int> vector<int>
DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, DynamicModel::removeEquationsHelper(
bool excluded_vars_change_type, set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
vector<BinaryOpNode*>& all_equations, vector<BinaryOpNode*>& all_equations, vector<optional<int>>& all_equations_lineno,
vector<optional<int>>& all_equations_lineno, vector<optional<tuple<int, expr_t, expr_t>>>& all_complementarity_conditions,
EquationTags& all_equation_tags, bool static_equations) const EquationTags& all_equation_tags, bool static_equations) const
{ {
if (all_equations.empty()) if (all_equations.empty())
return {}; return {};
...@@ -686,6 +695,7 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, ...@@ -686,6 +695,7 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
// remove from equations, equations_lineno, equation_tags // remove from equations, equations_lineno, equation_tags
vector<BinaryOpNode*> new_equations; vector<BinaryOpNode*> new_equations;
vector<optional<int>> new_equations_lineno; vector<optional<int>> new_equations_lineno;
vector<optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
map<int, int> old_eqn_num_2_new; map<int, int> old_eqn_num_2_new;
vector<int> excluded_vars; vector<int> excluded_vars;
for (size_t i = 0; i < all_equations.size(); i++) for (size_t i = 0; i < all_equations.size(); i++)
...@@ -717,11 +727,13 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, ...@@ -717,11 +727,13 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
new_equations.emplace_back(all_equations[i]); new_equations.emplace_back(all_equations[i]);
old_eqn_num_2_new[i] = new_equations.size() - 1; old_eqn_num_2_new[i] = new_equations.size() - 1;
new_equations_lineno.emplace_back(all_equations_lineno[i]); new_equations_lineno.emplace_back(all_equations_lineno[i]);
new_complementarity_conditions.emplace_back(all_complementarity_conditions[i]);
} }
int n_excl = all_equations.size() - new_equations.size(); int n_excl = all_equations.size() - new_equations.size();
all_equations = new_equations; all_equations = new_equations;
all_equations_lineno = new_equations_lineno; all_equations_lineno = new_equations_lineno;
all_complementarity_conditions = new_complementarity_conditions;
all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new); all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new);
...@@ -754,12 +766,13 @@ DynamicModel::removeEquations(const vector<map<string, string>>& listed_eqs_by_t ...@@ -754,12 +766,13 @@ DynamicModel::removeEquations(const vector<map<string, string>>& listed_eqs_by_t
vector<int> excluded_vars vector<int> excluded_vars
= removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type, equations, = removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type, equations,
equations_lineno, equation_tags, false); equations_lineno, complementarity_conditions, equation_tags, false);
// Ignore output because variables are not excluded when equations marked 'static' are excluded // Ignore output because variables are not excluded when equations marked 'static' are excluded
removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type, removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type,
static_only_equations, static_only_equations_lineno, static_only_equations, static_only_equations_lineno,
static_only_equations_equation_tags, true); static_only_complementarity_conditions, static_only_equations_equation_tags,
true);
if (!listed_eqs_by_tag2.empty()) if (!listed_eqs_by_tag2.empty())
{ {
...@@ -1134,6 +1147,11 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const ...@@ -1134,6 +1147,11 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
output << "'; " << endl; output << "'; " << endl;
} }
output << "};" << endl; output << "};" << endl;
output << "M_.dynamic_mcp_equations_reordering = [";
for (auto i : mcp_equations_reordering)
output << i + 1 << "; ";
output << "];" << endl;
} }
void void
...@@ -2507,6 +2525,8 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder, ...@@ -2507,6 +2525,8 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder,
cerr << "ERROR: Block decomposition requested but failed." << endl; cerr << "ERROR: Block decomposition requested but failed." << endl;
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
} }
computeMCPEquationsReordering();
} }
void void
...@@ -2810,6 +2830,8 @@ DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const strin ...@@ -2810,6 +2830,8 @@ DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const strin
writeSetAuxiliaryVariablesFile<true>(basename, julia); writeSetAuxiliaryVariablesFile<true>(basename, julia);
writeComplementarityConditionsFile<true>(basename);
// Support for model debugging // Support for model debugging
if (!julia) if (!julia)
writeDebugModelMFiles<true>(basename); writeDebugModelMFiles<true>(basename);
...@@ -2821,6 +2843,7 @@ DynamicModel::clearEquations() ...@@ -2821,6 +2843,7 @@ DynamicModel::clearEquations()
equations.clear(); equations.clear();
equations_lineno.clear(); equations_lineno.clear();
equation_tags.clear(); equation_tags.clear();
complementarity_conditions.clear();
} }
void void
...@@ -2828,14 +2851,26 @@ DynamicModel::replaceMyEquations(DynamicModel& dynamic_model) const ...@@ -2828,14 +2851,26 @@ DynamicModel::replaceMyEquations(DynamicModel& dynamic_model) const
{ {
dynamic_model.clearEquations(); dynamic_model.clearEquations();
auto clone_if_not_null = [&](expr_t e) { return e ? e->clone(dynamic_model) : nullptr; };
for (size_t i = 0; i < equations.size(); i++) for (size_t i = 0; i < equations.size(); i++)
dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i]); {
optional<tuple<int, expr_t, expr_t>> cc_clone;
if (complementarity_conditions[i])
{
auto& [symb_id, lower_bound, upper_bound] = *complementarity_conditions[i];
cc_clone = {symb_id, clone_if_not_null(lower_bound), clone_if_not_null(upper_bound)};
}
dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i],
move(cc_clone));
}
dynamic_model.equation_tags = equation_tags; dynamic_model.equation_tags = equation_tags;
} }
int int
DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model) DynamicModel::computeRamseyPolicyFOCs(const StaticModel& planner_objective,
map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints)
{ {
cout << "Ramsey Problem: added " << equations.size() << " multipliers." << endl; cout << "Ramsey Problem: added " << equations.size() << " multipliers." << endl;
...@@ -2849,8 +2884,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model) ...@@ -2849,8 +2884,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
} }
// Add Planner Objective to equations so that it appears in Lagrangian // Add Planner Objective to equations so that it appears in Lagrangian
assert(static_model.equations.size() == 1); assert(planner_objective.equations.size() == 1);
addEquation(static_model.equations[0]->clone(*this), nullopt); addEquation(planner_objective.equations[0]->clone(*this), nullopt);
// Get max endo lead and max endo lag // Get max endo lead and max endo lag
set<pair<int, int>> dynvars; set<pair<int, int>> dynvars;
...@@ -2892,9 +2927,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model) ...@@ -2892,9 +2927,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
lagrangian); lagrangian);
} }
// Save line numbers and tags, see below // Save line numbers, tags and complementarity conditions, see below
auto old_equations_lineno = equations_lineno; auto old_equations_lineno = equations_lineno;
auto old_equation_tags = equation_tags; auto old_equation_tags = equation_tags;
auto old_complementarity_conditions = complementarity_conditions;
// Prepare derivation of the Lagrangian // Prepare derivation of the Lagrangian
clearEquations(); clearEquations();
...@@ -2911,6 +2947,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model) ...@@ -2911,6 +2947,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
vector<expr_t> neweqs; vector<expr_t> neweqs;
vector<optional<int>> neweqs_lineno; vector<optional<int>> neweqs_lineno;
map<int, map<string, string>> neweqs_tags; map<int, map<string, string>> neweqs_tags;
map<int, optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
int orig_endo_nbr {0}; int orig_endo_nbr {0};
for (auto& [symb_id_and_lag, deriv_id] : deriv_id_table) for (auto& [symb_id_and_lag, deriv_id] : deriv_id_table)
{ {
...@@ -2923,12 +2960,20 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model) ...@@ -2923,12 +2960,20 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
{ {
// This is a derivative w.r.t. a Lagrange multiplier // This is a derivative w.r.t. a Lagrange multiplier
neweqs_lineno.push_back(old_equations_lineno[*i]); neweqs_lineno.push_back(old_equations_lineno[*i]);
neweqs_tags[neweqs.size() - 1] = old_equation_tags.getTagsByEqn(*i); neweqs_tags.emplace(neweqs.size() - 1, old_equation_tags.getTagsByEqn(*i));
new_complementarity_conditions.emplace(neweqs.size() - 1,
old_complementarity_conditions.at(*i));
} }
else else
{ {
orig_endo_nbr++; orig_endo_nbr++;
neweqs_lineno.emplace_back(nullopt); neweqs_lineno.emplace_back(nullopt);
if (cloned_ramsey_constraints.contains(symb_id))
{
auto& [lower_bound, upper_bound] = cloned_ramsey_constraints.at(symb_id);
new_complementarity_conditions.emplace(neweqs.size() - 1,
tuple {symb_id, lower_bound, upper_bound});
}
} }
} }
} }
...@@ -2936,7 +2981,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model) ...@@ -2936,7 +2981,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
// Overwrite equations with the Lagrangian derivatives // Overwrite equations with the Lagrangian derivatives
clearEquations(); clearEquations();
for (size_t i = 0; i < neweqs.size(); i++) for (size_t i = 0; i < neweqs.size(); i++)
addEquation(neweqs[i], neweqs_lineno[i], neweqs_tags[i]); addEquation(neweqs[i], neweqs_lineno[i], new_complementarity_conditions[i], neweqs_tags[i]);
return orig_endo_nbr; return orig_endo_nbr;
} }
...@@ -3773,6 +3818,7 @@ DynamicModel::fillEvalContext(eval_context_t& eval_context) const ...@@ -3773,6 +3818,7 @@ DynamicModel::fillEvalContext(eval_context_t& eval_context) const
void void
DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno, DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition,
map<string, string> eq_tags) map<string, string> eq_tags)
{ {
auto beq = dynamic_cast<BinaryOpNode*>(eq); auto beq = dynamic_cast<BinaryOpNode*>(eq);
...@@ -3781,6 +3827,7 @@ DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno, ...@@ -3781,6 +3827,7 @@ DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
static_only_equations_equation_tags.add(static_only_equations.size(), move(eq_tags)); static_only_equations_equation_tags.add(static_only_equations.size(), move(eq_tags));
static_only_equations.push_back(beq); static_only_equations.push_back(beq);
static_only_equations_lineno.push_back(lineno); static_only_equations_lineno.push_back(lineno);
static_only_complementarity_conditions.push_back(move(complementarity_condition));
} }
size_t size_t
...@@ -3834,7 +3881,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri ...@@ -3834,7 +3881,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
{ {
auto eq_tags_dynamic = eq_tags; auto eq_tags_dynamic = eq_tags;
eq_tags_dynamic["dynamic"] = ""; eq_tags_dynamic["dynamic"] = "";
addEquation(AddEqual(term, Zero), lineno, eq_tags_dynamic); addEquation(AddEqual(term, Zero), lineno, nullopt, eq_tags_dynamic);
} }
// Create or update the static equation (corresponding to the pure relax regime) // Create or update the static equation (corresponding to the pure relax regime)
...@@ -3856,7 +3903,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri ...@@ -3856,7 +3903,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
else else
{ {
eq_tags["static"] = ""; eq_tags["static"] = "";
addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, move(eq_tags)); addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, nullopt, move(eq_tags));
} }
} }
} }
...@@ -4171,8 +4218,8 @@ DynamicModel::simplifyEquations() ...@@ -4171,8 +4218,8 @@ DynamicModel::simplifyEquations()
{ {
size_t last_subst_table_size = 0; size_t last_subst_table_size = 0;
map<VariableNode*, NumConstNode*> subst_table; map<VariableNode*, NumConstNode*> subst_table;
// Equations with “mcp” tag are excluded, see dynare#1697 // Equations with a complementarity condition are excluded, see dynare#1697
findConstantEquationsWithoutMcpTag(subst_table); findConstantEquationsWithoutComplementarityCondition(subst_table);
while (subst_table.size() != last_subst_table_size) while (subst_table.size() != last_subst_table_size)
{ {
last_subst_table_size = subst_table.size(); last_subst_table_size = subst_table.size();
...@@ -4183,7 +4230,7 @@ DynamicModel::simplifyEquations() ...@@ -4183,7 +4230,7 @@ DynamicModel::simplifyEquations()
for (auto& equation : static_only_equations) for (auto& equation : static_only_equations)
equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table)); equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table));
subst_table.clear(); subst_table.clear();
findConstantEquationsWithoutMcpTag(subst_table); findConstantEquationsWithoutComplementarityCondition(subst_table);
} }
} }
......
src/DynamicModel.hh
View file @ 253c1866
...@@ -92,6 +92,9 @@ private: ...@@ -92,6 +92,9 @@ private:
//! Stores the equation tags of equations declared as [static] //! Stores the equation tags of equations declared as [static]
EquationTags static_only_equations_equation_tags; EquationTags static_only_equations_equation_tags;
// Complementarity conditions of equations declared as [static]
vector<optional<tuple<int, expr_t, expr_t>>> static_only_complementarity_conditions;
using deriv_id_table_t = map<pair<int, int>, int>; using deriv_id_table_t = map<pair<int, int>, int>;
//! Maps a pair (symbol_id, lag) to a deriv ID //! Maps a pair (symbol_id, lag) to a deriv ID
deriv_id_table_t deriv_id_table; deriv_id_table_t deriv_id_table;
...@@ -279,11 +282,11 @@ private: ...@@ -279,11 +282,11 @@ private:
Returns a list of excluded variables (empty if Returns a list of excluded variables (empty if
excluded_vars_change_type=false) */ excluded_vars_change_type=false) */
vector<int> removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, vector<int> removeEquationsHelper(
bool excluded_vars_change_type, set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
vector<BinaryOpNode*>& all_equations, vector<BinaryOpNode*>& all_equations, vector<optional<int>>& all_equations_lineno,
vector<optional<int>>& all_equations_lineno, vector<optional<tuple<int, expr_t, expr_t>>>& all_complementarity_conditions,
EquationTags& all_equation_tags, bool static_equations) const; EquationTags& all_equation_tags, bool static_equations) const;
//! Compute autoregressive matrices of trend component models //! Compute autoregressive matrices of trend component models
/* The algorithm uses matching rules over expression trees. It cannot handle /* The algorithm uses matching rules over expression trees. It cannot handle
...@@ -455,7 +458,8 @@ public: ...@@ -455,7 +458,8 @@ public:
Returns the number of optimality FOCs, which is by construction equal to Returns the number of optimality FOCs, which is by construction equal to
the number of endogenous before adding the Lagrange multipliers the number of endogenous before adding the Lagrange multipliers
(internally called ramsey_endo_nbr). */ (internally called ramsey_endo_nbr). */
int computeRamseyPolicyFOCs(const StaticModel& static_model); int computeRamseyPolicyFOCs(const StaticModel& planner_objective,
map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints);
//! Clears all equations //! Clears all equations
void clearEquations(); void clearEquations();
...@@ -464,7 +468,9 @@ public: ...@@ -464,7 +468,9 @@ public:
void replaceMyEquations(DynamicModel& dynamic_model) const; void replaceMyEquations(DynamicModel& dynamic_model) const;
//! Adds an equation marked as [static] //! Adds an equation marked as [static]
void addStaticOnlyEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags); void addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition,
map<string, string> eq_tags);
//! Returns number of static only equations //! Returns number of static only equations
size_t staticOnlyEquationsNbr() const; size_t staticOnlyEquationsNbr() const;
...@@ -708,7 +714,7 @@ public: ...@@ -708,7 +714,7 @@ public:
getStaticOnlyEquationsInfo() const getStaticOnlyEquationsInfo() const
{ {
return tuple {static_only_equations, static_only_equations_lineno, return tuple {static_only_equations, static_only_equations_lineno,
static_only_equations_equation_tags}; static_only_complementarity_conditions, static_only_equations_equation_tags};
}; };
//! Returns true if a parameter was used in the model block with a lead or lag //! Returns true if a parameter was used in the model block with a lead or lag
......
src/DynareBison.yy
View file @ 253c1866
...@@ -242,7 +242,7 @@ str_tolower(string s) ...@@ -242,7 +242,7 @@ str_tolower(string s)
%type <map<string, string>> tag_pair_list %type <map<string, string>> tag_pair_list
%type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt %type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
%type <vector<pair<int, int>>> period_list %type <vector<pair<int, int>>> period_list
%type <vector<expr_t>> matched_moments_list value_list %type <vector<expr_t>> matched_moments_list value_list ramsey_constraints_list
%type <tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> occbin_constraints_regime %type <tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> occbin_constraints_regime
%type <vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>>> occbin_constraints_regimes_list %type <vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>>> occbin_constraints_regimes_list
%type <map<string, expr_t>> occbin_constraints_regime_options_list %type <map<string, expr_t>> occbin_constraints_regime_options_list
...@@ -2616,23 +2616,20 @@ ramsey_policy : RAMSEY_POLICY ';' ...@@ -2616,23 +2616,20 @@ ramsey_policy : RAMSEY_POLICY ';'
{ driver.ramsey_policy($5); } { driver.ramsey_policy($5); }
; ;
ramsey_constraints : RAMSEY_CONSTRAINTS ';' ramsey_constraints_list END ';' ramsey_constraints : RAMSEY_CONSTRAINTS ';'
{ driver.add_ramsey_constraints_statement(); } { driver.begin_ramsey_constraints(); }
ramsey_constraints_list END ';'
{ driver.end_ramsey_constraints($4); }
; ;
ramsey_constraints_list : ramsey_constraints_list ramsey_constraint ramsey_constraints_list : ramsey_constraints_list hand_side ';'
| ramsey_constraint {
; $$ = $1;
$$.push_back($2);
ramsey_constraint : NAME LESS expression ';' }
{ driver.ramsey_constraint_add_less($1, $3); } | hand_side ';'
| NAME GREATER expression ';' { $$ = { $1 }; }
{ driver.ramsey_constraint_add_greater($1, $3); } ;
| NAME LESS_EQUAL expression ';'
{ driver.ramsey_constraint_add_less_equal($1, $3); }
| NAME GREATER_EQUAL expression ';'
{ driver.ramsey_constraint_add_greater_equal($1, $3); }
;
evaluate_planner_objective : EVALUATE_PLANNER_OBJECTIVE ';' evaluate_planner_objective : EVALUATE_PLANNER_OBJECTIVE ';'
{ driver.evaluate_planner_objective(); } { driver.evaluate_planner_objective(); }
......
src/ExprNode.cc
View file @ 253c1866
/* /*
* Copyright © 2007-2023 Dynare Team * Copyright © 2007-2024 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -422,7 +422,7 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs, ...@@ -422,7 +422,7 @@ ExprNode::fillErrorCorrectionRow(int eqn, const vector<int>& nontarget_lhs,
<< endl; << endl;
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
} }
if (*param_id) if (param_id)
{ {
cerr cerr
<< "ERROR in trend component model: spurious parameter in error correction term" << "ERROR in trend component model: spurious parameter in error correction term"
...@@ -9370,3 +9370,76 @@ ExprNode::toString() const ...@@ -9370,3 +9370,76 @@ ExprNode::toString() const
writeJsonOutput(ss, {}, {}); writeJsonOutput(ss, {}, {});
return ss.str(); return ss.str();
} }
tuple<int, expr_t, expr_t>
ExprNode::matchComplementarityCondition() const
{
throw MatchFailureException {"This expression is not an inequality"};
}
tuple<int, expr_t, expr_t>
BinaryOpNode::matchComplementarityCondition() 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 && datatree.symbol_table.getType(ve->symb_id) == SymbolType::endogenous)
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 {"Complementarity condition does not have the right form"};
auto id = match_contemporaneous_endogenous(barg1->arg2);
if (!id)
throw MatchFailureException {"Complementarity condition does not have the right form"};
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 @ 253c1866
...@@ -942,6 +942,11 @@ public: ...@@ -942,6 +942,11 @@ public:
// Substitutes orig_symb_id(±l) with exp(aux_symb_id(±l)) (used for “var(log)”) // 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; [[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;
}; };
//! Object used to compare two nodes (using their indexes) //! Object used to compare two nodes (using their indexes)
...@@ -1499,6 +1504,7 @@ public: ...@@ -1499,6 +1504,7 @@ public:
vector<int>& powers) const override; vector<int>& powers) const override;
[[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() 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 substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
[[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
}; };
//! Trinary operator node //! Trinary operator node
......
src/ModFile.cc
View file @ 253c1866
/* /*
* Copyright © 2006-2023 Dynare Team * Copyright © 2006-2024 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -167,7 +167,7 @@ ModFile::checkPass(bool nostrict, bool stochastic) ...@@ -167,7 +167,7 @@ ModFile::checkPass(bool nostrict, bool stochastic)
exit(EXIT_FAILURE); 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 " cerr << "ERROR: A ramsey_constraints block requires the presence of a ramsey_model or "
"ramsey_policy statement" "ramsey_policy statement"
...@@ -535,8 +535,15 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool ...@@ -535,8 +535,15 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
symbol_table, num_constants, external_functions_table, trend_component_model_table, symbol_table, num_constants, external_functions_table, trend_component_model_table,
var_model_table}; var_model_table};
ramsey_FOC_equations_dynamic_model = dynamic_model; 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 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); ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
} }
......
src/ModFile.hh
View file @ 253c1866
/* /*
* Copyright © 2006-2023 Dynare Team * Copyright © 2006-2024 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -117,6 +117,11 @@ public: ...@@ -117,6 +117,11 @@ public:
/*! (i.e. option parallel_local_files of model block) */ /*! (i.e. option parallel_local_files of model block) */
vector<string> parallel_local_files; 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: private:
//! List of statements //! List of statements
vector<unique_ptr<Statement>> statements; vector<unique_ptr<Statement>> statements;
......
src/ModelEquationBlock.cc
View file @ 253c1866
...@@ -377,7 +377,7 @@ Epilogue::checkPass(ModFileStructure& mod_file_struct) const ...@@ -377,7 +377,7 @@ Epilogue::checkPass(ModFileStructure& mod_file_struct) const
for (const auto& [symb_id, expr] : dynamic_def_table) for (const auto& [symb_id, expr] : dynamic_def_table)
if (so_far_defined.contains(symb_id)) 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; << "' is declared twice" << endl;
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
} }
......
src/ModelTree.cc
View file @ 253c1866
/* /*
* Copyright © 2003-2023 Dynare Team * Copyright © 2003-2024 Dynare Team
* *
* This file is part of Dynare. * This file is part of Dynare.
* *
...@@ -34,6 +34,7 @@ ...@@ -34,6 +34,7 @@
#endif #endif
#include <algorithm> #include <algorithm>
#include <numeric>
#include <regex> #include <regex>
#include <utility> #include <utility>
...@@ -51,13 +52,21 @@ vector<jthread> ModelTree::mex_compilation_workers; ...@@ -51,13 +52,21 @@ vector<jthread> ModelTree::mex_compilation_workers;
void void
ModelTree::copyHelper(const ModelTree& m) 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 // Equations
for (const auto& it : m.equations) for (const auto& it : m.equations)
equations.push_back(dynamic_cast<BinaryOpNode*>(f(it))); equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
for (const auto& it : m.aux_equations) for (const auto& it : m.aux_equations)
aux_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it))); 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) { auto convert_deriv_map = [f](const map<vector<int>, expr_t>& dm) {
map<vector<int>, expr_t> dm2; map<vector<int>, expr_t> dm2;
...@@ -184,6 +193,8 @@ ModelTree::operator=(const ModelTree& m) ...@@ -184,6 +193,8 @@ ModelTree::operator=(const ModelTree& m)
equations_lineno = m.equations_lineno; equations_lineno = m.equations_lineno;
aux_equations.clear(); aux_equations.clear();
equation_tags = m.equation_tags; equation_tags = m.equation_tags;
complementarity_conditions.clear();
computed_derivs_order = m.computed_derivs_order; computed_derivs_order = m.computed_derivs_order;
NNZDerivatives = m.NNZDerivatives; NNZDerivatives = m.NNZDerivatives;
...@@ -1394,28 +1405,33 @@ ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_b ...@@ -1394,28 +1405,33 @@ ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_b
} }
void 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); auto beq = dynamic_cast<BinaryOpNode*>(eq);
assert(beq && beq->op_code == BinaryOpcode::equal); assert(beq && beq->op_code == BinaryOpcode::equal);
equations.push_back(beq); equations.push_back(beq);
equations_lineno.push_back(lineno); equations_lineno.push_back(lineno);
complementarity_conditions.push_back(move(complementarity_condition));
} }
void 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++) for (size_t i = 0; i < equations.size(); i++)
if (!equation_tags.exists(i, "mcp")) if (!complementarity_conditions[i])
equations[i]->findConstantEquations(subst_table); equations[i]->findConstantEquations(subst_table);
} }
void 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)); equation_tags.add(equations.size(), move(eq_tags));
addEquation(eq, lineno); addEquation(eq, lineno, move(complementarity_condition));
} }
void void
...@@ -1593,6 +1609,27 @@ ModelTree::writeJsonModelEquations(ostream& output, bool residuals) const ...@@ -1593,6 +1609,27 @@ ModelTree::writeJsonModelEquations(ostream& output, bool residuals) const
eqtags.clear(); 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;
} }
output << endl << "]" << endl; output << endl << "]" << endl;
...@@ -2078,3 +2115,38 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o ...@@ -2078,3 +2115,38 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o
output << ";" << endl; output << ";" << endl;
} }
} }
void
ModelTree::computeMCPEquationsReordering()
{
/* 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()) != 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);
}
}
src/ModelTree.hh
View file @ 253c1866
...@@ -91,6 +91,9 @@ protected: ...@@ -91,6 +91,9 @@ protected:
vector<optional<int>> equations_lineno; vector<optional<int>> equations_lineno;
//! Stores equation tags //! Stores equation tags
EquationTags 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 ************** * ************** END **************
*/ */
...@@ -271,6 +274,12 @@ protected: ...@@ -271,6 +274,12 @@ protected:
Same remark as above regarding blocks of type “evaluate”. */ Same remark as above regarding blocks of type “evaluate”. */
vector<vector<int>> blocks_jacobian_sparse_colptr; 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 //! Computes derivatives
/*! \param order the derivation order /*! \param order the derivation order
\param vars the derivation IDs w.r.t. which compute the derivatives */ \param vars the derivation IDs w.r.t. which compute the derivatives */
...@@ -284,6 +293,11 @@ protected: ...@@ -284,6 +293,11 @@ protected:
//! Computes temporary terms for the file containing parameters derivatives //! Computes temporary terms for the file containing parameters derivatives
void computeParamsDerivativesTemporaryTerms(); void computeParamsDerivativesTemporaryTerms();
/* Computes the mcp_equations_reordering vector.
Also checks that a variable does not appear as constrained in two different equations. */
void computeMCPEquationsReordering();
//! Writes temporary terms //! Writes temporary terms
template<ExprNodeOutputType output_type> template<ExprNodeOutputType output_type>
void writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union, void writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union,
...@@ -424,6 +438,9 @@ protected: ...@@ -424,6 +438,9 @@ protected:
template<bool dynamic> template<bool dynamic>
void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const; void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
template<bool dynamic>
void writeComplementarityConditionsFile(const string& basename) const;
private: private:
//! Sparse matrix of double to store the values of the static Jacobian //! Sparse matrix of double to store the values of the static Jacobian
/*! First index is equation number, second index is endogenous type specific ID */ /*! First index is equation number, second index is endogenous type specific ID */
...@@ -651,10 +668,14 @@ protected: ...@@ -651,10 +668,14 @@ protected:
public: public:
//! Absolute value under which a number is considered to be zero //! Absolute value under which a number is considered to be zero
double cutoff {1e-15}; double cutoff {1e-15};
//! Declare a node as an equation of the model; also give its line number /* Declare a node as an equation of the model; also give its line number and complementarity
void addEquation(expr_t eq, const optional<int>& lineno); 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 //! 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 //! Declare a node as an auxiliary equation of the model, adding it at the end of the list of
//! auxiliary equations //! auxiliary equations
void addAuxEquation(expr_t eq); void addAuxEquation(expr_t eq);
...@@ -671,9 +692,10 @@ public: ...@@ -671,9 +692,10 @@ public:
/*! Reorder auxiliary variables so that they appear in recursive order in /*! Reorder auxiliary variables so that they appear in recursive order in
set_auxiliary_variables.m and dynamic_set_auxiliary_series.m */ set_auxiliary_variables.m and dynamic_set_auxiliary_series.m */
void reorderAuxiliaryEquations(); void reorderAuxiliaryEquations();
//! Find equations of the form “variable=constant”, excluding equations with “mcp” tag (see /* Find equations of the form “variable=constant”, excluding equations with a complementarity
//! dynare#1697) condition (see dynare#1697) */
void findConstantEquationsWithoutMcpTag(map<VariableNode*, NumConstNode*>& subst_table) const; void findConstantEquationsWithoutComplementarityCondition(
map<VariableNode*, NumConstNode*>& subst_table) const;
/* Given an expression, searches for the first equation that has exactly this /* Given an expression, searches for the first equation that has exactly this
expression on the LHS, and returns the RHS of that equation. expression on the LHS, and returns the RHS of that equation.
If no such equation can be found, throws an ExprNode::MatchFailureExpression */ If no such equation can be found, throws an ExprNode::MatchFailureExpression */
...@@ -3133,4 +3155,50 @@ ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) co ...@@ -3133,4 +3155,50 @@ ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) co
} }
} }
template<bool dynamic>
void
ModelTree::writeComplementarityConditionsFile(const string& basename) const
{
// TODO: when C++20 support is complete, mark the following string constexpr
const string funcname {(dynamic ? "dynamic"s : "static"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 #endif
src/ParsingDriver.cc
View file @ 253c1866
...@@ -2635,6 +2635,9 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_ ...@@ -2635,6 +2635,9 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
if (key == "endogenous") if (key == "endogenous")
declare_or_change_type(SymbolType::endogenous, value); 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 (eq_tags.contains("mcp"))
{ {
if (complementarity_condition) if (complementarity_condition)
...@@ -2643,20 +2646,54 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_ ...@@ -2643,20 +2646,54 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
else else
warning("Specifying complementarity conditions with the 'mcp' tag is obsolete. Please " warning("Specifying complementarity conditions with the 'mcp' tag is obsolete. Please "
"consider switching to the new syntax using the perpendicular symbol."); "consider switching to the new syntax using the perpendicular symbol.");
}
if (complementarity_condition) auto [var_name, is_lower_bound, constant] {[&]() -> tuple<string, bool, string> {
{ auto tagval = eq_tags.at("mcp");
if (auto bcomp = dynamic_cast<BinaryOpNode*>(complementarity_condition); if (auto less_pos = tagval.find('<'); less_pos != string::npos)
!(bcomp return {tagval.substr(0, less_pos), false, tagval.substr(less_pos + 1)};
&& (bcomp->op_code == BinaryOpcode::less || bcomp->op_code == BinaryOpcode::lessEqual else if (auto greater_pos = tagval.find('>'); greater_pos != string::npos)
|| bcomp->op_code == BinaryOpcode::greater return {tagval.substr(0, greater_pos), true, tagval.substr(greater_pos + 1)};
|| bcomp->op_code == BinaryOpcode::greaterEqual))) else
error("The complementarity constraint must be an inequality."); error("'mcp' tag does not contain an inequality");
}()};
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 (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);
}()};
eq_tags.emplace("mcp", complementarity_condition->toString()); matched_complementarity_condition = {symb_id, is_lower_bound ? matched_constant : nullptr,
is_lower_bound ? nullptr : matched_constant};
} }
if (complementarity_condition)
try
{
matched_complementarity_condition
= complementarity_condition->matchComplementarityCondition();
}
catch (ExprNode::MatchFailureException& e)
{
error("Complementarity condition has an incorrect form: " + e.message);
}
if (eq_tags.contains("static")) if (eq_tags.contains("static"))
{ {
// If the equation is tagged [static] // If the equation is tagged [static]
...@@ -2664,7 +2701,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_ ...@@ -2664,7 +2701,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 " error("An equation tagged [static] cannot contain leads, lags, expectations or "
"STEADY_STATE operators"); "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")) else if (eq_tags.contains("bind") || eq_tags.contains("relax"))
{ {
...@@ -2740,7 +2778,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_ ...@@ -2740,7 +2778,8 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2, map<string, string> eq_
} }
} }
else // General case 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; return id;
} }
...@@ -3629,48 +3668,31 @@ ParsingDriver::prior_posterior_function(bool prior_func) ...@@ -3629,48 +3668,31 @@ ParsingDriver::prior_posterior_function(bool prior_func)
} }
void 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();
}
void
ParsingDriver::ramsey_constraint_add_less(const string& name, const expr_t rhs)
{
add_ramsey_constraint(name, BinaryOpcode::less, rhs);
}
void
ParsingDriver::ramsey_constraint_add_greater(const string& name, const expr_t rhs)
{ {
add_ramsey_constraint(name, BinaryOpcode::greater, rhs); set_current_data_tree(&mod_file->dynamic_model);
}
void
ParsingDriver::ramsey_constraint_add_less_equal(const string& name, const expr_t rhs)
{
add_ramsey_constraint(name, BinaryOpcode::lessEqual, rhs);
} }
void void
ParsingDriver::ramsey_constraint_add_greater_equal(const string& name, const expr_t rhs) ParsingDriver::end_ramsey_constraints(const vector<expr_t>& constraints)
{ {
add_ramsey_constraint(name, BinaryOpcode::greaterEqual, rhs); for (expr_t c : constraints)
} try
{
auto [symb_id, lower_bound, upper_bound] = c->matchComplementarityCondition();
void auto [it, success]
ParsingDriver::add_ramsey_constraint(const string& name, BinaryOpcode op_code, const expr_t rhs) = mod_file->ramsey_constraints.try_emplace(symb_id, lower_bound, upper_bound);
{ if (!success)
check_symbol_is_endogenous(name); error("The ramsey_constraints block contains two constraints for variable "
int symb_id = mod_file->symbol_table.getID(name); + mod_file->symbol_table.getName(symb_id));
}
catch (ExprNode::MatchFailureException& e)
{
error("Ramsey constraint has an incorrect form: " + e.message);
}
RamseyConstraintsStatement::Constraint C; reset_data_tree();
C.endo = symb_id;
C.code = op_code;
C.expression = rhs;
ramsey_constraints.push_back(C);
} }
void void
......
src/ParsingDriver.hh
View file @ 253c1866
...@@ -186,8 +186,6 @@ private: ...@@ -186,8 +186,6 @@ private:
MomentCalibration::constraints_t moment_calibration_constraints; MomentCalibration::constraints_t moment_calibration_constraints;
//! Temporary storage for irf_calibration //! Temporary storage for irf_calibration
IrfCalibration::constraints_t irf_calibration_constraints; IrfCalibration::constraints_t irf_calibration_constraints;
//! Temporary storage for ramsey_constraints
RamseyConstraintsStatement::constraints_t ramsey_constraints;
//! Temporary storage for svar_identification blocks //! Temporary storage for svar_identification blocks
SvarIdentificationStatement::svar_identification_restrictions_t svar_ident_restrictions; SvarIdentificationStatement::svar_identification_restrictions_t svar_ident_restrictions;
//! Temporary storage for mapping the equation number to the restrictions within an //! Temporary storage for mapping the equation number to the restrictions within an
...@@ -667,18 +665,9 @@ public: ...@@ -667,18 +665,9 @@ public:
void end_planner_objective(expr_t expr); void end_planner_objective(expr_t expr);
//! Ramsey model statement //! Ramsey model statement
void ramsey_model(); void ramsey_model();
//! Ramsey constraints statement // Ramsey constraints statement
void add_ramsey_constraints_statement(); void begin_ramsey_constraints();
//! Ramsey less constraint void end_ramsey_constraints(const vector<expr_t>& constraints);
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 policy statement //! Ramsey policy statement
void ramsey_policy(vector<string> symbol_list); void ramsey_policy(vector<string> symbol_list);
//! Evaluate Planner Objective //! Evaluate Planner Objective
......
src/Statement.hh
View file @ 253c1866
...@@ -169,8 +169,6 @@ struct ModFileStructure ...@@ -169,8 +169,6 @@ struct ModFileStructure
bool endval_learnt_in_present {false}; bool endval_learnt_in_present {false};
// Whether an occbin_constraints block appears // Whether an occbin_constraints block appears
bool occbin_constraints_present {false}; bool occbin_constraints_present {false};
// Whether a ramsey_constraints block appears
bool ramsey_constraints_present {false};
}; };
class Statement class Statement
......
src/StaticModel.cc
View file @ 253c1866
...@@ -87,17 +87,18 @@ StaticModel::StaticModel(const DynamicModel& m) : ...@@ -87,17 +87,18 @@ StaticModel::StaticModel(const DynamicModel& m) :
if (dynamic_equations.contains(i)) if (dynamic_equations.contains(i))
{ {
auto [static_only_equations, static_only_equations_lineno, auto [static_only_equations, static_only_equations_lineno,
static_only_equations_equation_tags] static_only_complementarity_conditions, static_only_equations_equation_tags]
= m.getStaticOnlyEquationsInfo(); = m.getStaticOnlyEquationsInfo();
addEquation(static_only_equations[static_only_index]->toStatic(*this), addEquation(static_only_equations[static_only_index]->toStatic(*this),
static_only_equations_lineno[static_only_index], static_only_equations_lineno[static_only_index],
static_only_complementarity_conditions[static_only_index],
static_only_equations_equation_tags.getTagsByEqn(static_only_index)); static_only_equations_equation_tags.getTagsByEqn(static_only_index));
static_only_index++; static_only_index++;
} }
else else
addEquation(m.equations[i]->toStatic(*this), m.equations_lineno[i], addEquation(m.equations[i]->toStatic(*this), m.equations_lineno[i],
m.equation_tags.getTagsByEqn(i)); m.complementarity_conditions[i], m.equation_tags.getTagsByEqn(i));
} }
catch (DataTree::DivisionByZeroException) catch (DataTree::DivisionByZeroException)
{ {
...@@ -250,6 +251,8 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, ...@@ -250,6 +251,8 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder,
<< endl; << endl;
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
} }
computeMCPEquationsReordering();
} }
void void
...@@ -305,6 +308,8 @@ StaticModel::writeStaticFile(const string& basename, bool use_dll, const string& ...@@ -305,6 +308,8 @@ StaticModel::writeStaticFile(const string& basename, bool use_dll, const string&
writeSetAuxiliaryVariablesFile<false>(basename, julia); writeSetAuxiliaryVariablesFile<false>(basename, julia);
writeComplementarityConditionsFile<false>(basename);
// Support for model debugging // Support for model debugging
if (!julia) if (!julia)
writeDebugModelMFiles<false>(basename); writeDebugModelMFiles<false>(basename);
...@@ -331,6 +336,11 @@ StaticModel::writeDriverOutput(ostream& output) const ...@@ -331,6 +336,11 @@ StaticModel::writeDriverOutput(ostream& output) const
writeBlockDriverOutput(output); writeBlockDriverOutput(output);
writeDriverSparseIndicesHelper<false, false>(output); writeDriverSparseIndicesHelper<false, false>(output);
output << "M_.static_mcp_equations_reordering = [";
for (auto i : mcp_equations_reordering)
output << i + 1 << "; ";
output << "];" << endl;
} }
void void
......