From b4ac9aac8da9df9b8da41935b7821c1e2c9fa4cc Mon Sep 17 00:00:00 2001
From: Houtan Bastani <houtan@dynare.org>
Date: Mon, 6 Mar 2017 16:34:08 +0100
Subject: [PATCH] preprocessor: add lag info to equation cross references.
 closes #1398

---
 DynamicModel.cc | 110 ++++++++++++++++++++++++++++++++++++++++++++++--
 DynamicModel.hh |  23 +++++++++-
 DynareMain2.cc  |   4 +-
 ExprNode.cc     |  10 ++---
 ExprNode.hh     |  10 ++---
 ModFile.cc      |  15 ++++---
 ModFile.hh      |   6 +--
 ModelTree.cc    |  97 ------------------------------------------
 ModelTree.hh    |  15 -------
 9 files changed, 152 insertions(+), 138 deletions(-)

diff --git a/DynamicModel.cc b/DynamicModel.cc
index f3c90c25..b914eff6 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3164,7 +3164,7 @@ DynamicModel::runTrendTest(const eval_context_t &eval_context)
 void
 DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
-                            bool bytecode, bool compute_xrefs)
+                            bool bytecode)
 {
   assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder));
 
@@ -3272,9 +3272,113 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
         if (bytecode)
           computeTemporaryTermsMapping();
       }
+}
 
-  if (compute_xrefs)
-    computeXrefs();
+void
+DynamicModel::computeXrefs()
+{
+  int i = 0;
+  for (vector<BinaryOpNode *>::iterator it = equations.begin();
+       it != equations.end(); it++)
+    {
+      ExprNode::EquationInfo ei;
+      (*it)->computeXrefs(ei);
+      xrefs[i++] = ei;
+    }
+
+  i = 0;
+  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
+       it != xrefs.end(); it++, i++)
+    {
+      computeRevXref(xref_param, it->second.param, i);
+      computeRevXref(xref_endo, it->second.endo, i);
+      computeRevXref(xref_exo, it->second.exo, i);
+      computeRevXref(xref_exo_det, it->second.exo_det, i);
+    }
+}
+
+void
+DynamicModel::computeRevXref(map<pair<int, int>, set<int> > &xrefset, const set<pair<int, int> > &eiref, int eqn)
+{
+  for (set<pair<int, int> >::const_iterator it = eiref.begin();
+       it != eiref.end(); it++)
+    {
+      set<int> eq;
+      if (xrefset.find(*it) != xrefset.end())
+        eq = xrefset[*it];
+      eq.insert(eqn);
+      xrefset[*it] = eq;
+    }
+}
+
+void
+DynamicModel::writeXrefs(ostream &output) const
+{
+  output << "M_.xref1.param = cell(1, M_.eq_nbr);" << endl
+         << "M_.xref1.endo = cell(1, M_.eq_nbr);" << endl
+         << "M_.xref1.exo = cell(1, M_.eq_nbr);" << endl
+         << "M_.xref1.exo_det = cell(1, M_.eq_nbr);" << endl;
+  int i = 1;
+  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
+       it != xrefs.end(); it++, i++)
+    {
+      output << "M_.xref1.param{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.param.begin();
+           it1 != it->second.param.end(); it1++)
+        output << symbol_table.getTypeSpecificID(it1->first) + 1 << " ";
+      output << "];" << endl;
+
+      output << "M_.xref1.endo{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.endo.begin();
+           it1 != it->second.endo.end(); it1++)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      output << "];" << endl;
+
+      output << "M_.xref1.exo{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.exo.begin();
+           it1 != it->second.exo.end(); it1++)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      output << "];" << endl;
+
+      output << "M_.xref1.exo_det{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.exo_det.begin();
+           it1 != it->second.exo_det.end(); it1++)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      output << "];" << endl;
+    }
+
+  output << "M_.xref2.param = cell(1, M_.param_nbr);" << endl
+         << "M_.xref2.endo = cell(1, M_.endo_nbr);" << endl
+         << "M_.xref2.exo = cell(1, M_.exo_nbr);" << endl
+         << "M_.xref2.exo_det = cell(1, M_.exo_det_nbr);" << endl;
+  writeRevXrefs(output, xref_param, "param");
+  writeRevXrefs(output, xref_endo, "endo");
+  writeRevXrefs(output, xref_exo, "exo");
+  writeRevXrefs(output, xref_exo_det, "exo_det");
+}
+
+void
+DynamicModel::writeRevXrefs(ostream &output, const map<pair<int, int>, set<int> > &xrefmap, const string &type) const
+{
+  int last_tsid = -1;
+  for (map<pair<int, int>, set<int> >::const_iterator it = xrefmap.begin();
+       it != xrefmap.end(); it++)
+    {
+      int tsid = symbol_table.getTypeSpecificID(it->first.first) + 1;
+      output << "M_.xref2." << type << "{" << tsid << "} = [ ";
+      if (last_tsid == tsid)
+        output << "M_.xref2." << type << "{" << tsid << "}; ";
+      else
+        last_tsid = tsid;
+
+      for (set<int>::const_iterator it1 = it->second.begin();
+           it1 != it->second.end(); it1++)
+        if (type == "param")
+          output << *it1 + 1 << " ";
+        else
+          output << "struct('shift', " << it->first.second << ", 'eq', " << *it1+1 << ");";
+      output << "];" << endl;
+    }
 }
 
 map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 6b99423c..7f0177fb 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -62,6 +62,13 @@ private:
   /*! Set by computeDerivIDs() */
   int max_exo_det_lag, max_exo_det_lead;
 
