preprocessor: rework temporary terms

preprocessor/DynamicModel.cc

@@ -2116,21 +2116,17 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
   ostringstream model_output;             // Used for storing model temp vars and equations
   ostringstream jacobian_output;          // Used for storing jacobian equations
   ostringstream hessian_output;           // Used for storing Hessian equations
-  ostringstream third_derivatives_output;
+  ostringstream third_derivatives_output; // Used for storing third order derivatives equations
 
   ExprNodeOutputType output_type = (use_dll ? oCDynamicModel :
                                     julia ? oJuliaDynamicModel : oMatlabDynamicModel);
 
   deriv_node_temp_terms_t tef_terms;
-  temporary_terms_t temp_terms;
+  temporary_terms_t temp_term_union = temporary_terms_res;
-  if (julia)
-    temp_terms = temporary_terms_res;
-  else
-    temp_terms = temporary_terms;
 
   writeModelLocalVariables(model_local_vars_output, output_type, tef_terms);
 
-  writeTemporaryTerms(temp_terms, model_output, output_type, tef_terms);
+  writeTemporaryTerms(temporary_terms_res, model_output, output_type, tef_terms);
 
   writeModelEquations(model_output, output_type);
 
@@ -2138,12 +2134,12 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
   int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
 
   // Writing Jacobian
-  if (julia)
+  temp_term_union.insert(temporary_terms_g1.cbegin(), temporary_terms_g1.cend());
-    {
+  if (!first_derivatives.empty())
-      temp_terms = temporary_terms_g1;
+    if (julia)
-      if (!first_derivatives.empty())
+      writeTemporaryTerms(temp_term_union, jacobian_output, output_type, tef_terms);
-      writeTemporaryTerms(temp_terms, jacobian_output, output_type, tef_terms);
+    else
-    }
+      writeTemporaryTerms(temporary_terms_g1, jacobian_output, output_type, tef_terms);
   for (first_derivatives_t::const_iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     {
@@ -2153,17 +2149,17 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
 
       jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
       jacobian_output << "=";
-      d1->writeOutput(jacobian_output, output_type, temp_terms, tef_terms);
+      d1->writeOutput(jacobian_output, output_type, temp_term_union, tef_terms);
       jacobian_output << ";" << endl;
     }
 
   // Writing Hessian
-  if (julia)
+  temp_term_union.insert(temporary_terms_g2.cbegin(), temporary_terms_g2.cend());
-    {
+  if (!second_derivatives.empty())
-      temp_terms = temporary_terms_g2;
+    if (julia)
-      if (!second_derivatives.empty())
+      writeTemporaryTerms(temp_term_union, hessian_output, output_type, tef_terms);
-        writeTemporaryTerms(temp_terms, hessian_output, output_type, tef_terms);
+    else
-    }
+      writeTemporaryTerms(temporary_terms_g2, hessian_output, output_type, tef_terms);
   int k = 0; // Keep the line of a 2nd derivative in v2
   for (second_derivatives_t::const_iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
@@ -2184,7 +2180,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
         {
           for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
           hessian_output << "  @inbounds " << for_sym.str() << " = ";
-          d2->writeOutput(hessian_output, output_type, temp_terms, tef_terms);
+          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
           hessian_output << endl;
         }
       else
@@ -2197,7 +2193,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
 
           sparseHelper(2, hessian_output, k, 2, output_type);
           hessian_output << "=";
-          d2->writeOutput(hessian_output, output_type, temp_terms, tef_terms);
+          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
           hessian_output << ";" << endl;
 
           k++;
@@ -2226,12 +2222,12 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
     }
 
   // Writing third derivatives
-  if (julia)
+  temp_term_union.insert(temporary_terms_g3.cbegin(), temporary_terms_g3.cend());
-    {
+  if (!third_derivatives.empty())
-      temp_terms = temporary_terms_g3;
+    if (julia)
-      if (!third_derivatives.empty())
+      writeTemporaryTerms(temp_term_union, third_derivatives_output, output_type, tef_terms);
-        writeTemporaryTerms(temp_terms, third_derivatives_output, output_type, tef_terms);
+    else
-    }
+      writeTemporaryTerms(temporary_terms_g3, third_derivatives_output, output_type, tef_terms);
   k = 0; // Keep the line of a 3rd derivative in v3
   for (third_derivatives_t::const_iterator it = third_derivatives.begin();
        it != third_derivatives.end(); it++)
@@ -2254,7 +2250,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
         {
           for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
           third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
-          d3->writeOutput(third_derivatives_output, output_type, temp_terms, tef_terms);
+          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
           third_derivatives_output << endl;
         }
       else
@@ -2267,7 +2263,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
 
           sparseHelper(3, third_derivatives_output, k, 2, output_type);
           third_derivatives_output << "=";
-          d3->writeOutput(third_derivatives_output, output_type, temp_terms, tef_terms);
+          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
           third_derivatives_output << ";" << endl;
         }
 

