Compare revisions

Normann Rion · Sébastien Villemot · Johannes Pfeifer · Sébastien Villemot · Sébastien Villemot · Sébastien Villemot
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
 variables:
  TERM: linux
-  MINGW32_BOOST_VERSION: 1.84.0-1
-  MINGW64_BOOST_VERSION: 1.84.0-1
+  MINGW64_BOOST_VERSION: 1.85.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'
@@ -48,6 +47,9 @@ build_macos_x86_64:
  tags:
    - macOS
  script:
+    # Workaround for bug in Xcode 15.3 which does not include m4
+    # See https://github.com/Homebrew/homebrew-core/issues/165388 and https://trac.macports.org/ticket/69639
+    - export PATH="/usr/local/opt/m4/bin/:$PATH"
    - arch -x86_64 meson setup -D buildtype=release --native-file scripts/homebrew-native-x86_64.ini build
    - arch -x86_64 meson compile -C build -v
  artifacts:
@@ -60,6 +62,9 @@ build_macos_arm64:
    - macOS
  script:
    - export PATH="/opt/homebrew/bin:$PATH"
+    # Workaround for bug in Xcode 15.3 which does not include m4
+    # See https://github.com/Homebrew/homebrew-core/issues/165388 and https://trac.macports.org/ticket/69639
+    - export PATH="/opt/homebrew/opt/m4/bin/:$PATH"
    - arch -arm64 meson setup -D buildtype=release --native-file scripts/homebrew-native-arm64.ini build
    - arch -arm64 meson compile -C build -v
  artifacts:

--- a/doc/macroprocessor/macroprocessor.tex
+++ b/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}

--- a/scripts/homebrew-native-arm64.ini
+++ b/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'


--- a/scripts/homebrew-native-x86_64.ini
+++ b/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'


--- a/scripts/wasm.ini
+++ b/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'
--- a/src/Bytecode.cc
+++ b/src/Bytecode.cc
@@ -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;
 }

--- a/src/CommonEnums.hh
+++ b/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

--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -920,11 +920,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 ("
@@ -1183,9 +1181,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 ";
@@ -1766,21 +1763,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
@@ -1788,20 +1780,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 << "}"
         << "}";
 }
@@ -1817,21 +1806,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
@@ -1839,20 +1823,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 << "}"
         << "}";
 }
@@ -4231,18 +4212,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
@@ -4601,10 +4577,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
@@ -4812,10 +4785,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

--- a/src/Configuration.cc
+++ b/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

--- a/src/DataTree.cc
+++ b/src/DataTree.cc
 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -47,10 +47,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,6 +62,7 @@ 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}
 {
@@ -81,6 +84,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,6 +795,12 @@ 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
 {
@@ -842,8 +852,8 @@ 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(unary_op_node_map,
+                        [=](const auto& it) { return get<1>(it.first) == opcode; });
 }

 bool
@@ -863,8 +873,8 @@ 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(binary_op_node_map,
+                        [=](const auto& it) { return get<2>(it.first) == opcode; });
 }

 bool
@@ -959,8 +969,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 +982,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();
 }
--- a/src/DataTree.hh
+++ b/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,13 @@ 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;
+    return it->second->decreaseLeadsLags(-lead_lag);
  }

  static void

--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -50,9 +50,11 @@ DynamicModel::copyHelper(const DynamicModel& m)

 DynamicModel::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) :
-    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
+    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+               heterogeneity_table_arg, true},
    trend_component_model_table {trend_component_model_table_arg},
    var_model_table {var_model_table_arg}
 {
@@ -180,31 +182,29 @@ DynamicModel::writeDynamicBytecode(const string& basename) const
                                                 + temporary_terms_derivatives[1].size())};

  // Declare the (single) block
-  vector<int> exo(symbol_table.exo_nbr()), exo_det(symbol_table.exo_det_nbr());
-  iota(exo.begin(), exo.end(), 0);
-  iota(exo_det.begin(), exo_det.end(), 0);
-
-  int jacobian_ncols_endo {static_cast<int>(count_if(
-      dyn_jacobian_cols_table.begin(), dyn_jacobian_cols_table.end(),
-      [this](const auto& v) { return getTypeByDerivID(v.first) == SymbolType::endogenous; }))};
-  vector<int> eq_idx(equations.size());
-  iota(eq_idx.begin(), eq_idx.end(), 0);
-  vector<int> endo_idx(symbol_table.endo_nbr());
-  iota(endo_idx.begin(), endo_idx.end(), 0);
+  auto exo = views::iota(0, symbol_table.exo_nbr());
+  auto exo_det = views::iota(0, symbol_table.exo_det_nbr());
+  auto eq_idx = views::iota(0, static_cast<int>(equations.size()));
+  auto endo_idx = views::iota(0, symbol_table.endo_nbr());
+
+  int jacobian_ncols_endo {
+      static_cast<int>(ranges::count_if(dyn_jacobian_cols_table, [this](const auto& v) {
+        return getTypeByDerivID(v.first) == SymbolType::endogenous;
+      }))};

  code_file << Bytecode::FBEGINBLOCK {symbol_table.endo_nbr(),
                                      simulation_type,
                                      0,
                                      symbol_table.endo_nbr(),
-                                      endo_idx,
-                                      eq_idx,
+                                      {endo_idx.begin(), endo_idx.end()},
+                                      {eq_idx.begin(), eq_idx.end()},
                                      false,
                                      u_count_int,
                                      jacobian_ncols_endo,
                                      symbol_table.exo_det_nbr(),
                                      symbol_table.exo_nbr(),
-                                      exo_det,
-                                      exo};
+                                      {exo_det.begin(), exo_det.end()},
+                                      {exo.begin(), exo.end()}};

  writeBytecodeHelper<true>(code_file);
 }
@@ -1060,8 +1060,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
      try
        {
          getDerivID(symbol_table.getID(SymbolType::endogenous, endo_idx_block2orig[endoID]), lag);
-          if (find(state_var.begin(), state_var.end(), endo_idx_block2orig[endoID])
-              == state_var.end())
+          if (ranges::find(state_var, endo_idx_block2orig[endoID]) == state_var.end())
            state_var.push_back(endo_idx_block2orig[endoID]);
        }
      catch (UnknownDerivIDException& e)
