.clang-format

@@ -11,9 +11,12 @@ BreakInheritanceList: AfterColon
 Cpp11BracedListStyle: true
 DeriveLineEnding: false
 IndentPPDirectives: AfterHash
+InsertNewlineAtEOF: true
 PackConstructorInitializers: NextLine
 PPIndentWidth: 1
 PointerAlignment: Left
+# RemoveParentheses: ReturnStatement
+# RemoveSemicolon: true
 SpaceAfterTemplateKeyword: false
 SpaceBeforeParens: ControlStatements
 SpaceBeforeCpp11BracedList: true

.gitlab-ci.yml

 variables:
   TERM: linux
-  MINGW64_BOOST_VERSION: 1.85.0-2
+  MINGW64_BOOST_VERSION: 1.86.0-7
   WGET_OPTIONS: '--no-verbose --no-use-server-timestamps --retry-connrefused --retry-on-host-error'
   # To ensure that "false && true" fails, see https://gitlab.com/gitlab-org/gitlab-runner/-/issues/25394#note_412609647
   FF_ENABLE_BASH_EXIT_CODE_CHECK: 'true'
@@ -77,3 +77,16 @@ test_clang_format:
     - meson setup build-clang-format
     - ninja -C build-clang-format clang-format-check
   needs: []
+
+test_clang_tidy:
+  stage: test
+  script:
+    # Hack needed for meson < 1.6.0 which only looks for unversioned clang-tidy
+    - mkdir -p ~/.local/bin && ln -s /usr/bin/clang-tidy-19 ~/.local/bin/clang-tidy
+    - export PATH="$HOME/.local/bin:$PATH"
+    - meson setup build-clang-tidy
+    # Generate Flex and Bison files
+    - meson compile -C build-clang-tidy
+    - ninja -C build-clang-tidy clang-tidy
+  needs: []
+  when: manual

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/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

@@ -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
@@ -3713,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)}
 {
 }
 
@@ -3722,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) :
@@ -5347,8 +5342,8 @@ void
 ResidStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                             [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output, "options_resid_");
-  output << "display_static_residuals(M_, options_, oo_, options_resid_);" << endl;
+  options_list.writeOutput(output);
+  output << "display_static_residuals(M_, options_, oo_);" << endl;
 }
 
 void

src/ComputingTasks.hh

@@ -524,7 +524,7 @@ public:
   blockName() const override
   {
     return "estimated_params";
-  };
+  }
   void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
   void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream& output) const override;
@@ -542,7 +542,7 @@ public:
   blockName() const override
   {
     return "estimated_params_init";
-  };
+  }
   void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
   void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream& output) const override;
@@ -557,7 +557,7 @@ public:
   blockName() const override
   {
     return "estimated_params_bounds";
-  };
+  }
   void checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings) override;
   void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream& output) const override;
@@ -1007,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

@@ -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 ranges::any_of(unary_op_node_map,
-                        [=](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 ranges::any_of(binary_op_node_map,
-                        [=](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;
       }

src/DataTree.hh

@@ -33,6 +33,7 @@
 
 #include "ExprNode.hh"
 #include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
 #include "NumericalConstants.hh"
 #include "SubModel.hh"
 #include "SymbolTable.hh"
@@ -48,6 +49,8 @@ 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. */
@@ -139,7 +142,8 @@ private:
 
 public:
   DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-           ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg = false);
+           ExternalFunctionsTable& external_functions_table_arg,
+           HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);
 
   virtual ~DataTree() = default;
 
@@ -275,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
@@ -360,7 +367,13 @@ public:
     if (it == local_variables_table.end())
       throw UnknownLocalVariableException {symb_id};
 
-    return it->second->decreaseLeadsLags(-lead_lag);
+    /* In the following, the case without lead/lag is optimized. It makes a difference on models
+       with many nested model-local variables, see e.g.
+       https://forum.dynare.org/t/pre-processing-takes-very-long/26865 */
+    if (lead_lag == 0)
+      return it->second;
+    else
+      return it->second->decreaseLeadsLags(-lead_lag);
   }
 
   static void
@@ -389,7 +402,7 @@ DataTree::AddPossiblyNegativeConstant(double v)
   if (isnan(v))
     return NaN;
   if (isinf(v))