+  //! Cross reference information
+  map<int, ExprNode::EquationInfo> xrefs;
+  map<pair<int, int>, set<int> > xref_param;
+  map<pair<int, int>, set<int> > xref_endo;
+  map<pair<int, int>, set<int> > xref_exo;
+  map<pair<int, int>, set<int> > xref_exo_det;
+
   //! Number of columns of dynamic jacobian
   /*! Set by computeDerivID()s and computeDynJacobianCols() */
   int dynJacobianColsNbr;
@@ -189,6 +196,12 @@ private:
   /*! pair< pair<static, forward>, pair<backward,mixed> > */
   vector<pair< pair<int, int>, pair<int, int> > > block_col_type;
 
+  //! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
+  void computeRevXref(map<pair<int, int>, set<int> > &xrefset, const set<pair<int, int> > &eiref, int eqn);
+
+  //! Write reverse cross references
+  void writeRevXrefs(ostream &output, const map<pair<int, int>, set<int> > &xrefmap, const string &type) const;
+
   //! List for each variable its block number and its maximum lag and lead inside the block
   vector<pair<int, pair<int, int> > > variable_block_lead_lag;
   //! List for each equation its block number
@@ -202,7 +215,13 @@ public:
   //! Adds a variable node
   /*! This implementation allows for non-zero lag */
   virtual VariableNode *AddVariable(int symb_id, int lag = 0);