@@ -1136,7 +1135,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
    output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
  output << "];" << endl;

-  writeDriverSparseIndicesHelper<true, false>(output);
+  writeDriverSparseIndicesHelper("dynamic", output);

  // Write LHS of each equation in text form
  output << "M_.lhs = {" << endl;
@@ -1221,8 +1220,7 @@ DynamicModel::updateVarAndTrendModel() const
                        catch (...)
                          {
                          }
-                      if (find(trend_lhs.begin(), trend_lhs.end(), *trend_var_symb_id)
-                          == trend_lhs.end())
+                      if (ranges::find(trend_lhs, *trend_var_symb_id) == trend_lhs.end())
                        {
                          cerr << "ERROR: trend found in trend_component equation #" << eqn << " ("
                               << symbol_table.getName(*trend_var_symb_id)
@@ -1650,13 +1648,13 @@ DynamicModel::computeErrorComponentMatrices(const ExprNode::subst_table_t& diff_

      for (int i {0}; auto eqn : eqns)
        {
-          if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
+          if (ranges::find(nontarget_eqnums, eqn) != nontarget_eqnums.end())
            parsed_undiff_nontarget_lhs.push_back(undiff_nontarget_lhs.at(i));
          i++;
        }

      for (int i {0}; auto eqn : eqns)
-        if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
+        if (ranges::find(nontarget_eqnums, eqn) != nontarget_eqnums.end())
          equations[eqn]->arg2->fillErrorCorrectionRow(i++, parsed_undiff_nontarget_lhs, target_lhs,
                                                       A0, A0star);
      A0r[model_name] = A0;
@@ -1774,7 +1772,7 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,

  for (auto eqn : nontrend_eqnums)
    {
-      auto i = distance(eqnumber.begin(), find(eqnumber.begin(), eqnumber.end(), eqn));
+      auto i = distance(eqnumber.begin(), ranges::find(eqnumber, eqn));

      if (eqnumber[i] != eqn)
        {
@@ -2830,11 +2828,12 @@ DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const strin

  writeSetAuxiliaryVariablesFile<true>(basename, julia);

-  writeComplementarityConditionsFile<true>(basename);
-
-  // Support for model debugging
  if (!julia)
-    writeDebugModelMFiles<true>(basename);
+    {
+      writeComplementarityConditionsFile<true>(basename);
+      // Support for model debugging
+      writeDebugModelMFiles<true>(basename);
+    }
 }

 void
@@ -2968,6 +2967,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& planner_objective,
            {
              orig_endo_nbr++;
              neweqs_lineno.emplace_back(nullopt);
+              if (string eqname {"Ramsey FOC w.r.t. "s + symbol_table.getName(symb_id)};
+                  !equation_tags.exists("name", eqname))
+                neweqs_tags.emplace(neweqs.size() - 1, map<string, string> {{"name", eqname}});
+
              if (cloned_ramsey_constraints.contains(symb_id))
                {
                  auto& [lower_bound, upper_bound] = cloned_ramsey_constraints.at(symb_id);
@@ -3013,8 +3016,8 @@ DynamicModel::expandEqTags()
    if (!existing_tags.contains(eq))
      {
        if (auto lhs_expr = dynamic_cast<VariableNode*>(equations[eq]->arg1);
-            lhs_expr && !equation_tags.exists("name", symbol_table.getName(lhs_expr->symb_id)))
-          equation_tags.add(eq, "name", symbol_table.getName(lhs_expr->symb_id));
+            lhs_expr && !equation_tags.exists("name", lhs_expr->getName()))
+          equation_tags.add(eq, "name", lhs_expr->getName());
        else if (!equation_tags.exists("name", to_string(eq + 1)))
          equation_tags.add(eq, "name", to_string(eq + 1));
        else
@@ -3036,8 +3039,7 @@ DynamicModel::findUnusedEndogenous()
  for (auto& equation : static_only_equations)
    equation->collectVariables(SymbolType::endogenous, usedEndo);
  set<int> allEndo = symbol_table.getEndogenous();
-  set_difference(allEndo.begin(), allEndo.end(), usedEndo.begin(), usedEndo.end(),
-                 inserter(unusedEndo, unusedEndo.begin()));
+  ranges::set_difference(allEndo, usedEndo, inserter(unusedEndo, unusedEndo.begin()));
  return unusedEndo;
 }

@@ -3051,10 +3053,8 @@ DynamicModel::findUnusedExogenous()
    equation->collectVariables(SymbolType::exogenous, usedExo);
  set<int> observedExo = symbol_table.getObservedExogenous();
  set<int> allExo = symbol_table.getExogenous();
-  set_difference(allExo.begin(), allExo.end(), observedExo.begin(), observedExo.end(),
-                 inserter(unobservedExo, unobservedExo.begin()));
-  set_difference(unobservedExo.begin(), unobservedExo.end(), usedExo.begin(), usedExo.end(),
-                 inserter(unusedExo, unusedExo.begin()));
+  ranges::set_difference(allExo, observedExo, inserter(unobservedExo, unobservedExo.begin()));
+  ranges::set_difference(unobservedExo, usedExo, inserter(unusedExo, unusedExo.begin()));
  return unusedExo;
 }

@@ -3494,10 +3494,8 @@ pair<lag_equivalence_table_t, ExprNode::subst_table_t>
 DynamicModel::substituteUnaryOps(VarExpectationModelTable& var_expectation_model_table,
                                 PacModelTable& pac_model_table)
 {
-  vector<int> eqnumbers(equations.size());
-  iota(eqnumbers.begin(), eqnumbers.end(), 0);
-  return substituteUnaryOps({eqnumbers.begin(), eqnumbers.end()}, var_expectation_model_table,
-                            pac_model_table);
+  auto v = views::iota(0, static_cast<int>(equations.size()));
+  return substituteUnaryOps({v.begin(), v.end()}, var_expectation_model_table, pac_model_table);
 }

 pair<lag_equivalence_table_t, ExprNode::subst_table_t>
@@ -3701,6 +3699,36 @@ DynamicModel::substituteLogTransform()
    }
 }

+void
+DynamicModel::substituteAggregationOperators()
+{
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode*> neweqs;
+
+  for (auto& [symb_id, expr] : local_variables_table)
+    expr = expr->substituteAggregationOperators(subst_table, neweqs);
+
+  for (auto& equation : equations)
+    {
+      equation = dynamic_cast<BinaryOpNode*>(
+          equation->substituteAggregationOperators(subst_table, neweqs));
+      assert(equation);
+    }
+
+  for (auto& equation : static_only_equations)
+    {
+      equation = dynamic_cast<BinaryOpNode*>(
+          equation->substituteAggregationOperators(subst_table, neweqs));
+      assert(equation);
+    }
+
+  for (auto neweq : neweqs)
+    {
+      addEquation(neweq, nullopt);
+      addAuxEquation(neweq);
+    }
+}
+
 void
 DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
 {
@@ -3711,8 +3739,7 @@ DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
  const set<int>& with_log_transform = symbol_table.getVariablesWithLogTransform();

  vector<int> intersect;
-  set_intersection(endos.begin(), endos.end(), with_log_transform.begin(), with_log_transform.end(),
-                   back_inserter(intersect));
+  ranges::set_intersection(endos, with_log_transform, back_inserter(intersect));
  if (!intersect.empty())
    {
      cerr << "ERROR: the following variables are declared with var(log) and therefore cannot "
@@ -4394,8 +4421,9 @@ DynamicModel::OccbinRegimeTracker::checkAllRegimesPresent() const
      if (it == r.end())
        break;
      *it = true;
-      if (it != r.begin())
-        fill(r.begin(), prev(it), false);
+      /* NB: cannot use ranges::fill since vector<bool>::iterator does not satisfy
+         indirectly_writable concept. Should be fixed in C++23: https://wg21.link/p2321 */
+      fill(r.begin(), it, false);
    }
  while (true);
 }

--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -325,6 +325,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);

