Compare revisions

Sébastien Villemot · Sébastien Villemot · Sébastien Villemot · Sébastien Villemot · Johannes Pfeifer · Sébastien Villemot
--- a/.clang-format
+++ b/.clang-format
@@ -11,9 +11,12 @@ BreakInheritanceList: AfterColon
 Cpp11BracedListStyle: true
 DeriveLineEnding: false
 IndentPPDirectives: AfterHash
+InsertNewlineAtEOF: true
 PackConstructorInitializers: NextLine
 PPIndentWidth: 1
 PointerAlignment: Left
+# RemoveParentheses: ReturnStatement
+# RemoveSemicolon: true
 SpaceAfterTemplateKeyword: false
 SpaceBeforeParens: ControlStatements
 SpaceBeforeCpp11BracedList: true
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
 variables:
  TERM: linux
-  MINGW64_BOOST_VERSION: 1.85.0-2
+  MINGW64_BOOST_VERSION: 1.86.0-7
  WGET_OPTIONS: '--no-verbose --no-use-server-timestamps --retry-connrefused --retry-on-host-error'
  # To ensure that "false && true" fails, see https://gitlab.com/gitlab-org/gitlab-runner/-/issues/25394#note_412609647
  FF_ENABLE_BASH_EXIT_CODE_CHECK: 'true'
@@ -77,3 +77,16 @@ test_clang_format:
    - meson setup build-clang-format
    - ninja -C build-clang-format clang-format-check
  needs: []
+test_clang_tidy:
+  stage: test
+  script:
+    # Hack needed for meson < 1.6.0 which only looks for unversioned clang-tidy
+    - mkdir -p ~/.local/bin && ln -s /usr/bin/clang-tidy-19 ~/.local/bin/clang-tidy
+    - export PATH="$HOME/.local/bin:$PATH"
+    - meson setup build-clang-tidy
+    # Generate Flex and Bison files
+    - meson compile -C build-clang-tidy
+    - ninja -C build-clang-tidy clang-tidy
+  needs: []
+  when: manual
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -283,7 +283,9 @@ PerfectForesightWithExpectationErrorsSolverStatement::writeOutput(
    [[maybe_unused]] bool minimal_workspace) const
 {
  options_list.writeOutput(output);
-  output << "oo_ = perfect_foresight_with_expectation_errors_solver(M_, options_, oo_);" << endl;
+  output << "[oo_, Simulated_time_series] = perfect_foresight_with_expectation_errors_solver(M_, "
+            "options_, oo_);"
+         << endl;
 }
 void
@@ -5340,8 +5342,8 @@ void
 ResidStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                            [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output, "options_resid_");
+  options_list.writeOutput(output);
-  output << "display_static_residuals(M_, options_, oo_, options_resid_);" << endl;
+  output << "display_static_residuals(M_, options_, oo_);" << endl;
 }
 void

--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -524,7 +524,7 @@ public:
  blockName() const override
  {
    return "estimated_params";
-  };
+  }
  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;
@@ -542,7 +542,7 @@ public:
  blockName() const override
  {
    return "estimated_params_init";
-  };
+  }
  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;
@@ -557,7 +557,7 @@ public:
  blockName() const override
  {
    return "estimated_params_bounds";
-  };
+  }
  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;

--- a/src/Configuration.hh
+++ b/src/Configuration.hh
@@ -47,7 +47,7 @@ private:
    get_paths() const
    {
      return paths;
-    };
+    }
  private:
    map<string, vector<string>> paths;

--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -70,13 +70,14 @@ DataTree::DataTree(const DataTree& d) :
  // Constants must be initialized first because they are used in some Add* methods
  initConstants();
+  // See commment in DataTree::operator=() for the rationale
+  for (int symb_id : d.local_variables_vector)
+    local_variables_table[symb_id] = d.local_variables_table.at(symb_id)->clone(*this);
  for (const auto& it : d.node_list)
    it->clone(*this);
  assert(node_list.size() == d.node_list.size());
-  for (const auto& [symb_id, value] : d.local_variables_table)
-    local_variables_table[symb_id] = value->clone(*this);
 }
 DataTree&

--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -367,7 +367,13 @@ public:
    if (it == local_variables_table.end())
      throw UnknownLocalVariableException {symb_id};