-  
+
+  //! Compute cross references
+  void computeXrefs();
+
+  //! Write cross references
+  void writeXrefs(ostream &output) const;
+
   //! Execute computations (variable sorting + derivation)
   /*!
     \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
@@ -213,7 +232,7 @@ public:
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
   void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
-                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool compute_xrefs);
+                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
 
diff --git a/DynareMain2.cc b/DynareMain2.cc
index bc5065b2..b5d1e006 100644
--- a/DynareMain2.cc
+++ b/DynareMain2.cc
@@ -45,13 +45,13 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
   mod_file->checkPass(nostrict);
 
   // Perform transformations on the model (creation of auxiliary vars and equations)
-  mod_file->transformPass(nostrict);
+  mod_file->transformPass(nostrict, compute_xrefs);
 
   // Evaluate parameters initialization, initval, endval and pounds
   mod_file->evalAllExpressions(warn_uninit);
 
   // Do computations
-  mod_file->computingPass(no_tmp_terms, output_mode, compute_xrefs, params_derivs_order);
+  mod_file->computingPass(no_tmp_terms, output_mode, params_derivs_order);
 
   // Write outputs
   if (output_mode != none)
diff --git a/ExprNode.cc b/ExprNode.cc
index d7d8374c..4d723fa7 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2016 Dynare Team
+ * Copyright (C) 2007-2017 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -1069,16 +1069,16 @@ VariableNode::computeXrefs(EquationInfo &ei) const
   switch (type)
     {
     case eEndogenous:
-      ei.endo.insert(symb_id);
+      ei.endo.insert(make_pair(symb_id, lag));
       break;
     case eExogenous:
-      ei.exo.insert(symb_id);
+      ei.exo.insert(make_pair(symb_id, lag));
       break;
     case eExogenousDet:
-      ei.exo_det.insert(symb_id);
+      ei.exo_det.insert(make_pair(symb_id, lag));
       break;
     case eParameter:
-      ei.param.insert(symb_id);
+      ei.param.insert(make_pair(symb_id, 0));
       break;
     case eTrend:
     case eLogTrend:
diff --git a/ExprNode.hh b/ExprNode.hh
index 0440d807..e784c60d 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2016 Dynare Team
+ * Copyright (C) 2007-2017 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -165,10 +165,10 @@ protected:
   //! For creating equation cross references
   struct EquationInfo
   {
-    set<int> param;
-    set<int> endo;
-    set<int> exo;
-    set<int> exo_det;
+    set<pair<int, int> > param;
+    set<pair<int, int> > endo;
+    set<pair<int, int> > exo;
+    set<pair<int, int> > exo_det;
   };
 
 public:
diff --git a/ModFile.cc b/ModFile.cc
index 794f1e46..5719d019 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -319,7 +319,7 @@ ModFile::checkPass(bool nostrict)
 }
 
 void
-ModFile::transformPass(bool nostrict)
+ModFile::transformPass(bool nostrict, bool compute_xrefs)
 {
   // Save the original model (must be done before any model transformations by preprocessor)
   dynamic_model.cloneDynamic(original_model);
@@ -414,6 +414,9 @@ ModFile::transformPass(bool nostrict)
   // Freeze the symbol table
   symbol_table.freeze();
 
+  if (compute_xrefs)
+    dynamic_model.computeXrefs();
+
   /*
     Enforce the same number of equations and endogenous, except in three cases:
     - ramsey_model, ramsey_policy or discretionary_policy is used
@@ -483,7 +486,7 @@ ModFile::transformPass(bool nostrict)
 }
 
 void
-ModFile::computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xrefs, int params_derivs_order)
+ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order)
 {
   // Mod file may have no equation (for example in a standalone BVAR estimation)
   if (dynamic_model.equation_number() > 0)
@@ -517,7 +520,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xr
 	  || mod_file_struct.calib_smoother_present)
 	{
 	  if (mod_file_struct.perfect_foresight_solver_present)
-	    dynamic_model.computingPass(true, false, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+	    dynamic_model.computingPass(true, false, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
 	      else
 		{
 		  if (mod_file_struct.stoch_simul_present
@@ -542,13 +545,13 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xr
                   int paramsDerivsOrder = 0;
                   if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
                     paramsDerivsOrder = params_derivs_order;
-		  dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+		  dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
                   if (linear && mod_file_struct.ramsey_model_present)
-                    orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+                    orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
 		}
 	    }
 	  else // No computing task requested, compute derivatives up to 2nd order by default
-	    dynamic_model.computingPass(true, true, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+	    dynamic_model.computingPass(true, true, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
 
       if ((linear && !mod_file_struct.ramsey_model_present && !dynamic_model.checkHessianZero()) ||
           (linear && mod_file_struct.ramsey_model_present && !orig_ramsey_dynamic_model.checkHessianZero()))
diff --git a/ModFile.hh b/ModFile.hh
index 6ed24d33..68df1adf 100644
--- a/ModFile.hh
+++ b/ModFile.hh
@@ -130,12 +130,12 @@ public:
   /*! \todo add check for number of equations and endogenous if ramsey_policy is present */
   void checkPass(bool nostrict);
   //! Perform some transformations on the model (creation of auxiliary vars and equations)
-  void transformPass(bool nostrict);
+  /*! \param compute_xrefs if true, equation cross references will be computed */
+  void transformPass(bool nostrict, bool compute_xrefs);
   //! Execute computations
   /*! \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
-  /*! \param compute_xrefs if true, equation cross references will be computed */
   /*! \param params_derivs_order compute this order of derivs wrt parameters */