@@ -539,13 +540,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);

@@ -573,6 +567,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);


--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -87,7 +87,7 @@ string
 str_tolower(string s)
 {
  // Converting to unsigned char is needed, see https://en.cppreference.com/w/cpp/string/byte/tolower
-  transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return std::tolower(c); });
+  ranges::transform(s, s.begin(), [](unsigned char c){ return std::tolower(c); });
  return s;
 }
 }
@@ -100,7 +100,7 @@ str_tolower(string s)
 %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED PROPOSAL_DISTRIBUTION REALTIME VINTAGE
 %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
 %token COMMA CONSIDER_ALL_ENDOGENOUS CONSIDER_ALL_ENDOGENOUS_AND_AUXILIARY CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION
-%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
+%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_CYCLE_REDUCTION_MAXITER DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR EXPRESSION
 %token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS FIRST_SIMULATION_PERIOD LAST_SIMULATION_PERIOD LAST_OBS
 %token SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED
@@ -216,12 +216,13 @@ str_tolower(string s)
 %token HOMOTOPY_MAX_COMPLETION_SHARE HOMOTOPY_MIN_STEP_SIZE HOMOTOPY_INITIAL_STEP_SIZE HOMOTOPY_STEP_SIZE_INCREASE_SUCCESS_COUNT
 %token HOMOTOPY_LINEARIZATION_FALLBACK HOMOTOPY_MARGINAL_LINEARIZATION_FALLBACK HOMOTOPY_EXCLUDE_VAREXO FROM_INITVAL_TO_ENDVAL
 %token STATIC_MFS RELATIVE_TO_INITVAL MATCHED_IRFS MATCHED_IRFS_WEIGHTS WEIGHTS PERPENDICULAR
+%token HETEROGENEITY HETEROGENEITY_DIMENSION SUM

 %token <vector<string>> SYMBOL_VEC

 %type <expr_t> expression expression_or_empty
 %type <expr_t> equation hand_side
-%type <string> non_negative_number signed_number signed_integer date_str
+%type <string> non_negative_number signed_number signed_integer
 %type <string> filename symbol namespace_qualified_filename namespace_qualified_symbol
 %type <string> date_expr signed_inf signed_number_w_inf range
 %type <string> integer_range signed_integer_range boolean
@@ -272,6 +273,7 @@ statement : parameters
          | predetermined_variables
          | model_local_variable
          | change_type
+          | heterogeneity_dimension
          | model
          | initval
          | initval_file
