Occbin: new syntax for constructing multi-regime equations

Ref. #68

src/ComputingTasks.cc

@@ -25,6 +25,7 @@ using namespace std;
 
 #include "ComputingTasks.hh"
 #include "Statement.hh"
+#include "ParsingDriver.hh"
 
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wold-style-cast"
@@ -5262,7 +5263,7 @@ OccbinConstraintsStatement::writeOutput(ostream &output, const string &basename,
   output << "M_.occbin.constraint_nbr = " << constraints.size() << ';' << endl
          << "M_.occbin.pswitch = [" << endl;
   for (const auto &[name, bind, relax, error_bind, error_relax] : constraints)
-    output << symbol_table.getTypeSpecificID(name) + 1 << ' ';
+    output << symbol_table.getTypeSpecificID(ParsingDriver::buildOccbinBindParamName(name)) + 1 << ' ';
   output << "];" << endl
          << "options_.occbin = struct();" << endl
          << "options_.occbin = occbin.set_default_options(options_.occbin, M_);" << endl

src/DynamicModel.cc

@@ -27,6 +27,7 @@
 #include <sstream>
 
 #include "DynamicModel.hh"
+#include "ParsingDriver.hh"
 
 void
 DynamicModel::copyHelper(const DynamicModel &m)
@@ -5931,6 +5932,67 @@ DynamicModel::dynamicOnlyEquationsNbr() const
   return equation_tags.getDynamicEqns().size();
 }
 
