diff --git a/ComputingTasks.cc b/ComputingTasks.cc index f01bec39..0c0e6578 100644 --- a/ComputingTasks.cc +++ b/ComputingTasks.cc @@ -905,7 +905,7 @@ PlannerObjectiveStatement::getPlannerObjective() const void PlannerObjectiveStatement::computingPass() { - model_tree->computingPass(eval_context_t(), false, true, false, false); + model_tree->computingPass(eval_context_t(), false, true, false, false, false); } void diff --git a/DynamicModel.cc b/DynamicModel.cc index a753abc9..fffe1d63 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -2986,7 +2986,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative if (paramsDerivatives) { - cout << " - order 2 (derivatives of Jacobian w.r. to parameters)" << endl; + cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl; computeParamsDerivatives(); if (!no_tmp_terms) @@ -3715,108 +3715,6 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context } } -void -DynamicModel::computeParamsDerivatives() -{ - for (deriv_id_table_t::const_iterator it = deriv_id_table.begin(); - it != deriv_id_table.end(); it++) - { - if (symbol_table.getType(it->first.first) != eParameter) - continue; - - int param = it->second; - - for (int eq = 0; eq < (int) equations.size(); eq++) - { - expr_t d1 = equations[eq]->getDerivative(param); - if (d1 == Zero) - continue; - residuals_params_derivatives[make_pair(eq, param)] = d1; - } - - for (first_derivatives_t::const_iterator it2 = residuals_params_derivatives.begin(); - it2 != residuals_params_derivatives.end(); it2++) - { - int eq = it2->first.first; - int param1 = it2->first.second; - expr_t d1 = it2->second; - - expr_t d2 = d1->getDerivative(param); - if (d2 == Zero) - continue; - residuals_params_second_derivatives[make_pair(eq, make_pair(param1, param))] = d2; - } - - for (first_derivatives_t::const_iterator it2 = first_derivatives.begin(); - it2 != first_derivatives.end(); it2++) - { - int eq = it2->first.first; - int var = it2->first.second; - expr_t d1 = it2->second; - - expr_t d2 = d1->getDerivative(param); - if (d2 == Zero) - continue; - jacobian_params_derivatives[make_pair(eq, make_pair(var, param))] = d2; - } - - for (second_derivatives_t::const_iterator it2 = jacobian_params_derivatives.begin(); - it2 != jacobian_params_derivatives.end(); it2++) - { - int eq = it2->first.first; - int var = it2->first.second.first; - int param1 = it2->first.second.second; - expr_t d1 = it2->second; - - expr_t d2 = d1->getDerivative(param); - if (d2 == Zero) - continue; - jacobian_params_second_derivatives[make_pair(eq, make_pair(var, make_pair(param1, param)))] = d2; - } - - for (second_derivatives_t::const_iterator it2 = second_derivatives.begin(); - it2 != second_derivatives.end(); it2++) - { - int eq = it2->first.first; - int var1 = it2->first.second.first; - int var2 = it2->first.second.second; - expr_t d1 = it2->second; - - expr_t d2 = d1->getDerivative(param); - if (d2 == Zero) - continue; - hessian_params_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, param)))] = d2; - } - } -} - -void -DynamicModel::computeParamsDerivativesTemporaryTerms() -{ - map reference_count; - params_derivs_temporary_terms.clear(); - - 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); - - 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); - - 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); - - 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); - - 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); -} - void DynamicModel::writeParamsDerivativesFile(const string &basename) const { diff --git a/DynamicModel.hh b/DynamicModel.hh index 6f416f93..8bfde507 100644 --- a/DynamicModel.hh +++ b/DynamicModel.hh @@ -57,45 +57,6 @@ private: //! Number of columns of dynamic jacobian /*! Set by computeDerivID()s and computeDynJacobianCols() */ int dynJacobianColsNbr; - - //! Derivatives of the residuals w.r. to parameters - /*! First index is equation number, second is parameter. - Only non-null derivatives are stored in the map. - Parameter indices are those of the getDerivID() method. - */ - first_derivatives_t residuals_params_derivatives; - - //! Second derivatives of the residuals w.r. to parameters - /*! First index is equation number, second and third indeces are parameters. - Only non-null derivatives are stored in the map. - Parameter indices are those of the getDerivID() method. - */ - second_derivatives_t residuals_params_second_derivatives; - - //! Derivatives of the jacobian w.r. to parameters - /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter. - Only non-null derivatives are stored in the map. - Variable and parameter indices are those of the getDerivID() method. - */ - second_derivatives_t jacobian_params_derivatives; - - //! Second derivatives of the jacobian w.r. to parameters - /*! First index is equation number, second is endo/exo/exo_det variable, and third and fourth are parameters. - Only non-null derivatives are stored in the map. - Variable and parameter indices are those of the getDerivID() method. - */ - third_derivatives_t jacobian_params_second_derivatives; - - //! Derivatives of the hessian w.r. to parameters - /*! First index is equation number, first and second are endo/exo/exo_det variable, and third is parameter. - Only non-null derivatives are stored in the map. - Variable and parameter indices are those of the getDerivID() method. - */ - third_derivatives_t hessian_params_derivatives; - - //! Temporary terms for the file containing parameters derivatives - temporary_terms_t params_derivs_temporary_terms; - //! Temporary terms for block decomposed models vector< vector > v_temporary_terms; @@ -155,12 +116,8 @@ private: virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException); //! Compute the column indices of the dynamic Jacobian void computeDynJacobianCols(bool jacobianExo); - //! Computes derivatives of the Jacobian w.r. to parameters - void computeParamsDerivatives(); //! Computes derivatives of the Jacobian w.r. to trend vars and tests that they are equal to zero void testTrendDerivativesEqualToZero(const eval_context_t &eval_context); - //! Computes temporary terms for the file containing parameters derivatives - void computeParamsDerivativesTemporaryTerms(); //! Collect only the first derivatives map >, expr_t> collect_first_order_derivatives_endogenous(); diff --git a/ModFile.cc b/ModFile.cc index bd85f766..32fb91dc 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -380,8 +380,10 @@ ModFile::computingPass(bool no_tmp_terms) const bool static_hessian = mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation; + const bool paramsDerivatives = mod_file_struct.identification_present + || mod_file_struct.estimation_analytic_derivation; static_model.computingPass(global_eval_context, no_tmp_terms, static_hessian, - block, byte_code); + paramsDerivatives, block, byte_code); } // Set things to compute for dynamic model if (mod_file_struct.simul_present || mod_file_struct.check_present @@ -640,7 +642,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool no_log, b if (dynamic_model.equation_number() > 0) { if (!no_static) - static_model.writeStaticFile(basename, block, byte_code, use_dll); + { + static_model.writeStaticFile(basename, block, byte_code, use_dll); + static_model.writeParamsDerivativesFile(basename); + } dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option); dynamic_model.writeParamsDerivativesFile(basename); diff --git a/ModelTree.cc b/ModelTree.cc index 6ada1a47..fb4de128 100644 --- a/ModelTree.cc +++ b/ModelTree.cc @@ -1460,3 +1460,105 @@ ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, Expr output << row_nb + col_nb * NNZDerivatives[order-1]; output << RIGHT_ARRAY_SUBSCRIPT(output_type); } + +void +ModelTree::computeParamsDerivatives() +{ + set deriv_id_set; + addAllParamDerivId(deriv_id_set); + + for (set::const_iterator it = deriv_id_set.begin(); + it != deriv_id_set.end(); it++) + { + const int param = *it; + + for (int eq = 0; eq < (int) equations.size(); eq++) + { + expr_t d1 = equations[eq]->getDerivative(param); + if (d1 == Zero) + continue; + residuals_params_derivatives[make_pair(eq, param)] = d1; + } + + for (first_derivatives_t::const_iterator it2 = residuals_params_derivatives.begin(); + it2 != residuals_params_derivatives.end(); it2++) + { + int eq = it2->first.first; + int param1 = it2->first.second; + expr_t d1 = it2->second; + + expr_t d2 = d1->getDerivative(param); + if (d2 == Zero) + continue; + residuals_params_second_derivatives[make_pair(eq, make_pair(param1, param))] = d2; + } + + for (first_derivatives_t::const_iterator it2 = first_derivatives.begin(); + it2 != first_derivatives.end(); it2++) + { + int eq = it2->first.first; + int var = it2->first.second; + expr_t d1 = it2->second; + + expr_t d2 = d1->getDerivative(param); + if (d2 == Zero) + continue; + jacobian_params_derivatives[make_pair(eq, make_pair(var, param))] = d2; + } + + for (second_derivatives_t::const_iterator it2 = jacobian_params_derivatives.begin(); + it2 != jacobian_params_derivatives.end(); it2++) + { + int eq = it2->first.first; + int var = it2->first.second.first; + int param1 = it2->first.second.second; + expr_t d1 = it2->second; + + expr_t d2 = d1->getDerivative(param); + if (d2 == Zero) + continue; + jacobian_params_second_derivatives[make_pair(eq, make_pair(var, make_pair(param1, param)))] = d2; + } + + for (second_derivatives_t::const_iterator it2 = second_derivatives.begin(); + it2 != second_derivatives.end(); it2++) + { + int eq = it2->first.first; + int var1 = it2->first.second.first; + int var2 = it2->first.second.second; + expr_t d1 = it2->second; + + expr_t d2 = d1->getDerivative(param); + if (d2 == Zero) + continue; + hessian_params_derivatives[make_pair(eq, make_pair(var1, make_pair(var2, param)))] = d2; + } + } +} + +void +ModelTree::computeParamsDerivativesTemporaryTerms() +{ + map reference_count; + params_derivs_temporary_terms.clear(); + + 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); + + 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); + + 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); + + 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); + + 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); +} diff --git a/ModelTree.hh b/ModelTree.hh index 37413ae6..079c5564 100644 --- a/ModelTree.hh +++ b/ModelTree.hh @@ -92,8 +92,49 @@ protected: */ third_derivatives_t third_derivatives; - //! Temporary terms (those which will be noted Txxxx) + //! Derivatives of the residuals w.r. to parameters + /*! First index is equation number, second is parameter. + Only non-null derivatives are stored in the map. + Parameter indices are those of the getDerivID() method. + */ + first_derivatives_t residuals_params_derivatives; + + //! Second derivatives of the residuals w.r. to parameters + /*! First index is equation number, second and third indeces are parameters. + Only non-null derivatives are stored in the map. + Parameter indices are those of the getDerivID() method. + */ + second_derivatives_t residuals_params_second_derivatives; + + //! Derivatives of the jacobian w.r. to parameters + /*! First index is equation number, second is endo/exo/exo_det variable, and third is parameter. + Only non-null derivatives are stored in the map. + Variable and parameter indices are those of the getDerivID() method. + */ + second_derivatives_t jacobian_params_derivatives; + + //! Second derivatives of the jacobian w.r. to parameters + /*! First index is equation number, second is endo/exo/exo_det variable, and third and fourth are parameters. + Only non-null derivatives are stored in the map. + Variable and parameter indices are those of the getDerivID() method. + */ + third_derivatives_t jacobian_params_second_derivatives; + + //! Derivatives of the hessian w.r. to parameters + /*! First index is equation number, first and second are endo/exo/exo_det variable, and third is parameter. + Only non-null derivatives are stored in the map. + Variable and parameter indices are those of the getDerivID() method. + */ + third_derivatives_t hessian_params_derivatives; + + + //! Temporary terms for the static/dynamic file (those which will be noted Txxxx) temporary_terms_t temporary_terms; + + //! Temporary terms for the file containing parameters derivatives + temporary_terms_t params_derivs_temporary_terms; + + //! Trend variables and their growth factors trend_symbols_map_t trend_symbols_map; //! Nonstationary variables and their deflators @@ -114,12 +155,16 @@ protected: //! Computes 3rd derivatives /*! \param vars the derivation IDs w.r. to which derive the 2nd derivatives */ void computeThirdDerivatives(const set &vars); + //! Computes derivatives of the Jacobian and Hessian w.r. to parameters + void computeParamsDerivatives(); //! Write derivative of an equation w.r. to a variable void writeDerivative(ostream &output, int eq, int symb_id, int lag, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const; //! Computes temporary terms (for all equations and derivatives) void computeTemporaryTerms(bool is_matlab); - //! Writes temporary terms + //! Computes temporary terms for the file containing parameters derivatives + void computeParamsDerivativesTemporaryTerms(); +//! Writes temporary terms void writeTemporaryTerms(const temporary_terms_t &tt, ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const; //! Compiles temporary terms void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const; diff --git a/StaticModel.cc b/StaticModel.cc index 66a4366c..181bfbc8 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -1047,7 +1047,7 @@ StaticModel::collect_first_order_derivatives_endogenous() } void -StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode) +StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool paramsDerivatives, bool block, bool bytecode) { initializeVariablesAndEquations(); @@ -1070,6 +1070,15 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms computeHessian(vars); } +if (paramsDerivatives) + { + cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl; + computeParamsDerivatives(); + + if (!no_tmp_terms) + computeParamsDerivativesTemporaryTerms(); + } + if (block) { jacob_map_t contemporaneous_jacobian, static_jacobian; @@ -1542,7 +1551,12 @@ StaticModel::writeOutput(ostream &output, bool block) const SymbolType StaticModel::getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException) { - return symbol_table.getType(getSymbIDByDerivID(deriv_id)); + if (deriv_id < symbol_table.endo_nbr()) + return eEndogenous; + else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr()) + return eParameter; + else + throw UnknownDerivIDException(); } int @@ -1554,18 +1568,32 @@ StaticModel::getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException) int StaticModel::getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException) { - return deriv_id; + if (deriv_id < symbol_table.endo_nbr()) + return symbol_table.getID(eEndogenous, deriv_id); + else if (deriv_id < symbol_table.endo_nbr() + symbol_table.param_nbr()) + return symbol_table.getID(eParameter, deriv_id - symbol_table.endo_nbr()); + else + throw UnknownDerivIDException(); } int StaticModel::getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException) { if (symbol_table.getType(symb_id) == eEndogenous) - return symb_id; + return symbol_table.getTypeSpecificID(symb_id); + else if (symbol_table.getType(symb_id) == eParameter) + return symbol_table.getTypeSpecificID(symb_id) + symbol_table.endo_nbr(); else return -1; } +void +StaticModel::addAllParamDerivId(set &deriv_id_set) +{ + for (int i = 0; i < symbol_table.param_nbr(); i++) + deriv_id_set.insert(i + symbol_table.endo_nbr()); +} + map >, pair >, int> StaticModel::get_Derivatives(int block) { @@ -1794,3 +1822,159 @@ void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename) const output << ";" << endl; } } + +void +StaticModel::writeParamsDerivativesFile(const string &basename) const +{ + if (!residuals_params_derivatives.size() + && !residuals_params_second_derivatives.size() + && !jacobian_params_derivatives.size() + && !jacobian_params_second_derivatives.size() + && !hessian_params_derivatives.size()) + return; + + string filename = basename + "_static_params_derivs.m"; + + ofstream paramsDerivsFile; + paramsDerivsFile.open(filename.c_str(), ios::out | ios::binary); + if (!paramsDerivsFile.is_open()) + { + cerr << "ERROR: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_static_params_derivs(y, x, params)" << endl + << "%" << endl + << "% Warning : this file is generated automatically by Dynare" << endl + << "% from model file (.mod)" << endl << endl; + + deriv_node_temp_terms_t tef_terms; + writeModelLocalVariables(paramsDerivsFile, oMatlabStaticModel, tef_terms); + + writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, oMatlabStaticModel, tef_terms); + + // Write parameter derivative + paramsDerivsFile << "rp = zeros(" << equation_number() << ", " + << symbol_table.param_nbr() << ");" << endl; + + for (first_derivatives_t::const_iterator it = residuals_params_derivatives.begin(); + it != residuals_params_derivatives.end(); it++) + { + int eq = it->first.first; + int param = it->first.second; + expr_t d1 = it->second; + + int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; + + paramsDerivsFile << "rp(" << eq+1 << ", " << param_col << ") = "; + d1->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms); + paramsDerivsFile << ";" << endl; + } + + // Write jacobian derivatives + paramsDerivsFile << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", " + << symbol_table.param_nbr() << ");" << endl; + + for (second_derivatives_t::const_iterator it = jacobian_params_derivatives.begin(); + it != jacobian_params_derivatives.end(); it++) + { + int eq = it->first.first; + int var = it->first.second.first; + int param = it->first.second.second; + expr_t d2 = it->second; + + int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1; + int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; + + paramsDerivsFile << "gp(" << eq+1 << ", " << var_col << ", " << param_col << ") = "; + d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms); + paramsDerivsFile << ";" << endl; + } + + // If nargout >= 3... + paramsDerivsFile << "if nargout >= 3" << endl; + + // Write parameter second derivatives (only if nargout >= 3) + paramsDerivsFile << "rpp = zeros(" << residuals_params_second_derivatives.size() + << ",4);" << endl; + + int i = 1; + for (second_derivatives_t::const_iterator it = residuals_params_second_derivatives.begin(); + it != residuals_params_second_derivatives.end(); ++it, i++) + { + int eq = it->first.first; + int param1 = it->first.second.first; + int param2 = it->first.second.second; + expr_t d2 = it->second; + + int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1; + int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1; + + paramsDerivsFile << "rpp(" << i << ",1)=" << eq+1 << ";" << endl + << "rpp(" << i << ",2)=" << param1_col << ";" << endl + << "rpp(" << i << ",3)=" << param2_col << ";" << endl + << "rpp(" << i << ",4)="; + d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms); + paramsDerivsFile << ";" << endl; + } + + // Write jacobian second derivatives (only if nargout >= 3) + paramsDerivsFile << "gpp = zeros(" << jacobian_params_second_derivatives.size() + << ",5);" << endl; + + i = 1; + for (third_derivatives_t::const_iterator it = jacobian_params_second_derivatives.begin(); + it != jacobian_params_second_derivatives.end(); ++it, i++) + { + int eq = it->first.first; + int var = it->first.second.first; + int param1 = it->first.second.second.first; + int param2 = it->first.second.second.second; + expr_t d2 = it->second; + + int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1; + int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1; + int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1; + + paramsDerivsFile << "gpp(" << i << ",1)=" << eq+1 << ";" << endl + << "gpp(" << i << ",2)=" << var_col << ";" << endl + << "gpp(" << i << ",3)=" << param1_col << ";" << endl + << "gpp(" << i << ",4)=" << param2_col << ";" << endl + << "gpp(" << i << ",5)="; + d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms); + paramsDerivsFile << ";" << endl; + } + + // If nargout >= 5... + paramsDerivsFile << "end" << endl + << "if nargout >= 5" << endl; + + // Write hessian derivatives (only if nargout >= 5) + paramsDerivsFile << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl; + + i = 1; + for (third_derivatives_t::const_iterator it = hessian_params_derivatives.begin(); + it != hessian_params_derivatives.end(); ++it, i++) + { + int eq = it->first.first; + int var1 = it->first.second.first; + int var2 = it->first.second.second.first; + int param = it->first.second.second.second; + expr_t d2 = it->second; + + int var1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var1)) + 1; + int var2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var2)) + 1; + int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; + + paramsDerivsFile << "hp(" << i << ",1)=" << eq+1 << ";" << endl + << "hp(" << i << ",2)=" << var1_col << ";" << endl + << "hp(" << i << ",3)=" << var2_col << ";" << endl + << "hp(" << i << ",4)=" << param_col << ";" << endl + << "hp(" << i << ",5)="; + d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms); + paramsDerivsFile << ";" << endl; + } + + paramsDerivsFile << "end" << endl + << "end" << endl; + paramsDerivsFile.close(); +} diff --git a/StaticModel.hh b/StaticModel.hh index de80b457..260d7977 100644 --- a/StaticModel.hh +++ b/StaticModel.hh @@ -165,8 +165,9 @@ public: \param eval_context evaluation context for normalization \param no_tmp_terms if true, no temporary terms will be computed in the static files \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed + \param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed */ - void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode); + void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian,bool paramsDerivatives, bool block, bool bytecode); //! Adds informations for simulation in a binary file for a block decomposed model void Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num, @@ -175,6 +176,9 @@ public: //! Writes static model file void writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll) const; + //! Writes file containing static parameters derivatives + void writeParamsDerivativesFile(const string &basename) const; + //! Writes LaTeX file with the equations of the static model void writeLatexFile(const string &basename) const; @@ -185,6 +189,7 @@ public: void writeAuxVarRecursiveDefinitions(const string &basename) const; virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); + virtual void addAllParamDerivId(set &deriv_id_set); //! Return the number of blocks virtual unsigned int