-    return (v < 0 ? MinusInfinity : Infinity);
+    return v < 0 ? MinusInfinity : Infinity;
 
   bool neg = false;
   if (v < 0)

src/DynamicModel.cc

 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -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,32 +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);
+  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;
       }))};
-  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);
 
   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);
 }
@@ -1044,14 +1043,8 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
          << (static_only_equations.size() > 0) << ";" << endl;
 
   // Say if model contains an external function call
-  bool has_external_function = false;
-  for (auto equation : equations)
-    if (equation->containsExternalFunction())
-      {
-        has_external_function = true;
-        break;
-      }
-  output << "M_.has_external_function = " << boolalpha << has_external_function << ';' << endl;
+  output << "M_.has_external_function = " << boolalpha
+         << ranges::any_of(equations, &ExprNode::containsExternalFunction) << ';' << endl;
 
   // Compute list of state variables, ordered in block-order
   vector<int> state_var;
@@ -1136,7 +1129,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;
@@ -1360,8 +1353,7 @@ DynamicModel::fillVarModelTableFromOrigModel() const
           // save lhs variables
           equations[eqn]->arg1->collectVARLHSVariable(lhs);
 
-          equations[eqn]->arg1->countDiffs() > 0 ? diff_vec.push_back(true)
-                                                 : diff_vec.push_back(false);
+          diff_vec.push_back(equations[eqn]->arg1->countDiffs() > 0);
           if (diff_vec.back())
             {
               set<pair<int, int>> diff_set;
@@ -1705,10 +1697,7 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel() const
           // save lhs variables
           equations[eqn]->arg1->collectVARLHSVariable(lhs);
 
-          if (equations[eqn]->arg1->countDiffs() > 0)
-            diff_vec.push_back(true);
-          else
-            diff_vec.push_back(false);
+          diff_vec.push_back(equations[eqn]->arg1->countDiffs() > 0);
           if (diff_vec.back())
             {
               set<pair<int, int>> diff_set;
@@ -1781,7 +1770,7 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
           exit(EXIT_FAILURE);
         }
 
-      if (diff.at(i) != true)
+      if (!diff.at(i))
         {
           cerr << "ERROR: the variable on the LHS of equation #" << eqn
                << " does not have the diff operator applied to it yet you are trying to undiff it."
@@ -1789,7 +1778,6 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
           exit(EXIT_FAILURE);
         }
 
-      bool printerr = false;
       expr_t node = nullptr;
       expr_t aux_var = lhs_expr_t.at(i);
       for (const auto& it : diff_subst_table)
@@ -1806,11 +1794,8 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
         }
 
       node = node->undiff();
-      auto it1 = diff_subst_table.find(node);
-      if (it1 == diff_subst_table.end())
-        printerr = true;
 
-      if (printerr)
+      if (auto it1 = diff_subst_table.find(node); it1 == diff_subst_table.end())
         { // we have undiffed something like diff(x), hence x is not in diff_subst_table
           lhs_expr_t.at(i) = node;
           lhs.at(i) = dynamic_cast<VariableNode*>(node)->symb_id;
@@ -1831,7 +1816,7 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
   for (auto& equation : equations)
     if (equation->containsPacExpectation(name))
       {
-        if (!pac_eq_name[name].empty())
+        if (pac_eq_name.contains(name))
           {
             cerr << "It is not possible to use 'pac_expectation(" << name
                  << ")' in several equations." << endl;
@@ -1930,6 +1915,13 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
                                    move(additive_vars_params_and_constants),
                                    move(optim_additive_vars_params_and_constants)};
       }
+
+  if (!pac_eq_name.contains(name))
+    {
+      cerr << "ERROR: the model does not contain the 'pac_expectation(" << name << ")' operator."
+           << endl;
+      exit(EXIT_FAILURE);
+    }
 }
 
 int
@@ -3040,8 +3032,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;
 }
 
@@ -3055,10 +3046,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;
 }
 