preprocessor/ExprNode.cc

@@ -74,7 +74,21 @@ ExprNode::precedence(ExprNodeOutputType output_type, const temporary_terms_t &te
 }
 
 int
-ExprNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
+ExprNode::cost(int cost, bool is_matlab) const
+{
+  // For a terminal node, the cost is null
+  return 0;
+}
+
+int
+ExprNode::cost(const temporary_terms_t &temp_terms_map, bool is_matlab) const
+{
+  // For a terminal node, the cost is null
+  return 0;
+}
+
+int
+ExprNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
 {
   // For a terminal node, the cost is null
   return 0;
@@ -111,10 +125,7 @@ ExprNode::collectExogenous(set<pair<int, int> > &result) const
 
 void
 ExprNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                temporary_terms_t &temporary_terms_res,
+                                map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                temporary_terms_t &temporary_terms_g1,
-                                temporary_terms_t &temporary_terms_g2,
-                                temporary_terms_t &temporary_terms_g3,
                                 bool is_matlab, NodeTreeReference tr) const
 {
   // Nothing to do for a terminal node
@@ -1592,16 +1603,31 @@ UnaryOpNode::computeDerivative(int deriv_id)
   return composeDerivatives(darg, deriv_id);
 }
 
+int
+UnaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
+{
+  // For a temporary term, the cost is null
+  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.cbegin();
+       it != temp_terms_map.cend(); it++)
+    if (it->second.find(const_cast<UnaryOpNode *>(this)) != it->second.end())
+      return 0;
+
+  return cost(arg->cost(temp_terms_map, is_matlab), is_matlab);
+}
+
 int
 UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
 {
   // For a temporary term, the cost is null
-  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
+  if (temporary_terms.find(const_cast<UnaryOpNode *>(this)) != temporary_terms.end())
-  if (it != temporary_terms.end())
     return 0;
 
-  int cost = arg->cost(temporary_terms, is_matlab);
+  return cost(arg->cost(temporary_terms, is_matlab), is_matlab);
+}
 
+int
+UnaryOpNode::cost(int cost, bool is_matlab) const
+{
   if (is_matlab)
     // Cost for Matlab files
     switch (op_code)
@@ -1689,16 +1715,12 @@ UnaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) cons
       case oExpectation:
         return cost;
       }
-  // Suppress GCC warning
   exit(EXIT_FAILURE);
 }
 
 void
 UnaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                   temporary_terms_t &temporary_terms_res,
+                                   map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                   temporary_terms_t &temporary_terms_g1,
-                                   temporary_terms_t &temporary_terms_g2,
-                                   temporary_terms_t &temporary_terms_g3,
                                    bool is_matlab, NodeTreeReference tr) const
 {
   expr_t this2 = const_cast<UnaryOpNode *>(this);
@@ -1707,32 +1729,13 @@ UnaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &r
   if (it == reference_count.end())
     {
       reference_count[this2] = make_pair(1, tr);
-      arg->computeTemporaryTerms(reference_count,
+      arg->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
-                                 temporary_terms_res, temporary_terms_g1,
-                                 temporary_terms_g2, temporary_terms_g3,
-                                 is_matlab, tr);
     }
   else
     {
-      reference_count[this2] = make_pair(it->second.first++, it->second.second);
+      reference_count[this2] = make_pair(it->second.first + 1, it->second.second);
-      if (it->second.first * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
+      if (reference_count[this2].first * cost(temp_terms_map, is_matlab) > MIN_COST(is_matlab))
-        switch (it->second.second)
+        temp_terms_map[reference_count[this2].second].insert(this2);
-          {
-          case eResiduals:
-            temporary_terms_res.insert(this2);
-          case eFirstDeriv:
-            temporary_terms_g1.insert(this2);
-          case eSecondDeriv:
-            temporary_terms_g2.insert(this2);
-          case eThirdDeriv:
-            temporary_terms_g3.insert(this2);
-          case eResidualsParamsDeriv:
-          case eJacobianParamsDeriv:
-          case eResidualsParamsSecondDeriv:
-          case eJacobianParamsSecondDeriv:
-          case eHessianParamsDeriv:
-            temporary_terms_res.insert(this2);
-          }
     }
 }
 
@@ -2717,17 +2720,35 @@ BinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_t
   exit(EXIT_FAILURE);
 }
 
