.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
-  MINGW32_BOOST_VERSION: 1.84.0-1
-  MINGW64_BOOST_VERSION: 1.84.0-1
+  MINGW64_BOOST_VERSION: 1.88.0-2
   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'
@@ -72,3 +71,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

doc/macroprocessor/macroprocessor.tex

@@ -15,7 +15,7 @@
 \pgfdeclareimage[height=0.8cm]{logo}{dlogo}
 \institute[Dynare Team]{\pgfuseimage{logo}}
 
-\date{23 May 2023}
+\date{31 May 2024}
 
 \AtBeginSection[]
 {
@@ -156,9 +156,9 @@
     \item comparison operators: \texttt{< > <= >= == !=}
     \item logical operators: \verb+&& || !+
     \item range with unit increment: \texttt{1:4} is equivalent to
-      real array \texttt{[1, 2, 3, 4]}. (NB: \texttt{[1:4]} is equivalent to an
+      real array \texttt{[1, 2, 3, 4]} \\ (NB: \texttt{[1:4]} is equivalent to an
       array containing an array of reals, \textit{i.e.} \texttt{[[1, 2, 3, 4]]})
-    \item range with user-defined increment: \texttt{4:-1.1:-1} is equivalent to real array \texttt{[4, 2.9, 1.8, 0.7, -0.4]}.
+    \item range with user-defined increment: \\ \texttt{4:-1.1:-1} is equivalent to real array \texttt{[4, 2.9, 1.8, 0.7, -0.4]}
     \end{itemize}
   \end{block}
 
@@ -314,7 +314,7 @@ Then \texttt{distance(3, 4)} will be equivalent to \texttt{5}.
 \end{frame}
 
 \begin{frame}[fragile=singleslide]
-  \frametitle{Defining macro-variables}
+  \frametitle{Defining macro-variables (1/2)}
 
   The value of a macro-variable can be defined with the \verb+@#define+
   directive.
@@ -335,9 +335,31 @@ Then \texttt{distance(3, 4)} will be equivalent to \texttt{5}.
 @#define t = ("US" in w)    // Equals true
 \end{verbatim}
   \end{block}
-  NB: You can define macro variables on the Dynare command line by using the \texttt{-D} option
 \end{frame}
 
+\begin{frame}[fragile=singleslide]
+  \frametitle{Defining macro-variables (2/2)}
+  Macro variables can also be defined on the Dynare command line by using the
+  \texttt{-D} option, for easily switching between different flavours of a model.
+
+  \begin{block}{Example 1}
+\begin{verbatim}
+dynare myfile.mod -Dx=5
+\end{verbatim}
+The macro-variable \texttt{x} will be equal to \texttt{5} when running \texttt{myfile.mod}.
+\end{block}
+
+\begin{block}{Example 2}
+  Use single quotes around the \texttt{-D} option when there are spaces or
+  special characters in the variable definition.
+\begin{verbatim}
+dynare myfile.mod '-DA=[ i in [1,2,3] when i > 1 ]'
+\end{verbatim}
+The macro-variable \texttt{A} will be equal to \texttt{[2,3]} when running \texttt{myfile.mod}.
+\end{block}
+\end{frame}
+
+
 \begin{frame}[fragile=singleslide]
   \frametitle{Expression substitution}
   \framesubtitle{Dummy example}
@@ -611,7 +633,7 @@ end;
 \begin{frame}
   \frametitle{Macro-related command line options}
   \begin{itemize}
-  \item \texttt{savemacro}: Useful for debugging or learning purposes, saves the output of the macro processor. If your \texttt{.mod} file is called \texttt{file.mod}, the output is saved to \texttt{file-macroexp.mod}.
+  \item \texttt{savemacro}: Useful for debugging or learning purposes, saves the output of the macro processor. If your \texttt{.mod} file is called \texttt{file.mod}, the output is saved to \texttt{file\_macroexp.mod}.
   \item NB: \texttt{savemacro=filename} allows a user-defined file name
   \item \texttt{linemacro}: In the output of \texttt{savemacro}, print line numbers where the macro directives were placed.
   \item \texttt{onlymacro}: Stops processing after the macro processing step.
@@ -866,7 +888,7 @@ rhos = [ 0.8, 0.9, 1];
     \ccbysa
     \column{0.71\textwidth}
     \tiny
-    Copyright © 2008-2023 Dynare Team \\
+    Copyright © 2008-2024 Dynare Team \\
     License: \href{http://creativecommons.org/licenses/by-sa/4.0/}{Creative
       Commons Attribution-ShareAlike 4.0}
   \end{columns}

meson.build

@@ -2,7 +2,7 @@
 # It is not used when building Dynare as a whole.
 
 project('dynare-preprocessor', 'cpp',
-        version : '6-unstable',
+        version : '7-unstable',
         # NB: update C++ standard in .clang-format whenever the following is modified
         default_options : [ 'cpp_std=gnu++20', 'warning_level=2' ],
         meson_version : '>=0.64.0')

scripts/homebrew-native-arm64.ini

 # Meson native file for compiling under Homebrew / arm64
 
 [binaries]
-cpp = '/opt/homebrew/bin/g++-13'
+cpp = '/opt/homebrew/bin/g++-14'
 flex = '/opt/homebrew/opt/flex/bin/flex'
 bison = '/opt/homebrew/opt/bison/bin/bison'
 

scripts/homebrew-native-x86_64.ini

 # Meson native file for compiling under Homebrew / x86-64
 
 [binaries]
-cpp = '/usr/local/bin/g++-13'
+cpp = '/usr/local/bin/g++-14'
 flex = '/usr/local/opt/flex/bin/flex'
 bison = '/usr/local/opt/bison/bin/bison'
 

scripts/wasm.ini

+# Meson cross file for creating a WebAssembly version of the preprocessor.
+#
+# Requires emscripten to be installed.
+# Was successfully tested with emscripten 3.1.69 installed through emsdk
+# tool, as described on: https://emscripten.org/docs/getting_started/downloads.html
+# Don’t forget to source script snippet in current shell before running meson.
+#
+# Compilation creates a .wasm and .js wrapper under <builddir>/src/
+#
+# Can be run locally with node.js using:
+#  node dynare-preprocessor.js file.mod
+# NB: a version of node.js is shipped with emscripten (under the node/
+# subdirectory), but another version should also work.
+
+[binaries]
+cpp = 'em++'
+
+[host_machine]
+system = 'emscripten'
+# Could be changed to wasm64 if 4GB memory constraint is hit
+# Some background: https://v8.dev/blog/4gb-wasm-memory
+cpu_family = 'wasm32'
+cpu = 'wasm32'
+endian = 'little'
+
+[built-in options]
+# Never do a debug build, because otherwise the lack of optimisations can
+# overflow the memory capacities.
+buildtype = 'release'
+# The -fexceptions flag (for both compilation and linking) is needed for an
+# unknown reason (C++ compilers are supposed to always add exception support).
+# The -Wno-unqualified-std-cast-call flag removes many warnings about “move”
+# not being qualified with “std::” namespace.
+# The -fexperimental-library flag is needed to get std::jthread support (it was
+# supposed to no longer be necessary for LLVM 20, but for some reason we still
+# need it).
+cpp_args = [ '-fexceptions', '-Wno-unqualified-std-cast-call', '-fexperimental-library' ]
+# NODERAWFS=1 is needed for accessing the local filesystem
+cpp_link_args = [ '-s', 'NODERAWFS=1', '-fexceptions' ]
+
+[properties]
+# It’s necessary to use a different copy of Boost than the one under
+# /usr/include, because otherwise GCC headers confuse Clang
+boost_root = '/tmp/boost_1_86_0'

src/Bytecode.cc

 /*
- * Copyright © 2022-2023 Dynare Team
+ * Copyright © 2022-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -82,8 +82,8 @@ operator<<(Writer& code_file, const FBEGINBLOCK& instr)
   write_member(instr.type);
   for (int i = 0; i < instr.size; i++)
     {
-      write_member(instr.variable[i]);
-      write_member(instr.equation[i]);
+      write_member(instr.variables[i]);
+      write_member(instr.equations[i]);
     }
   if (instr.type == BlockSimulationType::solveTwoBoundariesSimple
       || instr.type == BlockSimulationType::solveTwoBoundariesComplete
@@ -97,8 +97,8 @@ operator<<(Writer& code_file, const FBEGINBLOCK& instr)
   write_member(instr.det_exo_size);
   write_member(instr.exo_size);
 
-  for_each_n(instr.det_exogenous.begin(), instr.det_exo_size, write_member);
-  for_each_n(instr.exogenous.begin(), instr.exo_size, write_member);
+  ranges::for_each_n(instr.det_exogenous.begin(), instr.det_exo_size, write_member);
+  ranges::for_each_n(instr.exogenous.begin(), instr.exo_size, write_member);
 
   return code_file;
 }

src/Bytecode.hh

 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -20,6 +20,7 @@
 #ifndef BYTECODE_HH
 #define BYTECODE_HH
 
+#include <concepts>
 #include <filesystem>
 #include <fstream>
 #include <ios>
@@ -74,9 +75,8 @@ enum class Tag
   FLDR,  // Loads a residual onto the stack
   FSTPR, // Stores a residual from the stack
 
-  FSTPG,  // Stores a derivative from the stack
-  FSTPG2, // Stores a derivative matrix for a static model from the stack
-  FSTPG3, // Stores a derivative matrix for a dynamic model from the stack
+  FSTPG,  // Stores the derivative of a simple (single-equation) block in simulate mode
+  FSTPG2, // Stores the derivative matrix of a block in evaluate mode
 
   FUNARY,   // A unary operator
   FBINARY,  // A binary operator
@@ -122,11 +122,6 @@ enum class ExternalFunctionCallType
   numericalSecondDerivative
 };
 
-struct Block_contain_type
-{
-  int Equation, Variable, Own_Derivative;
-};
-
 class Writer;
 
 struct Instruction
@@ -145,6 +140,9 @@ protected:
   ~Instruction() = default;
 };
 
+template<typename T>
+concept IsInstruction = derived_from<T, Instruction>;
+
 struct FLDZ final : public Instruction
 {
   FLDZ() : Instruction {Tag::FLDZ}
@@ -279,8 +277,7 @@ struct FSTPSU final : public Instruction
 
 struct FSTPG final : public Instruction
 {
-  const int pos;
-  explicit FSTPG(int pos_arg) : Instruction {Tag::FSTPG}, pos {pos_arg}
+  explicit FSTPG() : Instruction {Tag::FSTPG}
   {
   }
 };
@@ -293,15 +290,6 @@ struct FSTPG2 final : public Instruction
   }
 };
 
-struct FSTPG3 final : public Instruction
-{
-  const int row, col, lag, col_pos;
-  FSTPG3(int row_arg, int col_arg, int lag_arg, int col_pos_arg) :
-      Instruction {Tag::FSTPG3}, row {row_arg}, col {col_arg}, lag {lag_arg}, col_pos {col_pos_arg}
-  {
-  }
-};
-
 struct FUNARY final : public Instruction
 {
   const UnaryOpcode op_code;
@@ -442,7 +430,7 @@ struct FSTPV final : public Instruction
 
 class FCALL final : public Instruction
 {
-  template<typename B>
+  template<IsInstruction B>
   friend Writer& operator<<(Writer& code_file, const B& instr);
 
 private:
@@ -618,18 +606,17 @@ public:
 
 class FBEGINBLOCK final : public Instruction
 {
-  template<typename B>
+  template<IsInstruction B>
   friend Writer& operator<<(Writer& code_file, const B& instr);
 
 private:
   int size {0};
   BlockSimulationType type;
-  vector<int> variable;
-  vector<int> equation;
+  vector<int> variables;
+  vector<int> equations;
   vector<int> exogenous;
   vector<int> det_exogenous;
   bool is_linear {false};
-  vector<Block_contain_type> Block_Contain_;
   int u_count_int {0};
   int nb_col_jacob {0};
   int det_exo_size, exo_size;
@@ -639,16 +626,16 @@ public:
      sense when there is no block-decomposition, since there is no provision for
      derivatives w.r.t. endogenous not belonging to the block) */
   FBEGINBLOCK(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
-              const vector<int>& variable_arg, const vector<int>& equation_arg, bool is_linear_arg,
-              int u_count_int_arg, int nb_col_jacob_arg, int det_exo_size_arg, int exo_size_arg,
-              vector<int> det_exogenous_arg, vector<int> exogenous_arg) :
+              const vector<int>& variables_arg, const vector<int>& equations_arg,
+              bool is_linear_arg, int u_count_int_arg, int nb_col_jacob_arg, int det_exo_size_arg,
+              int exo_size_arg, vector<int> det_exogenous_arg, vector<int> exogenous_arg) :
       Instruction {Tag::FBEGINBLOCK},
       size {size_arg},
       type {type_arg},
-      variable {variable_arg.begin() + first_element,
-                variable_arg.begin() + (first_element + block_size)},
-      equation {equation_arg.begin() + first_element,
-                equation_arg.begin() + (first_element + block_size)},
+      variables {variables_arg.begin() + first_element,
+                 variables_arg.begin() + (first_element + block_size)},
+      equations {equations_arg.begin() + first_element,
+                 equations_arg.begin() + (first_element + block_size)},
       exogenous {move(exogenous_arg)},
       det_exogenous {move(det_exogenous_arg)},
       is_linear {is_linear_arg},
@@ -660,15 +647,15 @@ public:
   }
   // Constructor when derivatives w.r.t. exogenous are absent
   FBEGINBLOCK(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
-              const vector<int>& variable_arg, const vector<int>& equation_arg, bool is_linear_arg,
-              int u_count_int_arg, int nb_col_jacob_arg) :
+              const vector<int>& variables_arg, const vector<int>& equations_arg,
+              bool is_linear_arg, int u_count_int_arg, int nb_col_jacob_arg) :
       Instruction {Tag::FBEGINBLOCK},
       size {size_arg},
       type {type_arg},
