v4 preprocessor:

* fixed bug related to column numbering of exogenous variables in dynamic jacobian
* gave more explicit names to some methods in VariableTable, and other cosmetic changes


git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1990 ac1d8469-bf42-47a9-8791-bf33cf982152

ExprNode.cc

@@ -276,7 +276,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         {
         case oMatlabDynamicModel:
         case oCDynamicModel:
-          i = datatree.variable_table.getSortID(var_id) + OFFSET(output_type);
+          i = datatree.variable_table.getDynJacobianCol(var_id) + OFFSET(output_type);
           output <<  "y" << LPAR(output_type) << i << RPAR(output_type);
           break;
         case oMatlabStaticModel:

ModelTree.cc

@@ -1694,14 +1694,10 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
                     << "  g1 = NULL;" << endl
                     << "  if (nlhs >= 2)" << endl
                     << "  {" << endl
-                    << "     /* Set the output pointer to the output matrix g1. */" << endl;
+                    << "     /* Set the output pointer to the output matrix g1. */" << endl
 
-  if (computeJacobianExo)
-    mDynamicModelFile << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr() << ", mxREAL);" << endl;
-  else if (computeJacobian)
-    mDynamicModelFile << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.var_endo_nbr << ", mxREAL);" << endl;
-
-  mDynamicModelFile << "     /* Create a C pointer to a copy of the output matrix g1. */" << endl
+                    << "     plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl
+                    << "     /* Create a C pointer to a copy of the output matrix g1. */" << endl
                     << "     g1 = mxGetPr(plhs[1]);" << endl
                     << "  }" << endl
                     << endl
@@ -1709,7 +1705,7 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
                     << " if (nlhs >= 3)" << endl
                     << "  {" << endl
                     << "     /* Set the output pointer to the output matrix g2. */" << endl
-                    << "     plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.get_dyn_var_nbr()*variable_table.get_dyn_var_nbr() << ", mxREAL);" << endl
+                    << "     plhs[2] = mxCreateDoubleMatrix(" << equations.size() << ", " << variable_table.getDynJacobianColsNbr(computeJacobianExo)*variable_table.getDynJacobianColsNbr(computeJacobianExo) << ", mxREAL);" << endl
                     << "     /* Create a C pointer to a copy of the output matrix g1. */" << endl
                     << "     g2 = mxGetPr(plhs[2]);" << endl
                     << "  }" << endl
@@ -3295,11 +3291,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
   writeModelEquations(model_output, output_type);
 
   int nrows = equations.size();
-  int nvars;
-  if (computeJacobianExo)
-    nvars = variable_table.get_dyn_var_nbr();
-  else
-    nvars = variable_table.var_endo_nbr;
+  int nvars = variable_table.getDynJacobianColsNbr(computeJacobianExo);
   int nvars_sq = nvars * nvars;
 
   // Writing Jacobian