@@ -531,9 +533,9 @@ log_trend_var : LOG_TREND_VAR '(' LOG_GROWTH_FACTOR EQUAL { driver.begin_model()
              ;

 var : VAR symbol_list_with_tex_and_partition ';'
-      { driver.var($2, false); }
+      { driver.var($2, {}, false); }
    | VAR '(' LOG ')' symbol_list_with_tex_and_partition ';'
-      { driver.var($5, true); }
+      { driver.var($5, {}, true); }
    | VAR '(' DEFLATOR EQUAL { driver.begin_model(); } hand_side ')' symbol_list_with_tex_and_partition ';'
      { driver.end_nonstationary_var(false, $6, $8, false); }
    | VAR '(' LOG COMMA DEFLATOR EQUAL { driver.begin_model(); } hand_side ')' symbol_list_with_tex_and_partition ';'
@@ -542,6 +544,8 @@ var : VAR symbol_list_with_tex_and_partition ';'
      { driver.end_nonstationary_var(true, $6, $8, false); }
    /* The case LOG + LOG_DEFLATOR is omitted, because it does not make much sense
       from an economic point of view (amounts to taking the log two times) */
+    | VAR '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+      { driver.var($7, $5, false); }
    ;

 var_remove : VAR_REMOVE symbol_list ';' { driver.var_remove($2); };
@@ -615,7 +619,9 @@ var_expectation_model_option : VARIABLE EQUAL symbol
                             ;

 varexo : VAREXO symbol_list_with_tex_and_partition ';'
-         { driver.varexo($2); }
+         { driver.varexo($2, {}); }
+       | VAREXO '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+         { driver.varexo($7, $5); }
       ;

 varexo_det : VAREXO_DET symbol_list_with_tex_and_partition ';'
@@ -627,7 +633,9 @@ predetermined_variables : PREDETERMINED_VARIABLES symbol_list ';'
                        ;

 parameters : PARAMETERS symbol_list_with_tex_and_partition ';'
-             { driver.parameters($2); }
+             { driver.parameters($2, {}); }
+           | PARAMETERS '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+             { driver.parameters($7, $5); }
           ;

 model_local_variable : MODEL_LOCAL_VARIABLE symbol_list_with_tex ';'
@@ -648,6 +656,10 @@ change_type_arg : PARAMETERS
                  { $$ = SymbolType::exogenousDet; }
                ;

+heterogeneity_dimension : HETEROGENEITY_DIMENSION symbol_list ';'
+                          { driver.heterogeneity_dimension($2); }
+                        ;
+
 init_param : symbol EQUAL expression ';' { driver.init_param($1, $3); };

 expression : '(' expression ')'
@@ -995,6 +1007,8 @@ model : MODEL ';' { driver.begin_model(); }
        equation_list END ';' { driver.end_model(); }
      | MODEL '(' model_options_list ')' ';' { driver.begin_model(); }
        equation_list END ';' { driver.end_model(); }
+      | MODEL '(' HETEROGENEITY EQUAL symbol ')' { driver.begin_heterogeneous_model($5); }';'
+        equation_list END ';' { driver.end_model(); }
      ;

 equation_list : equation_list equation
@@ -1152,6 +1166,8 @@ hand_side : '(' hand_side ')'
            { $$ = driver.add_erfc($3); }
          | STEADY_STATE '(' hand_side ')'
            { $$ = driver.add_steady_state($3); }
+          | SUM '(' hand_side ')'
+            { $$ = driver.add_sum($3); }
          ;

 comma_hand_side : hand_side
@@ -1205,6 +1221,12 @@ shocks : SHOCKS ';' shock_list END ';' { driver.end_shocks(false); }
       | SHOCKS '(' LEARNT_IN EQUAL INT_NUMBER ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, false); }
       | SHOCKS '(' LEARNT_IN EQUAL INT_NUMBER COMMA OVERWRITE ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, true); }
       | SHOCKS '(' OVERWRITE COMMA LEARNT_IN EQUAL INT_NUMBER ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($7, true); }
+       | SHOCKS '(' HETEROGENEITY EQUAL symbol ')' ';' stoch_shock_list END ';'
+         { driver.end_heterogeneous_shocks($5, false); }
+       | SHOCKS '(' HETEROGENEITY EQUAL symbol COMMA OVERWRITE ')' ';' stoch_shock_list END ';'
+         { driver.end_heterogeneous_shocks($5, true); }
+       | SHOCKS '(' OVERWRITE COMMA HETEROGENEITY EQUAL symbol ')' ';' stoch_shock_list END ';'
+         { driver.end_heterogeneous_shocks($7, true); }
       ;

 shock_list : shock_list shock_elem
@@ -1212,16 +1234,23 @@ shock_list : shock_list shock_elem
           ;

 shock_elem : det_shock_elem
-           | VAR symbol ';' STDERR expression ';'
-             { driver.add_stderr_shock($2, $5); }
-           | VAR symbol EQUAL expression ';'
-             { driver.add_var_shock($2, $4); }
-           | VAR symbol COMMA symbol EQUAL expression ';'
-             { driver.add_covar_shock($2, $4, $6); }
-           | CORR symbol COMMA symbol EQUAL expression ';'
-             { driver.add_correl_shock($2, $4, $6); }
+           | stoch_shock_elem
           ;

+stoch_shock_elem : VAR symbol ';' STDERR expression ';'
+                   { driver.add_stderr_shock($2, $5); }
+                 | VAR symbol EQUAL expression ';'
+                   { driver.add_var_shock($2, $4); }
+                 | VAR symbol COMMA symbol EQUAL expression ';'
+                   { driver.add_covar_shock($2, $4, $6); }
+                 | CORR symbol COMMA symbol EQUAL expression ';'
+                   { driver.add_correl_shock($2, $4, $6); }
+                 ;
+
+stoch_shock_list : stoch_shock_list stoch_shock_elem
+                 | stoch_shock_elem
+                 ;
+
 det_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
                 { driver.add_det_shock($2, $5, $8, ParsingDriver::DetShockType::standard); }
               | VAR symbol ';' PERIODS period_list ';' ADD value_list ';'
@@ -1465,6 +1494,7 @@ steady_options : o_solve_algo
               | o_markowitz
               | o_steady_maxit
               | o_nocheck
+               | o_noprint
               | o_steady_tolf
               | o_steady_tolx
               | o_fsolve_options
@@ -1602,6 +1632,7 @@ method_of_moments_option : o_add_tiny_number_to_cholesky
                         | o_datafile
                         | o_dirname
                         | o_dr
+                         | o_dr_cycle_reduction_maxiter
                         | o_dr_cycle_reduction_tol
                         | o_dr_logarithmic_reduction_maxiter
                         | o_dr_logarithmic_reduction_tol
@@ -1779,6 +1810,7 @@ stoch_simul_primary_options : o_solve_algo
                            | o_pruning
                            | o_dr
                            | o_dr_cycle_reduction_tol
+                            | o_dr_cycle_reduction_maxiter
                            | o_dr_logarithmic_reduction_tol
                            | o_dr_logarithmic_reduction_maxiter
                            | o_irf_plot_threshold
@@ -2061,9 +2093,7 @@ prior_pdf : BETA_PDF
            { $$ = PriorDistributions::weibull; }
          ;

-date_str : DATES
-
-date_expr : date_str
+date_expr : DATES
          | date_expr PLUS INT_NUMBER
            { $$ = $1 + '+' + $3; }
          ;
@@ -2315,6 +2345,7 @@ estimation_options : o_datafile
                   | o_lyapunov_doubling_tol
                   | o_dr
                   | o_dr_cycle_reduction_tol
+                   | o_dr_cycle_reduction_maxiter
                   | o_dr_logarithmic_reduction_tol
                   | o_dr_logarithmic_reduction_maxiter
                   | o_analytic_derivation