-      variable {variable_arg.begin() + first_element,
-                variable_arg.begin() + (first_element + block_size)},
-      equation {equation_arg.begin() + first_element,
-                equation_arg.begin() + (first_element + block_size)},
+      variables {variables_arg.begin() + first_element,
+                 variables_arg.begin() + (first_element + block_size)},
+      equations {equations_arg.begin() + first_element,
+                 equations_arg.begin() + (first_element + block_size)},
       is_linear {is_linear_arg},
       u_count_int {u_count_int_arg},
       nb_col_jacob {nb_col_jacob_arg},
@@ -689,12 +676,12 @@ public:
 
     read_member(size);
     read_member(type);
+    variables.resize(size);
+    equations.resize(size);
     for (int i {0}; i < size; i++)
       {
-        Block_contain_type bc;
-        read_member(bc.Variable);
-        read_member(bc.Equation);
-        Block_Contain_.push_back(bc);
+        read_member(variables[i]);
+        read_member(equations[i]);
       }
     if (type == BlockSimulationType::solveTwoBoundariesSimple
         || type == BlockSimulationType::solveTwoBoundariesComplete
@@ -742,10 +729,15 @@ public:
   {
     return u_count_int;
   }
-  vector<Block_contain_type>
-  get_Block_Contain()
+  vector<int>
+  get_variables()
+  {
+    return variables;
+  }
+  vector<int>
+  get_equations()
   {
-    return Block_Contain_;
+    return equations;
   }
   int
   get_nb_col_jacob()
@@ -763,11 +755,6 @@ public:
     return det_exo_size;
   }
   vector<int>