+int
+BinaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
+{
+  // For a temporary term, the cost is null
+  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.cbegin();
+       it != temp_terms_map.cend(); it++)
+    if (it->second.find(const_cast<BinaryOpNode *>(this)) != it->second.end())
+      return 0;
+
+  int arg_cost = arg1->cost(temp_terms_map, is_matlab) + arg2->cost(temp_terms_map, is_matlab);
+
+  return cost(arg_cost, is_matlab);
+}
+
 int
 BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
 {
-  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
   // For a temporary term, the cost is null
-  if (it != temporary_terms.end())
+  if (temporary_terms.find(const_cast<BinaryOpNode *>(this)) != temporary_terms.end())
     return 0;
 
-  int cost = arg1->cost(temporary_terms, is_matlab);
+  int arg_cost = arg1->cost(temporary_terms, is_matlab) + arg2->cost(temporary_terms, is_matlab);
-  cost += arg2->cost(temporary_terms, is_matlab);
 
+  return cost(arg_cost, is_matlab);
+}
+
+int
+BinaryOpNode::cost(int cost, bool is_matlab) const
+{
   if (is_matlab)
     // Cost for Matlab files
     switch (op_code)
@@ -2787,7 +2808,7 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con
 
 void
 BinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                    temporary_terms_t &temporary_terms,
+                                    map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                     bool is_matlab, NodeTreeReference tr) const
 {
   expr_t this2 = const_cast<BinaryOpNode *>(this);
@@ -2797,40 +2818,18 @@ BinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &
       // If this node has never been encountered, set its ref count to one,
       //  and travel through its children
       reference_count[this2] = make_pair(1, tr);
-      arg1->computeTemporaryTerms(reference_count,
+      arg1->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
-                                  temporary_terms_res, temporary_terms_g1,
+      arg2->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
-                                  temporary_terms_g2, temporary_terms_g3,
-                                  is_matlab, tr);
-      arg2->computeTemporaryTerms(reference_count,
-                                  temporary_terms_res, temporary_terms_g1,
-                                  temporary_terms_g2, temporary_terms_g3,
-                                  is_matlab, tr);
     }
   else
     {
       /* If the node has already been encountered, increment its ref count
          and declare it as a temporary term if it is too costly (except if it is
          an equal node: we don't want them as temporary terms) */
-      reference_count[this2] = make_pair(it->second.first++, it->second.second);;
+      reference_count[this2] = make_pair(it->second.first + 1, it->second.second);;
-      if (it->second.first * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)
+      if (reference_count[this2].first * cost(temp_terms_map, is_matlab) > MIN_COST(is_matlab)
           && op_code != oEqual)
-        switch (it->second.second)
+        temp_terms_map[reference_count[this2].second].insert(this2);
-          {
-          case eResiduals:
-            temporary_terms_res.insert(this2);
-          case eFirstDeriv:
-            temporary_terms_g1.insert(this2);
-          case eSecondDeriv:
-            temporary_terms_g2.insert(this2);
-          case eThirdDeriv:
-            temporary_terms_g3.insert(this2);
-          case eResidualsParamsDeriv:
-          case eJacobianParamsDeriv:
-          case eResidualsParamsSecondDeriv:
-          case eJacobianParamsSecondDeriv:
-          case eHessianParamsDeriv:
-            temporary_terms_res.insert(this2);
-          }
     }
 }
 
@@ -3909,17 +3908,39 @@ TrinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_
   exit(EXIT_FAILURE);
 }
 
