diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh index 2a4e755f7..cdd6f2781 100644 --- a/preprocessor/CodeInterpreter.hh +++ b/preprocessor/CodeInterpreter.hh @@ -252,6 +252,19 @@ enum PriorDistributions eWeibull = 8 }; +enum NodeTreeReference + { + eResiduals = 0, + eFirstDeriv = 1, + eSecondDeriv = 2, + eThirdDeriv = 3, + eResidualsParamsDeriv = 4, + eJacobianParamsDeriv = 5, + eResidualsParamsSecondDeriv = 6, + eJacobianParamsSecondDeriv = 7, + eHessianParamsDeriv = 8 + }; + struct Block_contain_type { int Equation, Variable, Own_Derivative; diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index c48ee5e97..11ff41f0e 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -2112,26 +2112,34 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri void DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const { - ostringstream model_output; // Used for storing model - ostringstream model_eq_output; // Used for storing model equations - ostringstream jacobian_output; // Used for storing jacobian equations - ostringstream hessian_output; // Used for storing Hessian equations - ostringstream third_derivatives_output; + ostringstream model_local_vars_output; // Used for storing model local vars + 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; // Used for storing third order derivatives equations ExprNodeOutputType output_type = (use_dll ? oCDynamicModel : julia ? oJuliaDynamicModel : oMatlabDynamicModel); 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 hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr; // Writing Jacobian + temp_term_union.insert(temporary_terms_g1.begin(), temporary_terms_g1.end()); + 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++) { @@ -2141,11 +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, temporary_terms, tef_terms); + d1->writeOutput(jacobian_output, output_type, temp_term_union, tef_terms); jacobian_output << ";" << endl; } // Writing Hessian + temp_term_union.insert(temporary_terms_g2.begin(), temporary_terms_g2.end()); + 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++) @@ -2166,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, temporary_terms, tef_terms); + d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms); hessian_output << endl; } else @@ -2179,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, temporary_terms, tef_terms); + d2->writeOutput(hessian_output, output_type, temp_term_union, tef_terms); hessian_output << ";" << endl; k++; @@ -2208,6 +2222,12 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia } // Writing third derivatives + temp_term_union.insert(temporary_terms_g3.begin(), temporary_terms_g3.end()); + 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++) @@ -2230,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, temporary_terms, tef_terms); + d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms); third_derivatives_output << endl; } else @@ -2243,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, temporary_terms, tef_terms); + d3->writeOutput(third_derivatives_output, output_type, temp_term_union, tef_terms); third_derivatives_output << ";" << endl; } @@ -2287,8 +2307,8 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << "%" << endl << endl << "residual = zeros(" << nrows << ", 1);" << endl + << model_local_vars_output.str() << model_output.str() - << model_eq_output.str() // Writing initialization instruction for matrix g1 << "if nargout >= 2," << endl << " g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");" << endl @@ -2298,10 +2318,10 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " %" << endl << endl << jacobian_output.str() - << "end" << endl; + << endl // Initialize g2 matrix - DynamicOutput << "if nargout >= 3," << endl + << "if nargout >= 3," << endl << " %" << endl << " % Hessian matrix" << endl << " %" << endl @@ -2312,7 +2332,6 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << nrows << "," << hessianColsNbr << ");" << endl; else // Either hessian is all zero, or we didn't compute it DynamicOutput << " g2 = sparse([],[],[]," << nrows << "," << hessianColsNbr << ");" << endl; - DynamicOutput << "end" << endl; // Initialize g3 matrix DynamicOutput << "if nargout >= 4," << endl @@ -2328,7 +2347,9 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia else // Either 3rd derivatives is all zero, or we didn't compute it DynamicOutput << " g3 = sparse([],[],[]," << nrows << "," << ncols << ");" << endl; - DynamicOutput << "end" << endl; + DynamicOutput << "end" << endl + << "end" << endl + << "end" << endl; } else if (output_type == oCDynamicModel) { @@ -2337,33 +2358,30 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " 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 - << " else" << endl - << " {" << endl + << endl << jacobian_output.str() - << " }" << endl; + << endl; if (second_derivatives.size()) DynamicOutput << " /* Hessian for endogenous and exogenous variables */" << endl << " if (v2 == NULL)" << endl << " return;" << endl - << " else" << endl - << " {" << endl + << endl << hessian_output.str() - << " }" << endl; + << endl; if (third_derivatives.size()) DynamicOutput << " /* Third derivatives for endogenous and exogenous variables */" << endl << " if (v3 == NULL)" << endl << " return;" << endl - << " else" << endl - << " {" << endl + << endl << third_derivatives_output.str() - << " }" << endl; + << endl; DynamicOutput << "}" << endl << endl; } @@ -2396,8 +2414,8 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " #" << endl << " # Model equations" << endl << " #" << endl + << model_local_vars_output.str() << model_output.str() - << model_eq_output.str() << "end" << endl << endl << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, " << "params::Vector{Float64}," << endl @@ -2413,7 +2431,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " @assert size(g1) == (" << nrows << ", " << dynJacobianColsNbr << ")" << endl << " fill!(g1, 0.0)" << endl << " dynamic!(y, x, params, steady_state, it_, residual)" << endl - << model_output.str() + << model_local_vars_output.str() << " #" << endl << " # Jacobian matrix" << endl << " #" << endl @@ -2433,7 +2451,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl << " dynamic!(y, x, params, steady_state, it_, residual, g1)" << endl; if (second_derivatives.size()) - DynamicOutput << model_output.str() + DynamicOutput << model_local_vars_output.str() << " #" << endl << " # Hessian matrix" << endl << " #" << endl @@ -2456,7 +2474,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia << " @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl << " dynamic!(y, x, params, steady_state, it_, residual, g1, g2)" << endl; if (third_derivatives.size()) - DynamicOutput << model_output.str() + DynamicOutput << model_local_vars_output.str() << " #" << endl << " # Third order derivatives" << endl << " #" << endl diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index 8e60eea1d..057b780cc 100644 --- a/preprocessor/ExprNode.cc +++ b/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 &temp_terms_map, bool is_matlab) const { // For a terminal node, the cost is null return 0; @@ -110,9 +124,9 @@ ExprNode::collectExogenous(set > &result) const } void -ExprNode::computeTemporaryTerms(map &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +ExprNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { // Nothing to do for a terminal node } @@ -1621,16 +1635,31 @@ UnaryOpNode::computeDerivative(int deriv_id) return composeDerivatives(darg, deriv_id); } +int +UnaryOpNode::cost(const map &temp_terms_map, bool is_matlab) const +{ + // For a temporary term, the cost is null + for (map::const_iterator it = temp_terms_map.begin(); + it != temp_terms_map.end(); it++) + if (it->second.find(const_cast(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(this)); - if (it != temporary_terms.end()) + if (temporary_terms.find(const_cast(this)) != 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) @@ -1718,28 +1747,27 @@ 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 &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +UnaryOpNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { expr_t this2 = const_cast(this); - map::iterator it = reference_count.find(this2); + map >::iterator it = reference_count.find(this2); if (it == reference_count.end()) { - reference_count[this2] = 1; - arg->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + reference_count[this2] = make_pair(1, tr); + arg->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr); } else { - reference_count[this2]++; - if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)) - temporary_terms.insert(this2); + reference_count[this2] = make_pair(it->second.first + 1, it->second.second); + if (reference_count[this2].first * cost(temp_terms_map, is_matlab) > MIN_COST(is_matlab)) + temp_terms_map[reference_count[this2].second].insert(this2); } } @@ -2730,17 +2758,35 @@ BinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_t exit(EXIT_FAILURE); } +int +BinaryOpNode::cost(const map &temp_terms_map, bool is_matlab) const +{ + // For a temporary term, the cost is null + for (map::const_iterator it = temp_terms_map.begin(); + it != temp_terms_map.end(); it++) + if (it->second.find(const_cast(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(this)); // For a temporary term, the cost is null - if (it != temporary_terms.end()) + if (temporary_terms.find(const_cast(this)) != temporary_terms.end()) return 0; - int cost = arg1->cost(temporary_terms, is_matlab); - cost += arg2->cost(temporary_terms, is_matlab); + int arg_cost = arg1->cost(temporary_terms, is_matlab) + 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) @@ -2763,7 +2809,7 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con return cost + 990; case oPower: case oPowerDeriv: - return cost + 1160; + return cost + (MIN_COST_MATLAB/2+1); case oEqual: return cost; } @@ -2788,8 +2834,9 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con case oDivide: return cost + 15; case oPower: - case oPowerDeriv: return cost + 520; + case oPowerDeriv: + return cost + (MIN_COST_C/2+1);; case oEqual: return cost; } @@ -2798,29 +2845,29 @@ BinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) con } void -BinaryOpNode::computeTemporaryTerms(map &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +BinaryOpNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { expr_t this2 = const_cast(this); - map::iterator it = reference_count.find(this2); + map >::iterator it = reference_count.find(this2); if (it == reference_count.end()) { // If this node has never been encountered, set its ref count to one, // and travel through its children - reference_count[this2] = 1; - arg1->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); - arg2->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + reference_count[this2] = make_pair(1, tr); + arg1->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr); + arg2->computeTemporaryTerms(reference_count, temp_terms_map, 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]++; - if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab) + reference_count[this2] = make_pair(it->second.first + 1, it->second.second);; + if (reference_count[this2].first * cost(temp_terms_map, is_matlab) > MIN_COST(is_matlab) && op_code != oEqual) - temporary_terms.insert(this2); + temp_terms_map[reference_count[this2].second].insert(this2); } } @@ -3906,18 +3953,39 @@ TrinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_ exit(EXIT_FAILURE); } +int +TrinaryOpNode::cost(const map &temp_terms_map, bool is_matlab) const +{ + // For a temporary term, the cost is null + for (map::const_iterator it = temp_terms_map.begin(); + it != temp_terms_map.end(); it++) + if (it->second.find(const_cast(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(this)); // For a temporary term, the cost is null - if (it != temporary_terms.end()) + if (temporary_terms.find(const_cast(this)) != temporary_terms.end()) return 0; - int cost = arg1->cost(temporary_terms, is_matlab); - cost += arg2->cost(temporary_terms, is_matlab); - cost += arg3->cost(temporary_terms, is_matlab); + int arg_cost = arg1->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) @@ -3939,28 +4007,28 @@ TrinaryOpNode::cost(const temporary_terms_t &temporary_terms, bool is_matlab) co } void -TrinaryOpNode::computeTemporaryTerms(map &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +TrinaryOpNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { expr_t this2 = const_cast(this); - map::iterator it = reference_count.find(this2); + map >::iterator it = reference_count.find(this2); if (it == reference_count.end()) { // If this node has never been encountered, set its ref count to one, // and travel through its children - reference_count[this2] = 1; - arg1->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); - arg2->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); - arg3->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + reference_count[this2] = make_pair(1, tr); + arg1->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr); + arg2->computeTemporaryTerms(reference_count, temp_terms_map, is_matlab, tr); + arg3->computeTemporaryTerms(reference_count, temp_terms_map, 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]++; - if (reference_count[this2] * cost(temporary_terms, is_matlab) > MIN_COST(is_matlab)) - temporary_terms.insert(this2); + reference_count[this2] = make_pair(it->second.first + 1, it->second.second);; + if (reference_count[this2].first * cost(temp_terms_map, is_matlab) > MIN_COST(is_matlab)) + temp_terms_map[reference_count[this2].second].insert(this2); } } @@ -4771,11 +4839,11 @@ ExternalFunctionNode::composeDerivatives(const vector &dargs) } void -ExternalFunctionNode::computeTemporaryTerms(map &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +ExternalFunctionNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { - temporary_terms.insert(const_cast(this)); + temp_terms_map[tr].insert(const_cast(this)); } void @@ -5035,11 +5103,11 @@ FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatre } void -FirstDerivExternalFunctionNode::computeTemporaryTerms(map &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +FirstDerivExternalFunctionNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { - temporary_terms.insert(const_cast(this)); + temp_terms_map[tr].insert(const_cast(this)); } void @@ -5352,11 +5420,11 @@ SecondDerivExternalFunctionNode::SecondDerivExternalFunctionNode(DataTree &datat } void -SecondDerivExternalFunctionNode::computeTemporaryTerms(map &reference_count, - temporary_terms_t &temporary_terms, - bool is_matlab) const +SecondDerivExternalFunctionNode::computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const { - temporary_terms.insert(const_cast(this)); + temp_terms_map[tr].insert(const_cast(this)); } void diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh index 98f59800f..e46f567e4 100644 --- a/preprocessor/ExprNode.hh +++ b/preprocessor/ExprNode.hh @@ -157,7 +157,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 &temp_terms_map, bool is_matlab) const; //! For creating equation cross references struct EquationInfo @@ -193,7 +195,9 @@ 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 &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + 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 /*! @@ -516,7 +520,7 @@ public: 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; virtual void collectDynamicVariables(SymbolType type_arg, set > &result) const; - virtual void computeTemporaryTerms(map &reference_count, + virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, map > &first_occurence, int Curr_block, @@ -574,13 +578,17 @@ 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 &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 &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + 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; virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type, @@ -649,7 +657,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 &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; @@ -660,7 +670,9 @@ public: BinaryOpcode op_code_arg, const expr_t arg2_arg, int powerDerivOrder); virtual void prepareForDerivation(); virtual int precedence(ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const; - virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + 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; virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type, @@ -748,7 +760,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 &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: @@ -756,7 +770,9 @@ public: 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 &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + 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; virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type, @@ -829,7 +845,9 @@ public: AbstractExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg, const vector &arguments_arg); virtual void prepareForDerivation(); - virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const = 0; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + 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; virtual void writeExternalFunctionOutput(ostream &output, ExprNodeOutputType output_type, @@ -890,7 +908,9 @@ private: public: ExternalFunctionNode(DataTree &datatree_arg, int symb_id_arg, const vector &arguments_arg); - virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + 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, const temporary_terms_t &temporary_terms, @@ -922,7 +942,9 @@ public: int top_level_symb_id_arg, const vector &arguments_arg, int inputIndex_arg); - virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const; virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, map > &first_occurence, @@ -959,7 +981,9 @@ public: const vector &arguments_arg, int inputIndex1_arg, int inputIndex2_arg); - virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, bool is_matlab) const; + virtual void computeTemporaryTerms(map > &reference_count, + map &temp_terms_map, + bool is_matlab, NodeTreeReference tr) const; virtual void computeTemporaryTerms(map &reference_count, temporary_terms_t &temporary_terms, map > &first_occurence, diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index 72aa04889..2a34ddc4f 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -1189,24 +1189,50 @@ ModelTree::computeThirdDerivatives(const set &vars) void ModelTree::computeTemporaryTerms(bool is_matlab) { - map reference_count; + map > reference_count; temporary_terms.clear(); + temporary_terms_res.clear(); + temporary_terms_g1.clear(); + temporary_terms_g2.clear(); + temporary_terms_g3.clear(); + map 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::iterator it = equations.begin(); it != equations.end(); it++) - (*it)->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + (*it)->computeTemporaryTerms(reference_count, + temp_terms_map, + is_matlab, eResiduals); for (first_derivatives_t::iterator it = first_derivatives.begin(); it != first_derivatives.end(); it++) - it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + is_matlab, eFirstDeriv); for (second_derivatives_t::iterator it = second_derivatives.begin(); it != second_derivatives.end(); it++) - it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + is_matlab, eSecondDeriv); for (third_derivatives_t::iterator it = third_derivatives.begin(); it != third_derivatives.end(); it++) - it->second->computeTemporaryTerms(reference_count, temporary_terms, is_matlab); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + is_matlab, eThirdDeriv); + + for (map::const_iterator it = temp_terms_map.begin(); + it != temp_terms_map.end(); it++) + temporary_terms.insert(it->second.begin(), it->second.end()); + + 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 @@ -1215,7 +1241,6 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, ostream &output, { // Local var used to keep track of temp nodes already written temporary_terms_t tt2; - for (temporary_terms_t::const_iterator it = tt.begin(); it != tt.end(); it++) { @@ -1310,6 +1335,12 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t void ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const { + temporary_terms_t temp_terms; + if (IS_JULIA(output_type)) + temp_terms = temporary_terms_res; + else + temp_terms = temporary_terms; + for (int eq = 0; eq < (int) equations.size(); eq++) { BinaryOpNode *eq_node = equations[eq]; @@ -1331,13 +1362,13 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) if (IS_JULIA(output_type)) output << " @inbounds "; output << "lhs ="; - lhs->writeOutput(output, output_type, temporary_terms); + lhs->writeOutput(output, output_type, temp_terms); output << ";" << endl; if (IS_JULIA(output_type)) output << " @inbounds "; output << "rhs ="; - rhs->writeOutput(output, output_type, temporary_terms); + rhs->writeOutput(output, output_type, temp_terms); output << ";" << endl; if (IS_JULIA(output_type)) @@ -1355,7 +1386,7 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = "; - lhs->writeOutput(output, output_type, temporary_terms); + lhs->writeOutput(output, output_type, temp_terms); output << ";" << endl; } } @@ -1677,28 +1708,54 @@ ModelTree::computeParamsDerivatives() void ModelTree::computeParamsDerivativesTemporaryTerms() { - map reference_count; + map > reference_count; params_derivs_temporary_terms.clear(); + map 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, true); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + 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, true); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + 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, true); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + 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, true); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + 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, true); + it->second->computeTemporaryTerms(reference_count, + temp_terms_map, + true, eHessianParamsDeriv); + + for (map::const_iterator it = temp_terms_map.begin(); + it != temp_terms_map.end(); it++) + params_derivs_temporary_terms.insert(it->second.begin(), it->second.end()); + + 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 diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh index b68e3b3a0..259c3962a 100644 --- a/preprocessor/ModelTree.hh +++ b/preprocessor/ModelTree.hh @@ -137,9 +137,18 @@ protected: //! Temporary terms for the static/dynamic file (those which will be noted Txxxx) temporary_terms_t temporary_terms; + temporary_terms_t temporary_terms_res; + temporary_terms_t temporary_terms_g1; + temporary_terms_t temporary_terms_g2; + temporary_terms_t temporary_terms_g3; //! 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 diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 4fea0335e..92aa43f64 100644 --- a/preprocessor/StaticModel.cc +++ b/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_eq_output; // Used for storing model equations - ostringstream jacobian_output; // Used for storing jacobian equations - ostringstream hessian_output; // Used for storing Hessian equations + ostringstream model_local_vars_output; // Used for storing model local vars + ostringstream model_output; // Used for storing model + ostringstream jacobian_output; // Used for storing jacobian 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.begin(), temporary_terms_g1.end()); + 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.begin(), temporary_terms_g2.end()); + 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.begin(), temporary_terms_g3.end()); + 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 @@ -1372,7 +1392,6 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c << " if ~isreal(g1)" << endl << " g1 = real(g1)+2*imag(g1);" << endl << " end" << endl - << "end" << endl << "if nargout >= 3," << endl << " %" << endl << " % Hessian matrix" << endl @@ -1385,7 +1404,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c << " g2 = sparse(v2(:,1),v2(:,2),v2(:,3)," << equations.size() << "," << g2ncols << ");" << endl; else StaticOutput << " g2 = sparse([],[],[]," << equations.size() << "," << g2ncols << ");" << endl; - StaticOutput << "end" << endl; + // Initialize g3 matrix StaticOutput << "if nargout >= 4," << endl << " %" << endl @@ -1399,6 +1418,9 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c << " g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");" << endl; else // Either 3rd derivatives is all zero, or we didn't compute it StaticOutput << " g3 = sparse([],[],[]," << nrows << "," << ncols << ");" << endl; + StaticOutput << "end" << endl + << "end" << endl + << "end" << endl; } else if (output_type == oCStaticModel) { @@ -1407,32 +1429,29 @@ 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 - << " else" << endl - << " {" << endl + << endl << jacobian_output.str() - << " }" << endl; + << endl; if (second_derivatives.size()) StaticOutput << " /* Hessian for endogenous and exogenous variables */" << endl << " if (v2 == NULL)" << endl << " return;" << endl - << " else" << endl - << " {" << endl + << endl << hessian_output.str() - << " }" << endl; + << endl; if (third_derivatives.size()) StaticOutput << " /* Third derivatives for endogenous and exogenous variables */" << endl - << " if (v3 == NULL)" << endl - << " return;" << endl - << " else" << endl - << " {" << endl - << third_derivatives_output.str() - << " }" << endl; + << " if (v3 == NULL)" << endl + << " return;" << endl + << endl + << third_derivatives_output.str() + << endl; } else { @@ -1460,8 +1479,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 +1498,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 +1520,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 +1543,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