@@ -3949,6 +3980,7 @@ o_dr : DR EQUAL CYCLE_REDUCTION { driver.option_num("dr_cycle_reduction", "true"
       | DR EQUAL LOGARITHMIC_REDUCTION { driver.option_num("dr_logarithmic_reduction", "true"); }
       | DR EQUAL DEFAULT { driver.option_num("dr_cycle_reduction", "false"); driver.option_num("dr_logarithmic_reduction", "false"); };
 o_dr_cycle_reduction_tol : DR_CYCLE_REDUCTION_TOL EQUAL non_negative_number { driver.option_num("dr_cycle_reduction_tol", $3); };
+o_dr_cycle_reduction_maxiter : DR_CYCLE_REDUCTION_MAXITER EQUAL INT_NUMBER { driver.option_num("dr_cycle_reduction_maxiter", $3); };
 o_dr_logarithmic_reduction_tol : DR_LOGARITHMIC_REDUCTION_TOL EQUAL non_negative_number { driver.option_num("dr_logarithmic_reduction_tol", $3); };
 o_dr_logarithmic_reduction_maxiter : DR_LOGARITHMIC_REDUCTION_MAXITER EQUAL INT_NUMBER { driver.option_num("dr_logarithmic_reduction_maxiter", $3); };
 o_psd_detail_plot : DETAIL_PLOT { driver.option_num("plot_shock_decomp.detail_plot", "true"); };

--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -111,6 +111,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>predetermined_variables {BEGIN DYNARE_STATEMENT; return token::PREDETERMINED_VARIABLES;}
 <INITIAL>parameters {BEGIN DYNARE_STATEMENT; return token::PARAMETERS;}
 <INITIAL>model_local_variable {BEGIN DYNARE_STATEMENT; return token::MODEL_LOCAL_VARIABLE;}
+<INITIAL>heterogeneity_dimension {BEGIN DYNARE_STATEMENT; return token::HETEROGENEITY_DIMENSION;}
 <INITIAL>model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;}
 <INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
 <INITIAL>set_time {BEGIN DYNARE_STATEMENT; return token::SET_TIME;}
@@ -257,6 +258,9 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])

 /* Inside  of a Dynare statement */
 <DYNARE_STATEMENT>{DATE} {
+                           /* If a date is found within a statement, substitute it with a call to
+                              the dates() constructor in the input character stream. Then it will
+                              be handled by the rule that follows the present one. */
                           char* yycopy = strdup(yytext);
                           char* uput = yycopy + yyleng;
                           unput(')');
@@ -272,13 +276,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
                           unput('d');
                           free( yycopy );
                         }
-<DYNARE_STATEMENT>${DATE} { yylloc->step();
-#if (YY_FLEX_MAJOR_VERSION > 2) || (YY_FLEX_MAJOR_VERSION == 2 && YY_FLEX_MINOR_VERSION >= 6)
-                            yyout << yytext + 1;
-#else
-                            *yyout << yytext + 1;
-#endif
-                          }
 <DYNARE_STATEMENT>dates  {dates_parens_nb=0; BEGIN DATES_STATEMENT; yylval->build<string>("dates");}
 <DYNARE_STATEMENT>file                  {return token::FILE;}
 <DYNARE_STATEMENT>datafile 		{return token::DATAFILE;}
@@ -811,6 +808,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>static_mfs {return token::STATIC_MFS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>balanced_growth_test_tol {return token::BALANCED_GROWTH_TEST_TOL;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>heterogeneity {return token::HETEROGENEITY;}
 <DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;}
 <DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;}
 <DYNARE_BLOCK>normal_pdf {return token::NORMAL_PDF;}
@@ -882,6 +880,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>lyapunov_fixed_point_tol {return token::LYAPUNOV_FIXED_POINT_TOL;}
 <DYNARE_STATEMENT>lyapunov_doubling_tol {return token::LYAPUNOV_DOUBLING_TOL;}
 <DYNARE_STATEMENT>dr_cycle_reduction_tol {return token::DR_CYCLE_REDUCTION_TOL;}
+<DYNARE_STATEMENT>dr_cycle_reduction_maxiter {return token::DR_CYCLE_REDUCTION_MAXITER;}
 <DYNARE_STATEMENT>dr_logarithmic_reduction_tol {return token::DR_LOGARITHMIC_REDUCTION_TOL;}
 <DYNARE_STATEMENT>dr_logarithmic_reduction_maxiter {return token::DR_LOGARITHMIC_REDUCTION_MAXITER;}
 <DYNARE_STATEMENT>replic {return token::REPLIC;}
@@ -1010,6 +1009,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
 <DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
 <DYNARE_BLOCK>pac_target_nonstationary {return token::PAC_TARGET_NONSTATIONARY;}
+<DYNARE_BLOCK>sum {return token::SUM;}
 <DYNARE_STATEMENT>discount {return token::DISCOUNT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}

--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -572,38 +572,35 @@ public:
  {
  };

-  //! Returns the maximum lead of endogenous in this expression
+  //! Returns the maximum lead of endogenous in this expression (not incl. heterogeneous endo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxEndoLead() const = 0;

-  //! Returns the maximum lead of exogenous in this expression
+  //! Returns the maximum lead of exogenous in this expression (not incl. heterogeneous exo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxExoLead() const = 0;

-  //! Returns the maximum lag of endogenous in this expression
+  //! Returns the maximum lag of endogenous in this expression (not incl. heterogeneous endo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxEndoLag() const = 0;

-  //! Returns the maximum lag of exogenous in this expression
+  //! Returns the maximum lag of exogenous in this expression (not incl. heterogeneous exo)
  /*! Always returns a non-negative value */
  [[nodiscard]] virtual int maxExoLag() const = 0;

-  //! Returns the maximum lead of endo/exo/exodet in this expression
-  /*! A negative value means that the expression contains only lagged
-    variables. A value of numeric_limits<int>::min() means that there is
-    no variable. */
+  /* Returns the maximum lead of endo/exo/exodet in this expression (including heterogeneous
+     endo/exo). A negative value means that the expression contains only lagged variables. A value
+     of numeric_limits<int>::min() means that there is no variable. */
  [[nodiscard]] virtual int maxLead() const = 0;

