diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index 200195c3ccf4d616e4830a17041e05a1d8119a23..d84153427eedb7c0c19ca6d6c303c5689a600994 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -647,7 +647,7 @@ CalibVarStatement::writeOutput(ostream &output, const string &basename) const
     {
       const string &name = it->first;
       const string &weight = it->second.first;
-      const NodeID expression = it->second.second;
+      const expr_t expression = it->second.second;
 
       int id = symbol_table.getTypeSpecificID(name) + 1;
       if (symbol_table.getType(name) == eEndogenous)
@@ -675,7 +675,7 @@ CalibVarStatement::writeOutput(ostream &output, const string &basename) const
       const string &name1 = it->first.first;
       const string &name2 = it->first.second;
       const string &weight = it->second.first;
-      const NodeID expression = it->second.second;
+      const expr_t expression = it->second.second;
 
       int id1 = symbol_table.getTypeSpecificID(name1) + 1;
       int id2 = symbol_table.getTypeSpecificID(name2) + 1;
@@ -706,7 +706,7 @@ CalibVarStatement::writeOutput(ostream &output, const string &basename) const
       const string &name = it->first.first;
       int iar = it->first.second + 3;
       const string &weight = it->second.first;
-      const NodeID expression = it->second.second;
+      const expr_t expression = it->second.second;
 
       int id = symbol_table.getTypeSpecificID(name) + 1;
 
@@ -810,7 +810,7 @@ OptimWeightsStatement::writeOutput(ostream &output, const string &basename) cons
        it != var_weights.end(); it++)
     {
       const string &name = it->first;
-      const NodeID value = it->second;
+      const expr_t value = it->second;
       int id = symbol_table.getTypeSpecificID(name) + 1;
       output <<  "optim_weights_(" << id << "," << id << ") = ";
       value->writeOutput(output);
@@ -823,7 +823,7 @@ OptimWeightsStatement::writeOutput(ostream &output, const string &basename) cons
     {
       const string &name1 = it->first.first;
       const string &name2 = it->first.second;
-      const NodeID value = it->second;
+      const expr_t value = it->second;
       int id1 = symbol_table.getTypeSpecificID(name1) + 1;
       int id2 = symbol_table.getTypeSpecificID(name2) + 1;
       output <<  "optim_weights_(" << id1 << "," << id2 << ") = ";
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index 79c24dd67d60d0469ede716b8eb6b15d4bbcbc35..cbde458ecb8232b8434944802fb4eb98edfb9078 100644
--- a/ComputingTasks.hh
+++ b/ComputingTasks.hh
@@ -169,7 +169,7 @@ public:
 class ObservationTrendsStatement : public Statement
 {
 public:
-  typedef map<string, NodeID> trend_elements_t;
+  typedef map<string, expr_t> trend_elements_t;
 private:
   const trend_elements_t trend_elements;
   const SymbolTable &symbol_table;
@@ -241,7 +241,7 @@ class EstimationParams
 public:
   int type;
   string name, name2, prior;
-  NodeID init_val, low_bound, up_bound, mean, std, p3, p4, jscale;
+  expr_t init_val, low_bound, up_bound, mean, std, p3, p4, jscale;
 
   void
   init(const DataTree &datatree)
@@ -298,8 +298,8 @@ public:
 class OptimWeightsStatement : public Statement
 {
 public:
-  typedef map<string, NodeID> var_weights_t;
-  typedef map<pair<string, string>, NodeID> covar_weights_t;
+  typedef map<string, expr_t> var_weights_t;
+  typedef map<pair<string, string>, expr_t> covar_weights_t;
 private:
   const var_weights_t var_weights;
   const covar_weights_t covar_weights;
@@ -324,11 +324,11 @@ class CalibVarStatement : public Statement
 {
 public:
   //! Maps a variable to a pair (weight, expression)
-  typedef map<string, pair<string, NodeID> > calib_var_t;
+  typedef map<string, pair<string, expr_t> > calib_var_t;
   //! Maps a pair of variables to a pair (weight, expression)
-  typedef map<pair<string, string>, pair<string, NodeID> > calib_covar_t;
+  typedef map<pair<string, string>, pair<string, expr_t> > calib_covar_t;
   //! Maps a pair (variable, autocorr) to a pair (weight, expression)
-  typedef map<pair<string, int>, pair<string, NodeID> > calib_ac_t;
+  typedef map<pair<string, int>, pair<string, expr_t> > calib_ac_t;
 private:
   const calib_var_t calib_var;
   const calib_covar_t calib_covar;
diff --git a/DataTree.cc b/DataTree.cc
index 21eb50cc8a57ba973f2704577ccc2a21d9a17baf..df7044dc8a1676d70cf40597841cc280b625c236 100644
--- a/DataTree.cc
+++ b/DataTree.cc
@@ -50,7 +50,7 @@ DataTree::~DataTree()
     delete *it;
 }
 
-NodeID
+expr_t
 DataTree::AddNumConstant(const string &value)
 {
   int id = num_constants.AddConstant(value);
@@ -79,8 +79,8 @@ DataTree::AddVariable(int symb_id, int lag)
   return AddVariableInternal(symb_id, lag);
 }
 
-NodeID
-DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
 {
   if (iArg1 != Zero && iArg2 != Zero)
     {
@@ -93,7 +93,7 @@ DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
       // Nodes iArg1 and iArg2 are sorted by index
       if (iArg1->idx > iArg2->idx)
         {
-          NodeID tmp = iArg1;
+          expr_t tmp = iArg1;
           iArg1 = iArg2;
           iArg2 = tmp;
         }
@@ -107,8 +107,8 @@ DataTree::AddPlus(NodeID iArg1, NodeID iArg2)
     return Zero;
 }
 
-NodeID
-DataTree::AddMinus(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddMinus(expr_t iArg1, expr_t iArg2)
 {
   if (iArg2 == Zero)
     return iArg1;
@@ -122,8 +122,8 @@ DataTree::AddMinus(NodeID iArg1, NodeID iArg2)
   return AddBinaryOp(iArg1, oMinus, iArg2);
 }
 
-NodeID
-DataTree::AddUMinus(NodeID iArg1)
+expr_t
+DataTree::AddUMinus(expr_t iArg1)
 {
   if (iArg1 != Zero)
     {
@@ -138,8 +138,8 @@ DataTree::AddUMinus(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddTimes(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddTimes(expr_t iArg1, expr_t iArg2)
 {
   if (iArg1 == MinusOne)
     return AddUMinus(iArg2);
@@ -151,7 +151,7 @@ DataTree::AddTimes(NodeID iArg1, NodeID iArg2)
       // Nodes iArg1 and iArg2 are sorted by index
       if (iArg1->idx > iArg2->idx)
         {
-          NodeID tmp = iArg1;
+          expr_t tmp = iArg1;
           iArg1 = iArg2;
           iArg2 = tmp;
         }
@@ -167,8 +167,8 @@ DataTree::AddTimes(NodeID iArg1, NodeID iArg2)
     return Zero;
 }
 
-NodeID
-DataTree::AddDivide(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddDivide(expr_t iArg1, expr_t iArg2)
 {
   if (iArg2 == One)
     return iArg1;
@@ -189,44 +189,44 @@ DataTree::AddDivide(NodeID iArg1, NodeID iArg2)
   return AddBinaryOp(iArg1, oDivide, iArg2);
 }
 
-NodeID
-DataTree::AddLess(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddLess(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oLess, iArg2);
 }
 
-NodeID
-DataTree::AddGreater(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddGreater(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oGreater, iArg2);
 }
 
-NodeID
-DataTree::AddLessEqual(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddLessEqual(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oLessEqual, iArg2);
 }
 
-NodeID
-DataTree::AddGreaterEqual(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddGreaterEqual(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oGreaterEqual, iArg2);
 }
 
-NodeID
-DataTree::AddEqualEqual(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddEqualEqual(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oEqualEqual, iArg2);
 }
 
-NodeID
-DataTree::AddDifferent(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddDifferent(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oDifferent, iArg2);
 }
 
-NodeID
-DataTree::AddPower(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddPower(expr_t iArg1, expr_t iArg2)
 {
   if (iArg1 != Zero && iArg2 != Zero && iArg2 != One)
     return AddBinaryOp(iArg1, oPower, iArg2);
@@ -238,8 +238,8 @@ DataTree::AddPower(NodeID iArg1, NodeID iArg2)
     return Zero;
 }
 
-NodeID
-DataTree::AddExp(NodeID iArg1)
+expr_t
+DataTree::AddExp(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oExp, iArg1);
@@ -247,8 +247,8 @@ DataTree::AddExp(NodeID iArg1)
     return One;
 }
 
-NodeID
-DataTree::AddLog(NodeID iArg1)
+expr_t
+DataTree::AddLog(expr_t iArg1)
 {
   if (iArg1 != Zero && iArg1 != One)
     return AddUnaryOp(oLog, iArg1);
@@ -261,8 +261,8 @@ DataTree::AddLog(NodeID iArg1)
     }
 }
 
-NodeID
-DataTree::AddLog10(NodeID iArg1)
+expr_t
+DataTree::AddLog10(expr_t iArg1)
 {
   if (iArg1 != Zero && iArg1 != One)
     return AddUnaryOp(oLog10, iArg1);
@@ -275,8 +275,8 @@ DataTree::AddLog10(NodeID iArg1)
     }
 }
 
-NodeID
-DataTree::AddCos(NodeID iArg1)
+expr_t
+DataTree::AddCos(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oCos, iArg1);
@@ -284,8 +284,8 @@ DataTree::AddCos(NodeID iArg1)
     return One;
 }
 
-NodeID
-DataTree::AddSin(NodeID iArg1)
+expr_t
+DataTree::AddSin(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oSin, iArg1);
@@ -293,8 +293,8 @@ DataTree::AddSin(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddTan(NodeID iArg1)
+expr_t
+DataTree::AddTan(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oTan, iArg1);
@@ -302,8 +302,8 @@ DataTree::AddTan(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddAcos(NodeID iArg1)
+expr_t
+DataTree::AddAcos(expr_t iArg1)
 {
   if (iArg1 != One)
     return AddUnaryOp(oAcos, iArg1);
@@ -311,8 +311,8 @@ DataTree::AddAcos(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddAsin(NodeID iArg1)
+expr_t
+DataTree::AddAsin(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oAsin, iArg1);
@@ -320,8 +320,8 @@ DataTree::AddAsin(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddAtan(NodeID iArg1)
+expr_t
+DataTree::AddAtan(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oAtan, iArg1);
@@ -329,8 +329,8 @@ DataTree::AddAtan(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddCosh(NodeID iArg1)
+expr_t
+DataTree::AddCosh(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oCosh, iArg1);
@@ -338,8 +338,8 @@ DataTree::AddCosh(NodeID iArg1)
     return One;
 }
 
-NodeID
-DataTree::AddSinh(NodeID iArg1)
+expr_t
+DataTree::AddSinh(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oSinh, iArg1);
@@ -347,8 +347,8 @@ DataTree::AddSinh(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddTanh(NodeID iArg1)
+expr_t
+DataTree::AddTanh(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oTanh, iArg1);
@@ -356,8 +356,8 @@ DataTree::AddTanh(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddAcosh(NodeID iArg1)
+expr_t
+DataTree::AddAcosh(expr_t iArg1)
 {
   if (iArg1 != One)
     return AddUnaryOp(oAcosh, iArg1);
@@ -365,8 +365,8 @@ DataTree::AddAcosh(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddAsinh(NodeID iArg1)
+expr_t
+DataTree::AddAsinh(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oAsinh, iArg1);
@@ -374,8 +374,8 @@ DataTree::AddAsinh(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddAtanh(NodeID iArg1)
+expr_t
+DataTree::AddAtanh(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oAtanh, iArg1);
@@ -383,8 +383,8 @@ DataTree::AddAtanh(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddSqrt(NodeID iArg1)
+expr_t
+DataTree::AddSqrt(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oSqrt, iArg1);
@@ -392,8 +392,8 @@ DataTree::AddSqrt(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddErf(NodeID iArg1)
+expr_t
+DataTree::AddErf(expr_t iArg1)
 {
   if (iArg1 != Zero)
     return AddUnaryOp(oErf, iArg1);
@@ -401,69 +401,69 @@ DataTree::AddErf(NodeID iArg1)
     return Zero;
 }
 
-NodeID
-DataTree::AddMax(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddMax(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oMax, iArg2);
 }
 
-NodeID
-DataTree::AddMin(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddMin(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oMin, iArg2);
 }
 
-NodeID
-DataTree::AddNormcdf(NodeID iArg1, NodeID iArg2, NodeID iArg3)
+expr_t
+DataTree::AddNormcdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
 {
   return AddTrinaryOp(iArg1, oNormcdf, iArg2, iArg3);
 }
 
-NodeID
-DataTree::AddNormpdf(NodeID iArg1, NodeID iArg2, NodeID iArg3)
+expr_t
+DataTree::AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3)
 {
   return AddTrinaryOp(iArg1, oNormpdf, iArg2, iArg3);
 }
 
-NodeID
-DataTree::AddSteadyState(NodeID iArg1)
+expr_t
+DataTree::AddSteadyState(expr_t iArg1)
 {
   return AddUnaryOp(oSteadyState, iArg1);
 }
 
-NodeID
-DataTree::AddExpectation(int iArg1, NodeID iArg2)
+expr_t
+DataTree::AddExpectation(int iArg1, expr_t iArg2)
 {
   return AddUnaryOp(oExpectation, iArg2, iArg1);
 }
 
-NodeID
-DataTree::AddExpectation(string *iArg1, NodeID iArg2)
+expr_t
+DataTree::AddExpectation(string *iArg1, expr_t iArg2)
 {
   return AddUnaryOp(oExpectation, iArg2, 0, *iArg1);
 }
 
-NodeID
-DataTree::AddEqual(NodeID iArg1, NodeID iArg2)
+expr_t
+DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
 {
   return AddBinaryOp(iArg1, oEqual, iArg2);
 }
 
 void
-DataTree::AddLocalVariable(int symb_id, NodeID value) throw (LocalVariableException)
+DataTree::AddLocalVariable(int symb_id, expr_t value) throw (LocalVariableException)
 {
   assert(symbol_table.getType(symb_id) == eModelLocalVariable);
 
   // Throw an exception if symbol already declared
-  map<int, NodeID>::iterator it = local_variables_table.find(symb_id);
+  map<int, expr_t>::iterator it = local_variables_table.find(symb_id);
   if (it != local_variables_table.end())
     throw LocalVariableException(symbol_table.getName(symb_id));
 
   local_variables_table[symb_id] = value;
 }
 
-NodeID
-DataTree::AddExternalFunction(int symb_id, const vector<NodeID> &arguments)
+expr_t
+DataTree::AddExternalFunction(int symb_id, const vector<expr_t> &arguments)
 {
   assert(symbol_table.getType(symb_id) == eExternalFunction);
 
@@ -474,8 +474,8 @@ DataTree::AddExternalFunction(int symb_id, const vector<NodeID> &arguments)
   return new ExternalFunctionNode(*this, symb_id, arguments);
 }
 
-NodeID
-DataTree::AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index)
+expr_t
+DataTree::AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<expr_t> &arguments, int input_index)
 {
   assert(symbol_table.getType(top_level_symb_id) == eExternalFunction);
 
@@ -488,8 +488,8 @@ DataTree::AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<
   return new FirstDerivExternalFunctionNode(*this, top_level_symb_id, arguments, input_index);
 }
 
-NodeID
-DataTree::AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index1, int input_index2)
+expr_t
+DataTree::AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<expr_t> &arguments, int input_index1, int input_index2)
 {
   assert(symbol_table.getType(top_level_symb_id) == eExternalFunction);
 
diff --git a/DataTree.hh b/DataTree.hh
index ee5c029b4b7d0695b4a052e9b4eaa596370ff46b..0bd3a813ab36bec9993629a27c7738c18f9a6aa8 100644
--- a/DataTree.hh
+++ b/DataTree.hh
@@ -59,43 +59,43 @@ protected:
   //! Pair (symbol_id, lag) used as key
   typedef map<pair<int, int>, VariableNode *> variable_node_map_t;
   variable_node_map_t variable_node_map;
-  typedef map<pair<NodeID, UnaryOpcode>, UnaryOpNode *> unary_op_node_map_t;
+  typedef map<pair<expr_t, UnaryOpcode>, UnaryOpNode *> unary_op_node_map_t;
   unary_op_node_map_t unary_op_node_map;
-  typedef map<pair<pair<NodeID, NodeID>, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t;
+  typedef map<pair<pair<expr_t, expr_t>, BinaryOpcode>, BinaryOpNode *> binary_op_node_map_t;
   binary_op_node_map_t binary_op_node_map;
-  typedef map<pair<pair<pair<NodeID, NodeID>, NodeID>, TrinaryOpcode>, TrinaryOpNode *> trinary_op_node_map_t;
+  typedef map<pair<pair<pair<expr_t, expr_t>, expr_t>, TrinaryOpcode>, TrinaryOpNode *> trinary_op_node_map_t;
   trinary_op_node_map_t trinary_op_node_map;
-  typedef map<pair<vector<NodeID>, int>, ExternalFunctionNode *> external_function_node_map_t;
+  typedef map<pair<vector<expr_t>, int>, ExternalFunctionNode *> external_function_node_map_t;
   external_function_node_map_t external_function_node_map;
-  typedef map<pair<pair<vector<NodeID>, int>, int>, FirstDerivExternalFunctionNode *> first_deriv_external_function_node_map_t;
+  typedef map<pair<pair<vector<expr_t>, int>, int>, FirstDerivExternalFunctionNode *> first_deriv_external_function_node_map_t;
   first_deriv_external_function_node_map_t first_deriv_external_function_node_map;
-  typedef map<pair<pair<vector<NodeID>, pair<int, int> >, int>, SecondDerivExternalFunctionNode *> second_deriv_external_function_node_map_t;
+  typedef map<pair<pair<vector<expr_t>, pair<int, int> >, int>, SecondDerivExternalFunctionNode *> second_deriv_external_function_node_map_t;
   second_deriv_external_function_node_map_t second_deriv_external_function_node_map;
 
   //! Stores local variables value (maps symbol ID to corresponding node)
-  map<int, NodeID> local_variables_table;
+  map<int, expr_t> local_variables_table;
 
   //! Internal implementation of AddVariable(), without the check on the lag
   VariableNode *AddVariableInternal(int symb_id, int lag);
 
 private:
-  typedef list<NodeID> node_list_t;
+  typedef list<expr_t> node_list_t;
   //! The list of nodes
   node_list_t node_list;
   //! A counter for filling ExprNode's idx field
   int node_counter;
 
-  inline NodeID AddPossiblyNegativeConstant(double val);
-  inline NodeID AddUnaryOp(UnaryOpcode op_code, NodeID arg, int arg_exp_info_set = 0, const string &arg_exp_info_set_name="");
-  inline NodeID AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2);
-  inline NodeID AddTrinaryOp(NodeID arg1, TrinaryOpcode op_code, NodeID arg2, NodeID arg3);
+  inline expr_t AddPossiblyNegativeConstant(double val);
+  inline expr_t AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set = 0, const string &arg_exp_info_set_name="");
+  inline expr_t AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2);
+  inline expr_t AddTrinaryOp(expr_t arg1, TrinaryOpcode op_code, expr_t arg2, expr_t arg3);
 
 public:
   DataTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_arg);
   virtual ~DataTree();
 
   //! Some predefined constants
-  NodeID Zero, One, Two, MinusOne, NaN, Infinity, MinusInfinity, Pi;
+  expr_t Zero, One, Two, MinusOne, NaN, Infinity, MinusInfinity, Pi;
 
   //! Raised when a local parameter is declared twice
   class LocalVariableException
@@ -108,92 +108,92 @@ public:
   };
 
   //! Adds a numerical constant
-  NodeID AddNumConstant(const string &value);
+  expr_t AddNumConstant(const string &value);
   //! Adds a variable
   /*! The default implementation of the method refuses any lag != 0 */
   virtual VariableNode *AddVariable(int symb_id, int lag = 0);
   //! Adds "arg1+arg2" to model tree
-  NodeID AddPlus(NodeID iArg1, NodeID iArg2);
+  expr_t AddPlus(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1-arg2" to model tree
-  NodeID AddMinus(NodeID iArg1, NodeID iArg2);
+  expr_t AddMinus(expr_t iArg1, expr_t iArg2);
   //! Adds "-arg" to model tree
-  NodeID AddUMinus(NodeID iArg1);
+  expr_t AddUMinus(expr_t iArg1);
   //! Adds "arg1*arg2" to model tree
-  NodeID AddTimes(NodeID iArg1, NodeID iArg2);
+  expr_t AddTimes(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1/arg2" to model tree
-  NodeID AddDivide(NodeID iArg1, NodeID iArg2);
+  expr_t AddDivide(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1<arg2" to model tree
-  NodeID AddLess(NodeID iArg1, NodeID iArg2);
+  expr_t AddLess(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1>arg2" to model tree
-  NodeID AddGreater(NodeID iArg1, NodeID iArg2);
+  expr_t AddGreater(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1<=arg2" to model tree
-  NodeID AddLessEqual(NodeID iArg1, NodeID iArg2);
+  expr_t AddLessEqual(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1>=arg2" to model tree
-  NodeID AddGreaterEqual(NodeID iArg1, NodeID iArg2);
+  expr_t AddGreaterEqual(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1==arg2" to model tree
-  NodeID AddEqualEqual(NodeID iArg1, NodeID iArg2);
+  expr_t AddEqualEqual(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1!=arg2" to model tree
-  NodeID AddDifferent(NodeID iArg1, NodeID iArg2);
+  expr_t AddDifferent(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1^arg2" to model tree
-  NodeID AddPower(NodeID iArg1, NodeID iArg2);
+  expr_t AddPower(expr_t iArg1, expr_t iArg2);
   //! Adds "E(arg1)(arg2)" to model tree
-  NodeID AddExpectation(int iArg1, NodeID iArg2);
+  expr_t AddExpectation(int iArg1, expr_t iArg2);
   //! Adds "E(arg1)(arg2)" to model tree
-  NodeID AddExpectation(string *iArg1, NodeID iArg2);
+  expr_t AddExpectation(string *iArg1, expr_t iArg2);
   //! Adds "exp(arg)" to model tree
-  NodeID AddExp(NodeID iArg1);
+  expr_t AddExp(expr_t iArg1);
   //! Adds "log(arg)" to model tree
-  NodeID AddLog(NodeID iArg1);
+  expr_t AddLog(expr_t iArg1);
   //! Adds "log10(arg)" to model tree
-  NodeID AddLog10(NodeID iArg1);
+  expr_t AddLog10(expr_t iArg1);
   //! Adds "cos(arg)" to model tree
-  NodeID AddCos(NodeID iArg1);
+  expr_t AddCos(expr_t iArg1);
   //! Adds "sin(arg)" to model tree
-  NodeID AddSin(NodeID iArg1);
+  expr_t AddSin(expr_t iArg1);
   //! Adds "tan(arg)" to model tree
-  NodeID AddTan(NodeID iArg1);
+  expr_t AddTan(expr_t iArg1);
   //! Adds "acos(arg)" to model tree
-  NodeID AddAcos(NodeID iArg1);
+  expr_t AddAcos(expr_t iArg1);
   //! Adds "asin(arg)" to model tree
-  NodeID AddAsin(NodeID iArg1);
+  expr_t AddAsin(expr_t iArg1);
   //! Adds "atan(arg)" to model tree
-  NodeID AddAtan(NodeID iArg1);
+  expr_t AddAtan(expr_t iArg1);
   //! Adds "cosh(arg)" to model tree
-  NodeID AddCosh(NodeID iArg1);
+  expr_t AddCosh(expr_t iArg1);
   //! Adds "sinh(arg)" to model tree
-  NodeID AddSinh(NodeID iArg1);
+  expr_t AddSinh(expr_t iArg1);
   //! Adds "tanh(arg)" to model tree
-  NodeID AddTanh(NodeID iArg1);
+  expr_t AddTanh(expr_t iArg1);
   //! Adds "acosh(arg)" to model tree
-  NodeID AddAcosh(NodeID iArg1);
+  expr_t AddAcosh(expr_t iArg1);
   //! Adds "asinh(arg)" to model tree
-  NodeID AddAsinh(NodeID iArg1);
+  expr_t AddAsinh(expr_t iArg1);
   //! Adds "atanh(args)" to model tree
-  NodeID AddAtanh(NodeID iArg1);
+  expr_t AddAtanh(expr_t iArg1);
   //! Adds "sqrt(arg)" to model tree
-  NodeID AddSqrt(NodeID iArg1);
+  expr_t AddSqrt(expr_t iArg1);
   //! Adds "erf(arg)" to model tree
-  NodeID AddErf(NodeID iArg1);
+  expr_t AddErf(expr_t iArg1);
   //! Adds "max(arg1,arg2)" to model tree
-  NodeID AddMax(NodeID iArg1, NodeID iArg2);
+  expr_t AddMax(expr_t iArg1, expr_t iArg2);
   //! Adds "min(arg1,arg2)" to model tree
-  NodeID AddMin(NodeID iArg1, NodeID iArg2);
+  expr_t AddMin(expr_t iArg1, expr_t iArg2);
   //! Adds "normcdf(arg1,arg2,arg3)" to model tree
-  NodeID AddNormcdf(NodeID iArg1, NodeID iArg2, NodeID iArg3);
+  expr_t AddNormcdf(expr_t iArg1, expr_t iArg2, expr_t iArg3);
   //! Adds "normpdf(arg1,arg2,arg3)" to model tree
-  NodeID AddNormpdf(NodeID iArg1, NodeID iArg2, NodeID iArg3);
+  expr_t AddNormpdf(expr_t iArg1, expr_t iArg2, expr_t iArg3);
   //! Adds "steadyState(arg)" to model tree
-  NodeID AddSteadyState(NodeID iArg1);
+  expr_t AddSteadyState(expr_t iArg1);
   //! Adds "arg1=arg2" to model tree
-  NodeID AddEqual(NodeID iArg1, NodeID iArg2);
+  expr_t AddEqual(expr_t iArg1, expr_t iArg2);
   //! Adds a model local variable with its value
-  void AddLocalVariable(int symb_id, NodeID value) throw (LocalVariableException);
+  void AddLocalVariable(int symb_id, expr_t value) throw (LocalVariableException);
   //! Adds an external function node
-  NodeID AddExternalFunction(int symb_id, const vector<NodeID> &arguments);
+  expr_t AddExternalFunction(int symb_id, const vector<expr_t> &arguments);
   //! Adds an external function node for the first derivative of an external function
-  NodeID AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index);
+  expr_t AddFirstDerivExternalFunctionNode(int top_level_symb_id, const vector<expr_t> &arguments, int input_index);
   //! Adds an external function node for the second derivative of an external function
-  NodeID AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<NodeID> &arguments, int input_index1, int input_index2);
+  expr_t AddSecondDerivExternalFunctionNode(int top_level_symb_id, const vector<expr_t> &arguments, int input_index1, int input_index2);
   //! Checks if a given symbol is used somewhere in the data tree
   bool isSymbolUsed(int symb_id) const;
   //! Checks if a given unary op is used somewhere in the data tree
@@ -226,7 +226,7 @@ public:
   };
 };
 
-inline NodeID
+inline expr_t
 DataTree::AddPossiblyNegativeConstant(double v)
 {
   bool neg = false;
@@ -238,7 +238,7 @@ DataTree::AddPossiblyNegativeConstant(double v)
   ostringstream ost;
   ost << setprecision(CONSTANTS_PRECISION) << v;
 
-  NodeID cnode = AddNumConstant(ost.str());
+  expr_t cnode = AddNumConstant(ost.str());
 
   if (neg)
     return AddUMinus(cnode);
@@ -246,8 +246,8 @@ DataTree::AddPossiblyNegativeConstant(double v)
     return cnode;
 }
 
-inline NodeID
-DataTree::AddUnaryOp(UnaryOpcode op_code, NodeID arg, int arg_exp_info_set, const string &arg_exp_info_set_name)
+inline expr_t
+DataTree::AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set, const string &arg_exp_info_set_name)
 {
   // If the node already exists in tree, share it
   unary_op_node_map_t::iterator it = unary_op_node_map.find(make_pair(arg, op_code));
@@ -272,8 +272,8 @@ DataTree::AddUnaryOp(UnaryOpcode op_code, NodeID arg, int arg_exp_info_set, cons
   return new UnaryOpNode(*this, op_code, arg, arg_exp_info_set, arg_exp_info_set_name);
 }
 
-inline NodeID
-DataTree::AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2)
+inline expr_t
+DataTree::AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2)
 {
   binary_op_node_map_t::iterator it = binary_op_node_map.find(make_pair(make_pair(arg1, arg2), op_code));
   if (it != binary_op_node_map.end())
@@ -293,8 +293,8 @@ DataTree::AddBinaryOp(NodeID arg1, BinaryOpcode op_code, NodeID arg2)
   return new BinaryOpNode(*this, arg1, op_code, arg2);
 }
 
-inline NodeID
-DataTree::AddTrinaryOp(NodeID arg1, TrinaryOpcode op_code, NodeID arg2, NodeID arg3)
+inline expr_t
+DataTree::AddTrinaryOp(expr_t arg1, TrinaryOpcode op_code, expr_t arg2, expr_t arg3)
 {
   trinary_op_node_map_t::iterator it = trinary_op_node_map.find(make_pair(make_pair(make_pair(arg1, arg2), arg3), op_code));
   if (it != trinary_op_node_map.end())
diff --git a/DynamicModel.cc b/DynamicModel.cc
index 1bfc143544a959a4878635133efd1cd9dd371f0e..4817e0c3d7b564e0112c5bcf10c1f50aa4214ddc 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -72,7 +72,7 @@ DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int la
 void
 DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, const map_idx_t &map_idx) const
 {
-  map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
+  map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
   if (it != first_chain_rule_derivatives.end())
     (it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
   else
@@ -97,8 +97,8 @@ DynamicModel::initializeVariablesAndEquations()
 void
 DynamicModel::computeTemporaryTermsOrdered()
 {
-  map<NodeID, pair<int, int> > first_occurence;
-  map<NodeID, int> reference_count;
+  map<expr_t, pair<int, int> > first_occurence;
+  map<expr_t, int> reference_count;
   BinaryOpNode *eq_node;
   first_derivatives_t::const_iterator it;
   first_chain_rule_derivatives_t::const_iterator it_chr;
@@ -124,16 +124,16 @@ DynamicModel::computeTemporaryTermsOrdered()
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedNodeID(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
+                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
               else
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
                 }
             }
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
             {
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
             }
           for (derivative_t::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
@@ -158,16 +158,16 @@ DynamicModel::computeTemporaryTermsOrdered()
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedNodeID(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
+                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
               else
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
                 }
             }
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
             {
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
             }
           for (derivative_t::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
@@ -185,16 +185,16 @@ DynamicModel::computeTemporaryTermsOrdered()
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedNodeID(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+                getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
               else
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
                 }
             }
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
             {
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
             }
           for (derivative_t::const_iterator it = derivative_endo[block].begin(); it != derivative_endo[block].end(); it++)
@@ -223,10 +223,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
 {
   string tmp_s, sps;
   ostringstream tmp_output, tmp1_output, global_output;
-  NodeID lhs = NULL, rhs = NULL;
+  expr_t lhs = NULL, rhs = NULL;
   BinaryOpNode *eq_node;
   ostringstream Uf[symbol_table.endo_nbr()];
-  map<NodeID, int> reference_count;
+  map<expr_t, int> reference_count;
   temporary_terms_t local_temporary_terms;
   ofstream  output;
   int nze, nze_exo, nze_other_endo;
@@ -420,7 +420,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
           int equation_ID = getBlockEquationID(block, i);
           EquationType equ_type = getBlockEquationType(block, i);
           string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, variable_ID));
-          eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
           lhs = eq_node->get_arg1();
           rhs = eq_node->get_arg2();
           tmp_output.str("");
@@ -448,7 +448,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
                       rhs->writeOutput(output, local_output_type, local_temporary_terms);
                       output << "\n    ";
                       tmp_output.str("");
-                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedNodeID(block, i);
+                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                       lhs = eq_node->get_arg1();
                       rhs = eq_node->get_arg2();
                       lhs->writeOutput(output, local_output_type, local_temporary_terms);
@@ -517,7 +517,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               int eqr = getBlockInitialEquationID(block, eq);
               int varr = getBlockInitialVariableID(block, var);
 
-              NodeID id = it->second;
+              expr_t id = it->second;
 
               output << "      g1(" << eqr+1 << ", " << varr+1+(lag+block_max_lag)*block_size << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
@@ -532,7 +532,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               int eq = it->first.second.first;
               int var = it->first.second.second;
               int eqr = getBlockInitialEquationID(block, eq);
-              NodeID id = it->second;
+              expr_t id = it->second;
 
               output << "      g1_o(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
@@ -555,7 +555,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               int lag = it->first.first;
               unsigned int eq = it->first.second.first;
               unsigned int var = it->first.second.second;
-              NodeID id = it->second;
+              expr_t id = it->second;
 
               output << "    g1(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
@@ -570,7 +570,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               int lag = it->first.first;
               unsigned int eq = it->first.second.first;
               unsigned int var = it->first.second.second;
-              NodeID id = it->second;
+              expr_t id = it->second;
 
               output << "    g1_o(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
@@ -588,7 +588,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               unsigned int var = it->first.second;
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               int lag = it->second.first;
               output << "    g1(" << eq+1 << ", " << var+1-block_recursive << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
@@ -609,7 +609,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
               ostringstream tmp_output;
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               int lag = it->second.first;
               if (eq >= block_recursive && var >= block_recursive)
                 {
@@ -687,7 +687,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               int lag = it->first.first;
               unsigned int eq = it->first.second.first;
               unsigned int var = it->first.second.second;
-              NodeID id = it->second;
+              expr_t id = it->second;
               output << "      g1(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
               output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
@@ -700,7 +700,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               int lag = it->first.first;
               unsigned int eq = it->first.second.first;
               unsigned int var = it->first.second.second;
-              NodeID id = it->second;
+              expr_t id = it->second;
 
               output << "      g1_o(" << eq+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
@@ -784,7 +784,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
       int deriv_id = it->first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          NodeID d1 = it->second;
+          expr_t d1 = it->second;
           unsigned int eq = it->first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           unsigned int var = symbol_table.getTypeSpecificID(symb);
@@ -852,10 +852,10 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
   string tmp_s;
   ostringstream tmp_output;
   ofstream code_file;
-  NodeID lhs = NULL, rhs = NULL;
+  expr_t lhs = NULL, rhs = NULL;
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
-  map<NodeID, int> reference_count;
+  map<expr_t, int> reference_count;
   vector<int> feedback_variables;
   bool file_open = false;
 
@@ -960,7 +960,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
               }
               if (equ_type == E_EVALUATE)
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   lhs = eq_node->get_arg1();
                   rhs = eq_node->get_arg2();
                   rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
@@ -968,7 +968,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
                 }
               else if (equ_type == E_EVALUATE_S)
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->get_arg1();
                   rhs = eq_node->get_arg2();
                   rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
@@ -990,7 +990,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
             end:
               FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
               fnumexpr.write(code_file);
-              eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
               lhs = eq_node->get_arg1();
               rhs = eq_node->get_arg2();
               lhs->compile(code_file, false, temporary_terms, map_idx, true, false);
@@ -1660,7 +1660,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
     {
       int eq = it->first.first;
       int var = it->first.second;
-      NodeID d1 = it->second;
+      expr_t d1 = it->second;
 
       jacobian_output << "g1";
       jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
@@ -1677,7 +1677,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
       int eq = it->first.first;
       int var1 = it->first.second.first;
       int var2 = it->first.second.second;
-      NodeID d2 = it->second;
+      expr_t d2 = it->second;
 
       int id1 = getDynJacobianCol(var1);
       int id2 = getDynJacobianCol(var2);
@@ -1725,7 +1725,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
       int var1 = it->first.second.first;
       int var2 = it->first.second.second.first;
       int var3 = it->first.second.second.second;
-      NodeID d3 = it->second;
+      expr_t d3 = it->second;
 
       int id1 = getDynJacobianCol(var1);
       int id2 = getDynJacobianCol(var2);
@@ -2074,10 +2074,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
 
 }
 
-map<pair<int, pair<int, int > >, NodeID>
+map<pair<int, pair<int, int > >, expr_t>
 DynamicModel::collect_first_order_derivatives_endogenous()
 {
-  map<pair<int, pair<int, int > >, NodeID> endo_derivatives;
+  map<pair<int, pair<int, int > >, expr_t> endo_derivatives;
   for (first_derivatives_t::iterator it2 = first_derivatives.begin();
        it2 != first_derivatives.end(); it2++)
     {
@@ -2153,7 +2153,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
       computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
 
-      map<pair<int, pair<int, int> >, NodeID> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
+      map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
 
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
@@ -2243,7 +2243,7 @@ DynamicModel::get_Derivatives(int block)
 void
 DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
 {
-  map<int, NodeID> recursive_variables;
+  map<int, expr_t> recursive_variables;
   unsigned int nb_blocks = getNbBlocks();
   blocks_derivatives = blocks_derivatives_t(nb_blocks);
   for (unsigned int block = 0; block < nb_blocks; block++)
@@ -2260,9 +2260,9 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
           for (int i = 0; i < block_nb_recursives; i++)
             {
               if (getBlockEquationType(block, i) == E_EVALUATE_S)
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
               else
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
             }
           map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> Derivatives = get_Derivatives(block);
           map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>::const_iterator it = Derivatives.begin();
@@ -2297,9 +2297,9 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
           for (int i = 0; i < block_nb_recursives; i++)
             {
               if (getBlockEquationType(block, i) == E_EVALUATE_S)
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
               else
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
             }
           for (int eq = block_nb_recursives; eq < block_size; eq++)
             {
@@ -2307,7 +2307,7 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
               for (int var = block_nb_recursives; var < block_size; var++)
                 {
                   int varr = getBlockVariableID(block, var);
-                  NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
+                  expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
                   if (d1 == Zero)
                     continue;
                   first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
@@ -2464,7 +2464,7 @@ DynamicModel::toStatic(StaticModel &static_model) const
   assert(&symbol_table == &static_model.symbol_table);
 
   // Convert model local variables (need to be done first)
-  for (map<int, NodeID>::const_iterator it = local_variables_table.begin();
+  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     static_model.AddLocalVariable(it->first, it->second->toStatic(static_model));
 
@@ -2649,7 +2649,7 @@ DynamicModel::computeParamsDerivatives()
 
       for (int eq = 0; eq < (int) equations.size(); eq++)
         {
-          NodeID d1 = equations[eq]->getDerivative(param);
+          expr_t d1 = equations[eq]->getDerivative(param);
           if (d1 == Zero)
             continue;
           residuals_params_derivatives[make_pair(eq, param)] = d1;
@@ -2660,9 +2660,9 @@ DynamicModel::computeParamsDerivatives()
         {
           int eq = it2->first.first;
           int param1 = it2->first.second;
-          NodeID d1 = it2->second;
+          expr_t d1 = it2->second;
 
-          NodeID d2 = d1->getDerivative(param);
+          expr_t d2 = d1->getDerivative(param);
           if (d2 == Zero)
             continue;
           residuals_params_second_derivatives[make_pair(eq, make_pair(param1, param))] = d2;
@@ -2673,9 +2673,9 @@ DynamicModel::computeParamsDerivatives()
         {
           int eq = it2->first.first;
           int var = it2->first.second;
-          NodeID d1 = it2->second;
+          expr_t d1 = it2->second;
 
-          NodeID d2 = d1->getDerivative(param);
+          expr_t d2 = d1->getDerivative(param);
           if (d2 == Zero)
             continue;
           jacobian_params_derivatives[make_pair(eq, make_pair(var, param))] = d2;
@@ -2687,9 +2687,9 @@ DynamicModel::computeParamsDerivatives()
           int eq = it2->first.first;
           int var = it2->first.second.first;
           int param1 = it2->first.second.second;
-          NodeID d1 = it2->second;
+          expr_t d1 = it2->second;
 
-          NodeID d2 = d1->getDerivative(param);
+          expr_t d2 = d1->getDerivative(param);
           if (d2 == Zero)
             continue;
           jacobian_params_second_derivatives[make_pair(eq, make_pair(var, make_pair(param1, param)))] = d2;
@@ -2701,9 +2701,9 @@ DynamicModel::computeParamsDerivatives()
           int eq = it2->first.first;
           int var1 = it2->first.second.first;
           int var2 = it2->first.second.second;
-          NodeID d1 = it2->second;
+          expr_t d1 = it2->second;
 
-          NodeID d2 = d1->getDerivative(param);
+          expr_t d2 = d1->getDerivative(param);
           if (d2 == Zero)
             continue;
           hessian_params_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, param)))] = d2;
@@ -2714,7 +2714,7 @@ DynamicModel::computeParamsDerivatives()
 void
 DynamicModel::computeParamsDerivativesTemporaryTerms()
 {
-  map<NodeID, int> reference_count;
+  map<expr_t, int> reference_count;
   params_derivs_temporary_terms.clear();
 
   for (first_derivatives_t::iterator it = residuals_params_derivatives.begin();
@@ -2764,7 +2764,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
     {
       int eq = it->first.first;
       int param = it->first.second;
-      NodeID d1 = it->second;
+      expr_t d1 = it->second;
 
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
@@ -2783,7 +2783,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int eq = it->first.first;
       int var = it->first.second.first;
       int param = it->first.second.second;
-      NodeID d2 = it->second;
+      expr_t d2 = it->second;
 
       int var_col = getDynJacobianCol(var) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
@@ -2807,7 +2807,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int eq = it->first.first;
       int param1 = it->first.second.first;
       int param2 = it->first.second.second;
-      NodeID d2 = it->second;
+      expr_t d2 = it->second;
 
       int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
       int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
@@ -2832,7 +2832,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int var = it->first.second.first;
       int param1 = it->first.second.second.first;
       int param2 = it->first.second.second.second;
-      NodeID d2 = it->second;
+      expr_t d2 = it->second;
 
       int var_col = getDynJacobianCol(var) + 1;
       int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
@@ -2862,7 +2862,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int var1 = it->first.second.first;
       int var2 = it->first.second.second.first;
       int param = it->first.second.second.second;
-      NodeID d2 = it->second;
+      expr_t d2 = it->second;
 
       int var1_col = getDynJacobianCol(var1) + 1;
       int var2_col = getDynJacobianCol(var2) + 1;
@@ -2887,7 +2887,7 @@ DynamicModel::writeChainRuleDerivative(ostream &output, int eqr, int varr, int l
                                        ExprNodeOutputType output_type,
                                        const temporary_terms_t &temporary_terms) const
 {
-  map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
+  map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
   if (it != first_chain_rule_derivatives.end())
     (it->second)->writeOutput(output, output_type, temporary_terms);
   else
@@ -2953,10 +2953,10 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type)
   vector<BinaryOpNode *> neweqs;
 
   // Substitute in model local variables
-  for (map<int, NodeID>::iterator it = local_variables_table.begin();
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     {
-      NodeID subst;
+      expr_t subst;
       switch (type)
         {
         case avEndoLead:
@@ -2981,7 +2981,7 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type)
   // Substitute in equations
   for (int i = 0; i < (int) equations.size(); i++)
     {
-      NodeID subst;
+      expr_t subst;
       switch (type)
         {
         case avEndoLead:
@@ -3047,7 +3047,7 @@ DynamicModel::substituteExpectation(bool partial_information_model)
   vector<BinaryOpNode *> neweqs;
 
   // Substitute in model local variables
-  for (map<int, NodeID>::iterator it = local_variables_table.begin();
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     it->second = it->second->substituteExpectation(subst_table, neweqs, partial_information_model);
 
@@ -3108,12 +3108,12 @@ DynamicModel::fillEvalContext(eval_context_t &eval_context) const
     }
 
   // Second, model local variables
-  for (map<int, NodeID>::const_iterator it = local_variables_table.begin();
+  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     {
       try
         {
-          const NodeID expression = it->second;
+          const expr_t expression = it->second;
           double val = expression->eval(eval_context);
           eval_context[it->first] = val;
         }
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 3f76368f7000b9ae320f7c7402addfdaa0eb72fb..32e2296af291a8a49000264129df739aea361cd6 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -100,8 +100,8 @@ private:
 
   vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
 
-  //! Store the derivatives or the chainrule derivatives:map<pair< equation, pair< variable, lead_lag >, NodeID>
-  typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_t;
+  //! Store the derivatives or the chainrule derivatives:map<pair< equation, pair< variable, lead_lag >, expr_t>
+  typedef map< pair< int, pair< int, int> >, expr_t> first_chain_rule_derivatives_t;
   first_chain_rule_derivatives_t first_chain_rule_derivatives;
 
   //! Writes dynamic model file (Matlab version)
@@ -159,7 +159,7 @@ private:
   //! Computes temporary terms for the file containing parameters derivatives
   void computeParamsDerivativesTemporaryTerms();
   //! Collect only the first derivatives
-  map<pair<int, pair<int, int> >, NodeID> collect_first_order_derivatives_endogenous();
+  map<pair<int, pair<int, int> >, expr_t> collect_first_order_derivatives_endogenous();
 
   //! Allocates the derivation IDs for all dynamic variables of the model
   /*! Also computes max_{endo,exo}_{lead_lag}, and initializes dynJacobianColsNbr to the number of dynamic endos */
@@ -191,7 +191,7 @@ private:
   //! vector of block reordered variables and equations
   vector<int> equation_reordered, variable_reordered, inv_equation_reordered, inv_variable_reordered;
 
-  //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
+  //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
   equation_type_and_normalized_equation_t equation_type_and_normalized_equation;
 
   //! for each block contains pair< Simulation_Type, pair < Block_Size, Recursive_part_Size > >
@@ -206,8 +206,8 @@ private:
   //! Vector indicating if the block is linear in endogenous variable (true) or not (false)
   vector<bool> blocks_linear;
 
-  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), NodeID>
-  typedef map<pair< int, pair<int, int> >, NodeID> derivative_t;
+  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
+  typedef map<pair< int, pair<int, int> >, expr_t> derivative_t;
   //! Vector of derivative for each blocks
   vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
 
@@ -353,15 +353,15 @@ public:
   {
     return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
   };
-  //! Return the NodeID of the equation equation_number belonging to the block block_number
-  virtual NodeID
-  getBlockEquationNodeID(int block_number, int equation_number) const
+  //! Return the expr_t of the equation equation_number belonging to the block block_number
+  virtual expr_t
+  getBlockEquationExpr(int block_number, int equation_number) const
   {
     return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
   };
-  //! Return the NodeID of the renormalized equation equation_number belonging to the block block_number
-  virtual NodeID
-  getBlockEquationRenormalizedNodeID(int block_number, int equation_number) const
+  //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
+  virtual expr_t
+  getBlockEquationRenormalizedExpr(int block_number, int equation_number) const
   {
     return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
   };
diff --git a/DynareBison.yy b/DynareBison.yy
index 08de7ab4bc2ee5b9793e2877ec5e5493f8c6c425..1cdeb212e09b15af353850cbfb348322596e116f 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -71,7 +71,7 @@ class ParsingDriver;
 %union
 {
   string *string_val;
-  NodeID node_val;
+  expr_t node_val;
   SymbolType symbol_type_val;
   vector<string *> *vector_string_val;
   vector<int> *vector_int_val;
diff --git a/ExprNode.cc b/ExprNode.cc
index 47e0d2d18aa1e9fd578803a2ee5d084d9bb8b353..6dedb0c67c0de6d6d6fc55b337ed661f0ef3cecf 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -47,7 +47,7 @@ ExprNode::~ExprNode()
 {
 }
 
-NodeID
+expr_t
 ExprNode::getDerivative(int deriv_id)
 {
   if (!preparedForDerivation)
@@ -59,12 +59,12 @@ ExprNode::getDerivative(int deriv_id)
     return datatree.Zero;
 
   // If derivative is stored in cache, use the cached value, otherwise compute it (and cache it)
-  map<int, NodeID>::const_iterator it2 = derivatives.find(deriv_id);
+  map<int, expr_t>::const_iterator it2 = derivatives.find(deriv_id);
   if (it2 != derivatives.end())
     return it2->second;
   else
     {
-      NodeID d = computeDerivative(deriv_id);
+      expr_t d = computeDerivative(deriv_id);
       derivatives[deriv_id] = d;
       return d;
     }
@@ -114,7 +114,7 @@ ExprNode::collectModelLocalVariables(set<int> &result) const
 }
 
 void
-ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                 temporary_terms_t &temporary_terms,
                                 bool is_matlab) const
 {
@@ -122,9 +122,9 @@ ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
 }
 
 void
-ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ExprNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                 temporary_terms_t &temporary_terms,
-                                map<NodeID, pair<int, int> > &first_occurence,
+                                map<expr_t, pair<int, int> > &first_occurence,
                                 int Curr_block,
                                 vector<vector<temporary_terms_t> > &v_temporary_terms,
                                 int equation) const
@@ -132,10 +132,10 @@ ExprNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
   // Nothing to do for a terminal node
 }
 
-pair<int, NodeID >
-ExprNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+pair<int, expr_t >
+ExprNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
-  return (make_pair(0, (NodeID) NULL));
+  return (make_pair(0, (expr_t) NULL));
 }
 
 void
@@ -175,14 +175,14 @@ ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector
   if (it != subst_table.end())
     return const_cast<VariableNode *>(it->second);
 
-  NodeID substexpr = decreaseLeadsLags(n-1);
+  expr_t substexpr = decreaseLeadsLags(n-1);
   int lag = n-2;
 
   // Each iteration tries to create an auxvar such that auxvar(+1)=expr(-lag)
   // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to expr(-lag-1) (resp. expr(-lag))
   while (lag >= 0)
     {
-      NodeID orig_expr = decreaseLeadsLags(lag);
+      expr_t orig_expr = decreaseLeadsLags(lag);
       it = subst_table.find(orig_expr);
       if (it == subst_table.end())
         {
@@ -211,14 +211,14 @@ ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<
   if (it != subst_table.end())
     return const_cast<VariableNode *>(it->second);
 
-  NodeID substexpr = decreaseLeadsLags(n);
+  expr_t substexpr = decreaseLeadsLags(n);
   int lag = n-1;
 
   // Each iteration tries to create an auxvar such that auxvar(+1)=expr(-lag)
   // At the beginning (resp. end) of each iteration, substexpr is an expression (possibly an auxvar) equivalent to expr(-lag-1) (resp. expr(-lag))
   while (lag >= 0)
     {
-      NodeID orig_expr = decreaseLeadsLags(lag);
+      expr_t orig_expr = decreaseLeadsLags(lag);
       it = subst_table.find(orig_expr);
       if (it == subst_table.end())
         {
@@ -264,7 +264,7 @@ NumConstNode::prepareForDerivation()
   // All derivatives are null, so non_null_derivatives is left empty
 }
 
-NodeID
+expr_t
 NumConstNode::computeDerivative(int deriv_id)
 {
   return datatree.Zero;
@@ -311,19 +311,19 @@ NumConstNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result
 {
 }
 
-pair<int, NodeID >
-NumConstNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+pair<int, expr_t >
+NumConstNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
   return (make_pair(0, datatree.AddNumConstant(datatree.num_constants.get(id))));
 }
 
-NodeID
-NumConstNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables)
+expr_t
+NumConstNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
   return datatree.Zero;
 }
 
-NodeID
+expr_t
 NumConstNode::toStatic(DataTree &static_datatree) const
 {
   return static_datatree.AddNumConstant(datatree.num_constants.get(id));
@@ -353,43 +353,43 @@ NumConstNode::maxExoLag() const
   return 0;
 }
 
-NodeID
+expr_t
 NumConstNode::decreaseLeadsLags(int n) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
-NodeID
+expr_t
 NumConstNode::decreaseLeadsLagsPredeterminedVariables() const
 {
   return const_cast<NumConstNode *>(this);
 }
 
-NodeID
+expr_t
 NumConstNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
-NodeID
+expr_t
 NumConstNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
-NodeID
+expr_t
 NumConstNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
-NodeID
+expr_t
 NumConstNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   return const_cast<NumConstNode *>(this);
 }
 
-NodeID
+expr_t
 NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
   return const_cast<NumConstNode *>(this);
@@ -456,7 +456,7 @@ VariableNode::prepareForDerivation()
     }
 }
 
-NodeID
+expr_t
 VariableNode::computeDerivative(int deriv_id)
 {
   switch (type)
@@ -760,9 +760,9 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
 }
 
 void
-VariableNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+VariableNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                     temporary_terms_t &temporary_terms,
-                                    map<NodeID, pair<int, int> > &first_occurence,
+                                    map<expr_t, pair<int, int> > &first_occurence,
                                     int Curr_block,
                                     vector<vector<temporary_terms_t> > &v_temporary_terms,
                                     int equation) const
@@ -780,13 +780,13 @@ VariableNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result
     datatree.local_variables_table[symb_id]->collectVariables(type_arg, result);
 }
 
-pair<int, NodeID>
-VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+pair<int, expr_t>
+VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
   if (type == eEndogenous)
     {
       if (datatree.symbol_table.getTypeSpecificID(symb_id) == var_endo && lag == 0)
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       else
         return (make_pair(0, datatree.AddVariableInternal(symb_id, lag)));
     }
@@ -799,8 +799,8 @@ VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, Node
     }
 }
 
-NodeID
-VariableNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables)
+expr_t
+VariableNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
   switch (type)
     {
@@ -813,18 +813,18 @@ VariableNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recur
       else
         {
           //if there is in the equation a recursive variable we could use a chaine rule derivation
-          map<int, NodeID>::const_iterator it = recursive_variables.find(datatree.getDerivID(symb_id, lag));
+          map<int, expr_t>::const_iterator it = recursive_variables.find(datatree.getDerivID(symb_id, lag));
           if (it != recursive_variables.end())
             {
-              map<int, NodeID>::const_iterator it2 = derivatives.find(deriv_id);
+              map<int, expr_t>::const_iterator it2 = derivatives.find(deriv_id);
               if (it2 != derivatives.end())
                 return it2->second;
               else
                 {
-                  map<int, NodeID> recursive_vars2(recursive_variables);
+                  map<int, expr_t> recursive_vars2(recursive_variables);
                   recursive_vars2.erase(it->first);
-                  //NodeID c = datatree.AddNumConstant("1");
-                  NodeID d = datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id, recursive_vars2));
+                  //expr_t c = datatree.AddNumConstant("1");
+                  expr_t d = datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id, recursive_vars2));
                   //d = datatree.AddTimes(c, d);
                   derivatives[deriv_id] = d;
                   return d;
@@ -846,7 +846,7 @@ VariableNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recur
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 VariableNode::toStatic(DataTree &static_datatree) const
 {
   return static_datatree.AddVariable(symb_id);
@@ -908,7 +908,7 @@ VariableNode::maxExoLag() const
     }
 }
 
-NodeID
+expr_t
 VariableNode::decreaseLeadsLags(int n) const
 {
   switch (type)
@@ -924,7 +924,7 @@ VariableNode::decreaseLeadsLags(int n) const
     }
 }
 
-NodeID
+expr_t
 VariableNode::decreaseLeadsLagsPredeterminedVariables() const
 {
   if (datatree.symbol_table.isPredetermined(symb_id))
@@ -933,10 +933,10 @@ VariableNode::decreaseLeadsLagsPredeterminedVariables() const
     return const_cast<VariableNode *>(this);
 }
 
-NodeID
+expr_t
 VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID value;
+  expr_t value;
   switch (type)
     {
     case eEndogenous:
@@ -955,11 +955,11 @@ VariableNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vecto
     }
 }
 
-NodeID
+expr_t
 VariableNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   VariableNode *substexpr;
-  NodeID value;
+  expr_t value;
   subst_table_t::const_iterator it;
   int cur_lag;
   switch (type)
@@ -1006,10 +1006,10 @@ VariableNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
     }
 }
 
-NodeID
+expr_t
 VariableNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID value;
+  expr_t value;
   switch (type)
     {
     case eExogenous:
@@ -1028,11 +1028,11 @@ VariableNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode
     }
 }
 
-NodeID
+expr_t
 VariableNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   VariableNode *substexpr;
-  NodeID value;
+  expr_t value;
   subst_table_t::const_iterator it;
   int cur_lag;
   switch (type)
@@ -1079,7 +1079,7 @@ VariableNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *
     }
 }
 
-NodeID
+expr_t
 VariableNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
   return const_cast<VariableNode *>(this);
@@ -1100,7 +1100,7 @@ VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
     return false;
 }
 
-UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg) :
+UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg) :
   ExprNode(datatree_arg),
   arg(arg_arg),
   expectation_information_set(expectation_information_set_arg),
@@ -1125,10 +1125,10 @@ UnaryOpNode::prepareForDerivation()
   non_null_derivatives = arg->non_null_derivatives;
 }
 
-NodeID
-UnaryOpNode::composeDerivatives(NodeID darg)
+expr_t
+UnaryOpNode::composeDerivatives(expr_t darg)
 {
-  NodeID t11, t12, t13;
+  expr_t t11, t12, t13;
 
   switch (op_code)
     {
@@ -1212,10 +1212,10 @@ UnaryOpNode::composeDerivatives(NodeID darg)
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 UnaryOpNode::computeDerivative(int deriv_id)
 {
-  NodeID darg = arg->getDerivative(deriv_id);
+  expr_t darg = arg->getDerivative(deriv_id);
   return composeDerivatives(darg);
 }
 
@@ -1313,13 +1313,13 @@ UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) cons
 }
 
 void
-UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+UnaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                    temporary_terms_t &temporary_terms,
                                    bool is_matlab) const
 {
-  NodeID this2 = const_cast<UnaryOpNode *>(this);
+  expr_t this2 = const_cast<UnaryOpNode *>(this);
 
-  map<NodeID, int>::iterator it = reference_count.find(this2);
+  map<expr_t, int>::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       reference_count[this2] = 1;
@@ -1334,15 +1334,15 @@ UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
 }
 
 void
-UnaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+UnaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                    temporary_terms_t &temporary_terms,
-                                   map<NodeID, pair<int, int> > &first_occurence,
+                                   map<expr_t, pair<int, int> > &first_occurence,
                                    int Curr_block,
                                    vector< vector<temporary_terms_t> > &v_temporary_terms,
                                    int equation) const
 {
-  NodeID this2 = const_cast<UnaryOpNode *>(this);
-  map<NodeID, int>::iterator it = reference_count.find(this2);
+  expr_t this2 = const_cast<UnaryOpNode *>(this);
+  map<expr_t, int>::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       reference_count[this2] = 1;
@@ -1602,64 +1602,64 @@ UnaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result)
   arg->collectVariables(type_arg, result);
 }
 
-pair<int, NodeID>
-UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+pair<int, expr_t>
+UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
-  pair<bool, NodeID > res = arg->normalizeEquation(var_endo, List_of_Op_RHS);
+  pair<bool, expr_t > res = arg->normalizeEquation(var_endo, List_of_Op_RHS);
   int is_endogenous_present = res.first;
-  NodeID New_NodeID = res.second;
+  expr_t New_expr_t = res.second;
   /*if(res.second.second)*/
   if (is_endogenous_present == 2)
-    return (make_pair(2, (NodeID) NULL));
+    return (make_pair(2, (expr_t) NULL));
   else if (is_endogenous_present)
     {
       switch (op_code)
         {
         case oUminus:
-          List_of_Op_RHS.push_back(make_pair(oUminus, make_pair((NodeID) NULL, (NodeID) NULL)));
-          return (make_pair(1, (NodeID) NULL));
+          List_of_Op_RHS.push_back(make_pair(oUminus, make_pair((expr_t) NULL, (expr_t) NULL)));
+          return (make_pair(1, (expr_t) NULL));
         case oExp:
-          List_of_Op_RHS.push_back(make_pair(oLog, make_pair((NodeID) NULL, (NodeID) NULL)));
-          return (make_pair(1, (NodeID) NULL));
+          List_of_Op_RHS.push_back(make_pair(oLog, make_pair((expr_t) NULL, (expr_t) NULL)));
+          return (make_pair(1, (expr_t) NULL));
         case oLog:
-          List_of_Op_RHS.push_back(make_pair(oExp, make_pair((NodeID) NULL, (NodeID) NULL)));
-          return (make_pair(1, (NodeID) NULL));
+          List_of_Op_RHS.push_back(make_pair(oExp, make_pair((expr_t) NULL, (expr_t) NULL)));
+          return (make_pair(1, (expr_t) NULL));
         case oLog10:
-          List_of_Op_RHS.push_back(make_pair(oPower, make_pair((NodeID) NULL, datatree.AddNumConstant("10"))));
-          return (make_pair(1, (NodeID) NULL));
+          List_of_Op_RHS.push_back(make_pair(oPower, make_pair((expr_t) NULL, datatree.AddNumConstant("10"))));
+          return (make_pair(1, (expr_t) NULL));
         case oCos:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oSin:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oTan:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oAcos:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oAsin:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oAtan:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oCosh:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oSinh:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oTanh:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oAcosh:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oAsinh:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oAtanh:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oSqrt:
-          List_of_Op_RHS.push_back(make_pair(oPower, make_pair((NodeID) NULL, datatree.Two)));
-          return (make_pair(1, (NodeID) NULL));
+          List_of_Op_RHS.push_back(make_pair(oPower, make_pair((expr_t) NULL, datatree.Two)));
+          return (make_pair(1, (expr_t) NULL));
         case oSteadyState:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         case oExpectation:
           assert(0);
         case oErf:
-          return (make_pair(1, (NodeID) NULL));
+          return (make_pair(1, (expr_t) NULL));
         }
     }
   else
@@ -1667,59 +1667,59 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeI
       switch (op_code)
         {
         case oUminus:
-          return (make_pair(0, datatree.AddUMinus(New_NodeID)));
+          return (make_pair(0, datatree.AddUMinus(New_expr_t)));
         case oExp:
-          return (make_pair(0, datatree.AddExp(New_NodeID)));
+          return (make_pair(0, datatree.AddExp(New_expr_t)));
         case oLog:
-          return (make_pair(0, datatree.AddLog(New_NodeID)));
+          return (make_pair(0, datatree.AddLog(New_expr_t)));
         case oLog10:
-          return (make_pair(0, datatree.AddLog10(New_NodeID)));
+          return (make_pair(0, datatree.AddLog10(New_expr_t)));
         case oCos:
-          return (make_pair(0, datatree.AddCos(New_NodeID)));
+          return (make_pair(0, datatree.AddCos(New_expr_t)));
         case oSin:
-          return (make_pair(0, datatree.AddSin(New_NodeID)));
+          return (make_pair(0, datatree.AddSin(New_expr_t)));
         case oTan:
-          return (make_pair(0, datatree.AddTan(New_NodeID)));
+          return (make_pair(0, datatree.AddTan(New_expr_t)));
         case oAcos:
-          return (make_pair(0, datatree.AddAcos(New_NodeID)));
+          return (make_pair(0, datatree.AddAcos(New_expr_t)));
         case oAsin:
-          return (make_pair(0, datatree.AddAsin(New_NodeID)));
+          return (make_pair(0, datatree.AddAsin(New_expr_t)));
         case oAtan:
-          return (make_pair(0, datatree.AddAtan(New_NodeID)));
+          return (make_pair(0, datatree.AddAtan(New_expr_t)));
         case oCosh:
-          return (make_pair(0, datatree.AddCosh(New_NodeID)));
+          return (make_pair(0, datatree.AddCosh(New_expr_t)));
         case oSinh:
-          return (make_pair(0, datatree.AddSinh(New_NodeID)));
+          return (make_pair(0, datatree.AddSinh(New_expr_t)));
         case oTanh:
-          return (make_pair(0, datatree.AddTanh(New_NodeID)));
+          return (make_pair(0, datatree.AddTanh(New_expr_t)));
         case oAcosh:
-          return (make_pair(0, datatree.AddAcosh(New_NodeID)));
+          return (make_pair(0, datatree.AddAcosh(New_expr_t)));
         case oAsinh:
-          return (make_pair(0, datatree.AddAsinh(New_NodeID)));
+          return (make_pair(0, datatree.AddAsinh(New_expr_t)));
         case oAtanh:
-          return (make_pair(0, datatree.AddAtanh(New_NodeID)));
+          return (make_pair(0, datatree.AddAtanh(New_expr_t)));
         case oSqrt:
-          return (make_pair(0, datatree.AddSqrt(New_NodeID)));
+          return (make_pair(0, datatree.AddSqrt(New_expr_t)));
         case oSteadyState:
-          return (make_pair(0, datatree.AddSteadyState(New_NodeID)));
+          return (make_pair(0, datatree.AddSteadyState(New_expr_t)));
         case oExpectation:
           assert(0);
         case oErf:
-          return (make_pair(0, datatree.AddErf(New_NodeID)));
+          return (make_pair(0, datatree.AddErf(New_expr_t)));
         }
     }
-  return (make_pair(1, (NodeID) NULL));
+  return (make_pair(1, (expr_t) NULL));
 }
 
