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
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,
OptionsList options_list_arg,
const SymbolTable& symbol_table_arg) :
......@@ -2014,7 +1930,7 @@ void
OsrParamsStatement::checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings)
{
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;
try
......
src/ComputingTasks.hh
View file @ 253c1866
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2024 Dynare Team
*
* This file is part of Dynare.
*
......@@ -179,28 +179,6 @@ public:
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
{
private:
......
src/DynamicModel.cc
View file @ 253c1866
......@@ -34,10 +34,18 @@
void
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)
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,
......@@ -104,6 +112,7 @@ DynamicModel::operator=(const DynamicModel& m)
static_only_equations_lineno = m.static_only_equations_lineno;
static_only_equations_equation_tags = m.static_only_equations_equation_tags;
static_only_complementarity_conditions.clear();
deriv_id_table = m.deriv_id_table;
inv_deriv_id_table = m.inv_deriv_id_table;
dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
......@@ -651,11 +660,11 @@ DynamicModel::parseIncludeExcludeEquations(const string& inc_exc_option_value, b
}
vector<int>
DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
bool excluded_vars_change_type,
vector<BinaryOpNode*>& all_equations,
vector<optional<int>>& all_equations_lineno,
EquationTags& all_equation_tags, bool static_equations) const
DynamicModel::removeEquationsHelper(
set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
vector<BinaryOpNode*>& all_equations, 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
{
if (all_equations.empty())
return {};
......@@ -686,6 +695,7 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
// remove from equations, equations_lineno, equation_tags
vector<BinaryOpNode*> new_equations;
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;
vector<int> excluded_vars;
for (size_t i = 0; i < all_equations.size(); i++)
......@@ -717,11 +727,13 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
new_equations.emplace_back(all_equations[i]);
old_eqn_num_2_new[i] = new_equations.size() - 1;
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();
all_equations = new_equations;
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);
......@@ -754,12 +766,13 @@ DynamicModel::removeEquations(const vector<map<string, string>>& listed_eqs_by_t
vector<int> excluded_vars
= 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
removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type,
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())
{
......@@ -1134,6 +1147,11 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
output << "'; " << endl;
}
output << "};" << endl;
output << "M_.dynamic_mcp_equations_reordering = [";
for (auto i : mcp_equations_reordering)
output << i + 1 << "; ";
output << "];" << endl;
}
void
......@@ -2507,6 +2525,8 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder,
cerr << "ERROR: Block decomposition requested but failed." << endl;
exit(EXIT_FAILURE);
}
computeMCPEquationsReordering();
}
void
......@@ -2810,6 +2830,8 @@ DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const strin
writeSetAuxiliaryVariablesFile<true>(basename, julia);
writeComplementarityConditionsFile<true>(basename);
// Support for model debugging
if (!julia)
writeDebugModelMFiles<true>(basename);
......@@ -2821,6 +2843,7 @@ DynamicModel::clearEquations()
equations.clear();
equations_lineno.clear();
equation_tags.clear();
complementarity_conditions.clear();
}
void
......@@ -2828,14 +2851,26 @@ DynamicModel::replaceMyEquations(DynamicModel& dynamic_model) const
{
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++)
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;
}
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;
......@@ -2849,8 +2884,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
}
// Add Planner Objective to equations so that it appears in Lagrangian
assert(static_model.equations.size() == 1);
addEquation(static_model.equations[0]->clone(*this), nullopt);
assert(planner_objective.equations.size() == 1);
addEquation(planner_objective.equations[0]->clone(*this), nullopt);
// Get max endo lead and max endo lag
set<pair<int, int>> dynvars;
......@@ -2892,9 +2927,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
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_equation_tags = equation_tags;
auto old_complementarity_conditions = complementarity_conditions;
// Prepare derivation of the Lagrangian
clearEquations();
......@@ -2911,6 +2947,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
vector<expr_t> neweqs;
vector<optional<int>> neweqs_lineno;
map<int, map<string, string>> neweqs_tags;
map<int, optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
int orig_endo_nbr {0};
for (auto& [symb_id_and_lag, deriv_id] : deriv_id_table)
{
......@@ -2923,12 +2960,20 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
{
// This is a derivative w.r.t. a Lagrange multiplier
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
{
orig_endo_nbr++;
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)
// Overwrite equations with the Lagrangian derivatives
clearEquations();
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;
}
......@@ -3773,6 +3818,7 @@ DynamicModel::fillEvalContext(eval_context_t& eval_context) const
void
DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
optional<tuple<int, expr_t, expr_t>> complementarity_condition,
map<string, string> eq_tags)
{
auto beq = dynamic_cast<BinaryOpNode*>(eq);
......@@ -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.push_back(beq);
static_only_equations_lineno.push_back(lineno);
static_only_complementarity_conditions.push_back(move(complementarity_condition));
}
size_t
......@@ -3834,7 +3881,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
{
auto eq_tags_dynamic = eq_tags;
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)
......@@ -3856,7 +3903,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
else
{
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()
{
size_t last_subst_table_size = 0;
map<VariableNode*, NumConstNode*> subst_table;
// Equations with “mcp” tag are excluded, see dynare#1697
findConstantEquationsWithoutMcpTag(subst_table);
// Equations with a complementarity condition are excluded, see dynare#1697
findConstantEquationsWithoutComplementarityCondition(subst_table);
while (subst_table.size() != last_subst_table_size)
{
last_subst_table_size = subst_table.size();
......@@ -4183,7 +4230,7 @@ DynamicModel::simplifyEquations()
for (auto& equation : static_only_equations)
equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table));
subst_table.clear();
findConstantEquationsWithoutMcpTag(subst_table);
findConstantEquationsWithoutComplementarityCondition(subst_table);
}
}
......
src/DynamicModel.hh
View file @ 253c1866
......@@ -92,6 +92,9 @@ private:
//! Stores the equation tags of equations declared as [static]
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>;
//! Maps a pair (symbol_id, lag) to a deriv ID
deriv_id_table_t deriv_id_table;
......@@ -279,11 +282,11 @@ private:
Returns a list of excluded variables (empty if
excluded_vars_change_type=false) */
vector<int> removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
bool excluded_vars_change_type,
vector<BinaryOpNode*>& all_equations,
vector<optional<int>>& all_equations_lineno,
EquationTags& all_equation_tags, bool static_equations) const;
vector<int> removeEquationsHelper(
set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
vector<BinaryOpNode*>& all_equations, 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;
//! Compute autoregressive matrices of trend component models
/* The algorithm uses matching rules over expression trees. It cannot handle
......@@ -455,7 +458,8 @@ public:
Returns the number of optimality FOCs, which is by construction equal to
the number of endogenous before adding the Lagrange multipliers
(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
void clearEquations();
......@@ -464,7 +468,9 @@ public:
void replaceMyEquations(DynamicModel& dynamic_model) const;
//! 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
size_t staticOnlyEquationsNbr() const;
......@@ -708,7 +714,7 @@ public:
getStaticOnlyEquationsInfo() const
{
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
......
src/DynareBison.yy
View file @ 253c1866
......@@ -242,7 +242,7 @@ str_tolower(string s)
%type <map<string, string>> tag_pair_list
%type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
%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 <vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>>> occbin_constraints_regimes_list
%type <map<string, expr_t>> occbin_constraints_regime_options_list
......@@ -2616,23 +2616,20 @@ ramsey_policy : RAMSEY_POLICY ';'
{ driver.ramsey_policy($5); }
;
ramsey_constraints : RAMSEY_CONSTRAINTS ';' ramsey_constraints_list END ';'
{ driver.add_ramsey_constraints_statement(); }
ramsey_constraints : RAMSEY_CONSTRAINTS ';'
{ driver.begin_ramsey_constraints(); }
ramsey_constraints_list END ';'
{ driver.end_ramsey_constraints($4); }
;
ramsey_constraints_list : ramsey_constraints_list ramsey_constraint
| ramsey_constraint
;
ramsey_constraint : NAME LESS expression ';'
{ driver.ramsey_constraint_add_less($1, $3); }
| NAME GREATER expression ';'
{ 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); }
;
ramsey_constraints_list : ramsey_constraints_list hand_side ';'
{
$$ = $1;
$$.push_back($2);
}
| hand_side ';'
{ $$ = { $1 }; }
;
evaluate_planner_objective : 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.
*
......@@ -422,7 +422,7 @@ 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"
......@@ -9370,3 +9370,76 @@ ExprNode::toString() const
writeJsonOutput(ss, {}, {});
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:
// 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;
};
//! Object used to compare two nodes (using their indexes)
......@@ -1499,6 +1504,7 @@ 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]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
};
//! 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.
*
......@@ -167,7 +167,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"
......@@ -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,
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);
}
......
src/ModFile.hh
View file @ 253c1866
/*
* Copyright © 2006-2023 Dynare Team
* Copyright © 2006-2024 Dynare Team
*
* This file is part of Dynare.
*
......@@ -117,6 +117,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 @ 253c1866
......@@ -377,7 +377,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);
}
......
src/ModelTree.cc
View file @ 253c1866
/*
* Copyright © 2003-2023 Dynare Team
* Copyright © 2003-2024 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;
......@@ -184,6 +193,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;
......@@ -1394,28 +1405,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
......@@ -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;
......@@ -2078,3 +2115,38 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o
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:
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 **************
*/
......@@ -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();
//! Writes temporary terms
template<ExprNodeOutputType output_type>
void writeTemporaryTerms(const temporary_terms_t& tt, temporary_terms_t& temp_term_union,
......@@ -424,6 +438,9 @@ protected:
template<bool dynamic>
void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
template<bool dynamic>
void writeComplementarityConditionsFile(const string& basename) 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 */
......@@ -651,10 +668,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 +692,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 */
......@@ -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
src/ParsingDriver.cc
View file @ 253c1866
......@@ -2635,6 +2635,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,20 +2646,54 @@ 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.");
}
if (complementarity_condition)
{
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.");
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");
}()};
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 the equation is tagged [static]
......@@ -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 "
"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 +2778,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;
}
......@@ -3629,48 +3668,31 @@ ParsingDriver::prior_posterior_function(bool prior_func)
}
void
ParsingDriver::add_ramsey_constraints_statement()
{
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)
ParsingDriver::begin_ramsey_constraints()
{
add_ramsey_constraint(name, BinaryOpcode::greater, rhs);
}
void
ParsingDriver::ramsey_constraint_add_less_equal(const string& name, const expr_t rhs)
{
add_ramsey_constraint(name, BinaryOpcode::lessEqual, rhs);
set_current_data_tree(&mod_file->dynamic_model);
}
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
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);
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));
}
catch (ExprNode::MatchFailureException& e)
{
error("Ramsey constraint has an incorrect form: " + e.message);
}
RamseyConstraintsStatement::Constraint C;
C.endo = symb_id;
C.code = op_code;
C.expression = rhs;
ramsey_constraints.push_back(C);
reset_data_tree();
}
void
......
src/ParsingDriver.hh
View file @ 253c1866
......@@ -186,8 +186,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
......@@ -667,18 +665,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
......
src/Statement.hh
View file @ 253c1866
......@@ -169,8 +169,6 @@ struct ModFileStructure
bool endval_learnt_in_present {false};
// Whether an occbin_constraints block appears
bool occbin_constraints_present {false};
// Whether a ramsey_constraints block appears
bool ramsey_constraints_present {false};
};
class Statement
......
src/StaticModel.cc
View file @ 253c1866
......@@ -87,17 +87,18 @@ StaticModel::StaticModel(const DynamicModel& m) :
if (dynamic_equations.contains(i))
{
auto [static_only_equations, static_only_equations_lineno,
static_only_equations_equation_tags]
static_only_complementarity_conditions, static_only_equations_equation_tags]
= m.getStaticOnlyEquationsInfo();
addEquation(static_only_equations[static_only_index]->toStatic(*this),
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_index++;
}
else
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)
{
......@@ -250,6 +251,8 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder,
<< endl;
exit(EXIT_FAILURE);
}
computeMCPEquationsReordering();
}
void
......@@ -305,6 +308,8 @@ StaticModel::writeStaticFile(const string& basename, bool use_dll, const string&
writeSetAuxiliaryVariablesFile<false>(basename, julia);
writeComplementarityConditionsFile<false>(basename);
// Support for model debugging
if (!julia)
writeDebugModelMFiles<false>(basename);
......@@ -331,6 +336,11 @@ StaticModel::writeDriverOutput(ostream& output) const
writeBlockDriverOutput(output);
writeDriverSparseIndicesHelper<false, false>(output);
output << "M_.static_mcp_equations_reordering = [";
for (auto i : mcp_equations_reordering)
output << i + 1 << "; ";
output << "];" << endl;
}
void
......