diff: substitute auxiliary variables starting from the least lagged

src/DynamicModel.cc

@@ -4899,17 +4899,27 @@ DynamicModel::substituteAdl()
 void
 DynamicModel::substituteDiff(StaticModel &static_model)
 {
-  // Substitute in model local variables
+  // Find diff Nodes
-  vector<BinaryOpNode *> neweqs;
+  diff_table_t diff_table;
-  diff_subst_table_t diff_subst_table;
   for (map<int, expr_t>::iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
-    it->second = it->second->substituteDiff(static_model, diff_subst_table, neweqs);
+    it->second->findDiffNodes(static_model, diff_table);
+
+  for (int i = 0; i < (int) equations.size(); i++)
+    equations[i]->findDiffNodes(static_model, diff_table);
+
+  // Substitute in model local variables
+  vector<BinaryOpNode *> neweqs;
+  ExprNode::subst_table_t diff_subst_table;
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second = it->second->substituteDiff(static_model, diff_table, diff_subst_table, neweqs);
 
   // Substitute in equations
   for (int i = 0; i < (int) equations.size(); i++)
     {
-      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteDiff(static_model, diff_subst_table, neweqs));
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->
+                                                           substituteDiff(static_model, diff_table, diff_subst_table, neweqs));
       assert(substeq != NULL);
       equations[i] = substeq;
     }

src/ExprNode.cc

@@ -505,8 +505,13 @@ NumConstNode::substituteAdl() const
   return const_cast<NumConstNode *>(this);
 }
 
+void
+NumConstNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+{
+}
+
 expr_t
-NumConstNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+NumConstNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -1310,8 +1315,13 @@ VariableNode::substituteAdl() const
   return const_cast<VariableNode *>(this);
 }
 
+void
+VariableNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+{
+}
+
 expr_t
-VariableNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+VariableNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VariableNode *>(this);
 }
@@ -2857,38 +2867,83 @@ UnaryOpNode::isDiffPresent() const
   return arg->isDiffPresent();
 }
 
+void
+UnaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+{
+  if (op_code != oDiff)
+    return;
+
+  arg->findDiffNodes(static_datatree, diff_table);
+
+  expr_t sthis = this->toStatic(static_datatree);
+  int arg_max_lag = -arg->maxLag();
+  diff_table_t::iterator it = diff_table.find(sthis);
+  if (it != diff_table.end())
+    {
+      for (map<int, expr_t>::const_iterator it1 = it->second.begin();
+           it1 != it->second.end(); it1++)
+        if (arg == it1->second)
+          return;
+      it->second[arg_max_lag] = const_cast<UnaryOpNode *>(this);
+    }
+  else
+    diff_table[sthis][arg_max_lag] = const_cast<UnaryOpNode *>(this);
+}
+
 expr_t
-UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   if (op_code != oDiff)
     {
-      expr_t argsubst = arg->substituteDiff(static_datatree, diff_subst_table, neweqs);
+      expr_t argsubst = arg->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
 
-  expr_t sarg = arg->toStatic(static_datatree);
+  subst_table_t::const_iterator sit = subst_table.find(this);
-  int arg_max_lag = -arg->maxLag();
+  if (sit != subst_table.end())
-  diff_subst_table_t::iterator it = diff_subst_table.find(sarg);
+    return const_cast<VariableNode *>(sit->second);
-  if (it != diff_subst_table.end())
-    if (it->second.first.first == arg)
-      return const_cast<VariableNode *>(it->second.second);
-    else
-      return it->second.second->decreaseLeadsLags(it->second.first.second - arg_max_lag);
 
-  int symb_id = -1;
+  expr_t sthis = dynamic_cast<UnaryOpNode *>(this->toStatic(static_datatree));
-  expr_t argsubst = arg->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  diff_table_t::iterator it = diff_table.find(sthis);
+  assert(it != diff_table.end() && it->second[-arg->maxLag()] == this);
+
+  int origin_arg_max_lag;
+  VariableNode *origin_aux_var;
+  for (map<int, expr_t>::reverse_iterator rit = it->second.rbegin();
+       rit != it->second.rend(); rit++)
+    {
+      expr_t argsubst = dynamic_cast<UnaryOpNode *>(rit->second)->
+          get_arg()->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+
+      int symb_id;
       VariableNode *vn = dynamic_cast<VariableNode *>(argsubst);
       if (vn != NULL)
         symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
       else
         symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
 
-  //AUX_DIFF_IDX = argsubst - argsubst(-1)
+      if (rit == it->second.rbegin())
+        {
+          // make originating aux var & equation
+          origin_arg_max_lag = rit->first;
+          origin_aux_var = datatree.AddVariable(symb_id, 0);
+          //ORIG_AUX_DIFF = argsubst - argsubst(-1)
+          neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(origin_aux_var,
+                                                                          datatree.AddMinus(argsubst,
+                                                                                            argsubst->decreaseLeadsLags(1)))));
+          subst_table[rit->second] = dynamic_cast<VariableNode *>(origin_aux_var);
+        }
+      else
+        {
+          // just add equation of form: AUX_DIFF = ORIG_AUX_DIFF(origin_arg_max_lag - rit->first)
           VariableNode *new_aux_var = datatree.AddVariable(symb_id, 0);
           neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(new_aux_var,
-                                                                  datatree.AddMinus(argsubst, argsubst->decreaseLeadsLags(1)))));
+                                                                          origin_aux_var->decreaseLeadsLags(origin_arg_max_lag
-  diff_subst_table[sarg] = make_pair(make_pair(arg, arg_max_lag), new_aux_var);
+                                                                                                                - rit->first))));
-  return new_aux_var;
+          subst_table[rit->second] = dynamic_cast<VariableNode *>(new_aux_var);
+        }
+    }
+  return const_cast<VariableNode *>(subst_table[this]);
 }
 
 expr_t
