.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

.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

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_");
-  output << "display_static_residuals(M_, options_, oo_, options_resid_);" << endl;
+  options_list.writeOutput(output);
+  output << "display_static_residuals(M_, options_, oo_);" << endl;
 }
 
 void

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;

src/Configuration.hh

@@ -47,7 +47,7 @@ private:
     get_paths() const
     {
       return paths;
-    };
+    }
 
   private:
     map<string, vector<string>> paths;

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&

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)

src/DynamicModel.cc

 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -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

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;

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>{DATE} {
+<DYNARE_STATEMENT,DYNARE_BLOCK>{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 {
-  yylval->emplace<string>(yytext);
-  return token::WEIBULL;
-}
+<DYNARE_STATEMENT>weibull {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  {
-  yylval->emplace<string>(yytext);
-  return token::HESSIAN;
-}
-<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>hessian  {return token::HESSIAN;}
+<DYNARE_STATEMENT>prior_variance  {return token::PRIOR_VARIANCE;}
+<DYNARE_STATEMENT>identity_matrix  {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;}

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 << "(";
-      arg->writeJsonOutput(output, temporary_terms, tef_terms, isdynamic);
-      output << ")";
-      return;
+      output << "STEADY_STATE";
+      break;
     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)
-          || arg2->containsPacExpectation(pac_model_name)
-          || arg3->containsPacExpectation(pac_model_name));
+  return arg1->containsPacExpectation(pac_model_name)
+         || arg2->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
-                && !datatree.symbol_table.isAuxiliaryVariable(*avi.orig_symb_id)
-                && target->lag + avi.orig_lead_lag.value() == -1);
+        return avi.type == AuxVarType::unaryOp && avi.unary_op == "log" && avi.orig_symb_id
+               && !datatree.symbol_table.isAuxiliaryVariable(*avi.orig_symb_id)
+               && target->lag + avi.orig_lead_lag.value() == -1;
       }
     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);
 

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

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;

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(); }))
-    {
-      // Implementing it requires modifications in ModelTree::computeMCPEquationsReordering()
-      cerr << "ERROR: Complementarity conditions are not yet implemented in "
-              "model(heterogeneity=...) blocks"
-           << endl;
-      exit(EXIT_FAILURE);
-    }
+  computeMCPEquationsReordering(heterogeneity_dimension);
 }
 
 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;
 }

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);

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

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

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"
-         << "}" << endl
+  output << R"(\include{)" << latex_basename + "_content" << "}" << 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());

src/ModelTree.hh

 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -211,7 +211,7 @@ protected:
     getRecursiveSize() const
     {
       return size - mfs_size;
-    };
+    }
   };
 
   // Whether block decomposition has been successfully computed
@@ -296,7 +296,7 @@ protected:
 
   /* Computes the mcp_equations_reordering vector.
      Also checks that a variable does not appear as constrained in two different equations. */
-  void computeMCPEquationsReordering();
+  void computeMCPEquationsReordering(const optional<int>& heterogeneous_dimension = nullopt);
 
   //! Writes temporary terms
   template<ExprNodeOutputType output_type>
@@ -439,7 +439,9 @@ protected:
   void writeSetAuxiliaryVariablesFile(const string& basename, bool julia) const;
 
   template<bool dynamic>
-  void writeComplementarityConditionsFile(const string& basename) const;
+  void writeComplementarityConditionsFile(const string& basename,
+                                          const optional<int>& heterogeneous_dimension
+                                          = nullopt) const;
 
 private:
   //! Sparse matrix of double to store the values of the static Jacobian
@@ -559,7 +561,7 @@ protected:
   {
     return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
         .first;
-  };
+  }
   //! Return true if the equation has been normalized
   bool
   isBlockEquationRenormalized(int blk, int eq) const
@@ -567,44 +569,44 @@ protected:
     return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
                .first
            == EquationType::evaluateRenormalized;