+int
+TrinaryOpNode::cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const
+{
+  // For a temporary term, the cost is null
+  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.cbegin();
+       it != temp_terms_map.cend(); it++)
+    if (it->second.find(const_cast<TrinaryOpNode *>(this)) != it->second.end())
+      return 0;
+
+  int arg_cost = arg1->cost(temp_terms_map, is_matlab)
+    + arg2->cost(temp_terms_map, is_matlab)
+    + arg3->cost(temp_terms_map, is_matlab);
+
+  return cost(arg_cost, is_matlab);
+}
+
 int
 TrinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) const
 {
-  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
   // For a temporary term, the cost is null
-  if (it != temporary_terms.end())
+  if (temporary_terms.find(const_cast<TrinaryOpNode *>(this)) != temporary_terms.end())
     return 0;
 
-  int cost = arg1->cost(temporary_terms, is_matlab);
+  int arg_cost = arg1->cost(temporary_terms, is_matlab)
-  cost += arg2->cost(temporary_terms, is_matlab);
+    + arg2->cost(temporary_terms, is_matlab)
+    + arg3->cost(temporary_terms, is_matlab);
 
+  return cost(arg_cost, is_matlab);
+}
+
+int
+TrinaryOpNode::cost(int cost, bool is_matlab) const
+{
   if (is_matlab)
     // Cost for Matlab files
     switch (op_code)
@@ -3942,7 +3963,7 @@ TrinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) co
 
 void
 TrinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                      bool is_matlab, NodeTreeReference tr) const
 {
   expr_t this2 = const_cast<TrinaryOpNode *>(this);
@@ -3952,42 +3973,17 @@ TrinaryOpNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> >
       // If this node has never been encountered, set its ref count to one,
       //  and travel through its children
       reference_count[this2] = make_pair(1, tr);
-      arg1->computeTemporaryTerms(reference_count,
+      arg1->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
-                                  temporary_terms_res, temporary_terms_g1,
+      arg2->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
-                                  temporary_terms_g2, temporary_terms_g3,
+      arg3->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr);
-                                  is_matlab, tr);
-      arg2->computeTemporaryTerms(reference_count,
-                                  temporary_terms_res, temporary_terms_g1,
-                                  temporary_terms_g2, temporary_terms_g3,
-                                  is_matlab, tr);
-      arg3->computeTemporaryTerms(reference_count,
-                                  temporary_terms_res, temporary_terms_g1,
-                                  temporary_terms_g2, temporary_terms_g3,
-                                  is_matlab, tr);
     }
   else
     {
       // If the node has already been encountered, increment its ref count
       //  and declare it as a temporary term if it is too costly
-      reference_count[this2] = make_pair(it->second.first++, it->second.second);;
+      reference_count[this2] = make_pair(it->second.first + 1, it->second.second);;
-      if (it->second.first * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab))
+      if (reference_count[this2].first * cost(temp_terms_map, is_matlab) > MIN_COST(is_matlab))
-        switch (it->second.second)
+        temp_terms_map[reference_count[this2].second].insert(this2);
-          {
-          case eResiduals:
-            temporary_terms_res.insert(this2);
-          case eFirstDeriv:
-            temporary_terms_g1.insert(this2);
-          case eSecondDeriv:
-            temporary_terms_g2.insert(this2);
-          case eThirdDeriv:
-            temporary_terms_g3.insert(this2);
-          case eResidualsParamsDeriv:
-          case eJacobianParamsDeriv:
-          case eResidualsParamsSecondDeriv:
-          case eJacobianParamsSecondDeriv:
-          case eHessianParamsDeriv:
-            temporary_terms_res.insert(this2);
-          }
     }
 }
 
@@ -4791,10 +4787,10 @@ ExternalFunctionNode::composeDerivatives(const vector<expr_t> &dargs)
 
 void
 ExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                            temporary_terms_t &temporary_terms,
+                                            map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                             bool is_matlab, NodeTreeReference tr) const
 {
-  temporary_terms.insert(const_cast<ExternalFunctionNode *>(this));
+  temp_terms_map[tr].insert(const_cast<ExternalFunctionNode *>(this));
 }
 
 void
@@ -5046,10 +5042,10 @@ FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatre
 
 void
 FirstDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                                      temporary_terms_t &temporary_terms,
+                                                      map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                                       bool is_matlab, NodeTreeReference tr) const
 {
-  temporary_terms.insert(const_cast<FirstDerivExternalFunctionNode *>(this));
+  temp_terms_map[tr].insert(const_cast<FirstDerivExternalFunctionNode *>(this));
 }
 
 void