@@ -3498,10 +3487,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>
@@ -3705,6 +3692,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)
 {
@@ -3731,18 +3748,18 @@ void
 DynamicModel::detrendEquations()
 {
   // We go backwards in the list of trend_vars, to deal correctly with I(2) processes
-  for (const auto& it : std::ranges::reverse_view(nonstationary_symbols_map))
+  for (const auto& [symb_id, deflator] : std::ranges::reverse_view(nonstationary_symbols_map))
     {
       for (auto& equation : equations)
         {
           equation = dynamic_cast<BinaryOpNode*>(
-              equation->detrend(it.first, it.second.first, it.second.second));
+              equation->detrend(symb_id, deflator.first, deflator.second));
           assert(equation);
         }
       for (auto& equation : static_only_equations)
         {
           equation = dynamic_cast<BinaryOpNode*>(
-              equation->detrend(it.first, it.second.first, it.second.second));
+              equation->detrend(symb_id, deflator.first, deflator.second));
           assert(equation);
         }
     }

src/DynamicModel.hh

@@ -73,7 +73,8 @@ public:
     void checkAllRegimesPresent() const noexcept(false);
 
   private:
-    pair<vector<string>, vector<string>> convertBitVectorToRegimes(const vector<bool>& r) const;
+    [[nodiscard]] pair<vector<string>, vector<string>>
+    convertBitVectorToRegimes(const vector<bool>& r) const;
   };
 
 private:
@@ -325,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);
 
@@ -566,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);
 
@@ -708,7 +714,7 @@ public:
   {
     return tuple {static_only_equations, static_only_equations_lineno,
                   static_only_complementarity_conditions, static_only_equations_equation_tags};
-  };
+  }
 
   //! Returns true if a parameter was used in the model block with a lead or lag
   bool ParamUsedWithLeadLag() const;

src/DynareFlex.ll

 /* -*- C++ -*- */
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -20,7 +20,6 @@
 
 
 %{
-#include <cstring>
 #include "ParsingDriver.hh"
 
 using namespace std;
@@ -60,7 +59,6 @@ string eofbuff;
 %x VERBATIM_BLOCK
 %x NATIVE
 %x NATIVE_COMMENT
-%x DATES_STATEMENT
 %x LINE1
 %x LINE2
 %x LINE3
@@ -70,7 +68,9 @@ string eofbuff;
 #define YY_USER_ACTION location_increment(yylloc, yytext);
 %}
 
-DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
+NAME [a-z_][a-z0-9_]*
+FLOAT_NUMBER ((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+)
+DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|[sh][12])
 
 %%
  /* Code put at the beginning of yylex() */
@@ -92,12 +92,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
                 }
 
  /* spaces, tabs and carriage returns are ignored */
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,DATES_STATEMENT,LINE1,LINE2,LINE3>[[:space:]]+  { yylloc->step(); }
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,LINE1,LINE2,LINE3>[[:space:]]+  { yylloc->step(); }
 
  /* Comments */
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>%.*
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>"//".*
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>"/*"   {comment_caller = YY_START; BEGIN COMMENT;}
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>%.*
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>"//".*
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>"/*"   {comment_caller = YY_START; BEGIN COMMENT;}
 
 <COMMENT>"*/"        {BEGIN comment_caller;}
 <COMMENT>.
@@ -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;}
@@ -235,6 +236,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>pac_target_info {BEGIN DYNARE_BLOCK; return token::PAC_TARGET_INFO;}
 <INITIAL>matched_irfs {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS;}
 <INITIAL>matched_irfs_weights {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS_WEIGHTS;}
+<INITIAL>perfect_foresight_controlled_paths {BEGIN DYNARE_BLOCK; return token::PERFECT_FORESIGHT_CONTROLLED_PATHS;}
 
  /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
@@ -245,7 +247,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>subsamples {return token::SUBSAMPLES;}
 <DYNARE_STATEMENT>options {return token::OPTIONS;}
 <DYNARE_STATEMENT>prior {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::PRIOR;
 }
 <INITIAL>std {BEGIN DYNARE_STATEMENT; return token::STD;}
@@ -255,31 +257,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>prior_function {BEGIN DYNARE_STATEMENT; return token::PRIOR_FUNCTION;}
 <INITIAL>posterior_function {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_FUNCTION;}
 