-  get_endogenous()
-  {
-    return variable;
-  }
-  vector<int>
   get_exogenous()
   {
     return exogenous;
@@ -777,7 +764,7 @@ public:
 // Superclass of std::ofstream for writing a sequence of bytecode instructions
 class Writer : private ofstream
 {
-  template<typename B>
+  template<IsInstruction B>
   friend Writer& operator<<(Writer& code_file, const B& instr);
 
 private:
@@ -795,9 +782,8 @@ public:
   /* Overwrites an existing instruction, given its number.
      It is the responsibility of the caller to ensure that the new instruction
      occupies exactly as many bytes as the former one. */
-  template<typename B>
   void
-  overwriteInstruction(int instruction_number, const B& new_instruction)
+  overwriteInstruction(int instruction_number, const IsInstruction auto& new_instruction)
   {
     seekp(instructions_positions.at(instruction_number));
     *this << new_instruction;
@@ -808,7 +794,7 @@ public:
 
 // Overloads of operator<< for writing bytecode instructions
 
-template<typename B>
+template<IsInstruction B>
 Writer&
 operator<<(Writer& code_file, const B& instr)
 {

src/CommonEnums.hh

 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -25,14 +25,17 @@
  * command */
 enum class SymbolType
 {
-  endogenous = 0,            //!< Endogenous
-  exogenous = 1,             //!< Exogenous
-  exogenousDet = 2,          //!< Exogenous deterministic
-  parameter = 4,             //!< Parameter
-  modelLocalVariable = 10,   //!< Local variable whose scope is model (pound expression)
-  modFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
-  externalFunction = 12,     //!< External (user-defined) function
-  trend = 13,                //!< Trend variable
+  endogenous = 0,              // Endogenous (non-heterogeneous)
+  exogenous = 1,               // Exogenous (non-heterogeneous)
+  exogenousDet = 2,            // Exogenous deterministic (non-heterogeneous)
+  parameter = 4,               // Parameter (non-heterogeneous)
+  heterogeneousEndogenous = 5, // Endogenous that is heterogeneous across some dimension
+  heterogeneousExogenous = 6,  // Exogenous that is heterogeneous across some dimension
+  heterogeneousParameter = 7,  // Parameter that is heterogeneous across some dimension
+  modelLocalVariable = 10,     // Local variable whose scope is model (pound expression)
+  modFileLocalVariable = 11,   // Local variable whose scope is mod file (model excluded)
+  externalFunction = 12,       // External (user-defined) function
+  trend = 13,                  // Trend variable
   statementDeclaredVariable
   = 14, //!< Local variable assigned within a Statement (see subsample statement for example)
   logTrend = 15, //!< Log-trend variable
@@ -45,6 +48,13 @@ enum class SymbolType
   excludedVariable = 19 //!< Variable excluded via model_remove/var_remove/include_eqs/exclude_eqs
 };
 
+constexpr bool
+isHeterogeneous(SymbolType type)
+{
+  return type == SymbolType::heterogeneousEndogenous || type == SymbolType::heterogeneousExogenous
+         || type == SymbolType::heterogeneousParameter;
+}
+
 enum class UnaryOpcode
 {
   uminus,
@@ -75,7 +85,8 @@ enum class UnaryOpcode
   erf,
   erfc,
   diff,
-  adl
+  adl,
+  sum
 };
 
 enum class BinaryOpcode

src/ComputingTasks.cc

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -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
@@ -508,90 +510,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) :
@@ -1004,11 +922,9 @@ EstimationStatement::checkPass(ModFileStructure& mod_file_struct, WarningConsoli
   /* Check that we are not trying to estimate a parameter appearing in the
      planner discount factor (see dynare#1173) */
   vector<int> estimated_params_in_planner_discount;
-  set_intersection(mod_file_struct.estimated_parameters.begin(),
-                   mod_file_struct.estimated_parameters.end(),
-                   mod_file_struct.parameters_in_planner_discount.begin(),
-                   mod_file_struct.parameters_in_planner_discount.end(),
-                   back_inserter(estimated_params_in_planner_discount));
+  ranges::set_intersection(mod_file_struct.estimated_parameters,
+                           mod_file_struct.parameters_in_planner_discount,
+                           back_inserter(estimated_params_in_planner_discount));
   if (!estimated_params_in_planner_discount.empty())
     {
       cerr << "ERROR: It is not possible to estimate a parameter ("
@@ -1267,9 +1183,8 @@ AbstractEstimatedParamsStatement::commonCheckPass() const
       it.p4->collectVariables(SymbolType::parameter, used_params);
       it.jscale->collectVariables(SymbolType::parameter, used_params);
       vector<int> intersect;
-      set_intersection(declared_params.begin(), declared_params.end(), used_params.begin(),
-                       used_params.end(), back_inserter(intersect));
-      if (intersect.size() > 0)
+      ranges::set_intersection(declared_params, used_params, back_inserter(intersect));
+      if (!intersect.empty())
         {
           cerr << "ERROR: in `" << blockName() << "' block, the value of estimated parameter "
                << symbol_table.getName(intersect[0]) << " is used in the declaration for ";
@@ -1850,21 +1765,16 @@ DeterministicTrendsStatement::writeOutput(ostream& output, [[maybe_unused]] cons
                                           [[maybe_unused]] bool minimal_workspace) const
 {
   output << "options_.trend_coeff = {};" << endl;
-  for (const auto& trend_element : trend_elements)
-    {
-      SymbolType type = symbol_table.getType(trend_element.first);
-      if (type == SymbolType::endogenous)
-        {
-          output << "tmp1 = strmatch('" << trend_element.first << "',M_.endogenous_names,'exact');"
-                 << endl;
-          output << "options_.deterministic_trend_coeffs{tmp1} = '";
-          trend_element.second->writeOutput(output);
-          output << "';" << endl;
-        }
-      else
-        cerr << "Warning : Non-variable symbol used in deterministic_trends: "
-             << trend_element.first << endl;
-    }
+  for (const auto& [name, val] : trend_elements)
+    if (symbol_table.getType(name) == SymbolType::endogenous)
+      {
+        output << "tmp1 = strmatch('" << name << "',M_.endogenous_names,'exact');" << endl
+               << "options_.deterministic_trend_coeffs{tmp1} = '";
+        val->writeOutput(output);
+        output << "';" << endl;
+      }
+    else
+      cerr << "Warning: Non-variable symbol used in deterministic_trends: " << name << endl;
 }
 
 void
@@ -1872,20 +1782,17 @@ DeterministicTrendsStatement::writeJsonOutput(ostream& output) const
 {
   output << R"({"statementName": "deterministic_trends", )"
          << R"("trends" : {)";
-  for (bool printed_something {false}; const auto& trend_element : trend_elements)
-    {
-      if (symbol_table.getType(trend_element.first) == SymbolType::endogenous)
-        {
-          if (exchange(printed_something, true))
-            output << ", ";
-          output << R"(")" << trend_element.first << R"(": ")";
-          trend_element.second->writeJsonOutput(output, {}, {});
-          output << R"(")" << endl;
-        }
-      else
-        cerr << "Warning : Non-variable symbol used in deterministic_trends: "
-             << trend_element.first << endl;
-    }
+  for (bool printed_something {false}; const auto& [name, val] : trend_elements)
+    if (symbol_table.getType(name) == SymbolType::endogenous)
+      {
+        if (exchange(printed_something, true))
+          output << ", ";
+        output << R"(")" << name << R"(": ")";
+        val->writeJsonOutput(output, {}, {});
+        output << R"(")" << endl;
+      }
+    else
+      cerr << "Warning: Non-variable symbol used in deterministic_trends: " << name << endl;
   output << "}"
          << "}";
 }
@@ -1901,21 +1808,16 @@ ObservationTrendsStatement::writeOutput(ostream& output, [[maybe_unused]] const
                                         [[maybe_unused]] bool minimal_workspace) const
 {
   output << "options_.trend_coeff = {};" << endl;
-  for (const auto& trend_element : trend_elements)
-    {
-      SymbolType type = symbol_table.getType(trend_element.first);
-      if (type == SymbolType::endogenous)
-        {
-          output << "tmp1 = strmatch('" << trend_element.first << "',options_.varobs,'exact');"
-                 << endl;
-          output << "options_.trend_coeffs{tmp1} = '";
-          trend_element.second->writeOutput(output);
-          output << "';" << endl;
-        }
-      else
-        cerr << "Warning : Non-variable symbol used in observation_trends: " << trend_element.first
-             << endl;
-    }
+  for (const auto& [name, val] : trend_elements)
+    if (symbol_table.getType(name) == SymbolType::endogenous)
+      {
+        output << "tmp1 = strmatch('" << name << "',options_.varobs,'exact');" << endl
+               << "options_.trend_coeffs{tmp1} = '";
+        val->writeOutput(output);
+        output << "';" << endl;
+      }
+    else
+      cerr << "Warning: Non-variable symbol used in observation_trends: " << name << endl;
 }
 
 void
@@ -1923,20 +1825,17 @@ ObservationTrendsStatement::writeJsonOutput(ostream& output) const
 {
   output << R"({"statementName": "observation_trends", )"
          << R"("trends" : {)";
-  for (bool printed_something {false}; const auto& trend_element : trend_elements)
-    {
-      if (symbol_table.getType(trend_element.first) == SymbolType::endogenous)
-        {
-          if (exchange(printed_something, true))
-            output << ", ";
-          output << R"(")" << trend_element.first << R"(": ")";
-          trend_element.second->writeJsonOutput(output, {}, {});
-          output << R"(")" << endl;
-        }
-      else
-        cerr << "Warning : Non-variable symbol used in observation_trends: " << trend_element.first
-             << endl;
-    }
+  for (bool printed_something {false}; const auto& [name, val] : trend_elements)
+    if (symbol_table.getType(name) == SymbolType::endogenous)
+      {
+        if (exchange(printed_something, true))
+          output << ", ";
+        output << R"(")" << name << R"(": ")";
+        val->writeJsonOutput(output, {}, {});
+        output << R"(")" << endl;
+      }
+    else
+      cerr << "Warning: Non-variable symbol used in observation_trends: " << name << endl;
   output << "}"
          << "}";
 }