@@ -5354,10 +5350,10 @@ SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datat
 
 void
 SecondDerivExternalFunctionNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                                       temporary_terms_t &temporary_terms,
+                                                       map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
                                                        bool is_matlab, NodeTreeReference tr) const
 {
-  temporary_terms.insert(const_cast<SecondDerivExternalFunctionNode *>(this));
+  temp_terms_map[tr].insert(const_cast<SecondDerivExternalFunctionNode *>(this));
 }
 
 void

preprocessor/ExprNode.hh

@@ -159,7 +159,9 @@ protected:
 
   //! Cost of computing current node
   /*! Nodes included in temporary_terms are considered having a null cost */
+  virtual int cost(int cost, bool is_matlab) const;
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual int cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const;
 
 public:
   ExprNode(DataTree &datatree_arg);
@@ -187,10 +189,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<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const;
 
   //! Writes output of node, using a Txxx notation for nodes in temporary_terms, and specifiying the set of already written external functions
@@ -564,17 +563,16 @@ private:
   const int param1_symb_id, param2_symb_id;
   const UnaryOpcode op_code;
   virtual expr_t computeDerivative(int deriv_id);
+  virtual int cost(int cost, bool is_matlab) const;
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual int cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const;
   //! Returns the derivative of this node if darg is the derivative of the argument
   expr_t composeDerivatives(expr_t darg, int deriv_id);
 public:
   UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, int param1_symb_id_arg, int param2_symb_id_arg);
   virtual void prepareForDerivation();
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
@@ -643,7 +641,9 @@ private:
   const expr_t arg1, arg2;
   const BinaryOpcode op_code;
   virtual expr_t computeDerivative(int deriv_id);
+  virtual int cost(int cost, bool is_matlab) const;
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual int cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const;
   //! Returns the derivative of this node if darg1 and darg2 are the derivatives of the arguments
   expr_t composeDerivatives(expr_t darg1, expr_t darg2);
   const int powerDerivOrder;
@@ -655,10 +655,7 @@ public:
   virtual void prepareForDerivation();
   virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
@@ -746,7 +743,9 @@ private:
   const expr_t arg1, arg2, arg3;
   const TrinaryOpcode op_code;
   virtual expr_t computeDerivative(int deriv_id);
+  virtual int cost(int cost, bool is_matlab) const;
   virtual int cost(const temporary_terms_t &temporary_terms, bool is_matlab) const;
+  virtual int cost(const map<NodeTreeReference, temporary_terms_t> &temp_terms_map, bool is_matlab) const;
   //! Returns the derivative of this node if darg1, darg2 and darg3 are the derivatives of the arguments
   expr_t composeDerivatives(expr_t darg1, expr_t darg2, expr_t darg3);
 public:
@@ -755,10 +754,7 @@ public:
   virtual void prepareForDerivation();
   virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
   virtual bool containsExternalFunction() const;
@@ -832,10 +828,7 @@ public:
                                const vector<expr_t> &arguments_arg);
   virtual void prepareForDerivation();
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const = 0;
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const = 0;
   virtual bool containsExternalFunction() const;
@@ -897,10 +890,7 @@ public:
   ExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg,
                        const vector<expr_t> &arguments_arg);
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) 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,
@@ -933,10 +923,7 @@ public:
                                  const vector<expr_t> &arguments_arg,
                                  int inputIndex_arg);
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,
@@ -974,10 +961,7 @@ public:
                                   int inputIndex1_arg,
                                   int inputIndex2_arg);
   virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
-                                     temporary_terms_t &temporary_terms_res,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
-                                     temporary_terms_t &temporary_terms_g1,
-                                     temporary_terms_t &temporary_terms_g2,
-                                     temporary_terms_t &temporary_terms_g3,
                                      bool is_matlab, NodeTreeReference tr) const;
   virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
                                      temporary_terms_t &temporary_terms,

preprocessor/ModelTree.cc

@@ -1098,34 +1098,44 @@ ModelTree::computeTemporaryTerms(bool is_matlab)
   temporary_terms_g1.clear();
   temporary_terms_g2.clear();
   temporary_terms_g3.clear();