@@ -3315,7 +3307,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
           {
             ostringstream g1;
             g1 << "  g1";
-            matrixHelper(g1, eq, variable_table.getSortID(var), output_type);
+            matrixHelper(g1, eq, variable_table.getDynJacobianCol(var), output_type);
 
             jacobian_output << g1.str() << "=" << g1.str() << "+";
             d1->writeOutput(jacobian_output, output_type, temporary_terms);
@@ -3333,8 +3325,8 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
         int var2 = it->first.second.second;
         NodeID d2 = it->second;
 
-        int id1 = variable_table.getSortID(var1);
-        int id2 = variable_table.getSortID(var2);
+        int id1 = variable_table.getDynJacobianCol(var1);
+        int id2 = variable_table.getDynJacobianCol(var2);
 
         int col_nb = id1*nvars+id2;
         int col_nb_sym = id2*nvars+id1;
@@ -3367,9 +3359,9 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
         int var3 = it->first.second.second.second;
         NodeID d3 = it->second;
 
-        int id1 = variable_table.getSortID(var1);
-        int id2 = variable_table.getSortID(var2);
-        int id3 = variable_table.getSortID(var3);
+        int id1 = variable_table.getDynJacobianCol(var1);
+        int id2 = variable_table.getDynJacobianCol(var2);
+        int id3 = variable_table.getDynJacobianCol(var3);
 
         // Reference column number for the g3 matrix
         int ref_col = id1 * nvars_sq + id2 * nvars + id3;
@@ -3507,7 +3499,7 @@ ModelTree::writeOutput(ostream &output) const
           try
             {
               int varID = variable_table.getID(eEndogenous, endoID, lag);
-              output << " " << variable_table.getSortID(varID) + 1;
+              output << " " << variable_table.getDynJacobianCol(varID) + 1;
             }
           catch(VariableTable::UnknownVariableKeyException &e)
             {
@@ -3689,8 +3681,8 @@ ModelTree::computingPass(const eval_context_type &eval_context)
 {
   cout << equations.size() << " equation(s) found" << endl;
 
-  // Sorting variable table
-  variable_table.sort();
+  // Computes dynamic jacobian columns
+  variable_table.computeDynJacobianCols();
 
   // Determine derivation order
   int order = 1;

VariableTable.cc

@@ -17,6 +17,8 @@
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <cstdlib>
+
 #include "VariableTable.hh"
 
 VariableTable::VariableTable(const SymbolTable &symbol_table_arg) :
@@ -31,10 +33,10 @@ VariableTable::VariableTable(const SymbolTable &symbol_table_arg) :
 }
 
 int
-VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySortedException)
+VariableTable::addVariable(Type type, int symb_id, int lag) throw (DynJacobianColsAlreadyComputedException)
 {
-  if (sorted_ids_table.size() != 0)
-    throw AlreadySortedException();
+  if (dyn_jacobian_cols_table.size() != 0)
+    throw DynJacobianColsAlreadyComputedException();
 
   var_key_type key = make_pair(make_pair(type, lag), symb_id);
 
@@ -48,13 +50,6 @@ VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySorted
   variable_table[key] = var_id;
   inv_variable_table[var_id] = key;
 
-  if (type == eEndogenous)
-    var_endo_nbr++;
-  if (type == eExogenous)
-    var_exo_nbr++;
-  if (type == eExogenousDet)
-    var_exo_det_nbr++;
-
   // Setting maximum and minimum lags
   if (max_lead < lag)
     max_lead = lag;
@@ -64,18 +59,21 @@ VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySorted
   switch(type)
     {
     case eEndogenous:
+      var_endo_nbr++;
       if (max_endo_lead < lag)
         max_endo_lead = lag;
       else if (-max_endo_lag > lag)
         max_endo_lag = -lag;
       break;
     case eExogenous:
+      var_exo_nbr++;
       if (max_exo_lead < lag)
         max_exo_lead = lag;
       else if (-max_exo_lag > lag)
         max_exo_lag = -lag;
       break;
     case eExogenousDet:
+      var_exo_det_nbr++;
       if (max_exo_det_lead < lag)
         max_exo_det_lead = lag;
       else if (-max_exo_det_lag > lag)
@@ -88,21 +86,39 @@ VariableTable::addVariable(Type type, int symb_id, int lag) throw (AlreadySorted
         max_recur_lag = -lag;
       break;
     default:
-      ;
+      cerr << "VariableTable::addVariable(): forbidden variable type" << endl;
+      exit(-1);
     }
   return var_id;
 }
 
 void
-VariableTable::sort() throw (AlreadySortedException)
+VariableTable::computeDynJacobianCols() throw (DynJacobianColsAlreadyComputedException)
 {
-  if (sorted_ids_table.size() != 0)
-    throw AlreadySortedException();
+  if (dyn_jacobian_cols_table.size() != 0)
+    throw DynJacobianColsAlreadyComputedException();
+
+  dyn_jacobian_cols_table.resize(size());
 
+  variable_table_type::const_iterator it = variable_table.begin();
+
+  // Assign the first columns to endogenous, using the lexicographic order over (lag, symbol_id) implemented in variable_table map
   int sorted_id = 0;
-  sorted_ids_table.resize(size());
+  while(it->first.first.first == eEndogenous)
+    {
+      dyn_jacobian_cols_table[it->second] = sorted_id++;
+      it++;
+    }
 
-  for(variable_table_type::const_iterator it = variable_table.begin();
-      it != variable_table.end(); it++)
-    sorted_ids_table[it->second] = sorted_id++;
+  // Assign subsequent columns to exogenous and then exogenous deterministic, using an offset + symbol_id
+  while(it->first.first.first == eExogenous)
+    {
+      dyn_jacobian_cols_table[it->second] = var_endo_nbr + it->first.second;
+      it++;
+    }
+  while(it->first.first.first == eExogenousDet)
+    {
+      dyn_jacobian_cols_table[it->second] = var_endo_nbr + symbol_table.exo_nbr + it->first.second;
+      it++;
+    }
 }

include/CodeInterpreter.hh

@@ -58,7 +58,7 @@ const int EVALUATE_FOREWARD_R=8;
 const int EVALUATE_BACKWARD_R=9;
 
 //! Enumeration of possible symbol types
-/*! Warning: do not to change existing values: the order matters for VariableTable (at least ensure that eEndogenous is the first one), and the values matter for homotopy_setup command */
+/*! Warning: do not to change existing values: the order matters for VariableTable (at least for endogenous and exogenous types), and the values matter for homotopy_setup command */
 enum Type
   {
     eEndogenous = 0,               //!< Endogenous

include/VariableTable.hh

@@ -35,7 +35,7 @@ private:
   //! A reference to the symbol table
   const SymbolTable &symbol_table;
   //! A variable is a tuple (type, lag, symbol_id)
-  /*! Warning: don't change the order of elements in the tuple, since this determines the lexicographic ordering in sort() */
+  /*! Warning: don't change the order of elements in the tuple, since this determines the lexicographic ordering in computeDynJacobianCols() */
   typedef pair<pair<Type, int>, int> var_key_type;
 
   typedef map<var_key_type, int> variable_table_type;
@@ -45,16 +45,18 @@ private:
   typedef map<int, var_key_type> inv_variable_table_type;
   //! Maps a variable ID to a tuple (type, lag, symbol_id)
   inv_variable_table_type inv_variable_table;
-  //! Contains the sorted IDs (indexed by variable IDs)
-  vector<int> sorted_ids_table;
-public:
-  VariableTable(const SymbolTable &symbol_table_arg);
+
   //! Number of dynamic endogenous variables inside the model block
   int var_endo_nbr;
   //! Number of dynamic exogenous variables inside the model block
   int var_exo_nbr;
   //! Number of dynamic deterministic exogenous variables inside the model block
   int var_exo_det_nbr;
+
+  //! Contains the columns indices for the dynamic jacobian (indexed by variable IDs)
+  vector<int> dyn_jacobian_cols_table;
+public:
+  VariableTable(const SymbolTable &symbol_table_arg);
   //! Maximum lag over all types of variables (positive value)
   int max_lag;
   //! Maximum lead over all types of variables
@@ -91,17 +93,17 @@ public:
     int id;
     UnknownVariableIDException(int id_arg) : id(id_arg) {}
   };
-  //! Thrown when getSortID() called before sort()
-  class NotYetSortedException
+  //! Thrown when getDynJacobianCol() called before computeDynJacobianCols()
+  class DynJacobianColsNotYetComputedException
   {
   };
-  //! Thrown when sort() or addVariable() called after sort()
-  class AlreadySortedException
+  //! Thrown when computeDynJacobianCols() or addVariable() called after computeDynJacobianCols()
+  class DynJacobianColsAlreadyComputedException
   {
   };
   //! Adds a variable in the table, and returns its (newly allocated) variable ID
   /*! Also works if the variable already exists */
-  int addVariable(Type type, int symb_id, int lag) throw (AlreadySortedException);
+  int addVariable(Type type, int symb_id, int lag) throw (DynJacobianColsAlreadyComputedException);
   //! Return variable ID
   inline int getID(Type type, int symb_id, int lag) const throw (UnknownVariableKeyException);
   //! Return lag of variable
@@ -112,26 +114,24 @@ public:
   inline Type getType(int var_id) const throw (UnknownVariableIDException);
   //! Get number of variables
   inline int size() const;
-  //! Get variable ID of sorted variable table
-  /*! In practice, only used for endogenous variables */
-  inline int getSortID(int var_id) const throw (NotYetSortedException, UnknownVariableIDException);
-  //! Sorts variable table
-  /*! The order used is a lexicographic order over the tuple (type, lag, symbol_id) */
-  void sort() throw (AlreadySortedException);
-  //! Get the number of dynamic variables 
-  inline int get_dyn_var_nbr() const;
+  //! Get column index in dynamic jacobian
+  inline int getDynJacobianCol(int var_id) const throw (DynJacobianColsNotYetComputedException, UnknownVariableIDException);
+  //! Computes column indices in dynamic jacobian
+  void computeDynJacobianCols() throw (DynJacobianColsAlreadyComputedException);
+  //! Get the number of columns of dynamic jacobian (depending on whether we compute derivatives w.r. to exogenous)
+  inline int getDynJacobianColsNbr(bool computeJacobianExo) const;
 };
 
 inline int
-VariableTable::getSortID(int var_id) const throw (NotYetSortedException, UnknownVariableIDException)
+VariableTable::getDynJacobianCol(int var_id) const throw (DynJacobianColsNotYetComputedException, UnknownVariableIDException)
 {
-  if (sorted_ids_table.size() == 0)
-    throw NotYetSortedException();
+  if (dyn_jacobian_cols_table.size() == 0)
+    throw DynJacobianColsNotYetComputedException();
 
   if (var_id < 0 || var_id >= size())
     throw UnknownVariableIDException(var_id);
 
-  return sorted_ids_table[var_id];
+  return dyn_jacobian_cols_table[var_id];
 }
 
 inline int
@@ -181,9 +181,12 @@ VariableTable::size() const
 }
 
 inline int 
-VariableTable::get_dyn_var_nbr() const
+VariableTable::getDynJacobianColsNbr(bool computeJacobianExo) const
 {
-  return var_endo_nbr + symbol_table.exo_nbr + symbol_table.exo_det_nbr;
+  if (computeJacobianExo)
+    return var_endo_nbr + symbol_table.exo_nbr + symbol_table.exo_det_nbr;
+  else
+    return var_endo_nbr;
 }
 
 #endif