-  void computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xrefs, int params_derivs_order);
+  void computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order);
   //! Writes Matlab/Octave output files
   /*!
     \param basename The base name used for writing output files. Should be the name of the mod file without its extension
diff --git a/ModelTree.cc b/ModelTree.cc
index f295423e..0b88c298 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -229,103 +229,6 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
     }
 }
 
-void
-ModelTree::computeXrefs()
-{
-  int i = 0;
-  for (vector<BinaryOpNode *>::iterator it = equations.begin();
-       it != equations.end(); it++)
-    {
-      ExprNode::EquationInfo ei;
-      (*it)->computeXrefs(ei);
-      xrefs[i++] = ei;
-    }
-
-  i = 0;
-  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
-       it != xrefs.end(); it++, i++)
-    {
-      computeRevXref(xref_param, it->second.param, i);
-      computeRevXref(xref_endo, it->second.endo, i);
-      computeRevXref(xref_exo, it->second.exo, i);
-      computeRevXref(xref_exo_det, it->second.exo_det, i);
-    }
-}
-
-void
-ModelTree::computeRevXref(map<int, set<int> > &xrefset, const set<int> &eiref, int eqn)
-{
-  for (set<int>::const_iterator it1 = eiref.begin();
-       it1 != eiref.end(); it1++)
-    {
-      set<int> eq;
-      if (xrefset.find(symbol_table.getTypeSpecificID(*it1)) != xrefset.end())
-        eq = xrefset[symbol_table.getTypeSpecificID(*it1)];
-      eq.insert(eqn);
-      xrefset[symbol_table.getTypeSpecificID(*it1)] = eq;
-    }
-}
-
-void
-ModelTree::writeXrefs(ostream &output) const
-{
-  output << "M_.xref1.param = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref1.endo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref1.exo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref1.exo_det = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.param = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.endo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.exo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.exo_det = cell(1, M_.eq_nbr);" << endl;
-  int i = 1;
-  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
-       it != xrefs.end(); it++, i++)
-    {
-      output << "M_.xref1.param{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.param.begin();
-           it1 != it->second.param.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.xref1.endo{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.endo.begin();
-           it1 != it->second.endo.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.xref1.exo{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.exo.begin();
-           it1 != it->second.exo.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.xref1.exo_det{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.exo_det.begin();
-           it1 != it->second.exo_det.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-    }
-
-  writeRevXrefs(output, xref_param, "param");
-  writeRevXrefs(output, xref_endo, "endo");
-  writeRevXrefs(output, xref_exo, "exo");
-  writeRevXrefs(output, xref_exo_det, "exo_det");
-}
-
-void
-ModelTree::writeRevXrefs(ostream &output, const map<int, set<int> > &xrefmap, const string &type) const
-{
-  for (map<int, set<int> >::const_iterator it = xrefmap.begin();
-       it != xrefmap.end(); it++)
-    {
-      output << "M_.xref2." << type << "{" << it->first + 1 << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        output << *it1 + 1 << " ";
-      output << "];" << endl;
-    }
-}
-
 void
 ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
 {
diff --git a/ModelTree.hh b/ModelTree.hh
index 9b28fce4..b87debf1 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -100,13 +100,6 @@ protected:
   */
   first_derivatives_t residuals_params_derivatives;
 
-  //! Cross reference information
-  map<int, ExprNode::EquationInfo> xrefs;
-  map<int, set<int> > xref_param;
-  map<int, set<int> > xref_endo;
-  map<int, set<int> > xref_exo;
-  map<int, set<int> > xref_exo_det;
-
   //! Second derivatives of the residuals w.r. to parameters
   /*! First index is equation number, second and third indeces are parameters.
     Only non-null derivatives are stored in the map.
@@ -240,14 +233,6 @@ protected:
 
   //! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
   void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
-  //! Compute cross references
-  void computeXrefs();
-  //! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
-  void computeRevXref(map<int, set<int> > &xrefset, const set<int> &eiref, int eqn);
-  //! Write cross references
-  void writeXrefs(ostream &output) const;
-  //! Write reverse cross references
-  void writeRevXrefs(ostream &output, const map<int, set<int> > &xrefmap, const string &type) const;
   //! Evaluate the jacobian and suppress all the elements below the cutoff
   void evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose);
   //! Search the equations and variables belonging to the prologue and the epilogue of the model
-- 
GitLab