-NodeID
-UnaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables)
+expr_t
+UnaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
-  NodeID darg = arg->getChainRuleDerivative(deriv_id, recursive_variables);
+  expr_t darg = arg->getChainRuleDerivative(deriv_id, recursive_variables);
   return composeDerivatives(darg);
 }
 
-NodeID
-UnaryOpNode::buildSimilarUnaryOpNode(NodeID alt_arg, DataTree &alt_datatree) const
+expr_t
+UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) const
 {
   switch (op_code)
     {
@@ -1768,10 +1768,10 @@ UnaryOpNode::buildSimilarUnaryOpNode(NodeID alt_arg, DataTree &alt_datatree) con
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 UnaryOpNode::toStatic(DataTree &static_datatree) const
 {
-  NodeID sarg = arg->toStatic(static_datatree);
+  expr_t sarg = arg->toStatic(static_datatree);
   return buildSimilarUnaryOpNode(sarg, static_datatree);
 }
 
@@ -1799,26 +1799,26 @@ UnaryOpNode::maxExoLag() const
   return arg->maxExoLag();
 }
 
-NodeID
+expr_t
 UnaryOpNode::decreaseLeadsLags(int n) const
 {
-  NodeID argsubst = arg->decreaseLeadsLags(n);
+  expr_t argsubst = arg->decreaseLeadsLags(n);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
-NodeID
+expr_t
 UnaryOpNode::decreaseLeadsLagsPredeterminedVariables() const
 {
-  NodeID argsubst = arg->decreaseLeadsLagsPredeterminedVariables();
+  expr_t argsubst = arg->decreaseLeadsLagsPredeterminedVariables();
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
-NodeID
+expr_t
 UnaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   if (op_code == oUminus)
     {
-      NodeID argsubst = arg->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
+      expr_t argsubst = arg->substituteEndoLeadGreaterThanTwo(subst_table, neweqs);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
   else
@@ -1830,19 +1830,19 @@ UnaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector
     }
 }
 
-NodeID
+expr_t
 UnaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID argsubst = arg->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+  expr_t argsubst = arg->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
-NodeID
+expr_t
 UnaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   if (op_code == oUminus)
     {
-      NodeID argsubst = arg->substituteExoLead(subst_table, neweqs);
+      expr_t argsubst = arg->substituteExoLead(subst_table, neweqs);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
   else
@@ -1854,14 +1854,14 @@ UnaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *
     }
 }
 
-NodeID
+expr_t
 UnaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID argsubst = arg->substituteExoLag(subst_table, neweqs);
+  expr_t argsubst = arg->substituteExoLag(subst_table, neweqs);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
-NodeID
+expr_t
 UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
   switch (op_code)
@@ -1876,7 +1876,7 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
         //AUX_EXPECT_(LEAD/LAG)_(period)_(arg.idx) OR
         //AUX_EXPECT_(info_set_name)_(arg.idx)
         int symb_id = datatree.symbol_table.addExpectationAuxiliaryVar(expectation_information_set, arg->idx, expectation_information_set_name);
-        NodeID newAuxE = datatree.AddVariable(symb_id, 0);
+        expr_t newAuxE = datatree.AddVariable(symb_id, 0);
 
         if (partial_information_model && expectation_information_set == 0)
           {
@@ -1917,7 +1917,7 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
           {
             //take care of any nested expectation operators by calling arg->substituteExpectation(.), then decreaseLeadsLags for this oExpectation operator
             //arg(lag-period) (holds entire subtree of arg(lag-period)
-            NodeID substexpr = (arg->substituteExpectation(subst_table, neweqs, partial_information_model))->decreaseLeadsLags(expectation_information_set);
+            expr_t substexpr = (arg->substituteExpectation(subst_table, neweqs, partial_information_model))->decreaseLeadsLags(expectation_information_set);
             assert(substexpr != NULL);
             neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(newAuxE, substexpr))); //AUXE_period_arg.idx = arg(lag-period)
             newAuxE = datatree.AddVariable(symb_id, expectation_information_set);
@@ -1927,7 +1927,7 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
         return newAuxE;
       }
     default:
-      NodeID argsubst = arg->substituteExpectation(subst_table, neweqs, partial_information_model);
+      expr_t argsubst = arg->substituteExpectation(subst_table, neweqs, partial_information_model);
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
 }
@@ -1944,8 +1944,8 @@ UnaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag
   return false;
 }
 
-BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
-                           BinaryOpcode op_code_arg, const NodeID arg2_arg) :
+BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
+                           BinaryOpcode op_code_arg, const expr_t arg2_arg) :
   ExprNode(datatree_arg),
   arg1(arg1_arg),
   arg2(arg2_arg),
@@ -1974,10 +1974,10 @@ BinaryOpNode::prepareForDerivation()
             inserter(non_null_derivatives, non_null_derivatives.begin()));
 }
 
-NodeID
-BinaryOpNode::composeDerivatives(NodeID darg1, NodeID darg2)
+expr_t
+BinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2)
 {
-  NodeID t11, t12, t13, t14, t15;
+  expr_t t11, t12, t13, t14, t15;
 
   switch (op_code)
     {
@@ -2048,11 +2048,11 @@ BinaryOpNode::composeDerivatives(NodeID darg1, NodeID darg2)
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 BinaryOpNode::computeDerivative(int deriv_id)
 {
-  NodeID darg1 = arg1->getDerivative(deriv_id);
-  NodeID darg2 = arg2->getDerivative(deriv_id);
+  expr_t darg1 = arg1->getDerivative(deriv_id);
+  expr_t darg2 = arg2->getDerivative(deriv_id);
   return composeDerivatives(darg1, darg2);
 }
 
@@ -2162,12 +2162,12 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con
 }
 
 void
-BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                     temporary_terms_t &temporary_terms,
                                     bool is_matlab) const
 {
-  NodeID this2 = const_cast<BinaryOpNode *>(this);
-  map<NodeID, int>::iterator it = reference_count.find(this2);
+  expr_t this2 = const_cast<BinaryOpNode *>(this);
+  map<expr_t, int>::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       // If this node has never been encountered, set its ref count to one,
@@ -2189,15 +2189,15 @@ BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
 }
 
 void
-BinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+BinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                     temporary_terms_t &temporary_terms,
-                                    map<NodeID, pair<int, int> > &first_occurence,
+                                    map<expr_t, pair<int, int> > &first_occurence,
                                     int Curr_block,
                                     vector<vector<temporary_terms_t> > &v_temporary_terms,
                                     int equation) const
 {
-  NodeID this2 = const_cast<BinaryOpNode *>(this);
-  map<NodeID, int>::iterator it = reference_count.find(this2);
+  expr_t this2 = const_cast<BinaryOpNode *>(this);
+  map<expr_t, int>::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       reference_count[this2] = 1;
@@ -2487,8 +2487,8 @@ BinaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result
   arg2->collectVariables(type_arg, result);
 }
 
-NodeID
-BinaryOpNode::Compute_RHS(NodeID arg1, NodeID arg2, int op, int op_type) const
+expr_t
+BinaryOpNode::Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const
 {
   temporary_terms_t temp;
   switch (op_type)
@@ -2531,45 +2531,45 @@ BinaryOpNode::Compute_RHS(NodeID arg1, NodeID arg2, int op, int op_type) const
         }
       break;
     }
-  return ((NodeID) NULL);
+  return ((expr_t) NULL);
 }
 
-pair<int, NodeID>
-BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+pair<int, expr_t>
+BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
-  vector<pair<int, pair<NodeID, NodeID> > > List_of_Op_RHS1, List_of_Op_RHS2;
+  vector<pair<int, pair<expr_t, expr_t> > > List_of_Op_RHS1, List_of_Op_RHS2;
   int is_endogenous_present_1, is_endogenous_present_2;
-  pair<int, NodeID> res;
-  NodeID NodeID_1, NodeID_2;
+  pair<int, expr_t> res;
+  expr_t expr_t_1, expr_t_2;
   res = arg1->normalizeEquation(var_endo, List_of_Op_RHS1);
   is_endogenous_present_1 = res.first;
-  NodeID_1 = res.second;
+  expr_t_1 = res.second;
 
   res = arg2->normalizeEquation(var_endo, List_of_Op_RHS2);
   is_endogenous_present_2 = res.first;
-  NodeID_2 = res.second;
+  expr_t_2 = res.second;
   if (is_endogenous_present_1 == 2 || is_endogenous_present_2 == 2)
-    return (make_pair(2, (NodeID) NULL));
+    return (make_pair(2, (expr_t) NULL));
   else if (is_endogenous_present_1 && is_endogenous_present_2)
-    return (make_pair(2, (NodeID) NULL));
+    return (make_pair(2, (expr_t) NULL));
   else if (is_endogenous_present_1)
     {
       if (op_code == oEqual)
         {
-          pair<int, pair<NodeID, NodeID> > it;
+          pair<int, pair<expr_t, expr_t> > it;
           int oo = List_of_Op_RHS1.size();
           for (int i = 0; i < oo; i++)
             {
               it = List_of_Op_RHS1.back();
               List_of_Op_RHS1.pop_back();
               if (it.second.first && !it.second.second) /*Binary operator*/
-                NodeID_2 = Compute_RHS(NodeID_2, (BinaryOpNode *) it.second.first, it.first, 1);
+                expr_t_2 = Compute_RHS(expr_t_2, (BinaryOpNode *) it.second.first, it.first, 1);
               else if (it.second.second && !it.second.first) /*Binary operator*/
-                NodeID_2 = Compute_RHS(it.second.second, NodeID_2, it.first, 1);
+                expr_t_2 = Compute_RHS(it.second.second, expr_t_2, it.first, 1);
               else if (it.second.second && it.second.first) /*Binary operator*/
-                NodeID_2 = Compute_RHS(it.second.first, it.second.second, it.first, 1);
+                expr_t_2 = Compute_RHS(it.second.first, it.second.second, it.first, 1);
               else                                          /*Unary operator*/
-                NodeID_2 = Compute_RHS((UnaryOpNode *) NodeID_2, (UnaryOpNode *) it.second.first, it.first, 0);
+                expr_t_2 = Compute_RHS((UnaryOpNode *) expr_t_2, (UnaryOpNode *) it.second.first, it.first, 0);
             }
         }
       else
@@ -2582,17 +2582,17 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, Node
           int oo = List_of_Op_RHS2.size();
           for (int i = 0; i < oo; i++)
             {
-              pair<int, pair<NodeID, NodeID> > it;
+              pair<int, pair<expr_t, expr_t> > it;
               it = List_of_Op_RHS2.back();
               List_of_Op_RHS2.pop_back();
               if (it.second.first && !it.second.second) /*Binary operator*/
-                NodeID_1 = Compute_RHS((BinaryOpNode *) NodeID_1, (BinaryOpNode *) it.second.first, it.first, 1);
+                expr_t_1 = Compute_RHS((BinaryOpNode *) expr_t_1, (BinaryOpNode *) it.second.first, it.first, 1);
               else if (it.second.second && !it.second.first) /*Binary operator*/
-                NodeID_1 = Compute_RHS((BinaryOpNode *) it.second.second, (BinaryOpNode *) NodeID_1, it.first, 1);
+                expr_t_1 = Compute_RHS((BinaryOpNode *) it.second.second, (BinaryOpNode *) expr_t_1, it.first, 1);
               else if (it.second.second && it.second.first) /*Binary operator*/
-                NodeID_1 = Compute_RHS(it.second.first, it.second.second, it.first, 1);
+                expr_t_1 = Compute_RHS(it.second.first, it.second.second, it.first, 1);
               else
-                NodeID_1 = Compute_RHS((UnaryOpNode *) NodeID_1, (UnaryOpNode *) it.second.first, it.first, 0);
+                expr_t_1 = Compute_RHS((UnaryOpNode *) expr_t_1, (UnaryOpNode *) it.second.first, it.first, 0);
             }
         }
       else
@@ -2603,88 +2603,88 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, Node
     case oPlus:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(datatree.AddPlus(NodeID_1, NodeID_2), (NodeID) NULL)));
-          return (make_pair(0, datatree.AddPlus(NodeID_1, NodeID_2)));
+          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(datatree.AddPlus(expr_t_1, expr_t_2), (expr_t) NULL)));
+          return (make_pair(0, datatree.AddPlus(expr_t_1, expr_t_2)));
         }
       else if (is_endogenous_present_1 && is_endogenous_present_2)
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(NodeID_1, (NodeID) NULL)));
-          return (make_pair(1, NodeID_1));
+          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(expr_t_1, (expr_t) NULL)));
+          return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(NodeID_2, (NodeID) NULL)));
-          return (make_pair(1, NodeID_2));
+          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(expr_t_2, (expr_t) NULL)));
+          return (make_pair(1, expr_t_2));
         }
       break;
     case oMinus:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(datatree.AddMinus(NodeID_1, NodeID_2), (NodeID) NULL)));
-          return (make_pair(0, datatree.AddMinus(NodeID_1, NodeID_2)));
+          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(datatree.AddMinus(expr_t_1, expr_t_2), (expr_t) NULL)));
+          return (make_pair(0, datatree.AddMinus(expr_t_1, expr_t_2)));
         }
       else if (is_endogenous_present_1 && is_endogenous_present_2)
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oUminus, make_pair((NodeID) NULL, (NodeID) NULL)));
-          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(NodeID_1, (NodeID) NULL)));
-          return (make_pair(1, NodeID_1));
+          List_of_Op_RHS.push_back(make_pair(oUminus, make_pair((expr_t) NULL, (expr_t) NULL)));
+          List_of_Op_RHS.push_back(make_pair(oMinus, make_pair(expr_t_1, (expr_t) NULL)));
+          return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oPlus, make_pair(NodeID_2, (NodeID) NULL)));
-          return (make_pair(1, datatree.AddUMinus(NodeID_2)));
+          List_of_Op_RHS.push_back(make_pair(oPlus, make_pair(expr_t_2, (expr_t) NULL)));
+          return (make_pair(1, datatree.AddUMinus(expr_t_2)));
         }
       break;
     case oTimes:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddTimes(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddTimes(expr_t_1, expr_t_2)));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oDivide, make_pair(NodeID_1, (NodeID) NULL)));
-          return (make_pair(1, NodeID_1));
+          List_of_Op_RHS.push_back(make_pair(oDivide, make_pair(expr_t_1, (expr_t) NULL)));
+          return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oDivide, make_pair(NodeID_2, (NodeID) NULL)));
-          return (make_pair(1, NodeID_2));
+          List_of_Op_RHS.push_back(make_pair(oDivide, make_pair(expr_t_2, (expr_t) NULL)));
+          return (make_pair(1, expr_t_2));
         }
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oDivide:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddDivide(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddDivide(expr_t_1, expr_t_2)));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oDivide, make_pair((NodeID) NULL, NodeID_1)));
-          return (make_pair(1, NodeID_1));
+          List_of_Op_RHS.push_back(make_pair(oDivide, make_pair((expr_t) NULL, expr_t_1)));
+          return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oTimes, make_pair(NodeID_2, (NodeID) NULL)));
-          return (make_pair(1, NodeID_2));
+          List_of_Op_RHS.push_back(make_pair(oTimes, make_pair(expr_t_2, (expr_t) NULL)));
+          return (make_pair(1, expr_t_2));
         }
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oPower:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddPower(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddPower(expr_t_1, expr_t_2)));
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.push_back(make_pair(oPower, make_pair(datatree.AddDivide(datatree.One, NodeID_2), (NodeID) NULL)));
-          return (make_pair(1, (NodeID) NULL));
+          List_of_Op_RHS.push_back(make_pair(oPower, make_pair(datatree.AddDivide(datatree.One, expr_t_2), (expr_t) NULL)));
+          return (make_pair(1, (expr_t) NULL));
         }
       break;
     case oEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         {
           return (make_pair(0,
-                            datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getID(eEndogenous, var_endo), 0), datatree.AddMinus(NodeID_2, NodeID_1))
+                            datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getID(eEndogenous, var_endo), 0), datatree.AddMinus(expr_t_2, expr_t_1))
                             ));
         }
       else if (is_endogenous_present_1 && is_endogenous_present_2)