-    return it->second->decreaseLeadsLags(-lead_lag);
+    /* In the following, the case without lead/lag is optimized. It makes a difference on models
+       with many nested model-local variables, see e.g.
+       https://forum.dynare.org/t/pre-processing-takes-very-long/26865 */
+    if (lead_lag == 0)
+      return it->second;
+    else
+      return it->second->decreaseLeadsLags(-lead_lag);
  }
  static void
@@ -396,7 +402,7 @@ DataTree::AddPossiblyNegativeConstant(double v)
  if (isnan(v))
    return NaN;
  if (isinf(v))
-    return (v < 0 ? MinusInfinity : Infinity);
+    return v < 0 ? MinusInfinity : Infinity;
  bool neg = false;
  if (v < 0)

--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -731,9 +731,9 @@ DynamicModel::removeEquationsHelper(
      }
  int n_excl = all_equations.size() - new_equations.size();
-  all_equations = new_equations;
+  all_equations = move(new_equations);
-  all_equations_lineno = new_equations_lineno;
+  all_equations_lineno = move(new_equations_lineno);
-  all_complementarity_conditions = new_complementarity_conditions;
+  all_complementarity_conditions = move(new_complementarity_conditions);
  all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new);
@@ -1816,7 +1816,7 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
  for (auto& equation : equations)
    if (equation->containsPacExpectation(name))
      {
-        if (!pac_eq_name[name].empty())
+        if (pac_eq_name.contains(name))
          {
            cerr << "It is not possible to use 'pac_expectation(" << name
                 << ")' in several equations." << endl;
@@ -1915,6 +1915,13 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
                                   move(additive_vars_params_and_constants),
                                   move(optim_additive_vars_params_and_constants)};
      }
+  if (!pac_eq_name.contains(name))
+    {
+      cerr << "ERROR: the model does not contain the 'pac_expectation(" << name << ")' operator."
+           << endl;
+      exit(EXIT_FAILURE);
+    }
 }
 int

--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -73,7 +73,8 @@ public:
    void checkAllRegimesPresent() const noexcept(false);
  private:
-    pair<vector<string>, vector<string>> convertBitVectorToRegimes(const vector<bool>& r) const;
+    [[nodiscard]] pair<vector<string>, vector<string>>
+    convertBitVectorToRegimes(const vector<bool>& r) const;
  };
 private:
@@ -713,7 +714,7 @@ public:
  {
    return tuple {static_only_equations, static_only_equations_lineno,
                  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
  bool ParamUsedWithLeadLag() const;

--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
 /* -*- C++ -*- */
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -70,7 +70,7 @@ string eofbuff;
 NAME [a-z_][a-z0-9_]*
 FLOAT_NUMBER ((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+)
-DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
+DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|[sh][12])
 %%
 /* Code put at the beginning of yylex() */
@@ -236,6 +236,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>pac_target_info {BEGIN DYNARE_BLOCK; return token::PAC_TARGET_INFO;}
 <INITIAL>matched_irfs {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS;}
 <INITIAL>matched_irfs_weights {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS_WEIGHTS;}
+<INITIAL>perfect_foresight_controlled_paths {BEGIN DYNARE_BLOCK; return token::PERFECT_FORESIGHT_CONTROLLED_PATHS;}
 /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
@@ -256,11 +257,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>prior_function {BEGIN DYNARE_STATEMENT; return token::PRIOR_FUNCTION;}
 <INITIAL>posterior_function {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_FUNCTION;}
- /* Inside  of a Dynare statement */
+<DYNARE_STATEMENT,DYNARE_BLOCK>{DATE} {
-<DYNARE_STATEMENT>{DATE} {
  yylval->emplace<string>(yytext);
  return token::DATE;
 }
+ /* Inside a Dynare statement */
 <DYNARE_STATEMENT>file                  {return token::FILE;}
 <DYNARE_STATEMENT>datafile 		{return token::DATAFILE;}
 <DYNARE_STATEMENT>dirname       {return token::DIRNAME;}
@@ -404,6 +406,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>logarithmic_reduction {return token::LOGARITHMIC_REDUCTION;}
 <DYNARE_STATEMENT>use_univariate_filters_if_singularity_is_detected {return token::USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED;}
 <DYNARE_STATEMENT>hybrid {return token::HYBRID;}
+<DYNARE_STATEMENT>use_first_order_solution {return token::USE_FIRST_ORDER_SOLUTION;}
 <DYNARE_STATEMENT>default {return token::DEFAULT;}
 <DYNARE_STATEMENT>init2shocks {return token::INIT2SHOCKS;}
@@ -416,6 +419,10 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>kitagawa {return token::KITAGAWA;}
 <DYNARE_STATEMENT>smooth {return token::SMOOTH;}
 <DYNARE_STATEMENT>stratified {return token::STRATIFIED;}
+<DYNARE_STATEMENT>residual {
+  yylval->emplace<string>(yytext);
+  return token::RESIDUAL;
+}
 <DYNARE_STATEMENT>cpf_weights {return token::CPF_WEIGHTS;}
 <DYNARE_STATEMENT>amisanotristani {return token::AMISANOTRISTANI;}
 <DYNARE_STATEMENT>murrayjonesparslow {return token::MURRAYJONESPARSLOW;}
@@ -464,10 +471,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
  yylval->emplace<string>(yytext);
  return token::DIRICHLET;
 }
-<DYNARE_STATEMENT>weibull {
+<DYNARE_STATEMENT>weibull {return token::WEIBULL;}
-  yylval->emplace<string>(yytext);
-  return token::WEIBULL;
-}
 <DYNARE_STATEMENT>normal {
  yylval->emplace<string>(yytext);
  return token::NORMAL;
@@ -600,18 +604,9 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>substitute_libs {return token::SUBSTITUTE_LIBS;}
 <DYNARE_STATEMENT>compiler {return token::COMPILER;}
 <DYNARE_STATEMENT>instruments {return token::INSTRUMENTS;}
-<DYNARE_STATEMENT>hessian  {
+<DYNARE_STATEMENT>hessian  {return token::HESSIAN;}
-  yylval->emplace<string>(yytext);
+<DYNARE_STATEMENT>prior_variance  {return token::PRIOR_VARIANCE;}
-  return token::HESSIAN;
+<DYNARE_STATEMENT>identity_matrix  {return token::IDENTITY_MATRIX;}
-}
-<DYNARE_STATEMENT>prior_variance  {
-  yylval->emplace<string>(yytext);
-  return token::PRIOR_VARIANCE;
-}
-<DYNARE_STATEMENT>identity_matrix  {
-  yylval->emplace<string>(yytext);
-  return token::IDENTITY_MATRIX;
-}
 <DYNARE_STATEMENT>mcmc_jumping_covariance {return token::MCMC_JUMPING_COVARIANCE;}
 /* These four (var, varexo, varexo_det, parameters) are for change_type */
@@ -711,14 +706,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>lmmcp {return token::LMMCP;}
 <DYNARE_STATEMENT>additional_optimizer_steps	{return token::ADDITIONAL_OPTIMIZER_STEPS;}
 <DYNARE_STATEMENT>bartlett_kernel_lag {return token::BARTLETT_KERNEL_LAG; }
-<DYNARE_STATEMENT>optimal {
-  yylval->emplace<string>(yytext);
-  return token::OPTIMAL;
-}
-<DYNARE_STATEMENT>diagonal  {
-  yylval->emplace<string>(yytext);
-  return token::DIAGONAL;
-}
 <DYNARE_STATEMENT>gmm {return token::GMM;}
 <DYNARE_STATEMENT>smm {return token::SMM;}
 <DYNARE_STATEMENT>irf_matching {return token::IRF_MATCHING;}
@@ -845,6 +832,8 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
  return token::DD;
 }
 <DYNARE_BLOCK>weights {return token::WEIGHTS;}
+<DYNARE_BLOCK>exogenize {return token::EXOGENIZE;}
+<DYNARE_BLOCK>endogenize {return token::ENDOGENIZE;}
 /* Inside Dynare statement */
 <DYNARE_STATEMENT>solve_algo {return token::SOLVE_ALGO;}

--- a/src/EquationTags.cc
+++ b/src/EquationTags.cc
 /*
- * Copyright © 2020-2023 Dynare Team
+ * Copyright © 2020-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -74,13 +74,12 @@ EquationTags::erase(const set<int>& eqns, const map<int, int>& old_eqn_num_2_new
    eqn_tags.erase(eqn);
  for (const auto& [oldeqn, neweqn] : old_eqn_num_2_new)
-    for (auto& [eqn, tags] : eqn_tags)
+    if (eqn_tags.contains(oldeqn))
-      if (eqn == oldeqn)
+      {
-        {
+        auto tmp = eqn_tags.extract(oldeqn);
-          auto tmp = eqn_tags.extract(eqn);
+        tmp.key() = neweqn;
-          tmp.key() = neweqn;
+        eqn_tags.insert(move(tmp));
-          eqn_tags.insert(move(tmp));
+      }
-        }
 }
 void

--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -2978,10 +2978,8 @@ UnaryOpNode::writeJsonOutput(ostream& output, const temporary_terms_t& temporary
      output << "])";
      return;
    case UnaryOpcode::steadyState:
-      output << "(";
+      output << "STEADY_STATE";
-      arg->writeJsonOutput(output, temporary_terms, tef_terms, isdynamic);
+      break;
-      output << ")";
-      return;
    case UnaryOpcode::steadyStateParamDeriv:
      {
        auto varg = dynamic_cast<VariableNode*>(arg);
@@ -6458,9 +6456,9 @@ TrinaryOpNode::eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v
  switch (op_code)
    {
    case TrinaryOpcode::normcdf:
-      return (0.5 * (1 + erf((v1 - v2) / v3 / numbers::sqrt2)));
+      return 0.5 * (1 + erf((v1 - v2) / v3 / numbers::sqrt2));
    case TrinaryOpcode::normpdf:
-      return (1 / (v3 * sqrt(2 * numbers::pi) * exp(pow((v1 - v2) / v3, 2) / 2)));
+      return 1 / (v3 * sqrt(2 * numbers::pi) * exp(pow((v1 - v2) / v3, 2) / 2));
    }
  __builtin_unreachable(); // Silence GCC warning
 }
@@ -6989,9 +6987,9 @@ TrinaryOpNode::isVariableNodeEqualTo([[maybe_unused]] SymbolType type_arg,
 bool
 TrinaryOpNode::containsPacExpectation(const string& pac_model_name) const
 {
-  return (arg1->containsPacExpectation(pac_model_name)
+  return arg1->containsPacExpectation(pac_model_name)
-          || arg2->containsPacExpectation(pac_model_name)
+         || arg2->containsPacExpectation(pac_model_name)
-          || arg3->containsPacExpectation(pac_model_name));
+         || arg3->containsPacExpectation(pac_model_name);
 }
 bool
@@ -7907,7 +7905,7 @@ ExternalFunctionNode::sameTefTermPredicate() const
 {
  return [this](expr_t e) {
    auto e2 = dynamic_cast<ExternalFunctionNode*>(e);
-    return (e2 != nullptr && e2->symb_id == symb_id && e2->arguments == arguments);
+    return e2 != nullptr && e2->symb_id == symb_id && e2->arguments == arguments;
  };
 }
@@ -8245,12 +8243,12 @@ FirstDerivExternalFunctionNode::sameTefTermPredicate() const
  if (first_deriv_symb_id == symb_id)
    return [this](expr_t e) {
      auto e2 = dynamic_cast<ExternalFunctionNode*>(e);
-      return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
+      return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
    };
  else
    return [this](expr_t e) {
      auto e2 = dynamic_cast<FirstDerivExternalFunctionNode*>(e);
-      return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
+      return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
    };
 }
@@ -8600,12 +8598,12 @@ SecondDerivExternalFunctionNode::sameTefTermPredicate() const
  if (second_deriv_symb_id == symb_id)
    return [this](expr_t e) {
      auto e2 = dynamic_cast<ExternalFunctionNode*>(e);
-      return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
+      return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
    };
  else
    return [this](expr_t e) {
      auto e2 = dynamic_cast<SecondDerivExternalFunctionNode*>(e);
-      return (e2 && e2->symb_id == symb_id && e2->arguments == arguments);
+      return e2 && e2->symb_id == symb_id && e2->arguments == arguments;
    };
 }
@@ -9491,9 +9489,9 @@ ExprNode::matchParamTimesTargetMinusVariable(int symb_id) const
        auto& avi = datatree.symbol_table.getAuxVarInfo(target->symb_id);
        if (avi.type == AuxVarType::pacTargetNonstationary && target->lag == -1)
          return true;
-        return (avi.type == AuxVarType::unaryOp && avi.unary_op == "log" && avi.orig_symb_id
+        return avi.type == AuxVarType::unaryOp && avi.unary_op == "log" && avi.orig_symb_id
-                && !datatree.symbol_table.isAuxiliaryVariable(*avi.orig_symb_id)
+               && !datatree.symbol_table.isAuxiliaryVariable(*avi.orig_symb_id)
-                && target->lag + avi.orig_lead_lag.value() == -1);
+               && target->lag + avi.orig_lead_lag.value() == -1;
      }
    else
      return target->lag == -1;
@@ -9572,13 +9570,14 @@ ExprNode::toString() const
 }
 tuple<int, expr_t, expr_t>
-ExprNode::matchComplementarityCondition() const
+ExprNode::matchComplementarityCondition(
+    [[maybe_unused]] const optional<int>& heterogeneity_dimension) const
 {
  throw MatchFailureException {"This expression is not an inequality"};
 }
 tuple<int, expr_t, expr_t>
-BinaryOpNode::matchComplementarityCondition() const
+BinaryOpNode::matchComplementarityCondition(const optional<int>& heterogeneity_dimension) const
 {
  bool is_greater {[&] {
    switch (op_code)
@@ -9596,7 +9595,13 @@ BinaryOpNode::matchComplementarityCondition() const
  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)
+    if (ve && ve->lag == 0
+        && ((!heterogeneity_dimension
+             && datatree.symbol_table.getType(ve->symb_id) == SymbolType::endogenous)
+            || (heterogeneity_dimension
+                && datatree.symbol_table.getType(ve->symb_id) == SymbolType::heterogeneousEndogenous
+                && datatree.symbol_table.getHeterogeneityDimension(ve->symb_id)
+                       == *heterogeneity_dimension)))
      return ve->symb_id;
    else
      return nullopt;
@@ -9633,11 +9638,11 @@ BinaryOpNode::matchComplementarityCondition() const
            || (!is_greater
                && (barg1->op_code == BinaryOpcode::less
                    || barg1->op_code == BinaryOpcode::lessEqual)))))
-    throw MatchFailureException {"Complementarity condition does not have the right form"};
+    throw MatchFailureException {};
  auto id = match_contemporaneous_endogenous(barg1->arg2);
  if (!id)
-    throw MatchFailureException {"Complementarity condition does not have the right form"};
+    throw MatchFailureException {};
  check_bound_constant(barg1->arg1);
  check_bound_constant(arg2);

--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -283,7 +283,7 @@ protected:
  min_cost(bool is_matlab)
  {
    return is_matlab ? min_cost_matlab : min_cost_c;
-  };
+  }
  //! Initializes data member non_null_derivatives
  virtual void prepareForDerivation() = 0;
@@ -943,7 +943,8 @@ public:
  /* 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;
+  [[nodiscard]] virtual tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension = nullopt) const;
  /* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
     Also declares those aggregation operators in the HeterogeneityTable, so as to
@@ -1519,7 +1520,9 @@ public:
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
                                                      vector<BinaryOpNode*>& neweqs) const override;
-  [[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
+  [[nodiscard]] tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension
+                                = nullopt) const override;
 };
 //! Trinary operator node

--- a/src/HeterogeneityTable.hh
+++ b/src/HeterogeneityTable.hh
@@ -90,7 +90,7 @@ public:
  };
  void addSummedHeterogeneousEndogenous(int symb_id);
-  int getSummedHeterogenousEndogenousIndex(int symb_id) const;
+  [[nodiscard]] int getSummedHeterogenousEndogenousIndex(int symb_id) const;
  [[nodiscard]] int aggregateEndoSize() const;
  void writeOutput(ostream& output) const;

--- a/src/HeterogeneousModel.cc
+++ b/src/HeterogeneousModel.cc
 /*
- * Copyright © 2024 Dynare Team
+ * Copyright © 2024-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -33,14 +33,6 @@ HeterogeneousModel::HeterogeneousModel(SymbolTable& symbol_table_arg,
 {
 }
-HeterogeneousModel::HeterogeneousModel(const HeterogeneousModel& m) :
-    ModelTree {m},
-    heterogeneity_dimension {m.heterogeneity_dimension},
-    deriv_id_table {m.deriv_id_table},
-    inv_deriv_id_table {m.inv_deriv_id_table}
-{
-}
 HeterogeneousModel&
 HeterogeneousModel::operator=(const HeterogeneousModel& m)
 {
@@ -123,14 +115,7 @@ HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_d
  computeTemporaryTerms(!use_dll, no_tmp_terms);
-  if (ranges::any_of(complementarity_conditions, [](const auto& x) { return x.has_value(); }))
+  computeMCPEquationsReordering(heterogeneity_dimension);
-    {
-      // Implementing it requires modifications in ModelTree::computeMCPEquationsReordering()
-      cerr << "ERROR: Complementarity conditions are not yet implemented in "
-              "model(heterogeneity=...) blocks"
-           << endl;
-      exit(EXIT_FAILURE);
-    }
 }
 void
@@ -138,6 +123,7 @@ HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
 {
  assert(!julia); // Not yet implemented
  writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
+  writeComplementarityConditionsFile<true>(basename, heterogeneity_dimension);
 }
 int
@@ -245,4 +231,9 @@ HeterogeneousModel::writeDriverOutput(ostream& output) const
  output << "];" << endl;
  writeDriverSparseIndicesHelper(
      "heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").dynamic", output);
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1
+         << ").dynamic_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
 }
--- a/src/HeterogeneousModel.hh
+++ b/src/HeterogeneousModel.hh
@@ -35,7 +35,7 @@ public:
                     ExternalFunctionsTable& external_functions_table_arg,
                     HeterogeneityTable& heterogeneity_table_arg, int heterogeneity_dimension_arg);
-  HeterogeneousModel(const HeterogeneousModel& m);
+  HeterogeneousModel(const HeterogeneousModel& m) = default;
  HeterogeneousModel& operator=(const HeterogeneousModel& m);
  void computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll);

--- a/src/ModFile.cc
+++ b/src/ModFile.cc
 /*
- * Copyright © 2006-2024 Dynare Team
+ * Copyright © 2006-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -1137,7 +1137,8 @@ ModFile::writeMOutput(const string& basename, bool clear_all, bool clear_global,
              << "M_.heteroskedastic_shocks.Qvalue_orig = [];" << endl
              << "M_.heteroskedastic_shocks.Qscale_orig = [];" << endl
              << "M_.matched_irfs = {};" << endl
-              << "M_.matched_irfs_weights = {};" << endl;
+              << "M_.matched_irfs_weights = {};" << endl
+              << "M_.perfect_foresight_controlled_paths = [];" << endl;
  // NB: options_.{ramsey,discretionary}_policy should rather be fields of M_
  mOutputFile << boolalpha << "options_.linear = " << linear << ";" << endl

--- a/src/ModelEquationBlock.hh
+++ b/src/ModelEquationBlock.hh
 /*
- * Copyright © 2010-2024 Dynare Team
+ * Copyright © 2010-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -62,6 +62,17 @@ protected:
  {
    return "original Ramsey model";
  }
+  int
+  getJacobianCol(int deriv_id, [[maybe_unused]] bool sparse) const override
+  {
+    /* Override the DynamicModel method by returning a dummy Jacobian column number.
+       The override is necessary because the method from DynamicModel fails with
+       endos with lag/lead greater than 1 or exos with a lag/lead, while substitutions
+       are by definition not done for an original model.
+       In particular, this fixes dynare#1960 (equation derivatives are computed for models declared
+       as linear, to check whether they are truly linear). */
+    return deriv_id;
+  }
 };
 class SteadyStateModel : public DataTree