@@ -4525,11 +4580,18 @@ BinaryOpNode::substituteAdl() const
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
-expr_t
+void
-BinaryOpNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+BinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
-  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  arg1->findDiffNodes(static_datatree, diff_table);
-  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  arg2->findDiffNodes(static_datatree, diff_table);
+}
+
+expr_t
+BinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -5386,12 +5448,20 @@ TrinaryOpNode::substituteAdl() const
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
-expr_t
+void
-TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+TrinaryOpNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
 {
-  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  arg1->findDiffNodes(static_datatree, diff_table);
-  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  arg2->findDiffNodes(static_datatree, diff_table);
-  expr_t arg3subst = arg3->substituteDiff(static_datatree, diff_subst_table, neweqs);
+  arg3->findDiffNodes(static_datatree, diff_table);
+}
+
+expr_t
+TrinaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteDiff(static_datatree, diff_table, subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -5756,12 +5826,19 @@ AbstractExternalFunctionNode::substituteAdl() const
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+void
+AbstractExternalFunctionNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->findDiffNodes(static_datatree, diff_table);
+}
+
 expr_t
-AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+AbstractExternalFunctionNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteDiff(static_datatree, diff_subst_table, neweqs));
+    arguments_subst.push_back((*it)->substituteDiff(static_datatree, diff_table, subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -7285,8 +7362,13 @@ VarExpectationNode::substituteAdl() const
   return const_cast<VarExpectationNode *>(this);
 }
 
+void
+VarExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+{
+}
+
 expr_t
-VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+VarExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -7708,8 +7790,13 @@ PacExpectationNode::substituteAdl() const
   return const_cast<PacExpectationNode *>(this);
 }
 
+void
+PacExpectationNode::findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const
+{
+}
+
 expr_t
-PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const
+PacExpectationNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<PacExpectationNode *>(this);
 }

src/ExprNode.hh

@@ -34,6 +34,7 @@ using namespace std;
 
 class DataTree;
 class VariableNode;
+class UnaryOpNode;
 class BinaryOpNode;
 class PacExpectationNode;
 
@@ -58,7 +59,7 @@ typedef map<int, double> eval_context_t;
 typedef map<pair<int, vector<expr_t> >, int> deriv_node_temp_terms_t;
 
 //! Type for the substitution map used in the process of substitutitng diff expressions
-typedef map<expr_t, pair<pair<expr_t, int>, VariableNode *> > diff_subst_table_t;
+typedef map<expr_t, map<int, expr_t> > diff_table_t;
 
 //! Possible types of output when writing ExprNode(s)
 enum ExprNodeOutputType
@@ -472,7 +473,8 @@ class ExprNode
       virtual expr_t substituteAdl() const = 0;
 
       //! Substitute diff operator
-      virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+      virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const = 0;
+      virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
       virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
@@ -558,7 +560,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -640,7 +643,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -744,7 +748,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -864,7 +869,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -961,7 +967,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
@@ -1058,7 +1065,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -1246,7 +1254,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
@@ -1322,7 +1331,8 @@ public:
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
   virtual expr_t substituteAdl() const;
-  virtual expr_t substituteDiff(DataTree &static_datatree, diff_subst_table_t &diff_subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual void findDiffNodes(DataTree &static_datatree, diff_table_t &diff_table) const;
+  virtual expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,