-  //! Returns the maximum lag of endo/exo/exodet in this expression
-  /*! A negative value means that the expression contains only leaded
-    variables. A value of numeric_limits<int>::min() means that there is
-    no variable. */
+  /* Returns the maximum lag of endo/exo/exodet in this expression (including heterogeneous
+     endo/exo). A negative value means that the expression contains only leaded variables. A value
+     of numeric_limits<int>::min() means that there is no variable. */
  [[nodiscard]] virtual int maxLag() const = 0;

-  //! Returns the maximum lag of endo/exo/exodet, as if diffs were expanded
-  /*! This function behaves as maxLag(), except that it treats diff()
-    differently. For e.g., on diff(diff(x(-1))), maxLag() returns 1 while
-    maxLagWithDiffsExpanded() returns 3. */
+  /* Returns the maximum lag of endo/exo/exodet (including heterogeneous endo/exo), as if diffs were
+     expanded. This function behaves as maxLag(), except that it treats diff() differently. For
+     e.g., on diff(diff(x(-1))), maxLag() returns 1 while maxLagWithDiffsExpanded() returns 3. */
  [[nodiscard]] virtual int maxLagWithDiffsExpanded() const = 0;

  [[nodiscard]] virtual expr_t undiff() const = 0;
@@ -947,6 +944,13 @@ public:
     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;
+
+  /* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
+     Also declares those aggregation operators in the HeterogeneityTable, so as to
+     compute their index in the dedicated vector in argument of the dynamic/static files. */
+  [[nodiscard]] virtual expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                              vector<BinaryOpNode*>& neweqs) const
+      = 0;
 };

 //! Object used to compare two nodes (using their indexes)
@@ -1058,6 +1062,8 @@ public:
                                                    = "") const override;
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
  [[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;
 };

 //! Symbol or variable node
@@ -1090,6 +1096,10 @@ protected:

 public:
  VariableNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg, int lag_arg);
+  [[nodiscard]] SymbolType get_type() const;
+  [[nodiscard]] string getName() const;
+  [[nodiscard]] int getDerivID() const;
+  [[nodiscard]] int getTypeSpecificID() const;
  void writeOutput(ostream& output, ExprNodeOutputType output_type,
                   const temporary_terms_t& temporary_terms,
                   const temporary_terms_idxs_t& temporary_terms_idxs,
@@ -1107,7 +1117,6 @@ public:
                           const deriv_node_temp_terms_t& tef_terms) const override;
  expr_t toStatic(DataTree& static_datatree) const override;
  void computeXrefs(EquationInfo& ei) const override;
-  [[nodiscard]] SymbolType get_type() const;
  BinaryOpNode* normalizeEquationHelper(const set<expr_t>& contain_var, expr_t rhs) const override;
  [[nodiscard]] int maxEndoLead() const override;
  [[nodiscard]] int maxExoLead() const override;