@@ -2696,79 +2696,79 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, Node
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
           return (make_pair(0,
-                            datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getID(eEndogenous, var_endo), 0), /*datatree.AddUMinus(NodeID_1)*/ NodeID_1)
+                            datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getID(eEndogenous, var_endo), 0), /*datatree.AddUMinus(expr_t_1)*/ expr_t_1)
                             ));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
           return (make_pair(0,
-                            datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getID(eEndogenous, var_endo), 0), NodeID_2)
+                            datatree.AddEqual(datatree.AddVariable(datatree.symbol_table.getID(eEndogenous, var_endo), 0), expr_t_2)
                             ));
         }
       break;
     case oMax:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddMax(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddMax(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oMin:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddMin(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddMin(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oLess:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddLess(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddLess(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oGreater:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddGreater(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddGreater(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oLessEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddLessEqual(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddLessEqual(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oGreaterEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddGreaterEqual(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddGreaterEqual(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oEqualEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddEqualEqual(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddEqualEqual(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     case oDifferent:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        return (make_pair(0, datatree.AddDifferent(NodeID_1, NodeID_2)));
+        return (make_pair(0, datatree.AddDifferent(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (NodeID) NULL));
+        return (make_pair(1, (expr_t) NULL));
       break;
     }
   // Suppress GCC warning
   exit(EXIT_FAILURE);
 }
 
-NodeID
-BinaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables)
+expr_t
+BinaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
-  NodeID darg1 = arg1->getChainRuleDerivative(deriv_id, recursive_variables);
-  NodeID darg2 = arg2->getChainRuleDerivative(deriv_id, recursive_variables);
+  expr_t darg1 = arg1->getChainRuleDerivative(deriv_id, recursive_variables);
+  expr_t darg2 = arg2->getChainRuleDerivative(deriv_id, recursive_variables);
   return composeDerivatives(darg1, darg2);
 }
 
-NodeID
-BinaryOpNode::buildSimilarBinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, DataTree &alt_datatree) const
+expr_t
+BinaryOpNode::buildSimilarBinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, DataTree &alt_datatree) const
 {
   switch (op_code)
     {
@@ -2805,11 +2805,11 @@ BinaryOpNode::buildSimilarBinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, DataTre
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 BinaryOpNode::toStatic(DataTree &static_datatree) const
 {
-  NodeID sarg1 = arg1->toStatic(static_datatree);
-  NodeID sarg2 = arg2->toStatic(static_datatree);
+  expr_t sarg1 = arg1->toStatic(static_datatree);
+  expr_t sarg2 = arg2->toStatic(static_datatree);
   return buildSimilarBinaryOpNode(sarg1, sarg2, static_datatree);
 }
 
@@ -2837,26 +2837,26 @@ BinaryOpNode::maxExoLag() const
   return max(arg1->maxExoLag(), arg2->maxExoLag());
 }
 
-NodeID
+expr_t
 BinaryOpNode::decreaseLeadsLags(int n) const
 {
-  NodeID arg1subst = arg1->decreaseLeadsLags(n);
-  NodeID arg2subst = arg2->decreaseLeadsLags(n);
+  expr_t arg1subst = arg1->decreaseLeadsLags(n);
+  expr_t arg2subst = arg2->decreaseLeadsLags(n);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
-NodeID
+expr_t
 BinaryOpNode::decreaseLeadsLagsPredeterminedVariables() const
 {
-  NodeID arg1subst = arg1->decreaseLeadsLagsPredeterminedVariables();
-  NodeID arg2subst = arg2->decreaseLeadsLagsPredeterminedVariables();
+  expr_t arg1subst = arg1->decreaseLeadsLagsPredeterminedVariables();
+  expr_t arg2subst = arg2->decreaseLeadsLagsPredeterminedVariables();
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
-NodeID
+expr_t
 BinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst, arg2subst;
+  expr_t arg1subst, arg2subst;
   int maxendolead1 = arg1->maxEndoLead(), maxendolead2 = arg2->maxEndoLead();
 
   if (maxendolead1 < 2 && maxendolead2 < 2)
@@ -2889,18 +2889,18 @@ BinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vecto
     }
 }
 
-NodeID
+expr_t
 BinaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst = arg1->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
-  NodeID arg2subst = arg2->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
-NodeID
+expr_t
 BinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst, arg2subst;
+  expr_t arg1subst, arg2subst;
   int maxexolead1 = arg1->maxExoLead(), maxexolead2 = arg2->maxExoLead();
 
   if (maxexolead1 < 1 && maxexolead2 < 1)
@@ -2933,19 +2933,19 @@ BinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode
     }
 }
 
-NodeID
+expr_t
 BinaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst = arg1->substituteExoLag(subst_table, neweqs);
-  NodeID arg2subst = arg2->substituteExoLag(subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteExoLag(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteExoLag(subst_table, neweqs);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
-NodeID
+expr_t
 BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
-  NodeID arg1subst = arg1->substituteExpectation(subst_table, neweqs, partial_information_model);
-  NodeID arg2subst = arg2->substituteExpectation(subst_table, neweqs, partial_information_model);
+  expr_t arg1subst = arg1->substituteExpectation(subst_table, neweqs, partial_information_model);
+  expr_t arg2subst = arg2->substituteExpectation(subst_table, neweqs, partial_information_model);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -2961,8 +2961,8 @@ BinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
   return false;
 }
 
-TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
-                             TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) :
+TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
+                             TrinaryOpcode op_code_arg, const expr_t arg2_arg, const expr_t arg3_arg) :
   ExprNode(datatree_arg),
   arg1(arg1_arg),
   arg2(arg2_arg),
@@ -2999,17 +2999,17 @@ TrinaryOpNode::prepareForDerivation()
             inserter(non_null_derivatives, non_null_derivatives.begin()));
 }
 
-NodeID
-TrinaryOpNode::composeDerivatives(NodeID darg1, NodeID darg2, NodeID darg3)
+expr_t
+TrinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2, expr_t darg3)
 {
 
-  NodeID t11, t12, t13, t14, t15;
+  expr_t t11, t12, t13, t14, t15;
 
   switch (op_code)
     {
     case oNormcdf:
       // normal pdf is inlined in the tree
-      NodeID y;
+      expr_t y;
       // sqrt(2*pi)
       t14 = datatree.AddSqrt(datatree.AddTimes(datatree.Two, datatree.Pi));
       // x - mu
@@ -3067,12 +3067,12 @@ TrinaryOpNode::composeDerivatives(NodeID darg1, NodeID darg2, NodeID darg3)
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::computeDerivative(int deriv_id)
 {
-  NodeID darg1 = arg1->getDerivative(deriv_id);
-  NodeID darg2 = arg2->getDerivative(deriv_id);
-  NodeID darg3 = arg3->getDerivative(deriv_id);
+  expr_t darg1 = arg1->getDerivative(deriv_id);
+  expr_t darg2 = arg2->getDerivative(deriv_id);
+  expr_t darg3 = arg3->getDerivative(deriv_id);
   return composeDerivatives(darg1, darg2, darg3);
 }
 
@@ -3126,12 +3126,12 @@ TrinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) co
 }
 
 void
-TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+TrinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
                                      bool is_matlab) const
 {
-  NodeID this2 = const_cast<TrinaryOpNode *>(this);
-  map<NodeID, int>::iterator it = reference_count.find(this2);
+  expr_t this2 = const_cast<TrinaryOpNode *>(this);
+  map<expr_t, int>::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       // If this node has never been encountered, set its ref count to one,
@@ -3152,15 +3152,15 @@ TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
 }
 
 void
-TrinaryOpNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+TrinaryOpNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector<vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const
 {
-  NodeID this2 = const_cast<TrinaryOpNode *>(this);
-  map<NodeID, int>::iterator it = reference_count.find(this2);
+  expr_t this2 = const_cast<TrinaryOpNode *>(this);
+  map<expr_t, int>::iterator it = reference_count.find(this2);
   if (it == reference_count.end())
     {
       reference_count[this2] = 1;
@@ -3332,35 +3332,35 @@ TrinaryOpNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &resul
   arg3->collectVariables(type_arg, result);
 }
 
-pair<int, NodeID>
-TrinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const
+pair<int, expr_t>
+TrinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
-  pair<int, NodeID> res = arg1->normalizeEquation(var_endo, List_of_Op_RHS);
+  pair<int, expr_t> res = arg1->normalizeEquation(var_endo, List_of_Op_RHS);
   bool is_endogenous_present_1 = res.first;
-  NodeID NodeID_1 = res.second;
+  expr_t expr_t_1 = res.second;
   res = arg2->normalizeEquation(var_endo, List_of_Op_RHS);
   bool is_endogenous_present_2 = res.first;
-  NodeID NodeID_2 = res.second;
+  expr_t expr_t_2 = res.second;
   res = arg3->normalizeEquation(var_endo, List_of_Op_RHS);
   bool is_endogenous_present_3 = res.first;
-  NodeID NodeID_3 = res.second;
+  expr_t expr_t_3 = res.second;
   if (!is_endogenous_present_1 && !is_endogenous_present_2 && !is_endogenous_present_3)
-    return (make_pair(0, datatree.AddNormcdf(NodeID_1, NodeID_2, NodeID_3)));
+    return (make_pair(0, datatree.AddNormcdf(expr_t_1, expr_t_2, expr_t_3)));
   else
-    return (make_pair(1, (NodeID) NULL));
+    return (make_pair(1, (expr_t) NULL));
 }
 
-NodeID
-TrinaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables)
+expr_t
+TrinaryOpNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
-  NodeID darg1 = arg1->getChainRuleDerivative(deriv_id, recursive_variables);
-  NodeID darg2 = arg2->getChainRuleDerivative(deriv_id, recursive_variables);
-  NodeID darg3 = arg3->getChainRuleDerivative(deriv_id, recursive_variables);
+  expr_t darg1 = arg1->getChainRuleDerivative(deriv_id, recursive_variables);
+  expr_t darg2 = arg2->getChainRuleDerivative(deriv_id, recursive_variables);
+  expr_t darg3 = arg3->getChainRuleDerivative(deriv_id, recursive_variables);
   return composeDerivatives(darg1, darg2, darg3);
 }
 
-NodeID
-TrinaryOpNode::buildSimilarTrinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, NodeID alt_arg3, DataTree &alt_datatree) const
+expr_t
+TrinaryOpNode::buildSimilarTrinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, expr_t alt_arg3, DataTree &alt_datatree) const
 {
   switch (op_code)
     {
@@ -3373,12 +3373,12 @@ TrinaryOpNode::buildSimilarTrinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, NodeI
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::toStatic(DataTree &static_datatree) const
 {
-  NodeID sarg1 = arg1->toStatic(static_datatree);
-  NodeID sarg2 = arg2->toStatic(static_datatree);
-  NodeID sarg3 = arg3->toStatic(static_datatree);
+  expr_t sarg1 = arg1->toStatic(static_datatree);
+  expr_t sarg2 = arg2->toStatic(static_datatree);
+  expr_t sarg3 = arg3->toStatic(static_datatree);
   return buildSimilarTrinaryOpNode(sarg1, sarg2, sarg3, static_datatree);
 }
 
@@ -3406,25 +3406,25 @@ TrinaryOpNode::maxExoLag() const
   return max(arg1->maxExoLag(), max(arg2->maxExoLag(), arg3->maxExoLag()));
 }
 
-NodeID
+expr_t
 TrinaryOpNode::decreaseLeadsLags(int n) const
 {
-  NodeID arg1subst = arg1->decreaseLeadsLags(n);
-  NodeID arg2subst = arg2->decreaseLeadsLags(n);
-  NodeID arg3subst = arg3->decreaseLeadsLags(n);
+  expr_t arg1subst = arg1->decreaseLeadsLags(n);
+  expr_t arg2subst = arg2->decreaseLeadsLags(n);
+  expr_t arg3subst = arg3->decreaseLeadsLags(n);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::decreaseLeadsLagsPredeterminedVariables() const
 {
-  NodeID arg1subst = arg1->decreaseLeadsLagsPredeterminedVariables();
-  NodeID arg2subst = arg2->decreaseLeadsLagsPredeterminedVariables();
-  NodeID arg3subst = arg3->decreaseLeadsLagsPredeterminedVariables();
+  expr_t arg1subst = arg1->decreaseLeadsLagsPredeterminedVariables();
+  expr_t arg2subst = arg2->decreaseLeadsLagsPredeterminedVariables();
+  expr_t arg3subst = arg3->decreaseLeadsLagsPredeterminedVariables();
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   if (maxEndoLead() < 2)
@@ -3433,16 +3433,16 @@ TrinaryOpNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vect
     return createEndoLeadAuxiliaryVarForMyself(subst_table, neweqs);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst = arg1->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
-  NodeID arg2subst = arg2->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
-  NodeID arg3subst = arg3->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteEndoLagGreaterThanTwo(subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
   if (maxExoLead() == 0)
@@ -3451,21 +3451,21 @@ TrinaryOpNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode
     return createExoLeadAuxiliaryVarForMyself(subst_table, neweqs);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  NodeID arg1subst = arg1->substituteExoLag(subst_table, neweqs);
-  NodeID arg2subst = arg2->substituteExoLag(subst_table, neweqs);
-  NodeID arg3subst = arg3->substituteExoLag(subst_table, neweqs);
+  expr_t arg1subst = arg1->substituteExoLag(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteExoLag(subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteExoLag(subst_table, neweqs);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
-NodeID
+expr_t
 TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
-  NodeID arg1subst = arg1->substituteExpectation(subst_table, neweqs, partial_information_model);
-  NodeID arg2subst = arg2->substituteExpectation(subst_table, neweqs, partial_information_model);
-  NodeID arg3subst = arg3->substituteExpectation(subst_table, neweqs, partial_information_model);
+  expr_t arg1subst = arg1->substituteExpectation(subst_table, neweqs, partial_information_model);
+  expr_t arg2subst = arg2->substituteExpectation(subst_table, neweqs, partial_information_model);
+  expr_t arg3subst = arg3->substituteExpectation(subst_table, neweqs, partial_information_model);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -3483,7 +3483,7 @@ TrinaryOpNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int l
 
 ExternalFunctionNode::ExternalFunctionNode(DataTree &datatree_arg,
                                          int symb_id_arg,
-                                         const vector<NodeID> &arguments_arg) :
+                                         const vector<expr_t> &arguments_arg) :
   ExprNode(datatree_arg),
   symb_id(symb_id_arg),
   arguments(arguments_arg)
@@ -3498,7 +3498,7 @@ ExternalFunctionNode::prepareForDerivation()
   if (preparedForDerivation)
     return;
 
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     (*it)->prepareForDerivation();
 
   non_null_derivatives = arguments.at(0)->non_null_derivatives;
@@ -3512,40 +3512,40 @@ ExternalFunctionNode::prepareForDerivation()
   preparedForDerivation = true;
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::computeDerivative(int deriv_id)
 {
   assert(datatree.external_functions_table.getNargs(symb_id) > 0);
-  vector<NodeID> dargs;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> dargs;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     dargs.push_back((*it)->getDerivative(deriv_id));
   return composeDerivatives(dargs);
 }
 
-NodeID
-ExternalFunctionNode::composeDerivatives(const vector<NodeID> &dargs)
+expr_t
+ExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 {
-  vector<NodeID> dNodes;
+  vector<expr_t> dNodes;
   for (int i = 0; i < (int)dargs.size(); i++)
     if (dargs.at(i) != 0)
       dNodes.push_back(datatree.AddTimes(dargs.at(i),
                                          datatree.AddFirstDerivExternalFunctionNode(symb_id, arguments, i+1)));
 
-  NodeID theDeriv = datatree.Zero;
-  for (vector<NodeID>::const_iterator it = dNodes.begin(); it != dNodes.end(); it++)
+  expr_t theDeriv = datatree.Zero;
+  for (vector<expr_t>::const_iterator it = dNodes.begin(); it != dNodes.end(); it++)
     theDeriv = datatree.AddPlus(theDeriv, *it);
   return theDeriv;
 }
 
-NodeID
-ExternalFunctionNode::getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables)
+expr_t
+ExternalFunctionNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
 {
   cerr << "ExternalFunctionNode::getChainRuleDerivative: operation impossible!" << endl;
   exit(EXIT_FAILURE);
 }
 
 void
-ExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                            temporary_terms_t &temporary_terms,
                                            bool is_matlab) const
 {
@@ -3557,7 +3557,7 @@ ExternalFunctionNode::writeExternalFunctionArguments(ostream &output, ExprNodeOu
                                                      const temporary_terms_t &temporary_terms,
                                                      deriv_node_temp_terms_t &tef_terms) const
 {
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     {
       if (it != arguments.begin())
@@ -3601,7 +3601,7 @@ ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutpu
   int first_deriv_symb_id = datatree.external_functions_table.getFirstDerivSymbID(symb_id);
   assert(first_deriv_symb_id != eExtFunSetButNoNameProvided);
 
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     (*it)->writeExternalFunctionOutput(output, output_type, temporary_terms, tef_terms);
 
@@ -3627,9 +3627,9 @@ ExternalFunctionNode::writeExternalFunctionOutput(ostream &output, ExprNodeOutpu
 }
 
 void
-ExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+ExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                            temporary_terms_t &temporary_terms,
-                                           map<NodeID, pair<int, int> > &first_occurence,
+                                           map<expr_t, pair<int, int> > &first_occurence,
                                            int Curr_block,
                                            vector< vector<temporary_terms_t> > &v_temporary_terms,
                                            int equation) const
@@ -3641,7 +3641,7 @@ ExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
 void
 ExternalFunctionNode::collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const
 {
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     (*it)->collectVariables(type_arg, result);
 }
@@ -3671,30 +3671,30 @@ ExternalFunctionNode::compile(ostream &CompileCode, bool lhs_rhs, const temporar
   exit(EXIT_FAILURE);
 }
 
-pair<int, NodeID>
-ExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const
+pair<int, expr_t>
+ExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const
 {
-  vector<pair<bool, NodeID> > V_arguments;
-  vector<NodeID> V_NodeID;
+  vector<pair<bool, expr_t> > V_arguments;
+  vector<expr_t> V_expr_t;
   bool present = false;
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     {
       V_arguments.push_back((*it)->normalizeEquation(var_endo, List_of_Op_RHS));
       present = present || V_arguments[V_arguments.size()-1].first;
-      V_NodeID.push_back(V_arguments[V_arguments.size()-1].second);
+      V_expr_t.push_back(V_arguments[V_arguments.size()-1].second);
     }
   if (!present)
-    return (make_pair(0, datatree.AddExternalFunction(symb_id, V_NodeID)));
+    return (make_pair(0, datatree.AddExternalFunction(symb_id, V_expr_t)));
   else
-    return (make_pair(1, (NodeID) NULL));
+    return (make_pair(1, (expr_t) NULL));
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::toStatic(DataTree &static_datatree) const
 {
-  vector<NodeID> static_arguments;
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  vector<expr_t> static_arguments;
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     static_arguments.push_back((*it)->toStatic(static_datatree));
   return static_datatree.AddExternalFunction(symb_id, static_arguments);
@@ -3704,7 +3704,7 @@ int
 ExternalFunctionNode::maxEndoLead() const
 {
   int val = 0;
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     val = max(val, (*it)->maxEndoLead());
   return val;
@@ -3714,7 +3714,7 @@ int
 ExternalFunctionNode::maxExoLead() const
 {
   int val = 0;
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     val = max(val, (*it)->maxExoLead());
   return val;
@@ -3724,7 +3724,7 @@ int
 ExternalFunctionNode::maxEndoLag() const
 {
   int val = 0;
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     val = max(val, (*it)->maxEndoLag());
   return val;
@@ -3734,77 +3734,77 @@ int
 ExternalFunctionNode::maxExoLag() const
 {
   int val = 0;
-  for (vector<NodeID>::const_iterator it = arguments.begin();
+  for (vector<expr_t>::const_iterator it = arguments.begin();
        it != arguments.end(); it++)
     val = max(val, (*it)->maxExoLag());
   return val;
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::decreaseLeadsLags(int n) const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->decreaseLeadsLags(n));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::decreaseLeadsLagsPredeterminedVariables() const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->decreaseLeadsLagsPredeterminedVariables());
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->substituteEndoLeadGreaterThanTwo(subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->substituteEndoLagGreaterThanTwo(subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->substituteExoLead(subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->substituteExoLag(subst_table, neweqs));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
+expr_t
 ExternalFunctionNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
-  vector<NodeID> arguments_subst;
-  for (vector<NodeID>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     arguments_subst.push_back((*it)->substituteExpectation(subst_table, neweqs, partial_information_model));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
-NodeID
-ExternalFunctionNode::buildSimilarExternalFunctionNode(vector<NodeID> &alt_args, DataTree &alt_datatree) const
+expr_t
+ExternalFunctionNode::buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const
 {
   return alt_datatree.AddExternalFunction(symb_id, alt_args);
 }
@@ -3841,7 +3841,7 @@ ExternalFunctionNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id
 
 FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatree_arg,
                                                                int top_level_symb_id_arg,
-                                                               const vector<NodeID> &arguments_arg,
+                                                               const vector<expr_t> &arguments_arg,
                                                                int inputIndex_arg) :
   ExternalFunctionNode(datatree_arg, top_level_symb_id_arg, arguments_arg),
   inputIndex(inputIndex_arg)
@@ -3851,7 +3851,7 @@ FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatre
 }
 
 void
-FirstDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                                       temporary_terms_t &temporary_terms,
                                                       bool is_matlab) const
 {
@@ -3859,9 +3859,9 @@ FirstDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &referenc
 }
 
 void
-FirstDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                                       temporary_terms_t &temporary_terms,
-                                                      map<NodeID, pair<int, int> > &first_occurence,
+                                                      map<expr_t, pair<int, int> > &first_occurence,
                                                       int Curr_block,
                                                       vector< vector<temporary_terms_t> > &v_temporary_terms,
                                                       int equation) const
@@ -3870,16 +3870,16 @@ FirstDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &referenc
   exit(EXIT_FAILURE);
 }
 
-NodeID
-FirstDerivExternalFunctionNode::composeDerivatives(const vector<NodeID> &dargs)
+expr_t
+FirstDerivExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 {
-  vector<NodeID> dNodes;
+  vector<expr_t> dNodes;
   for (int i = 0; i < (int)dargs.size(); i++)
     if (dargs.at(i) != 0)
       dNodes.push_back(datatree.AddTimes(dargs.at(i),
                                          datatree.AddSecondDerivExternalFunctionNode(symb_id, arguments, inputIndex, i+1)));
-  NodeID theDeriv = datatree.Zero;
-  for (vector<NodeID>::const_iterator it = dNodes.begin(); it != dNodes.end(); it++)
+  expr_t theDeriv = datatree.Zero;
+  for (vector<expr_t>::const_iterator it = dNodes.begin(); it != dNodes.end(); it++)
     theDeriv = datatree.AddPlus(theDeriv, *it);
   return theDeriv;
 }
@@ -3947,7 +3947,7 @@ FirstDerivExternalFunctionNode::writeExternalFunctionOutput(ostream &output, Exp
 
 SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datatree_arg,
                                                                  int top_level_symb_id_arg,
-                                                                 const vector<NodeID> &arguments_arg,
+                                                                 const vector<expr_t> &arguments_arg,
                                                                  int inputIndex1_arg,
                                                                  int inputIndex2_arg) :
   ExternalFunctionNode(datatree_arg, top_level_symb_id_arg, arguments_arg),
@@ -3959,7 +3959,7 @@ SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datat
 }
 
 void
-SecondDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                                        temporary_terms_t &temporary_terms,
                                                        bool is_matlab) const
 {
@@ -3967,9 +3967,9 @@ SecondDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &referen
 }
 
 void
-SecondDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &reference_count,
+SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
                                                       temporary_terms_t &temporary_terms,
-                                                      map<NodeID, pair<int, int> > &first_occurence,
+                                                      map<expr_t, pair<int, int> > &first_occurence,
                                                       int Curr_block,
                                                       vector< vector<temporary_terms_t> > &v_temporary_terms,
                                                       int equation) const
@@ -3978,7 +3978,7 @@ SecondDerivExternalFunctionNode::computeTemporaryTerms(map<NodeID, int> &referen
   exit(EXIT_FAILURE);
 }
 
-NodeID
+expr_t
 SecondDerivExternalFunctionNode::computeDerivative(int deriv_id)
 {
   cerr << "ERROR: SecondDerivExternalFunctionNode::computeDerivative(). Not implemented" << endl;
diff --git a/ExprNode.hh b/ExprNode.hh
index 0f846ebbb5befdf9bcb273f2d1d98e606941e91d..79773a6be93967a8296c8fdc5e1c512722ae8d87 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -35,7 +35,7 @@ class DataTree;
 class VariableNode;
 class BinaryOpNode;
 
-typedef class ExprNode *NodeID;
+typedef class ExprNode *expr_t;
 
 struct Model_Block;
 
@@ -43,7 +43,7 @@ struct ExprNodeLess;
 
 //! Type for set of temporary terms
 /*! They are ordered by index number thanks to ExprNodeLess */
-typedef set<NodeID, ExprNodeLess> temporary_terms_t;
+typedef set<expr_t, ExprNodeLess> temporary_terms_t;
 
 //! set of temporary terms used in a block
 typedef set<int> temporary_terms_inuse_t;
@@ -55,7 +55,7 @@ typedef map<int, int> map_idx_t;
 typedef map<int, double> eval_context_t;
 
 //! Type for tracking first/second derivative functions that have already been written as temporary terms
-typedef map<pair<int, vector<NodeID> >, int> deriv_node_temp_terms_t;
+typedef map<pair<int, vector<expr_t> >, int> deriv_node_temp_terms_t;
 
 //! Possible types of output when writing ExprNode(s)
 enum ExprNodeOutputType
@@ -123,7 +123,7 @@ class ExprNode
 private:
   //! Computes derivative w.r. to a derivation ID (but doesn't store it in derivatives map)
   /*! You shoud use getDerivative() to get the benefit of symbolic a priori and of caching */
-  virtual NodeID computeDerivative(int deriv_id) = 0;
+  virtual expr_t computeDerivative(int deriv_id) = 0;
 
 protected:
   //! Reference to the enclosing DataTree
@@ -139,7 +139,7 @@ protected:
   set<int> non_null_derivatives;
 
   //! Used for caching of first order derivatives (when non-null)
-  map<int, NodeID> derivatives;
+  map<int, expr_t> derivatives;
 
   //! Cost of computing current node
   /*! Nodes included in temporary_terms are considered having a null cost */
@@ -155,14 +155,14 @@ public:
   //! Returns derivative w.r. to derivation ID
   /*! Uses a symbolic a priori to pre-detect null derivatives, and caches the result for other derivatives (to avoid computing it several times)
     For an equal node, returns the derivative of lhs minus rhs */
-  NodeID getDerivative(int deriv_id);
+  expr_t getDerivative(int deriv_id);
 
   //! Computes derivatives by applying the chain rule for some variables
   /*!
     \param deriv_id The derivation ID with respect to which we are derivating
     \param recursive_variables Contains the derivation ID for which chain rules must be applied. Keys are derivation IDs, values are equations of the form x=f(y) where x is the key variable and x doesn't appear in y
   */
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables) = 0;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables) = 0;
 
   //! Returns precedence of node
   /*! Equals 100 for constants, variables, unary ops, and temporary terms */
@@ -170,7 +170,7 @@ public:
 
   //! Fills temporary_terms set, using reference counts
   /*! A node will be marked as a temporary term if it is referenced at least two times (i.e. has at least two parents), and has a computing cost (multiplied by reference count) greater to datatree.min_cost */
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
 
   //! Writes output of node, using a Txxx notation for nodes in temporary_terms, and specifiying the set of already written external functions
   /*!
@@ -228,9 +228,9 @@ public:
   virtual void collectModelLocalVariables(set<int> &result) const;
 
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const = 0;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
@@ -247,9 +247,9 @@ public:
     This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
     adds the result in the static_datatree argument (and not in the original datatree), and returns it.
   */
-  virtual NodeID toStatic(DataTree &static_datatree) const = 0;
+  virtual expr_t toStatic(DataTree &static_datatree) const = 0;
   //! Try to normalize an equation linear in its endogenous variable
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > > &List_of_Op_RHS) const = 0;
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const = 0;
 
   //! Returns the maximum lead of endogenous in this expression
   /*! Always returns a non-negative value */
@@ -273,7 +273,7 @@ public:
     \param[in] n The number of lags by which to shift
     \return The same expression except that leads/lags have been shifted backwards
   */
-  virtual NodeID decreaseLeadsLags(int n) const = 0;
+  virtual expr_t decreaseLeadsLags(int n) const = 0;
 
   //! Type for the substitution map used in the process of creating auxiliary vars for leads >= 2
   typedef map<const ExprNode *, const VariableNode *> subst_table_t;
@@ -309,27 +309,27 @@ public:
 
     \return A new equivalent expression where sub-expressions with max endo lead >= 2 have been replaced by auxiliary variables
   */
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
   //! Constructs a new expression where endo variables with max endo lag >= 2 have been replaced by auxiliary variables
   /*!
     \param[in,out] subst_table Map used to store expressions that have already be substituted and their corresponding variable, in order to avoid creating two auxiliary variables for the same sub-expr.
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
   //! Constructs a new expression where exogenous variables with a lead have been replaced by auxiliary variables
   /*!
     \param[in,out] subst_table Map used to store expressions that have already be substituted and their corresponding variable, in order to avoid creating two auxiliary variables for the same sub-expr.
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
   //! Constructs a new expression where exogenous variables with a lag have been replaced by auxiliary variables
   /*!
     \param[in,out] subst_table Map used to store expressions that have already be substituted and their corresponding variable, in order to avoid creating two auxiliary variables for the same sub-expr.
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+  virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
   //! Constructs a new expression where the expectation operator has been replaced by auxiliary variables
   /*!
@@ -337,9 +337,9 @@ public:
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
     \param[in] partial_information_model Are we substituting in a partial information model?
   */
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const = 0;
+  virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const = 0;
 
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const = 0;
+  virtual expr_t decreaseLeadsLagsPredeterminedVariables() const = 0;
 
   //! Return true if the nodeID is a numerical constant equal to value and false otherwise
   /*!
@@ -361,7 +361,7 @@ public:
 struct ExprNodeLess
 {
   bool
-  operator()(NodeID arg1, NodeID arg2) const
+  operator()(expr_t arg1, expr_t arg2) const
   {
     return arg1->idx < arg2->idx;
   }
@@ -374,7 +374,7 @@ class NumConstNode : public ExprNode
 private:
   //! Id from numerical constants table
   const int id;
-  virtual NodeID computeDerivative(int deriv_id);
+  virtual expr_t computeDerivative(int deriv_id);
 public:
   NumConstNode(DataTree &datatree_arg, int id_arg);
   int
@@ -388,20 +388,20 @@ public:
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
   virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
-  virtual NodeID toStatic(DataTree &static_datatree) const;
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
+  virtual expr_t toStatic(DataTree &static_datatree) const;
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
   virtual int maxEndoLead() const;
   virtual int maxExoLead() const;
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
-  virtual NodeID decreaseLeadsLags(int n) const;
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  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 decreaseLeadsLagsPredeterminedVariables() const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
 };
@@ -415,41 +415,41 @@ private:
   const SymbolType type;
   //! A positive value is a lead, a negative is a lag
   const int lag;
-  virtual NodeID computeDerivative(int deriv_id);
+  virtual expr_t computeDerivative(int deriv_id);
 public:
   VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg);
   virtual void prepareForDerivation();
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
   virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
-  virtual NodeID toStatic(DataTree &static_datatree) const;
+  virtual expr_t toStatic(DataTree &static_datatree) const;
   int
   get_symb_id() const
   {
     return symb_id;
   };
   int get_lag() const { return lag; };
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
   virtual int maxEndoLead() const;
   virtual int maxExoLead() const;
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
-  virtual NodeID decreaseLeadsLags(int n) const;
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  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 decreaseLeadsLagsPredeterminedVariables() const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
 };
@@ -458,27 +458,27 @@ public:
 class UnaryOpNode : public ExprNode
 {
 private:
-  const NodeID arg;
+  const expr_t arg;
   //! Stores the information set. Only used for expectation operator
   const int expectation_information_set;
   //! Stores the information set name. Only used for expectation operator
   const string expectation_information_set_name;
   const UnaryOpcode op_code;
-  virtual NodeID computeDerivative(int deriv_id);
+  virtual expr_t computeDerivative(int deriv_id);
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
   //! Returns the derivative of this node if darg is the derivative of the argument
-  NodeID composeDerivatives(NodeID darg);
+  expr_t composeDerivatives(expr_t darg);
 public:
-  UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg);
+  UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, const int expectation_information_set_arg, const string &expectation_information_set_name_arg);
   virtual void prepareForDerivation();
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
@@ -488,7 +488,7 @@ public:
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
   virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
   //! Returns operand
-  NodeID
+  expr_t
   get_arg() const
   {
     return (arg);
@@ -499,22 +499,22 @@ public:
   {
     return (op_code);
   };
-  virtual NodeID toStatic(DataTree &static_datatree) const;
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
+  virtual expr_t toStatic(DataTree &static_datatree) const;
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
   virtual int maxEndoLead() const;
   virtual int maxExoLead() const;
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
-  virtual NodeID decreaseLeadsLags(int n) const;
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   //! Creates another UnaryOpNode with the same opcode, but with a possibly different datatree and argument
-  NodeID buildSimilarUnaryOpNode(NodeID alt_arg, DataTree &alt_datatree) const;
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
+  expr_t buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  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 decreaseLeadsLagsPredeterminedVariables() const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
 };
@@ -523,25 +523,25 @@ public:
 class BinaryOpNode : public ExprNode
 {
 private:
-  const NodeID arg1, arg2;
+  const expr_t arg1, arg2;
   const BinaryOpcode op_code;
-  virtual NodeID computeDerivative(int deriv_id);
+  virtual expr_t computeDerivative(int deriv_id);
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
   //! Returns the derivative of this node if darg1 and darg2 are the derivatives of the arguments
-  NodeID composeDerivatives(NodeID darg1, NodeID darg2);
+  expr_t composeDerivatives(expr_t darg1, expr_t darg2);
 public:
-  BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
-               BinaryOpcode op_code_arg, const NodeID arg2_arg);
+  BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
+               BinaryOpcode op_code_arg, const expr_t arg2_arg);
   virtual void prepareForDerivation();
   virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
@@ -550,15 +550,15 @@ public:
   static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
   virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
-  virtual NodeID Compute_RHS(NodeID arg1, NodeID arg2, int op, int op_type) const;
+  virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
   //! Returns first operand
-  NodeID
+  expr_t
   get_arg1() const
   {
     return (arg1);
   };
   //! Returns second operand
-  NodeID
+  expr_t
   get_arg2() const
   {
     return (arg2);
@@ -569,22 +569,22 @@ public:
   {
     return (op_code);
   };
-  virtual NodeID toStatic(DataTree &static_datatree) const;
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
+  virtual expr_t toStatic(DataTree &static_datatree) const;
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
   virtual int maxEndoLead() const;
   virtual int maxExoLead() const;
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
-  virtual NodeID decreaseLeadsLags(int n) const;
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   //! Creates another BinaryOpNode with the same opcode, but with a possibly different datatree and arguments
-  NodeID buildSimilarBinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, DataTree &alt_datatree) const;
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
+  expr_t buildSimilarBinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, DataTree &alt_datatree) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  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 decreaseLeadsLagsPredeterminedVariables() const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
 };
@@ -594,25 +594,25 @@ class TrinaryOpNode : public ExprNode
 {
   friend class ModelTree;
 private:
-  const NodeID arg1, arg2, arg3;
+  const expr_t arg1, arg2, arg3;
   const TrinaryOpcode op_code;
-  virtual NodeID computeDerivative(int deriv_id);
+  virtual expr_t computeDerivative(int deriv_id);
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
   //! Returns the derivative of this node if darg1, darg2 and darg3 are the derivatives of the arguments
-  NodeID composeDerivatives(NodeID darg1, NodeID darg2, NodeID darg3);
+  expr_t composeDerivatives(expr_t darg1, expr_t darg2, expr_t darg3);
 public:
-  TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
-                TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg);
+  TrinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
+                TrinaryOpcode op_code_arg, const expr_t arg2_arg, const expr_t arg3_arg);
   virtual void prepareForDerivation();
   virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
@@ -621,22 +621,22 @@ public:
   static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
   virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
-  virtual NodeID toStatic(DataTree &static_datatree) const;
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
+  virtual expr_t toStatic(DataTree &static_datatree) const;
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
   virtual int maxEndoLead() const;
   virtual int maxExoLead() const;
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
-  virtual NodeID decreaseLeadsLags(int n) const;
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   //! Creates another TrinaryOpNode with the same opcode, but with a possibly different datatree and arguments
-  NodeID buildSimilarTrinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, NodeID alt_arg3, DataTree &alt_datatree) const;
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
+  expr_t buildSimilarTrinaryOpNode(expr_t alt_arg1, expr_t alt_arg2, expr_t alt_arg3, DataTree &alt_datatree) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  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 decreaseLeadsLagsPredeterminedVariables() const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
 };
@@ -645,15 +645,15 @@ public:
 class ExternalFunctionNode : public ExprNode
 {
 private:
-  virtual NodeID computeDerivative(int deriv_id);
-  virtual NodeID composeDerivatives(const vector<NodeID> &dargs);
+  virtual expr_t computeDerivative(int deriv_id);
+  virtual expr_t composeDerivatives(const vector<expr_t> &dargs);
 protected:
   //! Thrown when trying to access an unknown entry in external_function_node_map
   class UnknownFunctionNameAndArgs
   {
   };
   const int symb_id;
-  const vector<NodeID> arguments;
+  const vector<expr_t> arguments;
   //! Returns true if the given external function has been written as a temporary term
   bool alreadyWrittenAsTefTerm(int the_symb_id, deriv_node_temp_terms_t &tef_terms) const;
   //! Returns the index in the tef_terms map of this external function
@@ -662,16 +662,16 @@ protected:
   void writeExternalFunctionArguments(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
 public:
   ExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg,
-                      const vector<NodeID> &arguments_arg);
+                      const vector<expr_t> &arguments_arg);
   virtual void prepareForDerivation();
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type,
                                            const temporary_terms_t &temporary_terms,
                                            deriv_node_temp_terms_t &tef_terms) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
@@ -679,21 +679,21 @@ public:
   virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
   virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
   virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
-  virtual NodeID toStatic(DataTree &static_datatree) const;
-  virtual pair<int, NodeID> normalizeEquation(int symb_id_endo, vector<pair<int, pair<NodeID, NodeID> > >  &List_of_Op_RHS) const;
-  virtual NodeID getChainRuleDerivative(int deriv_id, const map<int, NodeID> &recursive_variables);
+  virtual expr_t toStatic(DataTree &static_datatree) const;
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
   virtual int maxEndoLead() const;
   virtual int maxExoLead() const;
   virtual int maxEndoLag() const;
   virtual int maxExoLag() const;
-  virtual NodeID decreaseLeadsLags(int n) const;
-  virtual NodeID substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
-  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual NodeID buildSimilarExternalFunctionNode(vector<NodeID> &alt_args, DataTree &alt_datatree) const;
-  virtual NodeID decreaseLeadsLagsPredeterminedVariables() const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  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 buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const;
+  virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
 };
@@ -702,16 +702,16 @@ class FirstDerivExternalFunctionNode : public ExternalFunctionNode
 {
 private:
   const int inputIndex;
-  virtual NodeID composeDerivatives(const vector<NodeID> &dargs);
+  virtual expr_t composeDerivatives(const vector<expr_t> &dargs);
 public:
   FirstDerivExternalFunctionNode(DataTree &datatree_arg,
                                  int top_level_symb_id_arg,
-                                 const vector<NodeID> &arguments_arg,
+                                 const vector<expr_t> &arguments_arg,
                                  int inputIndex_arg);
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
@@ -726,17 +726,17 @@ class SecondDerivExternalFunctionNode : public ExternalFunctionNode
 private:
   const int inputIndex1;
   const int inputIndex2;
-  virtual NodeID computeDerivative(int deriv_id);
+  virtual expr_t computeDerivative(int deriv_id);
 public:
   SecondDerivExternalFunctionNode(DataTree &datatree_arg,
                                   int top_level_symb_id_arg,
-                                  const vector<NodeID> &arguments_arg,
+                                  const vector<expr_t> &arguments_arg,
                                   int inputIndex1_arg,
                                   int inputIndex2_arg);
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
-  virtual void computeTemporaryTerms(map<NodeID, int> &reference_count,
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
-                                     map<NodeID, pair<int, int> > &first_occurence,
+                                     map<expr_t, pair<int, int> > &first_occurence,
                                      int Curr_block,
                                      vector< vector<temporary_terms_t> > &v_temporary_terms,
                                      int equation) const;
diff --git a/ModelTree.cc b/ModelTree.cc
index 7bac60122d4e0f620fc5521fa86f41cdf0db871b..005c47937261b50cfff568b2faecf2e99f6337d4 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -252,7 +252,7 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
       int deriv_id = it->first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          NodeID Id = it->second;
+          expr_t Id = it->second;
           int eq = it->first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           int var = symbol_table.getTypeSpecificID(symb);
@@ -411,9 +411,9 @@ ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, ve
 }
 
 equation_type_and_normalized_equation_t
-ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
+ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int> >, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const
 {
-  NodeID lhs, rhs;
+  expr_t lhs, rhs;
   BinaryOpNode *eq_node;
   EquationType Equation_Simulation_Type;
   equation_type_and_normalized_equation_t V_Equation_Simulation_Type(equations.size());
@@ -425,8 +425,8 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int> >, NodeI
       lhs = eq_node->get_arg1();
       rhs = eq_node->get_arg2();
       Equation_Simulation_Type = E_SOLVE;
-      map<pair<int, pair<int, int> >, NodeID>::const_iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
-      pair<bool, NodeID> res;
+      map<pair<int, pair<int, int> >, expr_t>::const_iterator derivative = first_order_endo_derivatives.find(make_pair(eq, make_pair(var, 0)));
+      pair<bool, expr_t> res;
       if (derivative != first_order_endo_derivatives.end())
         {
           set<pair<int, int> > result;
@@ -439,7 +439,7 @@ ModelTree::equationTypeDetermination(const map<pair<int, pair<int, int> >, NodeI
             }
           else
             {
-              vector<pair<int, pair<NodeID, NodeID> > > List_of_Op_RHS;
+              vector<pair<int, pair<expr_t, expr_t> > > List_of_Op_RHS;
               res =  equations[eq]->normalizeEquation(var, List_of_Op_RHS);
               if (mfs == 2)
                 {
@@ -802,7 +802,7 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
               int lag = it->second.first;
               if (lag == 0)
                 {
-                  NodeID Id = it->second.second;
+                  expr_t Id = it->second.second;
                   set<pair<int, int> > endogenous;
                   Id->collectEndogenous(endogenous);
                   if (endogenous.size() > 0)
@@ -824,7 +824,7 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = derivatives_block.begin(); it != derivatives_block.end(); it++)
             {
               int lag = it->second.first;
-              NodeID Id = it->second.second; //
+              expr_t Id = it->second.second; //
               set<pair<int, int> > endogenous;
               Id->collectEndogenous(endogenous);
               if (endogenous.size() > 0)
@@ -880,7 +880,7 @@ ModelTree::computeJacobian(const set<int> &vars)
        it != vars.end(); it++)
     for (int eq = 0; eq < (int) equations.size(); eq++)
       {
-        NodeID d1 = equations[eq]->getDerivative(*it);
+        expr_t d1 = equations[eq]->getDerivative(*it);
         if (d1 == Zero)
           continue;
         first_derivatives[make_pair(eq, *it)] = d1;
@@ -896,7 +896,7 @@ ModelTree::computeHessian(const set<int> &vars)
     {
       int eq = it->first.first;
       int var1 = it->first.second;
-      NodeID d1 = it->second;
+      expr_t d1 = it->second;
 
       // Store only second derivatives with var2 <= var1
       for (set<int>::const_iterator it2 = vars.begin();
@@ -906,7 +906,7 @@ ModelTree::computeHessian(const set<int> &vars)
           if (var2 > var1)
             continue;
 
-          NodeID d2 = d1->getDerivative(var2);
+          expr_t d2 = d1->getDerivative(var2);
           if (d2 == Zero)
             continue;
           second_derivatives[make_pair(eq, make_pair(var1, var2))] = d2;
@@ -930,7 +930,7 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
       int var2 = it->first.second.second;
       // By construction, var2 <= var1
 
-      NodeID d2 = it->second;
+      expr_t d2 = it->second;
 
       // Store only third derivatives such that var3 <= var2 <= var1
       for (set<int>::const_iterator it2 = vars.begin();
@@ -940,7 +940,7 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
           if (var3 > var2)
             continue;
 
-          NodeID d3 = d2->getDerivative(var3);
+          expr_t d3 = d2->getDerivative(var3);
           if (d3 == Zero)
             continue;
           third_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, var3)))] = d3;
@@ -957,7 +957,7 @@ ModelTree::computeThirdDerivatives(const set<int> &vars)
 void
 ModelTree::computeTemporaryTerms(bool is_matlab)
 {
-  map<NodeID, int> reference_count;
+  map<expr_t, int> reference_count;
   temporary_terms.clear();
 
   for (vector<BinaryOpNode *>::iterator it = equations.begin();
@@ -1045,11 +1045,11 @@ ModelTree::compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt
 void
 ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const
 {
-  for (map<int, NodeID>::const_iterator it = local_variables_table.begin();
+  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     {
       int id = it->first;
-      NodeID value = it->second;
+      expr_t value = it->second;
 
       if (IS_C(output_type))
         output << "double ";
@@ -1067,8 +1067,8 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
   for (int eq = 0; eq < (int) equations.size(); eq++)
     {
       BinaryOpNode *eq_node = equations[eq];
-      NodeID lhs = eq_node->get_arg1();
-      NodeID rhs = eq_node->get_arg2();
+      expr_t lhs = eq_node->get_arg1();
+      expr_t rhs = eq_node->get_arg2();
 
       // Test if the right hand side of the equation is empty.
       double vrhs = 1.0;
@@ -1113,8 +1113,8 @@ ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt
   for (int eq = 0; eq < (int) equations.size(); eq++)
     {
       BinaryOpNode *eq_node = equations[eq];
-      NodeID lhs = eq_node->get_arg1();
-      NodeID rhs = eq_node->get_arg2();
+      expr_t lhs = eq_node->get_arg1();
+      expr_t rhs = eq_node->get_arg2();
       FNUMEXPR_ fnumexpr(ModelEquation, eq);
       fnumexpr.write(code_file);
       // Test if the right hand side of the equation is empty.
@@ -1209,11 +1209,11 @@ ModelTree::writeLatexModelFile(const string &filename, ExprNodeOutputType output
          << "\\footnotesize" << endl;
 
   // Write model local variables
-  for (map<int, NodeID>::const_iterator it = local_variables_table.begin();
+  for (map<int, expr_t>::const_iterator it = local_variables_table.begin();
        it != local_variables_table.end(); it++)
     {
       int id = it->first;
-      NodeID value = it->second;
+      expr_t value = it->second;
 
       output << "\\begin{equation*}" << endl
              << symbol_table.getName(id) << " = ";
@@ -1237,7 +1237,7 @@ ModelTree::writeLatexModelFile(const string &filename, ExprNodeOutputType output
 }
 
 void
-ModelTree::addEquation(NodeID eq)
+ModelTree::addEquation(expr_t eq)
 {
   BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
   assert(beq != NULL && beq->get_op_code() == oEqual);
@@ -1252,7 +1252,7 @@ ModelTree::addEquationTags(int i, const string &key, const string &value)
 }
 
 void
-ModelTree::addAuxEquation(NodeID eq)
+ModelTree::addAuxEquation(expr_t eq)
 {
   BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
   assert(beq != NULL && beq->get_op_code() == oEqual);
diff --git a/ModelTree.hh b/ModelTree.hh
index 9d0084e79eb522136af19e7f1b694fbe760e672d..9297a1b228c57308bb679fa0f091551205cc3f17 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -30,8 +30,8 @@ using namespace std;
 
 #include "DataTree.hh"
 
-//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
-typedef vector<pair<EquationType, NodeID > > equation_type_and_normalized_equation_t;
+//! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
+typedef vector<pair<EquationType, expr_t > > equation_type_and_normalized_equation_t;
 
 //! Vector describing variables: max_lag in the block, max_lead in the block
 typedef vector<pair< int, int> > lag_lead_vector_t;
@@ -39,8 +39,8 @@ typedef vector<pair< int, int> > lag_lead_vector_t;
 //! for each block contains pair< pair<Simulation_Type, first_equation>, pair < Block_Size, Recursive_part_Size > >
 typedef vector<pair< pair< BlockSimulationType, int>, pair<int, int> > > block_type_firstequation_size_mfs_t;
 
-//! for a block contains derivatives pair< pair<block_equation_number, block_variable_number> , pair<lead_lag, NodeID> >
-typedef vector< pair<pair<int, int>, pair< int, NodeID > > > block_derivatives_equation_variable_laglead_nodeid_t;
+//! for a block contains derivatives pair< pair<block_equation_number, block_variable_number> , pair<lead_lag, expr_t> >
+typedef vector< pair<pair<int, int>, pair< int, expr_t > > > block_derivatives_equation_variable_laglead_nodeid_t;
 
 //! for all blocks derivatives description
 typedef vector<block_derivatives_equation_variable_laglead_nodeid_t> blocks_derivatives_t;
@@ -64,7 +64,7 @@ protected:
   //! Number of non-zero derivatives
   int NNZDerivatives[3];
 
-  typedef map<pair<int, int>, NodeID> first_derivatives_t;
+  typedef map<pair<int, int>, expr_t> first_derivatives_t;
   //! First order derivatives
   /*! First index is equation number, second is variable w.r. to which is computed the derivative.
     Only non-null derivatives are stored in the map.
@@ -72,7 +72,7 @@ protected:
   */
   first_derivatives_t first_derivatives;
 
-  typedef map<pair<int, pair<int, int> >, NodeID> second_derivatives_t;
+  typedef map<pair<int, pair<int, int> >, expr_t> second_derivatives_t;
   //! Second order derivatives
   /*! First index is equation number, second and third are variables w.r. to which is computed the derivative.
     Only non-null derivatives are stored in the map.
@@ -81,7 +81,7 @@ protected:
   */
   second_derivatives_t second_derivatives;
 
-  typedef map<pair<int, pair<int, pair<int, int> > >, NodeID> third_derivatives_t;
+  typedef map<pair<int, pair<int, pair<int, int> > >, expr_t> third_derivatives_t;
   //! Third order derivatives
   /*! First index is equation number, second, third and fourth are variables w.r. to which is computed the derivative.
     Only non-null derivatives are stored in the map.
@@ -131,7 +131,7 @@ protected:
 
   //! Sparse matrix of double to store the values of the Jacobian
   /*! First index is lag, second index is equation number, third index is endogenous type specific ID */
-  typedef map<pair<int, pair<int, int> >, NodeID> dynamic_jacob_map_t;
+  typedef map<pair<int, pair<int, int> >, expr_t> dynamic_jacob_map_t;
 
   //! Normalization of equations
   /*! Maps endogenous type specific IDs to equation numbers */
@@ -165,7 +165,7 @@ protected:
   //! Search the equations and variables belonging to the prologue and the epilogue of the model
   void computePrologueAndEpilogue(const jacob_map_t &static_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered);
   //! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
-  equation_type_and_normalized_equation_t equationTypeDetermination(const map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
+  equation_type_and_normalized_equation_t equationTypeDetermination(const map<pair<int, pair<int, int> >, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
   //! Compute the block decomposition and for a non-recusive block find the minimum feedback set
   void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const;
   //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
@@ -199,10 +199,10 @@ protected:
   virtual EquationType getBlockEquationType(int block_number, int equation_number) const = 0;
   //! Return true if the equation has been normalized
   virtual bool isBlockEquationRenormalized(int block_number, int equation_number) const = 0;
-  //! Return the NodeID of the equation equation_number belonging to the block block_number
-  virtual NodeID getBlockEquationNodeID(int block_number, int equation_number) const = 0;
-  //! Return the NodeID of the renormalized equation equation_number belonging to the block block_number
-  virtual NodeID getBlockEquationRenormalizedNodeID(int block_number, int equation_number) const = 0;
+  //! Return the expr_t of the equation equation_number belonging to the block block_number
+  virtual expr_t getBlockEquationExpr(int block_number, int equation_number) const = 0;
+  //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
+  virtual expr_t getBlockEquationRenormalizedExpr(int block_number, int equation_number) const = 0;
   //! Return the original number of equation equation_number belonging to the block block_number
   virtual int getBlockEquationID(int block_number, int equation_number) const = 0;
   //! Return the original number of variable variable_number belonging to the block block_number
@@ -215,11 +215,11 @@ protected:
 public:
   ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_arg);
   //! Declare a node as an equation of the model
-  void addEquation(NodeID eq);
+  void addEquation(expr_t eq);
   //! Adds tags to equation number i
   void addEquationTags(int i, const string &key, const string &value);
   //! Declare a node as an auxiliary equation of the model, adding it at the end of the list of auxiliary equations
-  void addAuxEquation(NodeID eq);
+  void addAuxEquation(expr_t eq);
   //! Returns the number of equations in the model
   int equation_number() const;
 
diff --git a/NumericalInitialization.cc b/NumericalInitialization.cc
index 8f6da6ec34eba9e06286f7c54bdea5a7b88ff045..ceccb9733d8d79e81baac61fb2f5528d03dfefd4 100644
--- a/NumericalInitialization.cc
+++ b/NumericalInitialization.cc
@@ -24,7 +24,7 @@
 #include "NumericalInitialization.hh"
 
 InitParamStatement::InitParamStatement(int symb_id_arg,
-                                       const NodeID param_value_arg,
+                                       const expr_t param_value_arg,
                                        const SymbolTable &symbol_table_arg) :
   symb_id(symb_id_arg),
   param_value(param_value_arg),
@@ -93,7 +93,7 @@ InitOrEndValStatement::writeInitValues(ostream &output) const
        it != init_values.end(); it++)
     {
       const int symb_id = it->first;
-      const NodeID expression = it->second;
+      const expr_t expression = it->second;
 
       SymbolType type = symbol_table.getType(symb_id);
       int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
@@ -190,7 +190,7 @@ HistValStatement::writeOutput(ostream &output, const string &basename) const
     {
       int symb_id = it->first.first;
       int lag = it->first.second;
-      const NodeID expression = it->second;
+      const expr_t expression = it->second;
 
       SymbolType type = symbol_table.getType(symb_id);
       if (type == eEndogenous && lag < 0)
@@ -256,8 +256,8 @@ HomotopyStatement::writeOutput(ostream &output, const string &basename) const
        it != homotopy_values.end(); it++)
     {
       const int &symb_id = it->first;
-      const NodeID expression1 = it->second.first;
-      const NodeID expression2 = it->second.second;
+      const expr_t expression1 = it->second.first;
+      const expr_t expression2 = it->second.second;
 
       const SymbolType type = symbol_table.getType(symb_id);
       const int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
diff --git a/NumericalInitialization.hh b/NumericalInitialization.hh
index 85ae929c0334bd7246aa872757a08e95db67a576..bde6a7eceb4ae7558871068e9f8ab1a4c80f2e8a 100644
--- a/NumericalInitialization.hh
+++ b/NumericalInitialization.hh
@@ -34,10 +34,10 @@ class InitParamStatement : public Statement
 {
 private:
   const int symb_id;
-  const NodeID param_value;
+  const expr_t param_value;
   const SymbolTable &symbol_table;
 public:
-  InitParamStatement(int symb_id_arg, const NodeID param_value_arg,
+  InitParamStatement(int symb_id_arg, const expr_t param_value_arg,
                      const SymbolTable &symbol_table_arg);
   virtual void checkPass(ModFileStructure &mod_file_struct);
   virtual void writeOutput(ostream &output, const string &basename) const;
@@ -52,7 +52,7 @@ public:
     We use a vector instead of a map, since the order of declaration matters:
     an initialization can depend on a previously initialized variable inside the block
   */
-  typedef vector<pair<int, NodeID> > init_values_t;
+  typedef vector<pair<int, expr_t> > init_values_t;
 protected:
   const init_values_t init_values;
   const SymbolTable &symbol_table;
@@ -91,9 +91,9 @@ public:
   /*!
     Contrary to Initval and Endval, we use a map, since it is impossible to reuse
     a given initialization value in a second initialization inside the block.
-    Maps pairs (symbol_id, lag) to NodeID
+    Maps pairs (symbol_id, lag) to expr_t
   */
-  typedef map<pair<int, int>, NodeID> hist_values_t;
+  typedef map<pair<int, int>, expr_t> hist_values_t;
 private:
   const hist_values_t hist_values;
   const SymbolTable &symbol_table;
@@ -116,8 +116,8 @@ class HomotopyStatement : public Statement
 {
 public:
   //! Stores the declarations of homotopy_setup
-  /*! Order matter so we use a vector. First NodeID can be NULL if no initial value given. */
-  typedef vector<pair<int, pair<NodeID, NodeID> > > homotopy_values_t;
+  /*! Order matter so we use a vector. First expr_t can be NULL if no initial value given. */
+  typedef vector<pair<int, pair<expr_t, expr_t> > > homotopy_values_t;
 private:
   const homotopy_values_t homotopy_values;
   const SymbolTable &symbol_table;
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index b38486036e25ec3bf40f3f2a8679a73626632cf2..4b6a663df72e0aaff243ac38168b97541c989738 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -191,27 +191,27 @@ ParsingDriver::add_equation_tags(string *key, string *value)
   delete value;
 }
 
-NodeID
+expr_t
 ParsingDriver::add_constant(string *constant)
 {
-  NodeID id = data_tree->AddNumConstant(*constant);
+  expr_t id = data_tree->AddNumConstant(*constant);
   delete constant;
   return id;
 }
 
-NodeID
+expr_t
 ParsingDriver::add_nan_constant()
 {
   return data_tree->NaN;
 }
 
-NodeID
+expr_t
 ParsingDriver::add_inf_constant()
 {
   return data_tree->Infinity;
 }
 
-NodeID
+expr_t
 ParsingDriver::add_model_variable(string *name)
 {
   check_symbol_existence(*name);
@@ -220,7 +220,7 @@ ParsingDriver::add_model_variable(string *name)
   return add_model_variable(symb_id, 0);
 }
 
-NodeID
+expr_t
 ParsingDriver::add_model_variable(int symb_id, int lag)
 {
   assert(symb_id >= 0);
@@ -239,7 +239,7 @@ ParsingDriver::add_model_variable(int symb_id, int lag)
   return model_tree->AddVariable(symb_id, lag);
 }
 
-NodeID
+expr_t
 ParsingDriver::add_expression_variable(string *name)
 {
   // If symbol doesn't exist, then declare it as a mod file local variable
@@ -251,7 +251,7 @@ ParsingDriver::add_expression_variable(string *name)
     error("Variable " + *name + " not allowed outside model declaration. Its scope is only inside model.");
 
   int symb_id = mod_file->symbol_table.getID(*name);
-  NodeID id = data_tree->AddVariable(symb_id);
+  expr_t id = data_tree->AddVariable(symb_id);
 
   delete name;
   return id;
@@ -286,7 +286,7 @@ ParsingDriver::dsample(string *arg1, string *arg2)
 }
 
 void
-ParsingDriver::init_param(string *name, NodeID rhs)
+ParsingDriver::init_param(string *name, expr_t rhs)
 {
   check_symbol_existence(*name);
   int symb_id = mod_file->symbol_table.getID(*name);
@@ -299,7 +299,7 @@ ParsingDriver::init_param(string *name, NodeID rhs)
 }
 
 void
-ParsingDriver::init_val(string *name, NodeID rhs)
+ParsingDriver::init_val(string *name, expr_t rhs)
 {
   check_symbol_existence(*name);
   int symb_id = mod_file->symbol_table.getID(*name);
@@ -323,7 +323,7 @@ ParsingDriver::initval_file(string *filename)
 }
 
 void
-ParsingDriver::hist_val(string *name, string *lag, NodeID rhs)
+ParsingDriver::hist_val(string *name, string *lag, expr_t rhs)
 {
   check_symbol_existence(*name);
   int symb_id = mod_file->symbol_table.getID(*name);
@@ -347,7 +347,7 @@ ParsingDriver::hist_val(string *name, string *lag, NodeID rhs)
 }
 
 void
-ParsingDriver::homotopy_val(string *name, NodeID val1, NodeID val2)
+ParsingDriver::homotopy_val(string *name, expr_t val1, expr_t val2)
 {
   check_symbol_existence(*name);
   int symb_id = mod_file->symbol_table.getID(*name);
@@ -510,7 +510,7 @@ ParsingDriver::add_det_shock(string *var, bool conditional_forecast)
 }
 
 void
-ParsingDriver::add_stderr_shock(string *var, NodeID value)
+ParsingDriver::add_stderr_shock(string *var, expr_t value)
 {
   check_symbol_existence(*var);
   int symb_id = mod_file->symbol_table.getID(*var);
@@ -529,7 +529,7 @@ ParsingDriver::add_stderr_shock(string *var, NodeID value)
 }
 
 void
-ParsingDriver::add_var_shock(string *var, NodeID value)
+ParsingDriver::add_var_shock(string *var, expr_t value)
 {
   check_symbol_existence(*var);
   int symb_id = mod_file->symbol_table.getID(*var);
@@ -548,7 +548,7 @@ ParsingDriver::add_var_shock(string *var, NodeID value)
 }
 
 void
-ParsingDriver::add_covar_shock(string *var1, string *var2, NodeID value)
+ParsingDriver::add_covar_shock(string *var1, string *var2, expr_t value)
 {
   check_symbol_existence(*var1);
   check_symbol_existence(*var2);
@@ -577,7 +577,7 @@ ParsingDriver::add_covar_shock(string *var1, string *var2, NodeID value)
 }
 
 void
-ParsingDriver::add_correl_shock(string *var1, string *var2, NodeID value)
+ParsingDriver::add_correl_shock(string *var1, string *var2, expr_t value)
 {
   check_symbol_existence(*var1);
   check_symbol_existence(*var2);
@@ -626,7 +626,7 @@ ParsingDriver::add_period(string *p1)
 }
 
 void
-ParsingDriver::add_value(NodeID value)
+ParsingDriver::add_value(expr_t value)
 {
   det_shocks_values.push_back(value);
 }
@@ -746,7 +746,7 @@ ParsingDriver::add_to_row_const(string *s)
 }
 
 void
-ParsingDriver::add_to_row(NodeID v)
+ParsingDriver::add_to_row(expr_t v)
 {
   sigmae_row.push_back(v);
 }
@@ -990,7 +990,7 @@ ParsingDriver::set_trends()
 }
 
 void
-ParsingDriver::set_trend_element(string *arg1, NodeID arg2)
+ParsingDriver::set_trend_element(string *arg1, expr_t arg2)
 {
   check_symbol_existence(*arg1);
   if (trend_elements.find(*arg1) != trend_elements.end())
@@ -1000,7 +1000,7 @@ ParsingDriver::set_trend_element(string *arg1, NodeID arg2)
 }
 
 void
-ParsingDriver::set_optim_weights(string *name, NodeID value)
+ParsingDriver::set_optim_weights(string *name, expr_t value)
 {
   check_symbol_existence(*name);
   if (mod_file->symbol_table.getType(*name) != eEndogenous)
@@ -1012,7 +1012,7 @@ ParsingDriver::set_optim_weights(string *name, NodeID value)
 }
 
 void
-ParsingDriver::set_optim_weights(string *name1, string *name2, NodeID value)
+ParsingDriver::set_optim_weights(string *name1, string *name2, expr_t value)
 {
   check_symbol_existence(*name1);
   if (mod_file->symbol_table.getType(*name1) != eEndogenous)
@@ -1057,7 +1057,7 @@ ParsingDriver::run_osr()
 }
 
 void
-ParsingDriver::set_calib_var(string *name, string *weight, NodeID expression)
+ParsingDriver::set_calib_var(string *name, string *weight, expr_t expression)
 {
   check_symbol_existence(*name);
   if (mod_file->symbol_table.getType(*name) != eEndogenous
@@ -1075,7 +1075,7 @@ ParsingDriver::set_calib_var(string *name, string *weight, NodeID expression)
 
 void
 ParsingDriver::set_calib_covar(string *name1, string *name2,
-                               string *weight, NodeID expression)
+                               string *weight, expr_t expression)
 {
   check_symbol_existence(*name1);
   check_symbol_existence(*name2);
@@ -1100,7 +1100,7 @@ ParsingDriver::set_calib_covar(string *name1, string *name2,
 
 void
 ParsingDriver::set_calib_ac(string *name, string *ar,
-                            string *weight, NodeID expression)
+                            string *weight, expr_t expression)
 {
   check_symbol_existence(*name);
   if (mod_file->symbol_table.getType(*name) != eEndogenous)
@@ -1199,10 +1199,10 @@ ParsingDriver::begin_planner_objective()
 }
 
 void
-ParsingDriver::end_planner_objective(NodeID expr)
+ParsingDriver::end_planner_objective(expr_t expr)
 {
   // Add equation corresponding to expression
-  NodeID eq = model_tree->AddEqual(expr, model_tree->Zero);
+  expr_t eq = model_tree->AddEqual(expr, model_tree->Zero);
   model_tree->addEquation(eq);
 
   mod_file->addStatement(new PlannerObjectiveStatement(dynamic_cast<StaticModel *>(model_tree)));
@@ -1376,22 +1376,22 @@ ParsingDriver::conditional_forecast_paths()
   det_shocks.clear();
 }
 
-NodeID
-ParsingDriver::add_model_equal(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2)
 {
-  NodeID id = model_tree->AddEqual(arg1, arg2);
+  expr_t id = model_tree->AddEqual(arg1, arg2);
   model_tree->addEquation(id);
   return id;
 }
 
-NodeID
-ParsingDriver::add_model_equal_with_zero_rhs(NodeID arg)
+expr_t
+ParsingDriver::add_model_equal_with_zero_rhs(expr_t arg)
 {
   return add_model_equal(arg, model_tree->Zero);
 }
 
 void
-ParsingDriver::declare_and_init_model_local_variable(string *name, NodeID rhs)
+ParsingDriver::declare_and_init_model_local_variable(string *name, expr_t rhs)
 {
   int symb_id;
   try
@@ -1445,82 +1445,82 @@ ParsingDriver::change_type(SymbolType new_type, vector<string *> *var_list)
   delete var_list;
 }
 
-NodeID
-ParsingDriver::add_plus(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_plus(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddPlus(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_minus(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_minus(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddMinus(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_uminus(NodeID arg1)
+expr_t
+ParsingDriver::add_uminus(expr_t arg1)
 {
   return data_tree->AddUMinus(arg1);
 }
 
-NodeID
-ParsingDriver::add_times(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_times(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddTimes(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_divide(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_divide(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddDivide(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_less(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_less(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddLess(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_greater(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_greater(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddGreater(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_less_equal(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_less_equal(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddLessEqual(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_greater_equal(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_greater_equal(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddGreaterEqual(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_equal_equal(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_equal_equal(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddEqualEqual(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_different(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_different(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddDifferent(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_power(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_power(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddPower(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_expectation(string *arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_expectation(string *arg1, expr_t arg2)
 {
-  NodeID expectationNode;
+  expr_t expectationNode;
   if ("varobs"==*arg1 || "full"==*arg1)
     if (dynamic_cast<VariableNode *>(arg2) == NULL)
       error("EXPECTATION(" + *arg1  + ")(X) can only be used when X is a single variable.");
@@ -1536,146 +1536,146 @@ ParsingDriver::add_expectation(string *arg1, NodeID arg2)
   return expectationNode;
 }
 
-NodeID
-ParsingDriver::add_exp(NodeID arg1)
+expr_t
+ParsingDriver::add_exp(expr_t arg1)
 {
   return data_tree->AddExp(arg1);
 }
 
-NodeID
-ParsingDriver::add_log(NodeID arg1)
+expr_t
+ParsingDriver::add_log(expr_t arg1)
 {
   return data_tree->AddLog(arg1);
 }
 
-NodeID
-ParsingDriver::add_log10(NodeID arg1)
+expr_t
+ParsingDriver::add_log10(expr_t arg1)
 {
   return data_tree->AddLog10(arg1);
 }
 
-NodeID
-ParsingDriver::add_cos(NodeID arg1)
+expr_t
+ParsingDriver::add_cos(expr_t arg1)
 {
   return data_tree->AddCos(arg1);
 }
 
-NodeID
-ParsingDriver::add_sin(NodeID arg1)
+expr_t
+ParsingDriver::add_sin(expr_t arg1)
 {
   return data_tree->AddSin(arg1);
 }
 
-NodeID
-ParsingDriver::add_tan(NodeID arg1)
+expr_t
+ParsingDriver::add_tan(expr_t arg1)
 {
   return data_tree->AddTan(arg1);
 }
 
-NodeID
-ParsingDriver::add_acos(NodeID arg1)
+expr_t
+ParsingDriver::add_acos(expr_t arg1)
 {
   return data_tree->AddAcos(arg1);
 }
 
-NodeID
-ParsingDriver::add_asin(NodeID arg1)
+expr_t
+ParsingDriver::add_asin(expr_t arg1)
 {
   return data_tree->AddAsin(arg1);
 }
 
-NodeID
-ParsingDriver::add_atan(NodeID arg1)
+expr_t
+ParsingDriver::add_atan(expr_t arg1)
 {
   return data_tree->AddAtan(arg1);
 }
 
-NodeID
-ParsingDriver::add_cosh(NodeID arg1)
+expr_t
+ParsingDriver::add_cosh(expr_t arg1)
 {
   return data_tree->AddCosh(arg1);
 }
 
-NodeID
-ParsingDriver::add_sinh(NodeID arg1)
+expr_t
+ParsingDriver::add_sinh(expr_t arg1)
 {
   return data_tree->AddSinh(arg1);
 }
 
-NodeID
-ParsingDriver::add_tanh(NodeID arg1)
+expr_t
+ParsingDriver::add_tanh(expr_t arg1)
 {
   return data_tree->AddTanh(arg1);
 }
 
-NodeID
-ParsingDriver::add_acosh(NodeID arg1)
+expr_t
+ParsingDriver::add_acosh(expr_t arg1)
 {
   return data_tree->AddAcosh(arg1);
 }
 
-NodeID
-ParsingDriver::add_asinh(NodeID arg1)
+expr_t
+ParsingDriver::add_asinh(expr_t arg1)
 {
   return data_tree->AddAsinh(arg1);
 }
 
-NodeID
-ParsingDriver::add_atanh(NodeID arg1)
+expr_t
+ParsingDriver::add_atanh(expr_t arg1)
 {
   return data_tree->AddAtanh(arg1);
 }
 
-NodeID
-ParsingDriver::add_sqrt(NodeID arg1)
+expr_t
+ParsingDriver::add_sqrt(expr_t arg1)
 {
   return data_tree->AddSqrt(arg1);
 }
 
-NodeID
-ParsingDriver::add_max(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_max(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddMax(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_min(NodeID arg1, NodeID arg2)
+expr_t
+ParsingDriver::add_min(expr_t arg1, expr_t arg2)
 {
   return data_tree->AddMin(arg1, arg2);
 }
 
-NodeID
-ParsingDriver::add_normcdf(NodeID arg1, NodeID arg2, NodeID arg3)
+expr_t
+ParsingDriver::add_normcdf(expr_t arg1, expr_t arg2, expr_t arg3)
 {
   return data_tree->AddNormcdf(arg1, arg2, arg3);
 }
 
-NodeID
-ParsingDriver::add_normcdf(NodeID arg)
+expr_t
+ParsingDriver::add_normcdf(expr_t arg)
 {
   return add_normcdf(arg, data_tree->Zero, data_tree->One);
 }
 
-NodeID
-ParsingDriver::add_normpdf(NodeID arg1, NodeID arg2, NodeID arg3)
+expr_t
+ParsingDriver::add_normpdf(expr_t arg1, expr_t arg2, expr_t arg3)
 {
   return data_tree->AddNormpdf(arg1, arg2, arg3);
 }
 
-NodeID
-ParsingDriver::add_normpdf(NodeID arg)
+expr_t
+ParsingDriver::add_normpdf(expr_t arg)
 {
   return add_normpdf(arg, data_tree->Zero, data_tree->One);
 }
 
-NodeID
-ParsingDriver::add_erf(NodeID arg1)
+expr_t
+ParsingDriver::add_erf(expr_t arg1)
 {
   return data_tree->AddErf(arg1);
 }
 
-NodeID
-ParsingDriver::add_steady_state(NodeID arg1)
+expr_t
+ParsingDriver::add_steady_state(expr_t arg1)
 {
   return data_tree->AddSteadyState(arg1);
 }
@@ -1744,20 +1744,20 @@ ParsingDriver::external_function()
 void
 ParsingDriver::push_external_function_arg_vector_onto_stack()
 {
-  vector<NodeID> emptyvec;
+  vector<expr_t> emptyvec;
   stack_external_function_args.push(emptyvec);
 }
 
 void
-ParsingDriver::add_external_function_arg(NodeID arg)
+ParsingDriver::add_external_function_arg(expr_t arg)
 {
   stack_external_function_args.top().push_back(arg);
 }
 
-NodeID
+expr_t
 ParsingDriver::add_model_var_or_external_function(string *function_name, bool in_model_block)
 {
-  NodeID nid;
+  expr_t nid;
   if (mod_file->symbol_table.exists(*function_name))
     {
       if (mod_file->symbol_table.getType(*function_name) != eExternalFunction)
@@ -1865,7 +1865,7 @@ ParsingDriver::begin_steady_state_model()
 }
 
 void
-ParsingDriver::add_steady_state_model_equal(string *varname, NodeID expr)
+ParsingDriver::add_steady_state_model_equal(string *varname, expr_t expr)
 {
   int id;
   try
diff --git a/ParsingDriver.hh b/ParsingDriver.hh
index b347b7f84749e3b0a4ed542bcb2a56c4c2f6514e..597c735490323f3e6b66a88154e010e54be0db3c 100644
--- a/ParsingDriver.hh
+++ b/ParsingDriver.hh
@@ -127,7 +127,7 @@ private:
   //! Temporary storage for periods of deterministic shocks
   vector<pair<int, int> > det_shocks_periods;
   //! Temporary storage for values of deterministic shocks
-  vector<NodeID> det_shocks_values;
+  vector<expr_t> det_shocks_values;
   //! Temporary storage for variances of shocks
   ShocksStatement::var_and_std_shocks_t var_shocks;
   //! Temporary storage for standard errors of shocks
@@ -158,7 +158,7 @@ private:
   bool svar_lower_cholesky;
 
   //! Temporary storage for argument list of external function
-  stack<vector<NodeID> >  stack_external_function_args;
+  stack<vector<expr_t> >  stack_external_function_args;
   //! Temporary storage for the symb_id associated with the "name" symbol of the current external_function statement
   int current_external_function_id;
   //! Temporary storage for option list provided to external_function()
@@ -166,7 +166,7 @@ private:
   //! reset the values for temporary storage
   void reset_current_external_function_options();
   //! Adds a model lagged variable to ModelTree and VariableTable
-  NodeID add_model_variable(int symb_id, int lag);
+  expr_t add_model_variable(int symb_id, int lag);
 
   //! The mod file representation constructed by this ParsingDriver
   ModFile *mod_file;
@@ -218,21 +218,21 @@ public:
   //! Adds a predetermined_variable
   void add_predetermined_variable(string *name);
   //! Declares and initializes a local parameter
-  void declare_and_init_model_local_variable(string *name, NodeID rhs);
+  void declare_and_init_model_local_variable(string *name, expr_t rhs);
   //! Changes type of a symbol
   void change_type(SymbolType new_type, vector<string *> *var_list);
   //! Adds a list of tags for the current equation
   void add_equation_tags(string *key, string *value);
   //! Adds a constant to DataTree
-  NodeID add_constant(string *constant);
+  expr_t add_constant(string *constant);
   //! Adds a NaN constant to DataTree
-  NodeID add_nan_constant();
+  expr_t add_nan_constant();
   //! Adds an Inf constant to DataTree
-  NodeID add_inf_constant();
+  expr_t add_inf_constant();
   //! Adds a model variable to ModelTree and VariableTable
-  NodeID add_model_variable(string *name);
+  expr_t add_model_variable(string *name);
   //! Adds an Expression's variable
-  NodeID add_expression_variable(string *name);
+  expr_t add_expression_variable(string *name);
   //! Adds a "periods" statement
   void periods(string *periods);
   //! Adds a "dsample" statement
@@ -240,14 +240,14 @@ public:
   //! Adds a "dsample" statement
   void dsample(string *arg1, string *arg2);
   //! Writes parameter intitialisation expression
-  void init_param(string *name, NodeID rhs);
+  void init_param(string *name, expr_t rhs);
   //! Writes an initval block
-  void init_val(string *name, NodeID rhs);
+  void init_val(string *name, expr_t rhs);
   //! Writes an histval block
-  void hist_val(string *name, string *lag, NodeID rhs);
+  void hist_val(string *name, string *lag, expr_t rhs);
   //! Adds an entry in a homotopy_setup block
   /*! Second argument "val1" can be NULL if no initial value provided */
-  void homotopy_val(string *name, NodeID val1, NodeID val2);
+  void homotopy_val(string *name, expr_t val1, expr_t val2);
   //! Writes end of an initval block
   void end_initval();
   //! Writes end of an endval block
@@ -265,19 +265,19 @@ public:
   //! Adds a deterministic chock or a path element inside a conditional_forecast_paths block
   void add_det_shock(string *var, bool conditional_forecast);
   //! Adds a std error chock
-  void add_stderr_shock(string *var, NodeID value);
+  void add_stderr_shock(string *var, expr_t value);
   //! Adds a variance chock
-  void add_var_shock(string *var, NodeID value);
+  void add_var_shock(string *var, expr_t value);
   //! Adds a covariance chock
-  void add_covar_shock(string *var1, string *var2, NodeID value);
+  void add_covar_shock(string *var1, string *var2, expr_t value);
   //! Adds a correlated chock
-  void add_correl_shock(string *var1, string *var2, NodeID value);
+  void add_correl_shock(string *var1, string *var2, expr_t value);
   //! Adds a shock period range
   void add_period(string *p1, string *p2);
   //! Adds a shock period
   void add_period(string *p1);
   //! Adds a deterministic shock value
-  void add_value(NodeID value);
+  void add_value(expr_t value);
   //! Adds a deterministic shock value
   void add_value(string *p1);
   //! Writes a Sigma_e block
@@ -287,7 +287,7 @@ public:
   //! Adds a constant element to current row of Sigma_e
   void add_to_row_const(string *s);
   //! Adds an expression element to current row of Sigma_e
-  void add_to_row(NodeID v);
+  void add_to_row(expr_t v);
   //! Write a steady command
   void steady();
   //! Sets an option to a numerical value
@@ -359,17 +359,17 @@ public:
   //! Forecast Statement
   void forecast();
   void set_trends();
-  void set_trend_element(string *arg1, NodeID arg2);
+  void set_trend_element(string *arg1, expr_t arg2);
   void set_unit_root_vars();
   void optim_weights();
-  void set_optim_weights(string *name, NodeID value);
-  void set_optim_weights(string *name1, string *name2, NodeID value);
+  void set_optim_weights(string *name, expr_t value);
+  void set_optim_weights(string *name1, string *name2, expr_t value);
   void set_osr_params();
   void run_osr();
   void run_calib_var();
-  void set_calib_var(string *name, string *weight, NodeID expression);
-  void set_calib_covar(string *name1, string *name2, string *weight, NodeID expression);
-  void set_calib_ac(string *name, string *ar, string *weight, NodeID expression);
+  void set_calib_var(string *name, string *weight, expr_t expression);
+  void set_calib_covar(string *name1, string *name2, string *weight, expr_t expression);
+  void set_calib_ac(string *name, string *ar, string *weight, expr_t expression);
   void run_calib(int covar);
   void run_dynasave(string *filename);
   void run_dynatype(string *filename);
@@ -381,7 +381,7 @@ public:
   //! Begin a planner_objective statement
   void begin_planner_objective();
   //! End a planner objective statement
-  void end_planner_objective(NodeID expr);
+  void end_planner_objective(expr_t expr);
   //! ramsey policy statement
   void ramsey_policy();
   //! Adds a write_latex_dynamic_model statement
@@ -409,89 +409,89 @@ public:
   //! Plot conditional forecast statement
   void plot_conditional_forecast(string *periods = NULL);
   //! Writes token "arg1=arg2" to model tree
-  NodeID add_model_equal(NodeID arg1, NodeID arg2);
+  expr_t add_model_equal(expr_t arg1, expr_t arg2);
   //! Writes token "arg=0" to model tree
-  NodeID add_model_equal_with_zero_rhs(NodeID arg);
+  expr_t add_model_equal_with_zero_rhs(expr_t arg);
   //! Writes token "arg1+arg2" to model tree
-  NodeID add_plus(NodeID arg1, NodeID arg2);
+  expr_t add_plus(expr_t arg1, expr_t arg2);
   //! Writes token "arg1-arg2" to model tree
-  NodeID add_minus(NodeID arg1,  NodeID arg2);
+  expr_t add_minus(expr_t arg1,  expr_t arg2);
   //! Writes token "-arg1" to model tree
-  NodeID add_uminus(NodeID arg1);
+  expr_t add_uminus(expr_t arg1);
   //! Writes token "arg1*arg2" to model tree
-  NodeID add_times(NodeID arg1,  NodeID arg2);
+  expr_t add_times(expr_t arg1,  expr_t arg2);
   //! Writes token "arg1/arg2" to model tree
-  NodeID add_divide(NodeID arg1,  NodeID arg2);
+  expr_t add_divide(expr_t arg1,  expr_t arg2);
   //! Writes token "arg1<arg2" to model tree
-  NodeID add_less(NodeID arg1, NodeID arg2);
-  //! Writes token "arg1>arg2" to model treeNodeID
-  NodeID add_greater(NodeID arg1, NodeID arg2);
-  //! Writes token "arg1<=arg2" to model treeNodeID
-  NodeID add_less_equal(NodeID arg1, NodeID arg2);
-  //! Writes token "arg1>=arg2" to model treeNodeID
-  NodeID add_greater_equal(NodeID arg1, NodeID arg2);
-  //! Writes token "arg1==arg2" to model treeNodeIDNodeID
-  NodeID add_equal_equal(NodeID arg1, NodeID arg2);
-  //! Writes token "arg1!=arg2" to model treeNodeIDNodeID
-  NodeID add_different(NodeID arg1, NodeID arg2);
+  expr_t add_less(expr_t arg1, expr_t arg2);
+  //! Writes token "arg1>arg2" to model treeexpr_t
+  expr_t add_greater(expr_t arg1, expr_t arg2);
+  //! Writes token "arg1<=arg2" to model treeexpr_t
+  expr_t add_less_equal(expr_t arg1, expr_t arg2);
+  //! Writes token "arg1>=arg2" to model treeexpr_t
+  expr_t add_greater_equal(expr_t arg1, expr_t arg2);
+  //! Writes token "arg1==arg2" to model treeexpr_texpr_t
+  expr_t add_equal_equal(expr_t arg1, expr_t arg2);
+  //! Writes token "arg1!=arg2" to model treeexpr_texpr_t
+  expr_t add_different(expr_t arg1, expr_t arg2);
   //! Writes token "arg1^arg2" to model tree
-  NodeID add_power(NodeID arg1,  NodeID arg2);
+  expr_t add_power(expr_t arg1,  expr_t arg2);
   //! Writes token "E(arg1)(arg2)" to model tree
-  NodeID add_expectation(string *arg1,  NodeID arg2);
+  expr_t add_expectation(string *arg1,  expr_t arg2);
   //! Writes token "exp(arg1)" to model tree
-  NodeID add_exp(NodeID arg1);
+  expr_t add_exp(expr_t arg1);
   //! Writes token "log(arg1)" to model tree
-  NodeID add_log(NodeID arg1);
+  expr_t add_log(expr_t arg1);
   //! Writes token "log10(arg1)" to model tree
-  NodeID add_log10(NodeID arg1);
+  expr_t add_log10(expr_t arg1);
   //! Writes token "cos(arg1)" to model tree
-  NodeID add_cos(NodeID arg1);
+  expr_t add_cos(expr_t arg1);
   //! Writes token "sin(arg1)" to model tree
-  NodeID add_sin(NodeID arg1);
+  expr_t add_sin(expr_t arg1);
   //! Writes token "tan(arg1)" to model tree
-  NodeID add_tan(NodeID arg1);
+  expr_t add_tan(expr_t arg1);
   //! Writes token "acos(arg1)" to model tree
-  NodeID add_acos(NodeID arg1);
+  expr_t add_acos(expr_t arg1);
   //! Writes token "asin(arg1)" to model tree
-  NodeID add_asin(NodeID arg1);
+  expr_t add_asin(expr_t arg1);
   //! Writes token "atan(arg1)" to model tree
-  NodeID add_atan(NodeID arg1);
+  expr_t add_atan(expr_t arg1);
   //! Writes token "cosh(arg1)" to model tree
-  NodeID add_cosh(NodeID arg1);
+  expr_t add_cosh(expr_t arg1);
   //! Writes token "sinh(arg1)" to model tree
-  NodeID add_sinh(NodeID arg1);
+  expr_t add_sinh(expr_t arg1);
   //! Writes token "tanh(arg1)" to model tree
-  NodeID add_tanh(NodeID arg1);
+  expr_t add_tanh(expr_t arg1);
   //! Writes token "acosh(arg1)" to model tree
-  NodeID add_acosh(NodeID arg1);
+  expr_t add_acosh(expr_t arg1);
   //! Writes token "asin(arg1)" to model tree
-  NodeID add_asinh(NodeID arg1);
+  expr_t add_asinh(expr_t arg1);
   //! Writes token "atanh(arg1)" to model tree
-  NodeID add_atanh(NodeID arg1);
+  expr_t add_atanh(expr_t arg1);
   //! Writes token "sqrt(arg1)" to model tree
-  NodeID add_sqrt(NodeID arg1);
+  expr_t add_sqrt(expr_t arg1);
   //! Writes token "max(arg1,arg2)" to model tree
-  NodeID add_max(NodeID arg1, NodeID arg2);
+  expr_t add_max(expr_t arg1, expr_t arg2);
   //! Writes token "min(arg1,arg2)" to model tree
-  NodeID add_min(NodeID arg1, NodeID arg2);
+  expr_t add_min(expr_t arg1, expr_t arg2);
   //! Writes token "normcdf(arg1,arg2,arg3)" to model tree
-  NodeID add_normcdf(NodeID arg1, NodeID arg2, NodeID arg3);
+  expr_t add_normcdf(expr_t arg1, expr_t arg2, expr_t arg3);
   //! Writes token "normcdf(arg,0,1)" to model tree
-  NodeID add_normcdf(NodeID arg);
+  expr_t add_normcdf(expr_t arg);
   //! Writes token "normpdf(arg1,arg2,arg3)" to model tree
-  NodeID add_normpdf(NodeID arg1, NodeID arg2, NodeID arg3);
+  expr_t add_normpdf(expr_t arg1, expr_t arg2, expr_t arg3);
   //! Writes token "normpdf(arg,0,1)" to model tree
-  NodeID add_normpdf(NodeID arg);
+  expr_t add_normpdf(expr_t arg);
   //! Writes token "erf(arg)" to model tree
-  NodeID add_erf(NodeID arg);
+  expr_t add_erf(expr_t arg);
   //! Writes token "steadyState(arg1)" to model tree
-  NodeID add_steady_state(NodeID arg1);
+  expr_t add_steady_state(expr_t arg1);
   //! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)
   void push_external_function_arg_vector_onto_stack();
   //! Adds an external function argument
-  void add_external_function_arg(NodeID arg);
+  void add_external_function_arg(expr_t arg);
   //! Adds an external function call node
-  NodeID add_model_var_or_external_function(string *function_name, bool in_model_block);
+  expr_t add_model_var_or_external_function(string *function_name, bool in_model_block);
   //! Adds a native statement
   void add_native(const string &s);
   //! Adds a native statement, first removing the set of characters passed in token (and everything after)
@@ -501,7 +501,7 @@ public:
   //! Begin a steady_state_model block
   void begin_steady_state_model();
   //! Add an assignment equation in steady_state_model block
-  void add_steady_state_model_equal(string *varname, NodeID expr);
+  void add_steady_state_model_equal(string *varname, expr_t expr);
 };
 
 #endif // ! PARSING_DRIVER_HH
diff --git a/Shocks.cc b/Shocks.cc
index 3ccc06b3d0f0d448a19d074ca62c3e0ece52466f..dc66b7c192fce0c5e7087abed45024dfb3f64305 100644
--- a/Shocks.cc
+++ b/Shocks.cc
@@ -50,7 +50,7 @@ AbstractShocksStatement::writeDetShocks(ostream &output) const
         {
           const int &period1 = it->second[i].period1;
           const int &period2 = it->second[i].period2;
-          const NodeID value = it->second[i].value;
+          const expr_t value = it->second[i].value;
 
           if (period1 == period2)
             {
diff --git a/Shocks.hh b/Shocks.hh
index b90ac911e832e82764433e6553865f325d61a002..51475f71165851b16cae83d58d5b48f490be6195 100644
--- a/Shocks.hh
+++ b/Shocks.hh
@@ -37,7 +37,7 @@ public:
   {
     int period1;
     int period2;
-    NodeID value;
+    expr_t value;
   };
   typedef map<int, vector<DetShockElement> > det_shocks_t;
 protected:
@@ -55,8 +55,8 @@ protected:
 class ShocksStatement : public AbstractShocksStatement
 {
 public:
-  typedef map<int, NodeID> var_and_std_shocks_t;
-  typedef map<pair<int, int>, NodeID> covar_and_corr_shocks_t;
+  typedef map<int, expr_t> var_and_std_shocks_t;
+  typedef map<pair<int, int>, expr_t> covar_and_corr_shocks_t;
 private:
   const var_and_std_shocks_t var_shocks, std_shocks;
   const covar_and_corr_shocks_t covar_shocks, corr_shocks;
diff --git a/SigmaeInitialization.hh b/SigmaeInitialization.hh
index e619defedfe18442d825a14460367a20ea7aad78..c2decb3f831887b63284e7a489256d54c6dec636 100644
--- a/SigmaeInitialization.hh
+++ b/SigmaeInitialization.hh
@@ -39,7 +39,7 @@ public:
       eUpper = 1               //!< Upper triangular matrix
     };
   //! Type of a matrix row
-  typedef vector<NodeID> row_t;
+  typedef vector<expr_t> row_t;
   //! Type of a complete matrix
   typedef vector<row_t> matrix_t;
 
diff --git a/StaticModel.cc b/StaticModel.cc
index 4a35a46e56ff454d5e50936f905137b6c5f95c39..58a7328f520f70c2a3dbf23ce7b1e55e262dd1b4 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -61,7 +61,7 @@ StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx
 void
 StaticModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, map_idx_t &map_idx) const
 {
-  map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
+  map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
   if (it != first_chain_rule_derivatives.end())
     (it->second)->compile(code_file, false, temporary_terms, map_idx, false, false);
   else
@@ -84,8 +84,8 @@ StaticModel::initializeVariablesAndEquations()
 void
 StaticModel::computeTemporaryTermsOrdered()
 {
-  map<NodeID, pair<int, int> > first_occurence;
-  map<NodeID, int> reference_count;
+  map<expr_t, pair<int, int> > first_occurence;
+  map<expr_t, int> reference_count;
   BinaryOpNode *eq_node;
   first_derivatives_t::const_iterator it;
   first_chain_rule_derivatives_t::const_iterator it_chr;
@@ -113,16 +113,16 @@ StaticModel::computeTemporaryTermsOrdered()
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedNodeID(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
+                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
               else
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
                 }
             }
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
             {
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  block_size-1);
             }
           set<int> temporary_terms_in_use;
@@ -142,16 +142,16 @@ StaticModel::computeTemporaryTermsOrdered()
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedNodeID(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
+                getBlockEquationRenormalizedExpr(block, i)->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms,  i);
               else
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   eq_node->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, i);
                 }
             }
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
             {
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
             }
 
@@ -166,16 +166,16 @@ StaticModel::computeTemporaryTermsOrdered()
           for (unsigned int i = 0; i < block_size; i++)
             {
               if (i < block_nb_recursives && isBlockEquationRenormalized(block, i))
-                getBlockEquationRenormalizedNodeID(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
+                getBlockEquationRenormalizedExpr(block, i)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
               else
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   eq_node->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
                 }
             }
           for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
             {
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
             }
           for (int i = 0; i < (int) getBlockSize(block); i++)
@@ -203,9 +203,9 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
 {
   string tmp_s, sps;
   ostringstream tmp_output, tmp1_output, global_output;
-  NodeID lhs = NULL, rhs = NULL;
+  expr_t lhs = NULL, rhs = NULL;
   BinaryOpNode *eq_node;
-  map<NodeID, int> reference_count;
+  map<expr_t, int> reference_count;
   temporary_terms_t local_temporary_terms;
   ofstream  output;
   int nze;
@@ -306,7 +306,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
           int equation_ID = getBlockEquationID(block, i);
           EquationType equ_type = getBlockEquationType(block, i);
           string sModel = symbol_table.getName(symbol_table.getID(eEndogenous, variable_ID));
-          eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+          eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
           lhs = eq_node->get_arg1();
           rhs = eq_node->get_arg2();
           tmp_output.str("");
@@ -334,7 +334,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
                       rhs->writeOutput(output, local_output_type, local_temporary_terms);
                       output << "\n  ";
                       tmp_output.str("");
-                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedNodeID(block, i);
+                      eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                       lhs = eq_node->get_arg1();
                       rhs = eq_node->get_arg2();
                       lhs->writeOutput(output, local_output_type, local_temporary_terms);
@@ -384,7 +384,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
               unsigned int var = it->first.second;
               unsigned int eqr = getBlockEquationID(block, eq);
               unsigned int varr = getBlockVariableID(block, var);
-              NodeID id = it->second.second;
+              expr_t id = it->second.second;
               output << "    g1(" << eq+1-block_recursive << ", " << var+1-block_recursive << ") = ";
               id->writeOutput(output, local_output_type, local_temporary_terms);
               output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
@@ -462,7 +462,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
       int deriv_id = it->first.second;
       if (getTypeByDerivID(deriv_id) == eEndogenous)
         {
-          NodeID d1 = it->second;
+          expr_t d1 = it->second;
           unsigned int eq = it->first.first;
           int symb = getSymbIDByDerivID(deriv_id);
           unsigned int var = symbol_table.getTypeSpecificID(symb);
@@ -528,10 +528,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
   string tmp_s;
   ostringstream tmp_output;
   ofstream code_file;
-  NodeID lhs = NULL, rhs = NULL;
+  expr_t lhs = NULL, rhs = NULL;
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
-  map<NodeID, int> reference_count;
+  map<expr_t, int> reference_count;
   vector<int> feedback_variables;
   bool file_open = false;
 
@@ -619,7 +619,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               }
               if (equ_type == E_EVALUATE)
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
                   lhs = eq_node->get_arg1();
                   rhs = eq_node->get_arg2();
                   rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
@@ -627,7 +627,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
                 }
               else if (equ_type == E_EVALUATE_S)
                 {
-                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedNodeID(block, i);
+                  eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
                   lhs = eq_node->get_arg1();
                   rhs = eq_node->get_arg2();
                   rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
@@ -647,7 +647,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
             end:
               FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
               fnumexpr.write(code_file);
-              eq_node = (BinaryOpNode *) getBlockEquationNodeID(block, i);
+              eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
               lhs = eq_node->get_arg1();
               rhs = eq_node->get_arg2();
               lhs->compile(code_file, false, temporary_terms, map_idx, false, false);
@@ -815,10 +815,10 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const st
   SaveCode.close();
 }
 
-map<pair<int, pair<int, int > >, NodeID>
+map<pair<int, pair<int, int > >, expr_t>
 StaticModel::collect_first_order_derivatives_endogenous()
 {
-  map<pair<int, pair<int, int > >, NodeID> endo_derivatives;
+  map<pair<int, pair<int, int > >, expr_t> endo_derivatives;
   for (first_derivatives_t::iterator it2 = first_derivatives.begin();
        it2 != first_derivatives.end(); it2++)
     {
@@ -868,7 +868,7 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
 
       computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered);
 
-      map<pair<int, pair<int, int> >, NodeID> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
+      map<pair<int, pair<int, int> >, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous();
 
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
@@ -953,7 +953,7 @@ StaticModel::writeStaticMFile(const string &func_name) const
     {
       int eq = it->first.first;
       int symb_id = it->first.second;
-      NodeID d1 = it->second;
+      expr_t d1 = it->second;
 
       output << "  g1(" << eq+1 << "," << symbol_table.getTypeSpecificID(symb_id)+1 << ")=";
       d1->writeOutput(output, oMatlabStaticModel, temporary_terms);
@@ -983,7 +983,7 @@ StaticModel::writeStaticMFile(const string &func_name) const
           int eq = it->first.first;
           int symb_id1 = it->first.second.first;
           int symb_id2 = it->first.second.second;
-          NodeID d2 = it->second;
+          expr_t d2 = it->second;
 
           int tsid1 = symbol_table.getTypeSpecificID(symb_id1);
           int tsid2 = symbol_table.getTypeSpecificID(symb_id2);
@@ -1197,7 +1197,7 @@ StaticModel::get_Derivatives(int block)
 void
 StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
 {
-  map<int, NodeID> recursive_variables;
+  map<int, expr_t> recursive_variables;
   unsigned int nb_blocks = getNbBlocks();
   blocks_derivatives = blocks_derivatives_t(nb_blocks);
   for (unsigned int block = 0; block < nb_blocks; block++)
@@ -1214,9 +1214,9 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
           for (int i = 0; i < block_nb_recursives; i++)
             {
               if (getBlockEquationType(block, i) == E_EVALUATE_S)
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
               else
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
             }
           map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> Derivatives = get_Derivatives(block);
           map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>::const_iterator it = Derivatives.begin();
@@ -1251,9 +1251,9 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
           for (int i = 0; i < block_nb_recursives; i++)
             {
               if (getBlockEquationType(block, i) == E_EVALUATE_S)
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
               else
-                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationNodeID(block, i);
+                recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
             }
           for (int eq = block_nb_recursives; eq < block_size; eq++)
             {
@@ -1261,7 +1261,7 @@ StaticModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
               for (int var = block_nb_recursives; var < block_size; var++)
                 {
                   int varr = getBlockVariableID(block, var);
-                  NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
+                  expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
                   if (d1 == Zero)
                     continue;
                   first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
@@ -1322,7 +1322,7 @@ StaticModel::writeChainRuleDerivative(ostream &output, int eqr, int varr, int la
                                       ExprNodeOutputType output_type,
                                       const temporary_terms_t &temporary_terms) const
 {
-  map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
+  map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
   if (it != first_chain_rule_derivatives.end())
     (it->second)->writeOutput(output, output_type, temporary_terms);
   else
diff --git a/StaticModel.hh b/StaticModel.hh
index 0d6f6a5ea3676478e6ca5031c2003ed5894fcea6..108ed0afccfe42cf43f8820823f4d454e0bbe0ae 100644
--- a/StaticModel.hh
+++ b/StaticModel.hh
@@ -44,7 +44,7 @@ private:
 
   vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
 
-  typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_t;
+  typedef map< pair< int, pair< int, int> >, expr_t> first_chain_rule_derivatives_t;
   first_chain_rule_derivatives_t first_chain_rule_derivatives;
 
   //! Writes static model file (standard Matlab version)
@@ -99,7 +99,7 @@ private:
   //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
   void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
   //! Collect only the first derivatives
-  map<pair<int, pair<int, int> >, NodeID> collect_first_order_derivatives_endogenous();
+  map<pair<int, pair<int, int> >, expr_t> collect_first_order_derivatives_endogenous();
 
   //! Helper for writing the Jacobian elements in MATLAB and C
   /*! Writes either (i+1,j+1) or [i+j*no_eq] */
@@ -122,7 +122,7 @@ protected:
   //! vector of block reordered variables and equations
   vector<int> equation_reordered, variable_reordered, inv_equation_reordered, inv_variable_reordered;
 
-  //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a NodeID on the new normalized equation
+  //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
   equation_type_and_normalized_equation_t equation_type_and_normalized_equation;
 
   //! for each block contains pair< Simulation_Type, pair < Block_Size, Recursive_part_Size > >
@@ -137,8 +137,8 @@ protected:
   //! Vector indicating if the block is linear in endogenous variable (true) or not (false)
   vector<bool> blocks_linear;
 
-  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), NodeID>
-  typedef map<pair< int, pair<int, int> >, NodeID> derivative_t;
+  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
+  typedef map<pair< int, pair<int, int> >, expr_t> derivative_t;
   //! Vector of derivative for each blocks
   vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
 
@@ -249,15 +249,15 @@ public:
   {
     return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
   };
-  //! Return the NodeID of the equation equation_number belonging to the block block_number
-  virtual NodeID
-  getBlockEquationNodeID(int block_number, int equation_number) const
+  //! Return the expr_t of the equation equation_number belonging to the block block_number
+  virtual expr_t
+  getBlockEquationExpr(int block_number, int equation_number) const
   {
     return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
   };
-  //! Return the NodeID of the renormalized equation equation_number belonging to the block block_number
-  virtual NodeID
-  getBlockEquationRenormalizedNodeID(int block_number, int equation_number) const
+  //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
+  virtual expr_t
+  getBlockEquationRenormalizedExpr(int block_number, int equation_number) const
   {
     return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
   };
diff --git a/SteadyStateModel.cc b/SteadyStateModel.cc
index 15437b03f6c8f5ee227c1debda92c24babff5a7b..4a916fb933a9a71369c1946da62c2304d208d13c 100644
--- a/SteadyStateModel.cc
+++ b/SteadyStateModel.cc
@@ -28,7 +28,7 @@ SteadyStateModel::SteadyStateModel(SymbolTable &symbol_table_arg, NumericalConst
 }
 
 void
-SteadyStateModel::addDefinition(int symb_id, NodeID expr)
+SteadyStateModel::addDefinition(int symb_id, expr_t expr)
 {
   assert(symbol_table.getType(symb_id) == eEndogenous
          || symbol_table.getType(symb_id) == eModFileLocalVariable
@@ -56,7 +56,7 @@ SteadyStateModel::checkPass(bool ramsey_policy) const
       if (!ramsey_policy)
         {
           set<pair<int, int> > used_symbols;
-          NodeID expr = def_table.find(*it)->second;
+          expr_t expr = def_table.find(*it)->second;
           expr->collectVariables(eEndogenous, used_symbols);
           expr->collectVariables(eModFileLocalVariable, used_symbols);
           for(set<pair<int, int> >::const_iterator it2 = used_symbols.begin();
@@ -101,7 +101,7 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_polic
   for(size_t i = 0; i < recursive_order.size(); i++)
     {
       output << "    ";
-      map<int, NodeID>::const_iterator it = def_table.find(recursive_order[i]);
+      map<int, expr_t>::const_iterator it = def_table.find(recursive_order[i]);
       it->second->writeOutput(output, oSteadyStateFile);
       output << ";" << endl;
     }
diff --git a/SteadyStateModel.hh b/SteadyStateModel.hh
index 2cd25a47105fe51632cccdc9cef2ed986517e0b0..d5677565fef33e9c0ad8a102771542ee78e48b05 100644
--- a/SteadyStateModel.hh
+++ b/SteadyStateModel.hh
@@ -27,7 +27,7 @@ class SteadyStateModel : public DataTree
 {
 private:
   //! Associates a symbol ID to an expression of the form "var = expr"
-  map<int, NodeID> def_table;
+  map<int, expr_t> def_table;
   vector<int> recursive_order;
 
   //! Reference to static model (for writing auxiliary equations)
@@ -36,7 +36,7 @@ private:
 public:
   SteadyStateModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants, ExternalFunctionsTable &external_functions_table_arg, const StaticModel &static_model_arg);
   //! Add an expression of the form "var = expr;"
-  void addDefinition(int symb_id, NodeID expr);
+  void addDefinition(int symb_id, expr_t expr);
   //! Checks that definitions are in a recursive order, and that no variable is declared twice
   /*!
     \param[in] ramsey_policy Is there a ramsey_policy statement in the MOD file? If yes, then disable the check on the recursivity of the declarations