--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -895,7 +895,7 @@ ModelTree::determineLinearBlocks()
 int
 ModelTree::equation_number() const
 {
-  return (equations.size());
+  return equations.size();
 }
 void
@@ -1398,8 +1398,7 @@ ModelTree::writeLatexModelFile(const string& mod_basename, const string& latex_b
      content_output << endl << R"(\end{dmath})" << endl;
    }
-  output << R"(\include{)" << latex_basename + "_content"
+  output << R"(\include{)" << latex_basename + "_content" << "}" << endl
-         << "}" << endl
         << R"(\end{document})" << endl;
  output.close();
@@ -2119,11 +2118,13 @@ ModelTree::writeAuxVarRecursiveDefinitions(ostream& output, ExprNodeOutputType o
 }
 void
-ModelTree::computeMCPEquationsReordering()
+ModelTree::computeMCPEquationsReordering(const optional<int>& heterogeneous_dimension)
 {
  /* Optimal policy models (discretionary, or Ramsey before computing FOCs) do not have as many
     equations as variables. Do not even try to compute the reordering. */
-  if (static_cast<int>(equations.size()) != symbol_table.endo_nbr())
+  if (static_cast<int>(equations.size())
+      != (heterogeneous_dimension ? symbol_table.het_endo_nbr(*heterogeneous_dimension)
+                                  : symbol_table.endo_nbr()))
    return;
  assert(equations.size() == complementarity_conditions.size());
No results found