@@ -2014,7 +1913,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
@@ -2439,14 +2338,7 @@ void
 PlannerObjectiveStatement::writeOutput(ostream& output, const string& basename,
                                        [[maybe_unused]] bool minimal_workspace) const
 {
-  output << "M_.NNZDerivatives_objective = [";
-  for (int i = 1; i < static_cast<int>(model_tree->getNNZDerivatives().size()); i++)
-    output << (i > model_tree->getComputedDerivsOrder() ? -1 : model_tree->getNNZDerivatives()[i])
-           << ";";
-  output << "];" << endl << "M_.objective_tmp_nbr = [";
-  for (const auto& temporary_terms_derivative : model_tree->getTemporaryTermsDerivatives())
-    output << temporary_terms_derivative.size() << "; ";
-  output << "];" << endl;
+  model_tree->writeDriverOutput(output);
   model_tree->writeStaticFile(basename + ".objective", false, "", {}, false);
 }
 
@@ -3823,8 +3715,7 @@ SvarGlobalIdentificationCheckStatement::writeJsonOutput(ostream& output) const
   output << R"({"statementName": "svar_global_identification"})";
 }
 
-SetTimeStatement::SetTimeStatement(OptionsList options_list_arg) :
-    options_list {move(options_list_arg)}
+SetTimeStatement::SetTimeStatement(string period_arg) : period {move(period_arg)}
 {
 }
 