+  map<NodeTreeReference, temporary_terms_t> temp_terms_map;
+  temp_terms_map[eResiduals]=temporary_terms_res;
+  temp_terms_map[eFirstDeriv]=temporary_terms_g1;
+  temp_terms_map[eSecondDeriv]=temporary_terms_g2;
+  temp_terms_map[eThirdDeriv]=temporary_terms_g3;
 
   for (vector<BinaryOpNode *>::iterator it = equations.begin();
        it != equations.end(); it++)
     (*it)->computeTemporaryTerms(reference_count,
-                                 temporary_terms_res, temporary_terms_g1,
+                                 temp_terms_map,
-                                 temporary_terms_g2, temporary_terms_g3,
                                  is_matlab, eResiduals);
 
   for (first_derivatives_t::iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     it->second->computeTemporaryTerms(reference_count,
-                                      temporary_terms_res, temporary_terms_g1,
+                                      temp_terms_map,
-                                      temporary_terms_g2, temporary_terms_g3,
                                       is_matlab, eFirstDeriv);
 
   for (second_derivatives_t::iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
     it->second->computeTemporaryTerms(reference_count,
-                                     temporary_terms_res, temporary_terms_g1,
+                                      temp_terms_map,
-                                     temporary_terms_g2, temporary_terms_g3,
+                                      is_matlab, eSecondDeriv);
-                                     is_matlab, eSecondDeriv);
 
   for (third_derivatives_t::iterator it = third_derivatives.begin();
        it != third_derivatives.end(); it++)
     it->second->computeTemporaryTerms(reference_count,
-                                      temporary_terms_res, temporary_terms_g1,
+                                      temp_terms_map,
-                                      temporary_terms_g2, temporary_terms_g3,
                                       is_matlab, eThirdDeriv);
+
+  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
+       it != temp_terms_map.end(); it++)
+    temporary_terms.insert(it->second.cbegin(), it->second.cend());
+
+  temporary_terms_res = temp_terms_map[eResiduals];
+  temporary_terms_g1  = temp_terms_map[eFirstDeriv];
+  temporary_terms_g2  = temp_terms_map[eSecondDeriv];
+  temporary_terms_g3  = temp_terms_map[eThirdDeriv];
 }
 
 void
@@ -1603,41 +1613,52 @@ ModelTree::computeParamsDerivativesTemporaryTerms()
 {
   map<expr_t, pair<int, NodeTreeReference > > reference_count;
   params_derivs_temporary_terms.clear();
+  map<NodeTreeReference, temporary_terms_t> temp_terms_map;
+  temp_terms_map[eResidualsParamsDeriv]=params_derivs_temporary_terms_res;
+  temp_terms_map[eJacobianParamsDeriv]=params_derivs_temporary_terms_g1;
+  temp_terms_map[eResidualsParamsSecondDeriv]=params_derivs_temporary_terms_res2;
+  temp_terms_map[eJacobianParamsSecondDeriv]=params_derivs_temporary_terms_g12;
+  temp_terms_map[eHessianParamsDeriv]=params_derivs_temporary_terms_g2;
 
   for (first_derivatives_t::iterator it = residuals_params_derivatives.begin();
        it != residuals_params_derivatives.end(); it++)
     it->second->computeTemporaryTerms(reference_count,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      temp_terms_map,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
                                       true, eResidualsParamsDeriv);
 
   for (second_derivatives_t::iterator it = jacobian_params_derivatives.begin();
        it != jacobian_params_derivatives.end(); it++)
     it->second->computeTemporaryTerms(reference_count,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      temp_terms_map,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
                                       true, eJacobianParamsDeriv);
 
   for (second_derivatives_t::const_iterator it = residuals_params_second_derivatives.begin();
        it != residuals_params_second_derivatives.end(); ++it)
     it->second->computeTemporaryTerms(reference_count,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      temp_terms_map,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
                                       true, eResidualsParamsSecondDeriv);
 
   for (third_derivatives_t::const_iterator it = jacobian_params_second_derivatives.begin();
        it != jacobian_params_second_derivatives.end(); ++it)
     it->second->computeTemporaryTerms(reference_count,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      temp_terms_map,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
                                       true, eJacobianParamsSecondDeriv);
 
   for (third_derivatives_t::const_iterator it = hessian_params_derivatives.begin();
        it != hessian_params_derivatives.end(); ++it)
     it->second->computeTemporaryTerms(reference_count,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
+                                      temp_terms_map,
-                                      params_derivs_temporary_terms, params_derivs_temporary_terms,
                                       true, eHessianParamsDeriv);
+
+  for (map<NodeTreeReference, temporary_terms_t>::const_iterator it = temp_terms_map.begin();
+       it != temp_terms_map.end(); it++)
+    params_derivs_temporary_terms.insert(it->second.cbegin(), it->second.cend());
+
+  params_derivs_temporary_terms_res  = temp_terms_map[eResidualsParamsDeriv];
+  params_derivs_temporary_terms_g1   = temp_terms_map[eJacobianParamsDeriv];
+  params_derivs_temporary_terms_res2 = temp_terms_map[eResidualsParamsSecondDeriv];
+  params_derivs_temporary_terms_g12  = temp_terms_map[eJacobianParamsSecondDeriv];
+  params_derivs_temporary_terms_g2   = temp_terms_map[eHessianParamsDeriv];
 }
 
 bool ModelTree::isNonstationary(int symb_id) const