@@ -1162,6 +1171,8 @@ public:
                          vector<int>& powers) const override;
  [[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
 };

 //! Unary operator node
@@ -1309,6 +1320,8 @@ public:
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
  void decomposeAdditiveTerms(vector<pair<expr_t, int>>& terms, int current_sign) const override;
  [[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;
 };

 //! Binary operator node
@@ -1459,7 +1472,7 @@ public:
  [[nodiscard]] expr_t unpackPowerDeriv() const;
  //! Returns MULT_i*(lhs-rhs) = 0, creating multiplier MULT_i
  expr_t addMultipliersToConstraints(int i);
-  //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
+  //! Returns the non-zero hand side of an equation (that must have a hand side equal to zero)
  [[nodiscard]] expr_t getNonZeroPartofEquation() const;
  [[nodiscard]] bool isInStaticForm() const override;
  void fillAutoregressiveRow(int eqn, const vector<int>& lhs,
@@ -1504,6 +1517,8 @@ public:
                          vector<int>& powers) const override;
  [[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
  [[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
 };

@@ -1648,6 +1663,8 @@ public:
  [[nodiscard]] bool containsPacTargetNonstationary(const string& pac_model_name
                                                    = "") const override;
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
  [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
 };

@@ -1720,6 +1737,7 @@ protected:
 public:
  AbstractExternalFunctionNode(DataTree& datatree_arg, int idx_arg, int symb_id_arg,
                               vector<expr_t> arguments_arg);
+  [[nodiscard]] string getName() const;
  void computeTemporaryTerms(const pair<int, int>& derivOrder,
                             map<pair<int, int>, unordered_set<expr_t>>& temp_terms_map,
                             unordered_map<expr_t, pair<int, pair<int, int>>>& reference_count,
@@ -1822,6 +1840,8 @@ public:
                                                    = "") const override;
  [[nodiscard]] bool isParamTimesEndogExpr() const override;
  [[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;
 };

 class ExternalFunctionNode : public AbstractExternalFunctionNode
@@ -2024,6 +2044,8 @@ public:
  expr_t detrend(int symb_id, bool log_trend, expr_t trend) const override;
  [[nodiscard]] expr_t removeTrendLeadLag(const map<int, expr_t>& trend_symbols_map) const override;
  [[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;

 protected:
  void prepareForDerivation() override;

--- a/src/HeterogeneityTable.cc
+++ b/src/HeterogeneityTable.cc
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <cassert>
+#include <utility>
+
+#include "HeterogeneityTable.hh"
+#include "SymbolTable.hh"
+
+void
+HeterogeneityTable::setSymbolTable(SymbolTable* symbol_table_arg)
+{
+  symbol_table = symbol_table_arg;
+}
+
+int
+HeterogeneityTable::addDimension(string name)
+{
+  if (name_to_id.contains(name))
+    throw AlreadyDeclaredDimensionException {move(name)};
+
+  int id {static_cast<int>(id_to_name.size())};
+  name_to_id.emplace(name, id);
+  id_to_name.push_back(move(name));
+  return id;
+}
+
+bool
+HeterogeneityTable::exists(const string& name) const
+{
+  return name_to_id.contains(name);
+}
+
+int
+HeterogeneityTable::getID(const string& name) const
+{
+  if (auto it = name_to_id.find(name); it != name_to_id.end())
+    return it->second;
+  else
+    throw UnknownDimensionNameException {name};
+}
+
+string
+HeterogeneityTable::getName(int id) const
+{
+  if (id < 0 || id >= static_cast<int>(id_to_name.size()))
+    throw UnknownDimensionIDException {id};
+  else
+    return id_to_name[id];
+}
+
+bool
+HeterogeneityTable::empty() const
+{
+  return name_to_id.empty();
+}
+
+vector<string>
+HeterogeneityTable::getDimensions() const
+{
+  return id_to_name;
+}
+
+int
+HeterogeneityTable::size() const
+{
+  return static_cast<int>(name_to_id.size());
+}
+
+void
+HeterogeneityTable::addSummedHeterogeneousEndogenous(int symb_id)
+{
+  assert(symbol_table->getType(symb_id) == SymbolType::heterogeneousEndogenous);
+  if (summed_het_endo_to_index.contains(symb_id))
+    throw AlreadyDeclaredSummedHeterogeneousEndogenousException {symb_id};
+
+  int index {static_cast<int>(index_to_summed_het_endo.size())};
+  summed_het_endo_to_index.emplace(symb_id, index);
+  index_to_summed_het_endo.push_back(symb_id);
+}
+
+int
+HeterogeneityTable::getSummedHeterogenousEndogenousIndex(int symb_id) const
+{
+  if (auto it = summed_het_endo_to_index.find(symb_id); it != summed_het_endo_to_index.end())
+    return it->second;
+  else
+    throw UnknownSummedHeterogeneousEndogenousException {symb_id};
+}
+
+int
+HeterogeneityTable::aggregateEndoSize() const
+{
+  return index_to_summed_het_endo.size();
+}
+
+void
+HeterogeneityTable::writeOutput(ostream& output) const
+{
+  for (size_t id {0}; id < id_to_name.size(); id++)
+    output << "M_.heterogeneity(" << id + 1 << ").dimension_name = '" << id_to_name[id] << "';"
+           << endl;
+
+  output << "M_.heterogeneity_aggregates = {" << endl;
+  for (int symb_id : index_to_summed_het_endo)
+    output << "'sum', " << symbol_table->getHeterogeneityDimension(symb_id) + 1 << ", "
+           << symbol_table->getTypeSpecificID(symb_id) + 1 << ";" << endl;
+  output << "};" << endl;
+}
+
+void
+HeterogeneityTable::writeJsonOutput(ostream& output) const
+{
+  assert(!empty());
+  output << R"("heterogeneity_dimension": [)";
+  for (bool first_written {false}; const auto& dim : id_to_name)
+    {
+      if (exchange(first_written, true))
+        output << ", ";
+      output << '"' << dim << '"';
+    }
+  output << "]" << endl;
+}
--- a/src/HeterogeneityTable.hh
+++ b/src/HeterogeneityTable.hh
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef HETEROGENEITY_TABLE_HH
+#define HETEROGENEITY_TABLE_HH
+
+#include <map>
+#include <ostream>
+#include <string>
+#include <vector>
+
+using namespace std;
+
+class SymbolTable; // Forward declaration, to avoid circularity
+
+/*
+  There is a guarantee that heterogeneity IDs are increasing, i.e. if dimension A is added after
+  dimension B, then the ID of A is greater than the ID of B.
+  Moreover, the IDs form a contiguous interval starting at 0.
+*/
+class HeterogeneityTable
+{
+private:
+  // Maps dimension names to IDs
+  map<string, int> name_to_id;
+  // Maps dimension IDs to names
+  vector<string> id_to_name;
+
+  SymbolTable* symbol_table {nullptr}; // Cannot be a reference, because of circularity
+
+  /* Keeps track of the SUM() operator instances.
+     Maps a symbol ID that appears inside a SUM() operator into an index in
+     M_.heterogeneity_aggregates */
+  map<int, int> summed_het_endo_to_index;
+
+  // Maps an index in M_.heterogeneity_aggregates into a symbol ID
+  vector<int> index_to_summed_het_endo;
+
+public:
+  void setSymbolTable(SymbolTable* symbol_table_arg);
+
+  struct AlreadyDeclaredDimensionException
+  {
+    // Dimension name
+    const string name;
+  };
+  struct UnknownDimensionNameException
+  {
+    // Dimension name
+    const string name;
+  };
+  struct UnknownDimensionIDException
+  {
+    // Dimension ID
+    const int id;
+  };
+
+  // Returns the dimension ID
+  int addDimension(string name);
+  [[nodiscard]] bool exists(const string& name) const;
+  [[nodiscard]] int getID(const string& name) const;
+  [[nodiscard]] string getName(int id) const;
+  [[nodiscard]] bool empty() const;
+  [[nodiscard]] vector<string> getDimensions() const;
+  [[nodiscard]] int size() const;
+
+  struct AlreadyDeclaredSummedHeterogeneousEndogenousException
+  {
+    const int symb_id;
+  };
+  struct UnknownSummedHeterogeneousEndogenousException
+  {
+    const int symb_id;
+  };
+
+  void addSummedHeterogeneousEndogenous(int symb_id);
+  int getSummedHeterogenousEndogenousIndex(int symb_id) const;
+  [[nodiscard]] int aggregateEndoSize() const;
+
+  void writeOutput(ostream& output) const;
+  void writeJsonOutput(ostream& output) const;
+};
+
+#endif
--- a/src/HeterogeneousModel.cc
+++ b/src/HeterogeneousModel.cc
+/*
+ * Copyright © 2024 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <cstdlib>
+#include <iostream>
+
+#include "HeterogeneousModel.hh"
+
+HeterogeneousModel::HeterogeneousModel(SymbolTable& symbol_table_arg,
+                                       NumericalConstants& num_constants_arg,
+                                       ExternalFunctionsTable& external_functions_table_arg,
+                                       HeterogeneityTable& heterogeneity_table_arg,
+                                       int heterogeneity_dimension_arg) :
+    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+               heterogeneity_table_arg, true},
+    heterogeneity_dimension {heterogeneity_dimension_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)
+{
+  ModelTree::operator=(m);
+
+  assert(heterogeneity_dimension == m.heterogeneity_dimension);
+
+  deriv_id_table = m.deriv_id_table;
+  inv_deriv_id_table = m.inv_deriv_id_table;
+
+  return *this;
+}
+
+void
+HeterogeneousModel::computeChainRuleJacobian()
+{
+  cerr << "Heterogeneous::computeChainRuleJacobian(): unimplemented" << endl;
+  exit(EXIT_FAILURE);
+}
+
+int
+HeterogeneousModel::getBlockJacobianEndoCol([[maybe_unused]] int blk, [[maybe_unused]] int var,
+                                            [[maybe_unused]] int lead_lag) const
+{
+  cerr << "Heterogeneous::getBlockJacobianEndoCol(): unimplemented" << endl;
+  exit(EXIT_FAILURE);
+}
+
+int
+HeterogeneousModel::getMFS() const
+{
+  cerr << "Heterogeneous::getMFS(): unimplemented" << endl;
+  exit(EXIT_FAILURE);
+}
+
+void
+HeterogeneousModel::computeDerivIDs()
+{
+  set<pair<int, int>> dynvars;
+
+  for (auto& equation : equations)
+    {
+      equation->collectDynamicVariables(SymbolType::heterogeneousEndogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::heterogeneousExogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
+      equation->collectDynamicVariables(SymbolType::parameter, dynvars);
+    }
+
+  for (const auto& [symb_id, lead_lag] : dynvars)
+    {
+      auto type {symbol_table.getType(symb_id)};
+      if (isHeterogeneous(type))
+        assert(symbol_table.getHeterogeneityDimension(symb_id) == heterogeneity_dimension);
+      if (type == SymbolType::heterogeneousEndogenous || type == SymbolType::endogenous)
+        assert(abs(lead_lag) <= 1);
+      if (type == SymbolType::heterogeneousExogenous || type == SymbolType::exogenous)
+        assert(lead_lag == 0);
+      int deriv_id {static_cast<int>(deriv_id_table.size())};
+      deriv_id_table.emplace(pair {symb_id, lead_lag}, deriv_id);
+      inv_deriv_id_table.emplace_back(symb_id, lead_lag);
+    }
+}
+
+void
+HeterogeneousModel::computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll)
+{
+  assert(!use_dll); // Not yet implemented
+
+  computeDerivIDs();
+
+  set<int> vars;
+  for (auto& [symb_lag, deriv_id] : deriv_id_table)
+    if (symbol_table.getType(symb_lag.first) != SymbolType::parameter)
+      vars.insert(deriv_id);
+
+  cout << "Computing " << modelClassName() << " derivatives (order " << derivsOrder << ")." << endl;
+
+  computeDerivatives(derivsOrder, vars);
+
+  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);
+    }
+}
+
+void
+HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
+{
+  assert(!julia); // Not yet implemented
+  writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
+}
+
+int
+HeterogeneousModel::getJacobianCol(int deriv_id, bool sparse) const
+{
+  assert(sparse);
+
+  SymbolType type {getTypeByDerivID(deriv_id)};
+  int tsid {getTypeSpecificIDByDerivID(deriv_id)};
+  int lag {getLagByDerivID(deriv_id)};
+
+  if (type == SymbolType::heterogeneousEndogenous)
+    return tsid + (lag + 1) * symbol_table.het_endo_nbr(heterogeneity_dimension);
+
+  int shift {3 * symbol_table.het_endo_nbr(heterogeneity_dimension)};
+
+  if (type == SymbolType::heterogeneousExogenous)
+    return shift + tsid;
+
+  shift += symbol_table.het_exo_nbr(heterogeneity_dimension);
+
+  if (type == SymbolType::endogenous)
+    return shift + tsid + (lag + 1) * symbol_table.endo_nbr();
+
+  shift += symbol_table.endo_nbr();
+
+  if (type == SymbolType::exogenous)
+    return shift + tsid;
+
+  throw UnknownDerivIDException();
+}
+
+int
+HeterogeneousModel::getJacobianColsNbr(bool sparse) const
+{
+  assert(sparse);
+  return 3 * (symbol_table.het_endo_nbr(heterogeneity_dimension) + symbol_table.endo_nbr())
+         + symbol_table.het_exo_nbr(heterogeneity_dimension) + symbol_table.exo_nbr();
+}
+
+SymbolType
+HeterogeneousModel::getTypeByDerivID(int deriv_id) const noexcept(false)
+{
+  return symbol_table.getType(getSymbIDByDerivID(deriv_id));
+}
+
+int
+HeterogeneousModel::getLagByDerivID(int deriv_id) const noexcept(false)
+{
+  if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
+    throw UnknownDerivIDException();
+
+  return inv_deriv_id_table[deriv_id].second;
+}
+
+int
+HeterogeneousModel::getSymbIDByDerivID(int deriv_id) const noexcept(false)
+{
+  if (deriv_id < 0 || deriv_id >= static_cast<int>(inv_deriv_id_table.size()))
+    throw UnknownDerivIDException();
+
+  return inv_deriv_id_table[deriv_id].first;
+}
+
+int
+HeterogeneousModel::getTypeSpecificIDByDerivID(int deriv_id) const
+{
+  return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id));
+}
+
+int
+HeterogeneousModel::getDerivID(int symb_id, int lead_lag) const noexcept(false)
+{
+  if (auto it = deriv_id_table.find({symb_id, lead_lag}); it == deriv_id_table.end())
+    throw UnknownDerivIDException();
+  else
+    return it->second;
+}
+
+void
+HeterogeneousModel::writeDriverOutput(ostream& output) const
+{
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").dynamic_tmp_nbr = [";
+  for (const auto& it : temporary_terms_derivatives)
+    output << it.size() << "; ";
+  output << "];" << endl;
+  writeDriverSparseIndicesHelper(
+      "heterogeneity("s + to_string(heterogeneity_dimension + 1) + ").dynamic", output);
+}
No results found