@@ -3832,19 +3723,13 @@ void
 SetTimeStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                               [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output);
+  output << "options_.initial_period = " << period << endl;
 }
 
 void
 SetTimeStatement::writeJsonOutput(ostream& output) const
 {
-  output << R"({"statementName": "set_time")";
-  if (!options_list.empty())
-    {
-      output << ", ";
-      options_list.writeJsonOutput(output);
-    }
-  output << "}";
+  output << R"({"statementName": "set_time", "period": ")" << period << R"("})";
 }
 
 EstimationDataStatement::EstimationDataStatement(OptionsList options_list_arg) :
@@ -4322,18 +4207,13 @@ BasicPriorStatement::checkPass([[maybe_unused]] ModFileStructure& mod_file_struc
 bool
 BasicPriorStatement::is_structural_innovation(const SymbolType symb_type) const
 {
-  if (symb_type == SymbolType::exogenous)
-    return true;
-  return false;
+  return symb_type == SymbolType::exogenous;
 }
 
 void
 BasicPriorStatement::get_base_name(const SymbolType symb_type, string& lhs_field) const
 {
-  if (symb_type == SymbolType::exogenous)
-    lhs_field = "structural_innovation";
-  else
-    lhs_field = "measurement_error";
+  lhs_field = (symb_type == SymbolType::exogenous ? "structural_innovation" : "measurement_error");
 }
 
 void
@@ -4692,10 +4572,7 @@ BasicOptionsStatement::is_structural_innovation(const SymbolType symb_type) cons
 void
 BasicOptionsStatement::get_base_name(const SymbolType symb_type, string& lhs_field) const
 {
-  if (symb_type == SymbolType::exogenous)
-    lhs_field = "structural_innovation";
-  else
-    lhs_field = "measurement_error";
+  lhs_field = (symb_type == SymbolType::exogenous ? "structural_innovation" : "measurement_error");
 }
 
 void
@@ -4903,10 +4780,7 @@ OptionsEqualStatement::writeJsonOutput(ostream& output) const
 void
 OptionsEqualStatement::get_base_name(const SymbolType symb_type, string& lhs_field) const
 {
-  if (symb_type == SymbolType::exogenous)
-    lhs_field = "structural_innovation";
-  else
-    lhs_field = "measurement_error";
+  lhs_field = (symb_type == SymbolType::exogenous ? "structural_innovation" : "measurement_error");
 }
 
 void
@@ -5468,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

 /*
- * 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:
@@ -546,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;
@@ -564,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;
@@ -579,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;
@@ -1029,10 +1007,10 @@ public:
 class SetTimeStatement : public Statement
 {
 private:
-  const OptionsList options_list;
+  const string period;
 
 public:
-  explicit SetTimeStatement(OptionsList options_list_arg);
+  explicit SetTimeStatement(string period_arg);
   void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream& output) const override;
 };

src/Configuration.cc

 /*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -666,15 +666,14 @@ Configuration::transformPass()
       }
 #endif
 
-  auto cluster_it
-      = cluster_name.empty() ? clusters.find(firstClusterName) : clusters.find(cluster_name);
+  auto& cluster = cluster_name.empty() ? clusters.at(firstClusterName) : clusters.at(cluster_name);
 
   double weight_denominator {0.0};
-  for (const auto& it : cluster_it->second.member_nodes)
-    weight_denominator += it.second;
+  for (const auto& [name, weight] : cluster.member_nodes)
+    weight_denominator += weight;
 
-  for (auto& member_node : cluster_it->second.member_nodes)
-    member_node.second /= weight_denominator;
+  for (auto& [name, weight] : cluster.member_nodes)
+    weight /= weight_denominator;
 }
 
 vector<filesystem::path>
@@ -701,51 +700,35 @@ Configuration::writeCluster(ostream& output) const
   if (!parallel && !parallel_test)
     return;
 
-  auto cluster_it
-      = cluster_name.empty() ? clusters.find(firstClusterName) : clusters.find(cluster_name);
+  const auto& cluster
+      = cluster_name.empty() ? clusters.at(firstClusterName) : clusters.at(cluster_name);
 
-  for (int i {1}; const auto& follower_node : follower_nodes)
+  for (int i {1}; const auto& [name, node] : follower_nodes)
     {
-      bool follower_node_in_member_nodes = false;
-      for (const auto& itmn : cluster_it->second.member_nodes)
-        if (follower_node.first == itmn.first)
-          follower_node_in_member_nodes = true;
-
-      if (!follower_node_in_member_nodes)
-        continue;
+      if (!cluster.member_nodes.contains(name))
+        continue; // Skip nodes not in the selected cluster
 
       output << "options_.parallel";
       if (i > 1)
         output << "(" << i << ")";
       i++;
-      output << " = struct('Local', ";
-      if (follower_node.second.computerName == "localhost")
-        output << "1, ";
-      else
-        output << "0, ";
-
-      output << "'ComputerName', '" << follower_node.second.computerName << "', "
-             << "'Port', '" << follower_node.second.port << "', "
-             << "'CPUnbr', [" << follower_node.second.minCpuNbr << ":"
-             << follower_node.second.maxCpuNbr << "], "
-             << "'UserName', '" << follower_node.second.userName << "', "
-             << "'Password', '" << follower_node.second.password << "', "
-             << "'RemoteDrive', '" << follower_node.second.remoteDrive << "', "
-             << "'RemoteDirectory', '" << follower_node.second.remoteDirectory
+      output << " = struct('Local', " << noboolalpha << (node.computerName == "localhost") << ", "
+             << "'ComputerName', '" << node.computerName << "', "
+             << "'Port', '" << node.port << "', "
+             << "'CPUnbr', [" << node.minCpuNbr << ":" << node.maxCpuNbr << "], "
+             << "'UserName', '" << node.userName << "', "
+             << "'Password', '" << node.password << "', "
+             << "'RemoteDrive', '" << node.remoteDrive << "', "
+             << "'RemoteDirectory', '" << node.remoteDirectory
              << "', "
              // The following should be switched back to “ProgramPath” once we move to Dragonfly
-             << "'DynarePath', '" << follower_node.second.programPath << "', "
-             << "'ProgramConfig', '" << follower_node.second.programConfig << "', "
-             << "'MatlabOctavePath', '" << follower_node.second.matlabOctavePath << "', "
-             << "'OperatingSystem', '" << follower_node.second.operatingSystem << "', "
-             << "'NodeWeight', '" << cluster_it->second.member_nodes.at(follower_node.first)
-             << "', "
-             << "'NumberOfThreadsPerJob', " << follower_node.second.numberOfThreadsPerJob << ", ";
-
-      if (follower_node.second.singleCompThread)
-        output << "'SingleCompThread', 'true');" << endl;
-      else
-        output << "'SingleCompThread', 'false');" << endl;
+             << "'DynarePath', '" << node.programPath << "', "
+             << "'ProgramConfig', '" << node.programConfig << "', "
+             << "'MatlabOctavePath', '" << node.matlabOctavePath << "', "
+             << "'OperatingSystem', '" << node.operatingSystem << "', "
+             << "'NodeWeight', '" << cluster.member_nodes.at(name) << "', "
+             << "'NumberOfThreadsPerJob', " << node.numberOfThreadsPerJob << ", "
+             << "'SingleCompThread', '" << boolalpha << node.singleCompThread << "');" << endl;
     }
 
   // Default values for the following two are both in DynareMain.cc and

src/Configuration.hh

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

src/DataTree.cc

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -24,6 +24,7 @@
 #include <fstream>
 #include <iostream>
 #include <iterator>
+#include <ranges>
 
 #include "DataTree.hh"
 
@@ -47,10 +48,12 @@ DataTree::initConstants()
 }
 
 DataTree::DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-                   ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
+                   ExternalFunctionsTable& external_functions_table_arg,
+                   HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
     symbol_table {symbol_table_arg},
     num_constants {num_constants_arg},
     external_functions_table {external_functions_table_arg},
+    heterogeneity_table {heterogeneity_table_arg},
     is_dynamic {is_dynamic_arg}
 {
   initConstants();
@@ -60,19 +63,21 @@ DataTree::DataTree(const DataTree& d) :
     symbol_table {d.symbol_table},
     num_constants {d.num_constants},
     external_functions_table {d.external_functions_table},
+    heterogeneity_table {d.heterogeneity_table},
     is_dynamic {d.is_dynamic},
     local_variables_vector {d.local_variables_vector}
 {
   // 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&
@@ -81,6 +86,7 @@ DataTree::operator=(const DataTree& d)
   assert(&symbol_table == &d.symbol_table);
   assert(&num_constants == &d.num_constants);
   assert(&external_functions_table == &d.external_functions_table);
+  assert(&heterogeneity_table == &d.heterogeneity_table);
   assert(is_dynamic == d.is_dynamic);
 
   num_const_node_map.clear();
@@ -791,12 +797,18 @@ DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<exp
   return p;
 }
 
+expr_t
+DataTree::AddSum(expr_t arg)
+{
+  return AddUnaryOp(UnaryOpcode::sum, arg);
+}
+
 bool
 DataTree::isSymbolUsed(int symb_id) const
 {
-  for (const auto& [symb_lag, expr] : variable_node_map)
-    if (symb_lag.first == symb_id)
-      return true;
+  if (ranges::any_of(views::keys(variable_node_map),
+                     [=](const auto& symb_lag) { return symb_lag.first == symb_id; }))
+    return true;
 
   if (local_variables_table.contains(symb_id))
     return true;
@@ -842,18 +854,18 @@ DataTree::addAllParamDerivId([[maybe_unused]] set<int>& deriv_id_set)
 bool
 DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
 {
-  return any_of(unary_op_node_map.begin(), unary_op_node_map.end(),
-                [=](const auto& it) { return get<1>(it.first) == opcode; });
+  return ranges::any_of(views::keys(unary_op_node_map),
+                        [=](const auto& key) { return get<1>(key) == opcode; });
 }
 
 bool
 DataTree::isUnaryOpUsedOnType(SymbolType type, UnaryOpcode opcode) const
 {
   set<int> var;
-  for (const auto& it : unary_op_node_map)
-    if (get<1>(it.first) == opcode)
+  for (const auto& [key, value] : unary_op_node_map)
+    if (get<1>(key) == opcode)
       {
-        it.second->collectVariables(type, var);
+        value->collectVariables(type, var);
         if (!var.empty())
           return true;
       }
@@ -863,18 +875,18 @@ DataTree::isUnaryOpUsedOnType(SymbolType type, UnaryOpcode opcode) const
 bool
 DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
 {
-  return any_of(binary_op_node_map.begin(), binary_op_node_map.end(),
-                [=](const auto& it) { return get<2>(it.first) == opcode; });
+  return ranges::any_of(views::keys(binary_op_node_map),
+                        [=](const auto& key) { return get<2>(key) == opcode; });
 }
 
 bool
 DataTree::isBinaryOpUsedOnType(SymbolType type, BinaryOpcode opcode) const
 {
   set<int> var;
-  for (const auto& it : binary_op_node_map)
-    if (get<2>(it.first) == opcode)
+  for (const auto& [key, value] : binary_op_node_map)
+    if (get<2>(key) == opcode)
       {
-        it.second->collectVariables(type, var);
+        value->collectVariables(type, var);
         if (!var.empty())
           return true;
       }
@@ -900,7 +912,7 @@ DataTree::writeCHelpersDefinition(ostream& output) const
            << "getPowerDeriv(double x, double p, int k)" << endl
            << "{" << endl
            << "  if (fabs(x) < " << power_deriv_near_zero
-           << " && p > 0 && k > p && fabs(p-nearbyint(p)) < " << power_deriv_near_zero << ')'
+           << " && p >= 0 && k > p && fabs(p-nearbyint(p)) < " << power_deriv_near_zero << ')'
            << endl
            << "    return 0.0;" << endl
            << "  else" << endl
@@ -959,8 +971,8 @@ DataTree::writeToFileIfModified(stringstream& new_contents, const filesystem::pa
 {
   ifstream old_file {filename, ios::in | ios::binary};
   if (old_file.is_open()
-      && equal(istreambuf_iterator<char> {old_file}, istreambuf_iterator<char> {},
-               istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {}))
+      && ranges::equal(istreambuf_iterator<char> {old_file}, istreambuf_iterator<char> {},
+                       istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {}))
     return;
   old_file.close();
 
@@ -972,7 +984,7 @@ DataTree::writeToFileIfModified(stringstream& new_contents, const filesystem::pa
       cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  copy(istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {},
-       ostreambuf_iterator<char> {new_file});
+  ranges::copy(istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {},
+               ostreambuf_iterator<char> {new_file});
   new_file.close();
 }

src/DataTree.hh

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -33,6 +33,7 @@
 
 #include "ExprNode.hh"
 #include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
 #include "NumericalConstants.hh"
 #include "SubModel.hh"
 #include "SymbolTable.hh"
@@ -48,7 +49,11 @@ public:
   NumericalConstants& num_constants;
   //! A reference to the external functions table
   ExternalFunctionsTable& external_functions_table;
+  // A reference to the heterogeneity table
+  HeterogeneityTable& heterogeneity_table;
   //! Is it possible to use leads/lags on variable nodes?
+  /* NB: This data member cannot be replaced by a virtual method, because this information is needed
+     in AddVariable(), which itself can be called from the copy constructor. */
   const bool is_dynamic;
 
 private:
@@ -137,7 +142,8 @@ private:
 
 public:
   DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-           ExternalFunctionsTable& external_functions_table_arg, bool is_static_args = false);