+void
+DynamicModel::addOccbinEquation(expr_t eq, int lineno, const map<string, string> &eq_tags, const vector<string> &regimes_bind, const vector<string> &regimes_relax)
+{
+  auto beq = dynamic_cast<BinaryOpNode *>(eq);
+  assert(beq && beq->op_code == BinaryOpcode::equal);
+
+  // Construct the term to be added to the corresponding equation
+  expr_t basic_term = AddMinus(beq->arg1, beq->arg2);
+  expr_t term = basic_term;
+  for (auto &regime : regimes_bind)
+    {
+      int param_id = symbol_table.getID(ParsingDriver::buildOccbinBindParamName(regime));
+      term = AddTimes(term, AddVariable(param_id));
+    }
+  for (auto &regime : regimes_relax)
+    {
+      int param_id = symbol_table.getID(ParsingDriver::buildOccbinBindParamName(regime));
+      term = AddTimes(term, AddMinus(One, AddVariable(param_id)));
+    }
+
+  // Create or update the dynamic equation
+  try
+    {
+      int eqn = equation_tags.getEqnByTag("name", eq_tags.at("name"));
+      BinaryOpNode *orig_eq = equations[eqn];
+      /* In the following, we could have kept only orig_eq->arg1, but the
+         following adds a (somewhat bizarre) support for equation snippets
+         without “bind” nor “relax” */
+      equations[eqn] = AddEqual(AddPlus(AddMinus(orig_eq->arg1, orig_eq->arg2), term), Zero);
+      // It’s unclear how to update lineno and tags, so don’t do it
+    }
+  catch (EquationTags::TagNotFoundException &e)
+    {
+      auto eq_tags_dynamic = eq_tags;
+      eq_tags_dynamic["dynamic"] = "";
+      addEquation(AddEqual(term, Zero), lineno, eq_tags_dynamic);
+    }
+
+  // Create or update the static equation (corresponding to the pure relax regime)
+  if (regimes_bind.empty())
+    {
+      try
+        {
+          /* Similar remark as above. We could have entirely skipped this
+             equation updating, since normally there is only one such clause,
+             but the following adds a (somewhat bizarre) support for equation
+             snippets without “bind” nor “relax” */
+          int eqn = static_only_equations_equation_tags.getEqnByTag("name", eq_tags.at("name"));
+          BinaryOpNode *orig_eq = static_only_equations[eqn];
+          static_only_equations[eqn] = AddEqual(AddPlus(AddMinus(orig_eq->arg1, orig_eq->arg2), basic_term), Zero);
+          // It’s unclear how to update lineno and tags, so don’t do it
+        }
+      catch (EquationTags::TagNotFoundException &e)
+        {
+          auto eq_tags_static = eq_tags;
+          eq_tags_static["static"] = "";
+          addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, eq_tags_static);
+        }
+    }
+}
+
 bool
 DynamicModel::isChecksumMatching(const string &basename, bool block) const
 {

src/DynamicModel.hh

@@ -444,6 +444,13 @@ public:
   //! Returns number of dynamic only equations
   size_t dynamicOnlyEquationsNbr() const;
 
+  // Adds an occbin equation (with “bind” and/or “relax” tag)
+  /* This function assumes that there is a “name” tag, and that the relevant
+     auxiliary parameters have already been added to the symbol table.
+     It also assumes that the “bind” and “relax” tags have been cleared from
+     eq_tags. */
+  void addOccbinEquation(expr_t eq, int lineno, const map<string, string> &eq_tags, const vector<string> &regimes_bind, const vector<string> &regimes_relax);
+
   //! Writes LaTeX file with the equations of the dynamic model
   void writeLatexFile(const string &basename, bool write_equation_tags) const;
 

src/ParsingDriver.cc

@@ -2338,15 +2338,46 @@ ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2)
 {
   expr_t id = model_tree->AddEqual(arg1, arg2);
 
-  // Detect if the equation is tagged [static]
   if (eq_tags.find("static") != eq_tags.end())
     {
+      // If the equation is tagged [static]
       if (!id->isInStaticForm())
         error("An equation tagged [static] cannot contain leads, lags, expectations or STEADY_STATE operators");
 
       dynamic_model->addStaticOnlyEquation(id, location.begin.line, eq_tags);
     }
-  else
+  else if (eq_tags.find("bind") != eq_tags.end()
+           || eq_tags.find("relax") != eq_tags.end())
+    {
+      // If the equation has a “bind” or “relax” tag (occbin case)
+      if (eq_tags.find("name") == eq_tags.end())
+        error("An equation with a 'bind' or 'relax' tag must have a 'name' tag");
+      auto regimes_bind = strsplit(eq_tags["bind"], ',');
+      auto regimes_relax = strsplit(eq_tags["relax"], ',');
+      auto regimes_all = regimes_bind;
+      regimes_all.insert(regimes_all.end(), regimes_relax.begin(), regimes_relax.end()); // Concatenate the two vectors
+      for (const auto &regime : regimes_all)
+        {
+          if (!isSymbolIdentifier(regime))
+            error("The string '" + regime + "' is not a valid Occbin regime name (contains unauthorized characters)");
+          string param_name = buildOccbinBindParamName(regime);
+          try
+            {
+              if (mod_file->symbol_table.getType(param_name) != SymbolType::parameter)
+                error("The name '" + param_name + "' is already used. Please use another name for Occbin regime '" + regime + "'");
+            }
+          catch (SymbolTable::UnknownSymbolNameException &e)
+            {
+              // Declare and initialize the new parameter
+              int symb_id = mod_file->symbol_table.addSymbol(param_name, SymbolType::parameter);
+              mod_file->addStatement(make_unique<InitParamStatement>(symb_id, dynamic_model->Zero, mod_file->symbol_table));
+            }
+        }
+      eq_tags.erase("bind");
+      eq_tags.erase("relax");
+      dynamic_model->addOccbinEquation(id, location.begin.line, eq_tags, regimes_bind, regimes_relax);
+    }
+  else // General case
     model_tree->addEquation(id, location.begin.line, eq_tags);
 
   eq_tags.clear();
@@ -3428,7 +3459,9 @@ ParsingDriver::end_occbin_constraints(const vector<tuple<string, BinaryOpNode *,
   // Perform a few checks
   for (const auto &[name, bind, relax, error_bind, error_relax] : constraints)
     {
-      check_symbol_is_parameter(name);
+      string param_name = buildOccbinBindParamName(name);
+      if (!mod_file->symbol_table.exists(param_name))
+        error("No equation has been declared for regime '" + name + "'");
       if (!bind)
         error("The 'bind' expression is missing in constraint '" + name + "'");
       if (bind->hasExogenous())
@@ -3445,3 +3478,40 @@ ParsingDriver::end_occbin_constraints(const vector<tuple<string, BinaryOpNode *,
 
   reset_data_tree();
 }
+
+vector<string>
+ParsingDriver::strsplit(const string &str, char delim)
+{
+  vector<string> result;
+  size_t idx = 0;
+  while (idx < str.size())
+    {
+      size_t idx2 = str.find(delim, idx);
+      if (idx2 == string::npos)
+        {
+          result.push_back(str.substr(idx));
+          break;
+        }
+      result.push_back(str.substr(idx, idx2-idx));
+      idx = idx2 + 1;
+    }
+  return result;
+}
+
+bool
+ParsingDriver::isSymbolIdentifier(const string &str)
+{
+  if (str.empty())
+    return false;
+  auto myisalpha = [](char ch)
+  {
+    // We cannot use std::isalpha(), because it is locale-dependent
+    return (ch >= 'a' && ch <= 'z') || (ch >= 'A' && ch <= 'Z');
+  };
+  if (!(myisalpha(str[0]) || str[0] == '_'))
+    return false;
+  for (size_t i = 1; i < str.size(); i++)
+    if (!(myisalpha(str[i]) || isdigit(static_cast<unsigned char>(str[i])) || str[i] == '_'))
+      return false;
+  return true;
+}

src/ParsingDriver.hh

@@ -888,6 +888,15 @@ public:
   void begin_occbin_constraints();
   //! Add an occbin_constraints block
   void end_occbin_constraints(const vector<tuple<string, BinaryOpNode *, BinaryOpNode *, expr_t, expr_t>> &constraints);
+  // Equivalent of MATLAB’s strsplit. Returns an empty vector given an empty string.
+  static vector<string> strsplit(const string &str, char delim);
+  // Returns true iff the string is a legal symbol identifier (see NAME token in lexer)
+  static bool isSymbolIdentifier(const string &str);
+  // Given an Occbin regime name, returns the corresponding auxiliary parameter
+  static string buildOccbinBindParamName(const string &regime)
+  {
+    return "occbin_" + regime + "_bind";
+  }
 };
 
 #endif // ! PARSING_DRIVER_HH