- /* Inside  of a Dynare statement */
-<DYNARE_STATEMENT>{DATE} {
-                           char* yycopy = strdup(yytext);
-                           char* uput = yycopy + yyleng;
-                           unput(')');
-                           unput('\'');
-                           while (uput > yycopy)
-                             unput(*--uput);
-                           unput('\'');
-                           unput('(');
-                           unput('s');
-                           unput('e');
-                           unput('t');
-                           unput('a');
-                           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,DYNARE_BLOCK>{DATE} {
+  yylval->emplace<string>(yytext);
+  return token::DATE;
+}
+
+ /* Inside a Dynare statement */
 <DYNARE_STATEMENT>file                  {return token::FILE;}
 <DYNARE_STATEMENT>datafile 		{return token::DATAFILE;}
 <DYNARE_STATEMENT>dirname       {return token::DIRNAME;}
@@ -312,10 +295,10 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>posterior_nograph   		{return token::POSTERIOR_NOGRAPH;}
 <DYNARE_STATEMENT>nodisplay     {return token::NODISPLAY;}
 <DYNARE_STATEMENT>graph_format  {return token::GRAPH_FORMAT;}
-<DYNARE_STATEMENT>eps  {yylval->build<string>(yytext); return token::EPS;}
-<DYNARE_STATEMENT>pdf  {yylval->build<string>(yytext); return token::PDF;}
-<DYNARE_STATEMENT>fig  {yylval->build<string>(yytext); return token::FIG;}
-<DYNARE_STATEMENT>none  {yylval->build<string>(yytext); return token::NONE;}
+<DYNARE_STATEMENT>eps  {yylval->emplace<string>(yytext); return token::EPS;}
+<DYNARE_STATEMENT>pdf  {yylval->emplace<string>(yytext); return token::PDF;}
+<DYNARE_STATEMENT>fig  {yylval->emplace<string>(yytext); return token::FIG;}
+<DYNARE_STATEMENT>none  {yylval->emplace<string>(yytext); return token::NONE;}
 <DYNARE_STATEMENT>print   		{return token::PRINT;}
 <DYNARE_STATEMENT>noprint   		{return token::NOPRINT;}
 <DYNARE_STATEMENT>conf_sig  		{return token::CONF_SIG;}
@@ -423,6 +406,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>logarithmic_reduction {return token::LOGARITHMIC_REDUCTION;}
 <DYNARE_STATEMENT>use_univariate_filters_if_singularity_is_detected {return token::USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED;}
 <DYNARE_STATEMENT>hybrid {return token::HYBRID;}
+<DYNARE_STATEMENT>use_first_order_solution {return token::USE_FIRST_ORDER_SOLUTION;}
 <DYNARE_STATEMENT>default {return token::DEFAULT;}
 <DYNARE_STATEMENT>init2shocks {return token::INIT2SHOCKS;}
 
@@ -435,6 +419,10 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>kitagawa {return token::KITAGAWA;}
 <DYNARE_STATEMENT>smooth {return token::SMOOTH;}
 <DYNARE_STATEMENT>stratified {return token::STRATIFIED;}
+<DYNARE_STATEMENT>residual {
+  yylval->emplace<string>(yytext);
+  return token::RESIDUAL;
+}
 <DYNARE_STATEMENT>cpf_weights {return token::CPF_WEIGHTS;}
 <DYNARE_STATEMENT>amisanotristani {return token::AMISANOTRISTANI;}
 <DYNARE_STATEMENT>murrayjonesparslow {return token::MURRAYJONESPARSLOW;}
@@ -456,43 +444,40 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>fsolve_options {return token::FSOLVE_OPTIONS;}
 
 <DYNARE_STATEMENT>alpha {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::ALPHA;
 }
 <DYNARE_STATEMENT>beta {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::BETA;
 }
 <DYNARE_STATEMENT>gamma {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::GAMMA;
 }
 <DYNARE_STATEMENT>inv_gamma {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::INV_GAMMA;
 }
 <DYNARE_STATEMENT>inv_gamma1 {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::INV_GAMMA1;
 }
 <DYNARE_STATEMENT>inv_gamma2 {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::INV_GAMMA2;
 }
 <DYNARE_STATEMENT>dirichlet {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::DIRICHLET;
 }
-<DYNARE_STATEMENT>weibull {
-  yylval->build<string>(yytext);
-  return token::WEIBULL;
-}
+<DYNARE_STATEMENT>weibull {return token::WEIBULL;}
 <DYNARE_STATEMENT>normal {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::NORMAL;
 }
 <DYNARE_STATEMENT>uniform {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::UNIFORM;
 }
 <DYNARE_STATEMENT>gsig2_lmdm {return token::GSIG2_LMDM;}
@@ -503,13 +488,13 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>ncsk {return token::NCSK;}
 <DYNARE_STATEMENT>nstd {return token::NSTD;}
 <DYNARE_STATEMENT>ninv {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::NINV;
 }
 <DYNARE_STATEMENT>indxparr {return token::INDXPARR;}
 <DYNARE_STATEMENT>indxovr {return token::INDXOVR;}
 <DYNARE_STATEMENT>aband {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::ABAND;
 }
 <DYNARE_STATEMENT>write_equation_tags {return token::WRITE_EQUATION_TAGS;}
@@ -526,18 +511,18 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>indxgdls {return token::INDXGDLS;}
 <DYNARE_STATEMENT>eq_ms {return token::EQ_MS;}
 <DYNARE_STATEMENT>cms {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::CMS;
 }
 <DYNARE_STATEMENT>ncms {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::NCMS;
 }
 <DYNARE_STATEMENT>eq_cms {return token::EQ_CMS;}
 <DYNARE_STATEMENT>tlindx {return token::TLINDX;}
 <DYNARE_STATEMENT>tlnumber {return token::TLNUMBER;}
 <DYNARE_STATEMENT>cnum {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::CNUM;
 }
 <DYNARE_STATEMENT>nodecomposition {return token::NODECOMPOSITION;};
@@ -619,18 +604,9 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>substitute_libs {return token::SUBSTITUTE_LIBS;}
 <DYNARE_STATEMENT>compiler {return token::COMPILER;}
 <DYNARE_STATEMENT>instruments {return token::INSTRUMENTS;}
-<DYNARE_STATEMENT>hessian  {
-  yylval->build<string>(yytext);
-  return token::HESSIAN;
-}
-<DYNARE_STATEMENT>prior_variance  {
-  yylval->build<string>(yytext);
-  return token::PRIOR_VARIANCE;
-}
-<DYNARE_STATEMENT>identity_matrix  {
-  yylval->build<string>(yytext);
-  return token::IDENTITY_MATRIX;
-}
+<DYNARE_STATEMENT>hessian  {return token::HESSIAN;}
+<DYNARE_STATEMENT>prior_variance  {return token::PRIOR_VARIANCE;}
+<DYNARE_STATEMENT>identity_matrix  {return token::IDENTITY_MATRIX;}
 <DYNARE_STATEMENT>mcmc_jumping_covariance {return token::MCMC_JUMPING_COVARIANCE;}
 
  /* These four (var, varexo, varexo_det, parameters) are for change_type */
@@ -730,14 +706,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>lmmcp {return token::LMMCP;}
 <DYNARE_STATEMENT>additional_optimizer_steps	{return token::ADDITIONAL_OPTIMIZER_STEPS;}
 <DYNARE_STATEMENT>bartlett_kernel_lag {return token::BARTLETT_KERNEL_LAG; }
-<DYNARE_STATEMENT>optimal {
-  yylval->build<string>(yytext);
-  return token::OPTIMAL;
-}
-<DYNARE_STATEMENT>diagonal  {
-  yylval->build<string>(yytext);
-  return token::DIAGONAL;
-}
 <DYNARE_STATEMENT>gmm {return token::GMM;}
 <DYNARE_STATEMENT>smm {return token::SMM;}
 <DYNARE_STATEMENT>irf_matching {return token::IRF_MATCHING;}
@@ -783,8 +751,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>non_zero {return token::NON_ZERO;}
 
 <DYNARE_STATEMENT>\$[^$]*\$ {
-  strtok(yytext + 1, "$");
-  yylval->build<string>(yytext + 1);
+  yylval->emplace<string>(yytext + 1).pop_back();
   return token::TEX_NAME;
 }
 
@@ -797,20 +764,21 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>periods {return token::PERIODS;}
 <DYNARE_BLOCK>scales {return token::SCALES;}
 <DYNARE_BLOCK>add {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::ADD;
 }
 <DYNARE_BLOCK>multiply {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::MULTIPLY;
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>cutoff {return token::CUTOFF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>mfs {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::MFS;
 }
 <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;}
@@ -822,19 +790,19 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>dsge_prior_weight {return token::DSGE_PRIOR_WEIGHT;}
 <DYNARE_BLOCK>surprise {return token::SURPRISE;}
 <DYNARE_BLOCK>bind {
-    yylval->build<string>(yytext);
+    yylval->emplace<string>(yytext);
     return token::BIND;
 }
 <DYNARE_BLOCK>relax {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::RELAX;
 }
 <DYNARE_BLOCK>error_bind {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::ERROR_BIND;
 }
 <DYNARE_BLOCK>error_relax {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::ERROR_RELAX;
 }
 <DYNARE_BLOCK>relative_to_initval {return token::RELATIVE_TO_INITVAL;}
@@ -848,22 +816,24 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK,DYNARE_STATEMENT>auxname {return token::AUXNAME;}
 <DYNARE_BLOCK>auxname_target_nonstationary {return token::AUXNAME_TARGET_NONSTATIONARY;}
 <DYNARE_BLOCK,DYNARE_STATEMENT>kind {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::KIND;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>ll {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::LL;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>dl {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::DL;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>dd {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::DD;
 }
 <DYNARE_BLOCK>weights {return token::WEIGHTS;}
+<DYNARE_BLOCK>exogenize {return token::EXOGENIZE;}
+<DYNARE_BLOCK>endogenize {return token::ENDOGENIZE;}
 
  /* Inside Dynare statement */
 <DYNARE_STATEMENT>solve_algo {return token::SOLVE_ALGO;}
@@ -875,7 +845,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>order {return token::ORDER;}
 <DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
 <DYNARE_STATEMENT>dr {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::DR;
  }
 <DYNARE_STATEMENT>lyapunov_complex_threshold {return token::LYAPUNOV_COMPLEX_THRESHOLD;}
@@ -939,11 +909,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>time_shift {return token::TIME_SHIFT;}
 <DYNARE_STATEMENT>structural {return token::STRUCTURAL;}
 <DYNARE_STATEMENT>true {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::TRUE;
 }
 <DYNARE_STATEMENT>false {
-    yylval->build<string>(yytext);
+    yylval->emplace<string>(yytext);
     return token::FALSE;
 }
 
@@ -1011,6 +981,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;}
@@ -1066,38 +1037,27 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>use_shock_groups {return token::USE_SHOCK_GROUPS;}
 <DYNARE_STATEMENT>colormap {return token::COLORMAP;}
 
-<DYNARE_STATEMENT,DYNARE_BLOCK>[a-z_][a-z0-9_]* {
-  yylval->build<string>(yytext);
+<DYNARE_STATEMENT,DYNARE_BLOCK>{NAME} {
+  yylval->emplace<string>(yytext);
   return token::NAME;
 }
 
-<DYNARE_STATEMENT,DYNARE_BLOCK>((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+) {
-  yylval->build<string>(yytext);
+<DYNARE_STATEMENT,DYNARE_BLOCK>{FLOAT_NUMBER} {
+  yylval->emplace<string>(yytext);
   return token::FLOAT_NUMBER;
 }
 
 <DYNARE_STATEMENT,DYNARE_BLOCK>[0-9]+ {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
   return token::INT_NUMBER;
 }
 
-<DATES_STATEMENT>\( { yylval->as<string>().append(yytext); dates_parens_nb++; }
-<DATES_STATEMENT>\) {
-                      yylval->as<string>().append(yytext);
-                      if (--dates_parens_nb == 0)
-                      {
-                        BEGIN DYNARE_STATEMENT;
-                        return token::DATES;
-                      }
-                    }
-<DATES_STATEMENT>.  { yylval->as<string>().append(yytext); }
-
 <DYNARE_BLOCK>\|e { return token::PIPE_E; }
 <DYNARE_BLOCK>\|x { return token::PIPE_X; }
 <DYNARE_BLOCK>\|p { return token::PIPE_P; }
 
 <DYNARE_STATEMENT,DYNARE_BLOCK>\'[^\']*\' {
-  yylval->build<string>(yytext + 1).pop_back();
+  yylval->emplace<string>(yytext + 1).pop_back();
   return token::QUOTED_STRING;
 }
 
@@ -1129,11 +1089,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
     element in initval (in which case Dynare recognizes the matrix name as an external
     function symbol), and may want to modify the matrix later with Matlab statements.
  */
-<INITIAL>[a-z_][a-z0-9_]* {
+<INITIAL>{NAME} {
   if (driver.symbol_exists_and_is_not_modfile_local_or_external_function(yytext))
     {
       BEGIN DYNARE_STATEMENT;
-      yylval->build<string>(yytext);
+      yylval->emplace<string>(yytext);
       return token::NAME;
     }
   else
@@ -1154,7 +1114,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
     be able to back out of the statement if we realize it's a native statement
     and move to the NATIVE context
  */
-<INITIAL>\[([[:space:]]*[a-z_][a-z0-9_]*[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*[a-z_][a-z0-9_]*[[:space:]]*){1}\] {
+<INITIAL>\[([[:space:]]*{NAME}[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*{NAME}[[:space:]]*){1}\] {
   string yytextcpy{yytext};
   yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), '['), yytextcpy.end());
   yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), ']'), yytextcpy.end());
@@ -1178,7 +1138,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
   if (dynare_statement)
     {
       BEGIN DYNARE_STATEMENT;
-      yylval->build<vector<string>>(val);
+      yylval->emplace<vector<string>>(val);
       return token::SYMBOL_VEC;
     }
 }
@@ -1223,7 +1183,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <NATIVE_COMMENT>"*/"[[:space:]]*\n   { BEGIN NATIVE; }
 <NATIVE_COMMENT>.
 
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,DATES_STATEMENT,LINE1,LINE2,LINE3,NATIVE_COMMENT><<EOF>> { yyterminate(); }
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,LINE1,LINE2,LINE3,NATIVE_COMMENT><<EOF>> { yyterminate(); }
 
 <*>.      { driver.error(*yylloc, "character unrecognized by lexer"); }
 %%

src/ExprNode.hh

 /*
- * Copyright © 2007-2024 Dynare Team
+ * Copyright © 2007-2025 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -283,7 +283,7 @@ protected:
   min_cost(bool is_matlab)
   {
     return is_matlab ? min_cost_matlab : min_cost_c;
-  };
+  }
 
   //! Initializes data member non_null_derivatives
   virtual void prepareForDerivation() = 0;
@@ -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;
@@ -946,7 +943,15 @@ public:
   /* Matches an expression that constitutes a complementarity condition.
      If successful, returns a triplet (endo_symb_id, lower_bound, upper_bound).
      Otherwise, throws a MatchFailureException. */
-  [[nodiscard]] virtual tuple<int, expr_t, expr_t> matchComplementarityCondition() const;
+  [[nodiscard]] virtual tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension = nullopt) const;
+
+  /* Replaces aggregation operators (e.g. SUM()) by new auxiliary variables.
+     Also declares those aggregation operators in the HeterogeneityTable, so as to
+     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 +1063,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
@@ -1165,6 +1172,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
@@ -1312,6 +1321,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
@@ -1507,7 +1518,11 @@ public:
                           vector<int>& powers) const override;
   [[nodiscard]] pair<int, expr_t> matchEndogenousTimesConstant() const override;
   [[nodiscard]] expr_t substituteLogTransform(int orig_symb_id, int aux_symb_id) const override;
-  [[nodiscard]] tuple<int, expr_t, expr_t> matchComplementarityCondition() const override;
+  [[nodiscard]] expr_t substituteAggregationOperators(subst_table_t& subst_table,
+                                                      vector<BinaryOpNode*>& neweqs) const override;
+  [[nodiscard]] tuple<int, expr_t, expr_t>
+  matchComplementarityCondition(const optional<int>& heterogeneity_dimension
+                                = nullopt) const override;
 };
 
 //! Trinary operator node
@@ -1651,6 +1666,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;
 };
 
@@ -1826,6 +1843,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
@@ -2028,6 +2047,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;

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;
+}

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);
+  [[nodiscard]] int getSummedHeterogenousEndogenousIndex(int symb_id) const;
+  [[nodiscard]] int aggregateEndoSize() const;
+
+  void writeOutput(ostream& output) const;
+  void writeJsonOutput(ostream& output) const;
+};
+
+#endif

src/HeterogeneousModel.cc

+/*
+ * Copyright © 2024-2025 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::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);
+
+  computeMCPEquationsReordering(heterogeneity_dimension);
+}
+
+void
+HeterogeneousModel::writeModelFiles(const string& basename, bool julia) const
+{
+  assert(!julia); // Not yet implemented
+  writeSparseModelMFiles<true>(basename, heterogeneity_dimension);
+  writeComplementarityConditionsFile<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
+{
+  std::vector<int> state_var;
+  for (int endoID = 0; endoID < symbol_table.het_endo_nbr(heterogeneity_dimension); endoID++)
+    try
+      {
+        getDerivID(symbol_table.getID(SymbolType::heterogeneousEndogenous, endoID,
+                                      heterogeneity_dimension),
+                   -1);
+        if (ranges::find(state_var, endoID) == state_var.end())
+          state_var.push_back(endoID);
+      }
+    catch (UnknownDerivIDException& e)
+      {
+      }
+
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1 << ").state_var = [";
+  for (int it : state_var)
+    output << it + 1 << " ";
+  output << "];" << endl;
+
+  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);
+  output << "M_.heterogeneity(" << heterogeneity_dimension + 1
+         << ").dynamic_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
+}

src/HeterogeneousModel.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 HETEROGENEOUS_MODEL_HH
+#define HETEROGENEOUS_MODEL_HH
+
+#include <string>
+
+#include "ModelTree.hh"
+
+using namespace std;
+
+class HeterogeneousModel : public ModelTree
+{
+public:
+  const int heterogeneity_dimension;
+
+  HeterogeneousModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
+                     ExternalFunctionsTable& external_functions_table_arg,
+                     HeterogeneityTable& heterogeneity_table_arg, int heterogeneity_dimension_arg);
+
+  HeterogeneousModel(const HeterogeneousModel& m) = default;
+  HeterogeneousModel& operator=(const HeterogeneousModel& m);
+
+  void computingPass(int derivsOrder, bool no_tmp_terms, bool use_dll);
+
+  void writeModelFiles(const string& basename, bool julia) const;
+  void writeDriverOutput(ostream& output) const;
+
+  [[nodiscard]] int getJacobianCol(int deriv_id, bool sparse) const override;
+  [[nodiscard]] int getJacobianColsNbr(bool sparse) const override;
+
+#if 0
+  void substituteEndoLeadGreaterThanTwo();
+
+  //! Transforms the model by removing all lags greater or equal than 2 on endos
+  void substituteEndoLagGreaterThanTwo();
+
+  //! Transforms the model by removing all leads on exos
+  /*! Note that this can create new lags on endos and exos */
+  void substituteExoLead();
+
+  //! Transforms the model by removing all lags on exos
+  void substituteExoLag();
+
+  //! Transforms the model by removing all UnaryOpcode::expectation
+  void substituteExpectation(bool partial_information_model);
+
+  //! Transforms the model by decreasing the lead/lag of predetermined variables in model equations
+  //! by one
+  void transformPredeterminedVariables();
+
+  //! Substitutes out all model-local variables
+  void substituteModelLocalVariables();
+#endif
+
+  // FIXME: the following 5 functions are identical to those in DynamicModel. Factorization?
+  [[nodiscard]] int getDerivID(int symb_id, int lead_lag) const noexcept(false) override;
+  [[nodiscard]] SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override;
+  [[nodiscard]] int getLagByDerivID(int deriv_id) const noexcept(false) override;
+  [[nodiscard]] int getSymbIDByDerivID(int deriv_id) const noexcept(false) override;
+  [[nodiscard]] int getTypeSpecificIDByDerivID(int deriv_id) const override;
+
+protected:
+  void computeChainRuleJacobian() override;
+  int getBlockJacobianEndoCol(int blk, int var, int lead_lag) const override;
+  string
+  modelClassName() const override
+  {
+    return "dynamic model for heterogeneity dimension '"
+           + heterogeneity_table.getName(heterogeneity_dimension) + "'";
+  }
+  int getMFS() const override;
+
+private:
+  // Maps a pair (symbol ID, lead/lag) to a deriv ID
+  map<pair<int, int>, int> deriv_id_table;
+  // Maps a deriv ID to a pair (symbol ID, lead/lag)
+  vector<pair<int, int>> inv_deriv_id_table;
+
+  // Allocates the derivation IDs for all endogenous variables for this heterogeneity dimension
+  void computeDerivIDs();
+};
+
+#endif