+           ExternalFunctionsTable& external_functions_table_arg,
+           HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);
 
   virtual ~DataTree() = default;
 
@@ -273,6 +279,9 @@ public:
   //! Adds an external function node for the second derivative of an external function
   expr_t AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t>& arguments,
                                         int input_index1, int input_index2);
+  // Adds "SUM(arg)" to model tree
+  expr_t AddSum(expr_t arg);
+
   //! Checks if a given symbol is used somewhere in the data tree
   [[nodiscard]] bool isSymbolUsed(int symb_id) const;
   //! Checks if a given unary op is used somewhere in the data tree
@@ -345,14 +354,6 @@ public:
   //! Adds to the set all the deriv IDs corresponding to parameters
   virtual void addAllParamDerivId(set<int>& deriv_id_set);
 
-  //! Returns bool indicating whether DataTree represents a Dynamic Model (returns true in
-  //! DynamicModel.hh)
-  [[nodiscard]] virtual bool
-  isDynamic() const
-  {
-    return false;
-  };
-
   struct UnknownLocalVariableException
   {
     //! Symbol ID
@@ -360,13 +361,19 @@ public:
   };
 
   [[nodiscard]] expr_t
-  getLocalVariable(int symb_id) const
+  getLocalVariable(int symb_id, int lead_lag) const
   {
     auto it = local_variables_table.find(symb_id);
     if (it == local_variables_table.end())
       throw UnknownLocalVariableException {symb_id};
 
-    return it->second;
+    /* 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
@@ -395,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.hh

 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -43,6 +43,40 @@ public:
      “model” block. The default should not be too small (see dynare#1389). */
   double balanced_growth_test_tol {1e-6};
 
+  /* For a given equation, tracks all the regimes and the declared alternatives with combinations of
+     bind and relax tags */
+  class OccbinRegimeTracker
+  {
+  private:
+    // The list of constraints used for this equation
+    vector<string> constraints;
+    /* The list of regimes present for this equation; each regime is a vector of boolean, of same
+       length as “constraints”; each boolean represents a constraint (in the order of
+       “constraints”): false for relax, true for bind */
+    set<vector<bool>> regimes_present;
+
+  public:
+    struct ConstraintInBothBindAndRelaxException
+    {
+      const string constraint;
+    };
+    struct RegimeAlreadyPresentException
+    {
+      const vector<string> constraints_bind, constraints_relax;
+    };
+    void addRegime(const vector<string>& constraints_bind,
+                   const vector<string>& constraints_relax) noexcept(false);
+    struct MissingRegimeException
+    {
+      const vector<string> constraints_bind, constraints_relax;
+    };
+    void checkAllRegimesPresent() const noexcept(false);
+
+  private:
+    [[nodiscard]] pair<vector<string>, vector<string>>
+    convertBitVectorToRegimes(const vector<bool>& r) const;
+  };
+
 private:
   /* Used in the balanced growth test, for skipping equations where the test
      cannot be performed (i.e. when LHS=RHS at the initial values). Should not
@@ -59,6 +93,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;
@@ -246,11 +283,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
@@ -276,6 +313,9 @@ private:
     return blocks_jacob_cols_endo[blk].at({var, lag});
   }
 
+  // Used to check consistency of bind/relax tags; the keys are equation names
+  map<string, OccbinRegimeTracker> occbin_regime_trackers;
+
 protected:
   string
   modelClassName() const override
@@ -286,6 +326,7 @@ protected:
 public:
   DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
                ExternalFunctionsTable& external_functions_table_arg,
+               HeterogeneityTable& heterogeneity_table_arg,
                TrendComponentModelTable& trend_component_model_table_arg,
                VarModelTable& var_model_table_arg);
 
@@ -419,7 +460,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();
@@ -428,7 +470,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;
@@ -442,7 +486,8 @@ public:
      It also assumes that the “bind” and “relax” tags have been cleared from
      eq_tags. */
   void addOccbinEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags,
-                         const vector<string>& regimes_bind, const vector<string>& regimes_relax);
+                         const vector<string>& constraints_bind,
+                         const vector<string>& constraints_relax);
 
   //! Writes LaTeX file with the equations of the dynamic model
   void writeLatexFile(const string& basename, bool write_equation_tags) const;