preprocessor/ModelTree.hh

@@ -137,6 +137,11 @@ protected:
 
   //! Temporary terms for the file containing parameters derivatives
   temporary_terms_t params_derivs_temporary_terms;
+  temporary_terms_t params_derivs_temporary_terms_res;
+  temporary_terms_t params_derivs_temporary_terms_g1;
+  temporary_terms_t params_derivs_temporary_terms_res2;
+  temporary_terms_t params_derivs_temporary_terms_g12;
+  temporary_terms_t params_derivs_temporary_terms_g2;
 
 
   //! Trend variables and their growth factors

preprocessor/StaticModel.cc

@@ -1183,27 +1183,35 @@ StaticModel::writeStaticMFile(const string &func_name) const
 void
 StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const
 {
-  ostringstream model_output;    // Used for storing model
+  ostringstream model_local_vars_output;   // Used for storing model local vars
-  ostringstream model_eq_output; // Used for storing model equations
+  ostringstream model_output;              // Used for storing model
-  ostringstream jacobian_output; // Used for storing jacobian equations
+  ostringstream jacobian_output;           // Used for storing jacobian equations
-  ostringstream hessian_output;  // Used for storing Hessian equations
+  ostringstream hessian_output;            // Used for storing Hessian equations
   ostringstream third_derivatives_output;  // Used for storing third order derivatives equations
   ostringstream for_sym;
   ExprNodeOutputType output_type = (use_dll ? oCStaticModel :
                                     julia ? oJuliaStaticModel : oMatlabStaticModel);
 
   deriv_node_temp_terms_t tef_terms;
-  writeModelLocalVariables(model_output, output_type, tef_terms);
+  temporary_terms_t temp_term_union = temporary_terms_res;
 
-  writeTemporaryTerms(temporary_terms, model_output, output_type, tef_terms);
+  writeModelLocalVariables(model_local_vars_output, output_type, tef_terms);
 
-  writeModelEquations(model_eq_output, output_type);
+  writeTemporaryTerms(temporary_terms_res, model_output, output_type, tef_terms);
+
+  writeModelEquations(model_output, output_type);
 
   int nrows = equations.size();
   int JacobianColsNbr = symbol_table.endo_nbr();
   int hessianColsNbr = JacobianColsNbr*JacobianColsNbr;
 
   // Write Jacobian w.r. to endogenous only
+  temp_term_union.insert(temporary_terms_g1.cbegin(), temporary_terms_g1.cend());
+  if (!first_derivatives.empty())
+    if (julia)
+      writeTemporaryTerms(temp_term_union, jacobian_output, output_type, tef_terms);
+    else
+      writeTemporaryTerms(temporary_terms_g1, jacobian_output, output_type, tef_terms);
   for (first_derivatives_t::const_iterator it = first_derivatives.begin();
        it != first_derivatives.end(); it++)
     {
@@ -1213,12 +1221,18 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
       jacobianHelper(jacobian_output, eq, symbol_table.getTypeSpecificID(symb_id), output_type);
       jacobian_output << "=";
-      d1->writeOutput(jacobian_output, output_type, temporary_terms, tef_terms);
+      d1->writeOutput(jacobian_output, output_type, temp_term_union, tef_terms);
       jacobian_output << ";" << endl;
     }
 
   int g2ncols = symbol_table.endo_nbr() * symbol_table.endo_nbr();
   // Write Hessian w.r. to endogenous only (only if 2nd order derivatives have been computed)
+  temp_term_union.insert(temporary_terms_g2.cbegin(), temporary_terms_g2.cend());
+  if (!second_derivatives.empty())
+    if (julia)
+      writeTemporaryTerms(temp_term_union, hessian_output, output_type, tef_terms);
+    else
+      writeTemporaryTerms(temporary_terms_g2, hessian_output, output_type, tef_terms);
   int k = 0; // Keep the line of a 2nd derivative in v2
   for (second_derivatives_t::const_iterator it = second_derivatives.begin();
        it != second_derivatives.end(); it++)
@@ -1238,7 +1252,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
         {
           for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
           hessian_output << "  @inbounds " << for_sym.str() << " = ";
-          d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
+          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
           hessian_output << endl;
         }
       else
@@ -1251,7 +1265,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
           sparseHelper(2, hessian_output, k, 2, output_type);
           hessian_output << "=";
-          d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
+          d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms);
           hessian_output << ";" << endl;
 
           k++;