-  };
+  }
   //! Return the expr_t of equation belonging to the block
   BinaryOpNode*
   getBlockEquationExpr(int blk, int eq) const
   {
     return equations[eq_idx_block2orig[blocks[blk].first_equation + eq]];
-  };
+  }
   //! Return the expr_t of renormalized equation belonging to the block
   BinaryOpNode*
   getBlockEquationRenormalizedExpr(int blk, int eq) const
   {
     return equation_type_and_normalized_equation[eq_idx_block2orig[blocks[blk].first_equation + eq]]
         .second;
-  };
+  }
   //! Return the original number of equation belonging to the block
   int
   getBlockEquationID(int blk, int eq) const
   {
     return eq_idx_block2orig[blocks[blk].first_equation + eq];
-  };
+  }
   //! Return the original number of variable belonging to the block
   int
   getBlockVariableID(int blk, int var) const
   {
     return endo_idx_block2orig[blocks[blk].first_equation + var];
-  };
+  }
   //! Return the position of an equation (given by its original index) inside its block
   int
   getBlockInitialEquationID(int blk, int eq) const
   {
     return eq_idx_orig2block[eq] - blocks[blk].first_equation;
-  };
+  }
   //! Return the position of a variable (given by its original index) inside its block
   int
   getBlockInitialVariableID(int blk, int var) const
   {
     return endo_idx_orig2block[var] - blocks[blk].first_equation;
-  };
+  }
   //! Initialize equation_reordered & variable_reordered
   void initializeVariablesAndEquations();
 
@@ -1004,9 +1006,8 @@ ModelTree::writeModelCFile(const string& basename, const string& mexext,
 
   const filesystem::path model_src_dir {filesystem::path {basename} / "model" / "src"};
 
-  auto [d_output, tt_output] = writeModelFileHelper < dynamic
-                                   ? ExprNodeOutputType::CDynamicModel
-                                   : ExprNodeOutputType::CStaticModel > ();
+  auto [d_output, tt_output] = writeModelFileHelper<dynamic ? ExprNodeOutputType::CDynamicModel
+                                                            : ExprNodeOutputType::CStaticModel>();
   vector<filesystem::path> header_files, object_files;
 
   // TODO: when C++20 support is complete, mark the following strings constexpr
@@ -2376,9 +2377,9 @@ ModelTree::writeSparseModelJuliaFiles(const string& basename) const
 {
   assert(heterogeneity_table.empty());
 
-  auto [d_sparse_output, tt_sparse_output] = writeModelFileHelper < dynamic
-                                                 ? ExprNodeOutputType::juliaSparseDynamicModel
-                                                 : ExprNodeOutputType::juliaSparseStaticModel > ();
+  auto [d_sparse_output, tt_sparse_output]
+      = writeModelFileHelper<dynamic ? ExprNodeOutputType::juliaSparseDynamicModel
+                                     : ExprNodeOutputType::juliaSparseStaticModel>();
 
   filesystem::path julia_dir {filesystem::path {basename} / "model" / "julia"};
   // TODO: when C++20 support is complete, mark the following strings constexpr
@@ -3154,10 +3155,13 @@ ModelTree::writeSetAuxiliaryVariablesFile(const string& basename, bool julia) co
 
 template<bool dynamic>
 void
-ModelTree::writeComplementarityConditionsFile(const string& basename) const
+ModelTree::writeComplementarityConditionsFile(const string& basename,
+                                              const optional<int>& heterogeneous_dimension) const
 {
-  // TODO: when C++20 support is complete, mark the following string constexpr
-  const string funcname {(dynamic ? "dynamic"s : "static"s) + "_complementarity_conditions"};
+  const string funcname {
+      (dynamic ? "dynamic"s : "static"s)
+      + (heterogeneous_dimension ? "_het"s + to_string(*heterogeneous_dimension + 1) : ""s)
+      + "_complementarity_conditions"};
   const filesystem::path filename {packageDir(basename) / (funcname + ".m")};
   /* Can’t use matlabOutsideModel for output type, since it uses M_.
      Static is ok even for the dynamic model, since there are no leads/lags. */