@@ -496,13 +541,6 @@ public:
 
   void addAllParamDerivId(set<int>& deriv_id_set) override;
 
-  //! Returns true indicating that this is a dynamic model
-  bool
-  isDynamic() const override
-  {
-    return true;
-  };
-
   //! Drive test of detrended equations
   void runTrendTest(const eval_context_t& eval_context);
 
@@ -530,6 +568,10 @@ public:
   // Performs the transformations associated to variables declared with “var(log)”
   void substituteLogTransform();
 
+  /* Performs the transformations associated to aggregation operators in heterogeneous models, such
+     as SUM(…) */
+  void substituteAggregationOperators();
+
   // Check that no variable was declared with “var(log)” in the given equations
   void checkNoWithLogTransform(const set<int>& eqnumbers);
 
@@ -601,6 +643,18 @@ public:
       map<string, int>& pac_aux_var_symb_ids, map<string, vector<int>>& pac_aux_param_symb_ids,
       map<string, expr_t>& pac_expectation_substitution);
 
+  /* For a PAC MCE model with an associated pac_target_info, fill pac_expectation_substitution with
+     the expression that will be substituted for the pac_expectation operator. In the process, add
+     the variables and the equations defining Z₁ and Z₀ for each component. The new auxiliary
+     parameters are added to pac_mce_alpha_symb_ids. The routine also creates the auxiliary
+     variables for the components, and adds the corresponding equations.
+  */
+  void computePacModelConsistentExpectationSubstitutionWithComponents(
+      const string& name, int discount_symb_id, int pac_eq_max_lag,
+      ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+      vector<PacModelTable::target_component_t>& pac_target_components,
+      map<string, expr_t>& pac_expectation_substitution);
+
   /* For a PAC backward model, fill pac_expectation_substitution with the
      expression that will be substituted for the pac_expectation operator.
      The symbol IDs of the new parameters are also added to pac_aux_param_symb_ids.
@@ -659,8 +713,8 @@ 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
   bool ParamUsedWithLeadLag() const;
@@ -695,6 +749,9 @@ public:
   {
     static_mfs = static_mfs_arg;
   }
+
+  // Checks that all alternatives are declared for all Occbin regimes in all equations
+  void checkOccbinRegimes() const;
 };
 
 template<bool julia>