@@ -1280,6 +1294,12 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
     }
 
   // Writing third derivatives
+  temp_term_union.insert(temporary_terms_g3.cbegin(), temporary_terms_g3.cend());
+  if (!third_derivatives.empty())
+    if (julia)
+      writeTemporaryTerms(temp_term_union, third_derivatives_output, output_type, tef_terms);
+    else
+      writeTemporaryTerms(temporary_terms_g3, third_derivatives_output, output_type, tef_terms);
   k = 0; // Keep the line of a 3rd derivative in v3
   for (third_derivatives_t::const_iterator it = third_derivatives.begin();
        it != third_derivatives.end(); it++)
@@ -1302,7 +1322,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
         {
           for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
           third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
-          d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
+          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
           third_derivatives_output << endl;
         }
       else
@@ -1315,7 +1335,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
           sparseHelper(3, third_derivatives_output, k, 2, output_type);
           third_derivatives_output << "=";
-          d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
+          d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms);
           third_derivatives_output << ";" << endl;
         }
 
@@ -1358,8 +1378,8 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << "%" << endl
                    << "% Model equations" << endl
                    << "%" << endl << endl
+                   << model_local_vars_output.str()
                    << model_output.str()
-                   << model_eq_output.str()
                    << "if ~isreal(residual)" << endl
                    << "  residual = real(residual)+imag(residual).^2;" << endl
                    << "end" << endl
@@ -1407,8 +1427,8 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << "  double lhs, rhs;" << endl
                    << endl
                    << "  /* Residual equations */" << endl
+                   << model_local_vars_output.str()
                    << model_output.str()
-                   << model_eq_output.str()
                    << "  /* Jacobian  */" << endl
                    << "  if (g1 == NULL)" << endl
                    << "    return;" << endl
@@ -1460,8 +1480,8 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << "  #" << endl
                    << "  # Model equations" << endl
                    << "  #" << endl
+                   << model_local_vars_output.str()
                    << model_output.str()
-                   << model_eq_output.str()
                    << "if ~isreal(residual)" << endl
                    << "  residual = real(residual)+imag(residual).^2;" << endl
                    << "end" << endl
@@ -1479,7 +1499,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << ")" << endl
                    << "  fill!(g1, 0.0)" << endl
                    << "  static!(y, x, params, residual)" << endl
-                   << model_output.str()
+                   << model_local_vars_output.str()
                    << "  #" << endl
                    << "  # Jacobian matrix" << endl
                    << "  #" << endl
@@ -1501,7 +1521,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << "  @assert size(g2) == (" << equations.size() << ", " << g2ncols << ")" << endl
                    << "  static!(y, x, params, residual, g1)" << endl;
       if (second_derivatives.size())
-        StaticOutput << model_output.str()
+        StaticOutput << model_local_vars_output.str()
                      << "  #" << endl
                      << "  # Hessian matrix" << endl
                      << "  #" << endl
@@ -1524,7 +1544,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
                    << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
                    << "  static!(y, x, params, residual, g1, g2)" << endl;
       if (third_derivatives.size())
-        StaticOutput << model_output.str()
+        StaticOutput << model_local_vars_output.str()
                      << "  #" << endl
                      << "  # Third order derivatives" << endl
                      << "  #" << endl