diff --git a/src/DataTree.cc b/src/DataTree.cc index ec0590d4..363a1647 100644 --- a/src/DataTree.cc +++ b/src/DataTree.cc @@ -832,12 +832,6 @@ DataTree::addAllParamDerivId([[maybe_unused]] set &deriv_id_set) { } -int -DataTree::getDynJacobianCol([[maybe_unused]] int deriv_id) const noexcept(false) -{ - throw UnknownDerivIDException(); -} - bool DataTree::isUnaryOpUsed(UnaryOpcode opcode) const { diff --git a/src/DataTree.hh b/src/DataTree.hh index e4e36f04..d8e9439e 100644 --- a/src/DataTree.hh +++ b/src/DataTree.hh @@ -305,8 +305,21 @@ public: virtual SymbolType getTypeByDerivID(int deriv_id) const noexcept(false); virtual int getLagByDerivID(int deriv_id) const noexcept(false); virtual int getSymbIDByDerivID(int deriv_id) const noexcept(false); - //! Returns the column of the dynamic Jacobian associated to a derivation ID - virtual int getDynJacobianCol(int deriv_id) const noexcept(false); + + //! Returns the column of the Jacobian associated to a derivation ID + virtual int + getJacobianCol([[maybe_unused]] int deriv_id) const + { + throw UnknownDerivIDException(); + } + + //! Returns the number of columns of the Jacobian + virtual int + getJacobianColsNbr() const + { + throw UnknownDerivIDException(); + } + //! Adds to the set all the deriv IDs corresponding to parameters virtual void addAllParamDerivId(set &deriv_id_set); diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index d9a46304..585b82d2 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -96,7 +96,6 @@ DynamicModel::DynamicModel(const DynamicModel &m) : xref_exo{m.xref_exo}, xref_exo_det{m.xref_exo_det}, nonzero_hessian_eqs{m.nonzero_hessian_eqs}, - dynJacobianColsNbr{m.dynJacobianColsNbr}, variableMapping{m.variableMapping}, blocks_other_endo{m.blocks_other_endo}, blocks_exo{m.blocks_exo}, @@ -147,7 +146,6 @@ DynamicModel::operator=(const DynamicModel &m) xref_exo = m.xref_exo; xref_exo_det = m.xref_exo_det; nonzero_hessian_eqs = m.nonzero_hessian_eqs; - dynJacobianColsNbr = m.dynJacobianColsNbr; variableMapping = m.variableMapping; blocks_derivatives_other_endo.clear(); blocks_derivatives_exo.clear(); @@ -1268,13 +1266,280 @@ DynamicModel::writeDynamicBlockBytecode(const string &basename) const void DynamicModel::writeDynamicMFile(const string &basename) const { - writeDynamicModel(basename, false, false); + auto [d_output, tt_output] = writeModelFileHelper(); + + // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201 + map tmp_paren_vars; + bool message_printed{false}; + for (auto &it : tt_output) + fixNestedParenthesis(it, tmp_paren_vars, message_printed); + for (auto &it : d_output) + fixNestedParenthesis(it, tmp_paren_vars, message_printed); + + ostringstream init_output, end_output; + init_output << "residual = zeros(" << equations.size() << ", 1);"; + writeDynamicModelHelper(basename, "dynamic_resid", "residual", + "dynamic_resid_tt", + temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(), + "", init_output, end_output, + d_output[0], tt_output[0]); + + init_output.str(""); + init_output << "g1 = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ");"; + writeDynamicModelHelper(basename, "dynamic_g1", "g1", + "dynamic_g1_tt", + temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(), + "dynamic_resid_tt", + init_output, end_output, + d_output[1], tt_output[1]); + writeWrapperFunctions(basename, "g1"); + + // For order ≥ 2 + int ncols{getJacobianColsNbr()}; + int ntt{static_cast(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())}; + for (size_t i{2}; i < derivatives.size(); i++) + { + ncols *= getJacobianColsNbr(); + ntt += temporary_terms_derivatives[i].size(); + string gname{"g" + to_string(i)}; + string gprevname{"g" + to_string(i-1)}; + + init_output.str(""); + end_output.str(""); + if (derivatives[i].size()) + { + init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl + << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl + << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl; + end_output << gname << " = sparse(" + << gname << "_i," << gname << "_j," << gname << "_v," + << equations.size() << "," << ncols << ");"; + } + else + init_output << gname << " = sparse([],[],[]," << equations.size() << "," << ncols << ");"; + writeDynamicModelHelper(basename, "dynamic_" + gname, gname, + "dynamic_" + gname + "_tt", + ntt, + "dynamic_" + gprevname + "_tt", + init_output, end_output, + d_output[i], tt_output[i]); + if (i <= 3) + writeWrapperFunctions(basename, gname); + } + + writeDynamicMatlabCompatLayer(basename); } void DynamicModel::writeDynamicJuliaFile(const string &basename) const { - writeDynamicModel(basename, false, true); + auto [d_output, tt_output] = writeModelFileHelper(); + + stringstream output; + + output << "module " << basename << "Dynamic" << endl + << "#" << endl + << "# NB: this file was automatically generated by Dynare" << endl + << "# from " << basename << ".mod" << endl + << "#" << endl + << "using StatsFuns" << endl << endl + << "export tmp_nbr, dynamic!, dynamicResid!, dynamicG1!, dynamicG2!, dynamicG3!" << endl << endl + << "#=" << endl + << "# The comments below apply to all functions contained in this module #" << endl + << " NB: The arguments contained on the first line of the function" << endl + << " definition are those that are modified in place" << endl << endl + << "## Exported Functions ##" << endl + << " dynamic! : Wrapper function; computes residuals, Jacobian, Hessian," << endl + << " and third derivatives depending on the arguments provided" << endl + << " dynamicResid! : Computes the dynamic model residuals" << endl + << " dynamicG1! : Computes the dynamic model Jacobian" << endl + << " dynamicG2! : Computes the dynamic model Hessian" << endl + << " dynamicG3! : Computes the dynamic model third derivatives" << endl << endl + << "## Exported Variables ##" << endl + << " tmp_nbr : Vector{Int}(4) respectively the number of temporary variables" << endl + << " for the residuals, g1, g2 and g3." << endl << endl + << "## Local Functions ##" << endl + << " dynamicResidTT! : Computes the dynamic model temporary terms for the residuals" << endl + << " dynamicG1TT! : Computes the dynamic model temporary terms for the Jacobian" << endl + << " dynamicG2TT! : Computes the dynamic model temporary terms for the Hessian" << endl + << " dynamicG3TT! : Computes the dynamic model temporary terms for the third derivatives" << endl << endl + << "## Function Arguments ##" << endl + << " T : Vector{Float64}(num_temp_terms), temporary terms" << endl + << " y : Vector{Float64}(num_dynamic_vars), endogenous variables in the order stored model_.lead_lag_incidence; see the manual" << endl + << " x : Matrix{Float64}(nperiods,model_.exo_nbr), exogenous variables (in declaration order) for all simulation periods" << endl + << " params : Vector{Float64}(model_.param_nbr), parameter values in declaration order" << endl + << " steady_state : Vector{Float64}(model_endo_nbr)" << endl + << " it_ : Int, time period for exogenous variables for which to evaluate the model" << endl + << " residual : Vector{Float64}(model_.eq_nbr), residuals of the dynamic model equations in order of declaration of the equations." << endl + << " g1 : Matrix{Float64}(model_.eq_nbr, num_dynamic_vars), Jacobian matrix of the dynamic model equations" << endl + << " The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl + << " stored in model_.lead_lag_incidence" << endl + << " g2 : spzeros(model_.eq_nbr, (num_dynamic_vars)^2) Hessian matrix of the dynamic model equations" << endl + << " The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl + << " stored in model_.lead_lag_incidence" << endl + << " g3 : spzeros(model_.eq_nbr, (num_dynamic_vars)^3) Third order derivative matrix of the dynamic model equations;" << endl + << " The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl + << " stored in model_.lead_lag_incidence" << endl << endl + << "## Remarks ##" << endl + << " [1] `num_dynamic_vars` is the number of non zero entries in the lead lag incidence matrix, `model_.lead_lag_incidence.`" << endl + << " [2] The size of `T`, ie the value of `num_temp_terms`, depends on the version of the dynamic model called. The number of temporary variables" << endl + << " used for the different returned objects (residuals, jacobian, hessian or third order derivatives) is given by the elements in `tmp_nbr`" << endl + << " exported vector. The first element is the number of temporaries used for the computation of the residuals, the second element is the" << endl + << " number of temporaries used for the evaluation of the jacobian matrix, etc. If one calls the version of the dynamic model computing the" << endl + << " residuals, the jacobian and hessian matrices, then `T` must have at least `sum(tmp_nbr[1:3])` elements." << endl + << "=#" << endl << endl; + + // Write the number of temporary terms + output << "tmp_nbr = zeros(Int,4)" << endl + << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl + << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl + << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl + << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl; + + // dynamicResidTT! + output << "function dynamicResidTT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << tt_output[0].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamic! + output << "function dynamicResid!(T::Vector{Float64}, residual::AbstractVector{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl + << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl + << " @assert length(residual) == " << equations.size() << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T_flag" << endl + << " dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl + << " end" << endl + << d_output[0].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamicG1TT! + output << "function dynamicG1TT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl + << tt_output[1].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamicG1! + output << "function dynamicG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl + << " @assert length(T) >= " + << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl + << " @assert size(g1) == (" << equations.size() << ", " << getJacobianColsNbr() << ")" << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T_flag" << endl + << " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl + << " end" << endl + << " fill!(g1, 0.0)" << endl + << d_output[1].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamicG2TT! + output << "function dynamicG2TT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl + << tt_output[2].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamicG2! + int hessianColsNbr{getJacobianColsNbr() * getJacobianColsNbr()}; + output << "function dynamicG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl + << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl + << " @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T_flag" << endl + << " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl + << " end" << endl + << " fill!(g2, 0.0)" << endl + << d_output[2].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamicG3TT! + output << "function dynamicG3TT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl + << tt_output[3].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamicG3! + int ncols{hessianColsNbr * getJacobianColsNbr()}; + output << "function dynamicG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl + << " @assert length(T) >= " + << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl + << " @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl + << " @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T_flag" << endl + << " dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl + << " end" << endl + << " fill!(g3, 0.0)" << endl + << d_output[3].str() + << " return nothing" << endl + << "end" << endl << endl; + + // dynamic! + output << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicResid!(T, residual, y, x, params, steady_state, it_, true)" << endl + << " return nothing" << endl + << "end" << endl + << endl + << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicG1!(T, g1, y, x, params, steady_state, it_, true)" << endl + << " dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl + << " return nothing" << endl + << "end" << endl + << endl + << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicG2!(T, g2, y, x, params, steady_state, it_, true)" << endl + << " dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl + << " dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl + << " return nothing" << endl + << "end" << endl + << endl + << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Matrix{Float64}, " + << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl + << " dynamicG3!(T, g3, y, x, params, steady_state, it_, true)" << endl + << " dynamicG2!(T, g2, y, x, params, steady_state, it_, false)" << endl + << " dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl + << " dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl + << " return nothing" << endl + << "end" << endl + << endl; + + // Write function definition if BinaryOpcode::powerDeriv is used + writePowerDerivJulia(output); + + output << "end" << endl; + + writeToFileIfModified(output, basename + "Dynamic.jl"); } void @@ -1282,7 +1547,7 @@ DynamicModel::writeDynamicCFile(const string &basename) const { string filename = basename + "/model/src/dynamic.c"; - int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size(); + int ntt{static_cast(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size())}; ofstream output{filename, ios::out | ios::binary}; if (!output.is_open()) @@ -1307,7 +1572,31 @@ DynamicModel::writeDynamicCFile(const string &basename) const output << endl; - writeDynamicModel(output, true, false); + auto [d_output, tt_output] = writeModelFileHelper(); + + for (size_t i{0}; i < d_output.size(); i++) + { + string funcname{i == 0 ? "resid" : "g" + to_string(i)}; + output << "void dynamic_" << funcname << "_tt(const double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, int it_, double *restrict T)" << endl + << "{" << endl + << tt_output[i].str() + << "}" << endl + << endl + << "void dynamic_" << funcname << "(const double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, int it_, const double *restrict T, "; + if (i == 0) + output << "double *restrict residual"; + else if (i == 1) + output << "double *restrict g1"; + else + output << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v"; + output << ")" << endl + << "{" << endl; + if (i == 0) + output << " double lhs, rhs;" << endl; + output << d_output[i].str() + << "}" << endl + << endl; + } output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl << "{" << endl @@ -1335,7 +1624,7 @@ DynamicModel::writeDynamicCFile(const string &basename) const << endl << " if (nlhs >= 2)" << endl << " {" << endl - << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << dynJacobianColsNbr << ", mxREAL);" << endl + << " plhs[1] = mxCreateDoubleMatrix(" << equations.size() << ", " << getJacobianColsNbr() << ", mxREAL);" << endl << " double *g1 = mxGetPr(plhs[1]);" << endl << " dynamic_g1_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl << " dynamic_g1(y, x, nb_row_x, params, steady_state, it_, T, g1);" << endl @@ -1349,7 +1638,7 @@ DynamicModel::writeDynamicCFile(const string &basename) const << " dynamic_g2_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl << " dynamic_g2(y, x, nb_row_x, params, steady_state, it_, T, mxGetPr(g2_i), mxGetPr(g2_j), mxGetPr(g2_v));" << endl << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl - << " mxArray *n = mxCreateDoubleScalar(" << dynJacobianColsNbr*dynJacobianColsNbr << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g2_i, g2_j, g2_v, m, n };" << endl << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl << " plhs[2] = plhs_sparse[0];" << endl @@ -1368,7 +1657,7 @@ DynamicModel::writeDynamicCFile(const string &basename) const << " dynamic_g3_tt(y, x, nb_row_x, params, steady_state, it_, T);" << endl << " dynamic_g3(y, x, nb_row_x, params, steady_state, it_, T, mxGetPr(g3_i), mxGetPr(g3_j), mxGetPr(g3_v));" << endl << " mxArray *m = mxCreateDoubleScalar(" << equations.size() << ");" << endl - << " mxArray *n = mxCreateDoubleScalar(" << dynJacobianColsNbr*dynJacobianColsNbr*dynJacobianColsNbr << ");" << endl + << " mxArray *n = mxCreateDoubleScalar(" << getJacobianColsNbr()*getJacobianColsNbr()*getJacobianColsNbr() << ");" << endl << " mxArray *plhs_sparse[1], *prhs_sparse[5] = { g3_i, g3_j, g3_v, m, n };" << endl << R"( mexCallMATLAB(1, plhs_sparse, 5, prhs_sparse, "sparse");)" << endl << " plhs[3] = plhs_sparse[0];" << endl @@ -1757,447 +2046,6 @@ DynamicModel::writeDynamicMatlabCompatLayer(const string &basename) const output.close(); } -void -DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const -{ - writeDynamicModel("", DynamicOutput, use_dll, julia); -} - -void -DynamicModel::writeDynamicModel(const string &basename, bool use_dll, bool julia) const -{ - ofstream DynamicOutput; - writeDynamicModel(basename, DynamicOutput, use_dll, julia); -} - -void -DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput, bool use_dll, bool julia) const -{ - vector d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual) - vector tt_output(derivatives.size()); // Temp terms output (at all orders) - - ExprNodeOutputType output_type = (use_dll ? ExprNodeOutputType::CDynamicModel : - julia ? ExprNodeOutputType::juliaDynamicModel : ExprNodeOutputType::matlabDynamicModel); - - deriv_node_temp_terms_t tef_terms; - temporary_terms_t temp_term_union; - - writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs, - tt_output[0], output_type, tef_terms); - - writeTemporaryTerms(temporary_terms_derivatives[0], - temp_term_union, - temporary_terms_idxs, - tt_output[0], output_type, tef_terms); - - writeModelEquations(d_output[0], output_type, temp_term_union); - - int nrows = equations.size(); - int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr; - - // Writing Jacobian - if (!derivatives[1].empty()) - { - writeTemporaryTerms(temporary_terms_derivatives[1], - temp_term_union, - temporary_terms_idxs, - tt_output[1], output_type, tef_terms); - - for (const auto &first_derivative : derivatives[1]) - { - auto [eq, var] = vectorToTuple<2>(first_derivative.first); - expr_t d1 = first_derivative.second; - - jacobianHelper(d_output[1], eq, getDynJacobianCol(var), output_type); - d_output[1] << "="; - d1->writeOutput(d_output[1], output_type, - temp_term_union, temporary_terms_idxs, tef_terms); - d_output[1] << ";" << endl; - } - } - - // Write derivatives for order ≥ 2 - for (size_t i = 2; i < derivatives.size(); i++) - if (!derivatives[i].empty()) - { - writeTemporaryTerms(temporary_terms_derivatives[i], - temp_term_union, - temporary_terms_idxs, - tt_output[i], output_type, tef_terms); - - /* When creating the sparse matrix (in MATLAB or C mode), since storage - is in column-major order, output the first column, then the second, - then the third. This gives a significant performance boost in use_dll - mode (at both compilation and runtime), because it facilitates memory - accesses and expression reusage. */ - ostringstream i_output, j_output, v_output; - - for (int k{0}; // Current line index in the 3-column matrix - const auto &[vidx, d] : derivatives[i]) - { - int eq = vidx[0]; - - int col_idx = 0; - for (size_t j = 1; j < vidx.size(); j++) - { - col_idx *= dynJacobianColsNbr; - col_idx += getDynJacobianCol(vidx[j]); - } - - if (output_type == ExprNodeOutputType::juliaDynamicModel) - { - d_output[i] << " @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = "; - d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms); - d_output[i] << endl; - } - else - { - i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << eq + 1 << ";" << endl; - j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << col_idx + 1 << ";" << endl; - v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; - d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms); - v_output << ";" << endl; - - k++; - } - - // Output symetric elements at order 2 - if (i == 2 && vidx[1] != vidx[2]) - { - int col_idx_sym = getDynJacobianCol(vidx[2]) * dynJacobianColsNbr + getDynJacobianCol(vidx[1]); - - if (output_type == ExprNodeOutputType::juliaDynamicModel) - d_output[2] << " @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = " - << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl; - else - { - i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << eq + 1 << ";" << endl; - j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << col_idx_sym + 1 << ";" << endl; - v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" - << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl; - - k++; - } - } - } - if (output_type != ExprNodeOutputType::juliaDynamicModel) - d_output[i] << i_output.str() << j_output.str() << v_output.str(); - } - - if (output_type == ExprNodeOutputType::matlabDynamicModel) - { - // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201 - map tmp_paren_vars; - bool message_printed = false; - for (auto &it : tt_output) - fixNestedParenthesis(it, tmp_paren_vars, message_printed); - for (auto &it : d_output) - fixNestedParenthesis(it, tmp_paren_vars, message_printed); - - ostringstream init_output, end_output; - init_output << "residual = zeros(" << nrows << ", 1);"; - writeDynamicModelHelper(basename, "dynamic_resid", "residual", - "dynamic_resid_tt", - temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(), - "", init_output, end_output, - d_output[0], tt_output[0]); - - init_output.str(""); - init_output << "g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");"; - writeDynamicModelHelper(basename, "dynamic_g1", "g1", - "dynamic_g1_tt", - temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(), - "dynamic_resid_tt", - init_output, end_output, - d_output[1], tt_output[1]); - writeWrapperFunctions(basename, "g1"); - - // For order ≥ 2 - int ncols = dynJacobianColsNbr; - int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(); - for (size_t i = 2; i < derivatives.size(); i++) - { - ncols *= dynJacobianColsNbr; - ntt += temporary_terms_derivatives[i].size(); - string gname = "g" + to_string(i); - string gprevname = "g" + to_string(i-1); - - init_output.str(""); - end_output.str(""); - if (derivatives[i].size()) - { - init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl - << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl - << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl; - end_output << gname << " = sparse(" - << gname << "_i," << gname << "_j," << gname << "_v," - << nrows << "," << ncols << ");"; - } - else - init_output << gname << " = sparse([],[],[]," << nrows << "," << ncols << ");"; - writeDynamicModelHelper(basename, "dynamic_" + gname, gname, - "dynamic_" + gname + "_tt", - ntt, - "dynamic_" + gprevname + "_tt", - init_output, end_output, - d_output[i], tt_output[i]); - if (i <= 3) - writeWrapperFunctions(basename, gname); - } - - writeDynamicMatlabCompatLayer(basename); - } - else if (output_type == ExprNodeOutputType::CDynamicModel) - for (size_t i = 0; i < d_output.size(); i++) - { - string funcname = i == 0 ? "resid" : "g" + to_string(i); - DynamicOutput << "void dynamic_" << funcname << "_tt(const double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, int it_, double *restrict T)" << endl - << "{" << endl - << tt_output[i].str() - << "}" << endl - << endl - << "void dynamic_" << funcname << "(const double *restrict y, const double *restrict x, int nb_row_x, const double *restrict params, const double *restrict steady_state, int it_, const double *restrict T, "; - if (i == 0) - DynamicOutput << "double *restrict residual"; - else if (i == 1) - DynamicOutput << "double *restrict g1"; - else - DynamicOutput << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v"; - DynamicOutput << ")" << endl - << "{" << endl; - if (i == 0) - DynamicOutput << " double lhs, rhs;" << endl; - DynamicOutput << d_output[i].str() - << "}" << endl - << endl; - } - else - { - stringstream output; - - output << "module " << basename << "Dynamic" << endl - << "#" << endl - << "# NB: this file was automatically generated by Dynare" << endl - << "# from " << basename << ".mod" << endl - << "#" << endl - << "using StatsFuns" << endl << endl - << "export tmp_nbr, dynamic!, dynamicResid!, dynamicG1!, dynamicG2!, dynamicG3!" << endl << endl - << "#=" << endl - << "# The comments below apply to all functions contained in this module #" << endl - << " NB: The arguments contained on the first line of the function" << endl - << " definition are those that are modified in place" << endl << endl - << "## Exported Functions ##" << endl - << " dynamic! : Wrapper function; computes residuals, Jacobian, Hessian," << endl - << " and third derivatives depending on the arguments provided" << endl - << " dynamicResid! : Computes the dynamic model residuals" << endl - << " dynamicG1! : Computes the dynamic model Jacobian" << endl - << " dynamicG2! : Computes the dynamic model Hessian" << endl - << " dynamicG3! : Computes the dynamic model third derivatives" << endl << endl - << "## Exported Variables ##" << endl - << " tmp_nbr : Vector{Int}(4) respectively the number of temporary variables" << endl - << " for the residuals, g1, g2 and g3." << endl << endl - << "## Local Functions ##" << endl - << " dynamicResidTT! : Computes the dynamic model temporary terms for the residuals" << endl - << " dynamicG1TT! : Computes the dynamic model temporary terms for the Jacobian" << endl - << " dynamicG2TT! : Computes the dynamic model temporary terms for the Hessian" << endl - << " dynamicG3TT! : Computes the dynamic model temporary terms for the third derivatives" << endl << endl - << "## Function Arguments ##" << endl - << " T : Vector{Float64}(num_temp_terms), temporary terms" << endl - << " y : Vector{Float64}(num_dynamic_vars), endogenous variables in the order stored model_.lead_lag_incidence; see the manual" << endl - << " x : Matrix{Float64}(nperiods,model_.exo_nbr), exogenous variables (in declaration order) for all simulation periods" << endl - << " params : Vector{Float64}(model_.param_nbr), parameter values in declaration order" << endl - << " steady_state : Vector{Float64}(model_endo_nbr)" << endl - << " it_ : Int, time period for exogenous variables for which to evaluate the model" << endl - << " residual : Vector{Float64}(model_.eq_nbr), residuals of the dynamic model equations in order of declaration of the equations." << endl - << " g1 : Matrix{Float64}(model_.eq_nbr, num_dynamic_vars), Jacobian matrix of the dynamic model equations" << endl - << " The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl - << " stored in model_.lead_lag_incidence" << endl - << " g2 : spzeros(model_.eq_nbr, (num_dynamic_vars)^2) Hessian matrix of the dynamic model equations" << endl - << " The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl - << " stored in model_.lead_lag_incidence" << endl - << " g3 : spzeros(model_.eq_nbr, (num_dynamic_vars)^3) Third order derivative matrix of the dynamic model equations;" << endl - << " The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl - << " stored in model_.lead_lag_incidence" << endl << endl - << "## Remarks ##" << endl - << " [1] `num_dynamic_vars` is the number of non zero entries in the lead lag incidence matrix, `model_.lead_lag_incidence.`" << endl - << " [2] The size of `T`, ie the value of `num_temp_terms`, depends on the version of the dynamic model called. The number of temporary variables" << endl - << " used for the different returned objects (residuals, jacobian, hessian or third order derivatives) is given by the elements in `tmp_nbr`" << endl - << " exported vector. The first element is the number of temporaries used for the computation of the residuals, the second element is the" << endl - << " number of temporaries used for the evaluation of the jacobian matrix, etc. If one calls the version of the dynamic model computing the" << endl - << " residuals, the jacobian and hessian matrices, then `T` must have at least `sum(tmp_nbr[1:3])` elements." << endl - << "=#" << endl << endl; - - // Write the number of temporary terms - output << "tmp_nbr = zeros(Int,4)" << endl - << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl - << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl - << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl - << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl; - - // dynamicResidTT! - output << "function dynamicResidTT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << tt_output[0].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamic! - output << "function dynamicResid!(T::Vector{Float64}, residual::AbstractVector{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl - << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl - << " @assert length(residual) == " << nrows << endl - << " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T_flag" << endl - << " dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl - << " end" << endl - << d_output[0].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamicG1TT! - output << "function dynamicG1TT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl - << tt_output[1].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamicG1! - output << "function dynamicG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl - << " @assert length(T) >= " - << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl - << " @assert size(g1) == (" << nrows << ", " << dynJacobianColsNbr << ")" << endl - << " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T_flag" << endl - << " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl - << " end" << endl - << " fill!(g1, 0.0)" << endl - << d_output[1].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamicG2TT! - output << "function dynamicG2TT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl - << tt_output[2].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamicG2! - output << "function dynamicG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl - << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl - << " @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl - << " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T_flag" << endl - << " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl - << " end" << endl - << " fill!(g2, 0.0)" << endl - << d_output[2].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamicG3TT! - output << "function dynamicG3TT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl - << tt_output[3].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamicG3! - int ncols = hessianColsNbr * dynJacobianColsNbr; - output << "function dynamicG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl - << " @assert length(T) >= " - << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl - << " @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl - << " @assert length(y)+size(x, 2) == " << dynJacobianColsNbr << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T_flag" << endl - << " dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl - << " end" << endl - << " fill!(g3, 0.0)" << endl - << d_output[3].str() - << " return nothing" << endl - << "end" << endl << endl; - - // dynamic! - output << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicResid!(T, residual, y, x, params, steady_state, it_, true)" << endl - << " return nothing" << endl - << "end" << endl - << endl - << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicG1!(T, g1, y, x, params, steady_state, it_, true)" << endl - << " dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl - << " return nothing" << endl - << "end" << endl - << endl - << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicG2!(T, g2, y, x, params, steady_state, it_, true)" << endl - << " dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl - << " dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl - << " return nothing" << endl - << "end" << endl - << endl - << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Matrix{Float64}, " - << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl - << " dynamicG3!(T, g3, y, x, params, steady_state, it_, true)" << endl - << " dynamicG2!(T, g2, y, x, params, steady_state, it_, false)" << endl - << " dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl - << " dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl - << " return nothing" << endl - << "end" << endl - << endl; - - // Write function definition if BinaryOpcode::powerDeriv is used - writePowerDerivJulia(output); - - output << "end" << endl; - - writeToFileIfModified(output, basename + "Dynamic.jl"); - } -} - void DynamicModel::writeDynamicJacobianNonZeroElts(const string &basename) const { @@ -2808,7 +2656,7 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl try { int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag); - output << " " << getDynJacobianCol(varID) + 1; + output << " " << getJacobianCol(varID) + 1; if (lag == -1) { sstatic = 0; @@ -4068,8 +3916,8 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO indices, check that we will not overflow (see #89). Note that such a check is not needed for parameter derivatives, since tensors for those are not stored as matrices. This check cannot be done before since - dynJacobianColsNbr is not yet set.*/ - if (log2(dynJacobianColsNbr)*derivsOrder >= numeric_limits::digits) + getJacobianColsNbr() is not yet set.*/ + if (log2(getJacobianColsNbr())*derivsOrder >= numeric_limits::digits) { cerr << "ERROR: The dynamic derivatives matrix is too large. Please decrease the approximation order." << endl; exit(EXIT_FAILURE); @@ -4738,13 +4586,9 @@ DynamicModel::computeDerivIDs() { set> dynvars; - for (auto &equation : equations) - equation->collectDynamicVariables(SymbolType::endogenous, dynvars); - - dynJacobianColsNbr = dynvars.size(); - for (auto &equation : equations) { + equation->collectDynamicVariables(SymbolType::endogenous, dynvars); equation->collectDynamicVariables(SymbolType::exogenous, dynvars); equation->collectDynamicVariables(SymbolType::exogenousDet, dynvars); equation->collectDynamicVariables(SymbolType::parameter, dynvars); @@ -4838,61 +4682,30 @@ DynamicModel::addAllParamDerivId(set &deriv_id_set) void DynamicModel::computeDynJacobianCols(bool jacobianExo) { - /* Sort the dynamic endogenous variables by lexicographic order over (lag, type_specific_symbol_id) - and fill the dynamic columns for exogenous and exogenous deterministic */ + // Sort the dynamic endogenous variables by lexicographic order over (lag, type_specific_symbol_id) map, int> ordered_dyn_endo; + for (const auto &[symb_lag, deriv_id] : deriv_id_table) + if (const auto &[symb_id, lag] = symb_lag; + symbol_table.getType(symb_id) == SymbolType::endogenous) + ordered_dyn_endo[{ lag, symbol_table.getTypeSpecificID(symb_id) }] = deriv_id; - for (auto &[symb_lag, deriv_id] : deriv_id_table) - { - auto &[symb_id, lag] = symb_lag; - SymbolType type = symbol_table.getType(symb_id); - int tsid = symbol_table.getTypeSpecificID(symb_id); - - switch (type) - { - case SymbolType::endogenous: - ordered_dyn_endo[{ lag, tsid }] = deriv_id; - break; - case SymbolType::exogenous: - // At this point, dynJacobianColsNbr contains the number of dynamic endogenous - if (jacobianExo) - dyn_jacobian_cols_table[deriv_id] = dynJacobianColsNbr + tsid; - break; - case SymbolType::exogenousDet: - // At this point, dynJacobianColsNbr contains the number of dynamic endogenous - if (jacobianExo) - dyn_jacobian_cols_table[deriv_id] = dynJacobianColsNbr + symbol_table.exo_nbr() + tsid; - break; - case SymbolType::parameter: - case SymbolType::trend: - case SymbolType::logTrend: - // We don't assign a dynamic jacobian column to parameters or trend variables - break; - default: - // Shut up GCC - cerr << "DynamicModel::computeDynJacobianCols: impossible case" << endl; - exit(EXIT_FAILURE); - } - } - - // Fill in dynamic jacobian columns for endogenous + // Fill the dynamic jacobian columns for endogenous for (int sorted_id{0}; - auto &it : ordered_dyn_endo) - dyn_jacobian_cols_table[it.second] = sorted_id++; + const auto &[ignore, deriv_id] : ordered_dyn_endo) + dyn_jacobian_cols_table[deriv_id] = sorted_id++; - // Set final value for dynJacobianColsNbr + // Fill the dynamic columns for exogenous and exogenous deterministic if (jacobianExo) - dynJacobianColsNbr += symbol_table.exo_nbr() + symbol_table.exo_det_nbr(); -} - -int -DynamicModel::getDynJacobianCol(int deriv_id) const noexcept(false) -{ - if (auto it = dyn_jacobian_cols_table.find(deriv_id); - it == dyn_jacobian_cols_table.end()) - throw UnknownDerivIDException(); - else - return it->second; + for (const auto &[symb_lag, deriv_id] : deriv_id_table) + { + int symb_id{symb_lag.first}; + int tsid{symbol_table.getTypeSpecificID(symb_id)}; + if (SymbolType type{symbol_table.getType(symb_id)}; + type == SymbolType::exogenous) + dyn_jacobian_cols_table[deriv_id] = ordered_dyn_endo.size() + tsid; + else if (type == SymbolType::exogenousDet) + dyn_jacobian_cols_table[deriv_id] = ordered_dyn_endo.size() + symbol_table.exo_nbr() + tsid; + } } void @@ -4971,7 +4784,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con { auto [eq, var, param] = vectorToTuple<3>(indices); - int var_col = getDynJacobianCol(var) + 1; + int var_col = getJacobianCol(var) + 1; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; gp_output << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << var_col @@ -5023,7 +4836,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con { auto [eq, var, param1, param2] = vectorToTuple<4>(indices); - int var_col = getDynJacobianCol(var) + 1; + int var_col = getJacobianCol(var) + 1; int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1; int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1; @@ -5066,8 +4879,8 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con { auto [eq, var1, var2, param] = vectorToTuple<4>(indices); - int var1_col = getDynJacobianCol(var1) + 1; - int var2_col = getDynJacobianCol(var2) + 1; + int var1_col = getJacobianCol(var1) + 1; + int var2_col = getJacobianCol(var2) + 1; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; hp_output << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1" @@ -5109,9 +4922,9 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con { auto [eq, var1, var2, var3, param] = vectorToTuple<5>(indices); - int var1_col = getDynJacobianCol(var1) + 1; - int var2_col = getDynJacobianCol(var2) + 1; - int var3_col = getDynJacobianCol(var3) + 1; + int var1_col = getJacobianCol(var1) + 1; + int var2_col = getJacobianCol(var2) + 1; + int var3_col = getJacobianCol(var3) + 1; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; g3p_output << "g3p" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1" @@ -5209,7 +5022,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << rp_output.str() - << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl << gp_output.str() << "if nargout >= 3" << endl << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl @@ -5240,7 +5053,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << rp_output.str() - << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl << gp_output.str() << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl << rpp_output.str() @@ -6052,7 +5865,7 @@ DynamicModel::writeJsonDynamicModelInfo(ostream &output) const if (lag != -max_endo_lag) output << ","; int varID = getDerivID(symbol_table.getID(SymbolType::endogenous, endoID), lag); - output << " " << getDynJacobianCol(varID) + 1; + output << " " << getJacobianCol(varID) + 1; if (lag == -1) { sstatic = 0; @@ -6135,7 +5948,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c d_output[0] << ", "; writeJsonModelEquations(d_output[0], true); - int ncols = dynJacobianColsNbr; + int ncols{getJacobianColsNbr()}; for (size_t i = 1; i < derivatives.size(); i++) { string matrix_name = i == 1 ? "jacobian" : i == 2 ? "hessian" : i == 3 ? "third_derivative" : to_string(i) + "th_derivative"; @@ -6157,8 +5970,8 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c int col_idx = 0; for (size_t j = 1; j < vidx.size(); j++) { - col_idx *= dynJacobianColsNbr; - col_idx += getDynJacobianCol(vidx[j]); + col_idx *= getJacobianColsNbr(); + col_idx += getJacobianCol(vidx[j]); } if (writeDetails) @@ -6170,7 +5983,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c if (i == 2 && vidx[1] != vidx[2]) // Symmetric elements in hessian { - int col_idx_sym = getDynJacobianCol(vidx[2]) * dynJacobianColsNbr + getDynJacobianCol(vidx[1]); + int col_idx_sym = getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1]); d_output[i] << ", " << col_idx_sym + 1; } if (i > 1) @@ -6187,7 +6000,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c } d_output[i] << "]}"; - ncols *= dynJacobianColsNbr; + ncols *= getJacobianColsNbr(); } if (writeDetails) @@ -6254,7 +6067,7 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) gp_output << R"("deriv_jacobian_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvarcols": )" << dynJacobianColsNbr + << R"(, "nvarcols": )" << getJacobianColsNbr() << R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; for (bool printed_something{false}; @@ -6265,7 +6078,7 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) auto [eq, var, param] = vectorToTuple<3>(vidx); - int var_col = getDynJacobianCol(var) + 1; + int var_col = getJacobianCol(var) + 1; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; if (writeDetails) @@ -6322,7 +6135,7 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) gpp_output << R"("second_deriv_jacobian_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvarcols": )" << dynJacobianColsNbr + << R"(, "nvarcols": )" << getJacobianColsNbr() << R"(, "nparam1cols": )" << symbol_table.param_nbr() << R"(, "nparam2cols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; @@ -6334,7 +6147,7 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) auto [eq, var, param1, param2] = vectorToTuple<4>(vidx); - int var_col = getDynJacobianCol(var) + 1; + int var_col = getJacobianCol(var) + 1; int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1; int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1; @@ -6361,8 +6174,8 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) hp_output << R"("derivative_hessian_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvar1cols": )" << dynJacobianColsNbr - << R"(, "nvar2cols": )" << dynJacobianColsNbr + << R"(, "nvar1cols": )" << getJacobianColsNbr() + << R"(, "nvar2cols": )" << getJacobianColsNbr() << R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; for (bool printed_something{false}; @@ -6373,8 +6186,8 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) auto [eq, var1, var2, param] = vectorToTuple<4>(vidx); - int var1_col = getDynJacobianCol(var1) + 1; - int var2_col = getDynJacobianCol(var2) + 1; + int var1_col = getJacobianCol(var1) + 1; + int var2_col = getJacobianCol(var2) + 1; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; if (writeDetails) @@ -6401,9 +6214,9 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) g3p_output << R"("derivative_g3_wrt_params": {)" << R"( "neqs": )" << equations.size() - << R"(, "nvar1cols": )" << dynJacobianColsNbr - << R"(, "nvar2cols": )" << dynJacobianColsNbr - << R"(, "nvar3cols": )" << dynJacobianColsNbr + << R"(, "nvar1cols": )" << getJacobianColsNbr() + << R"(, "nvar2cols": )" << getJacobianColsNbr() + << R"(, "nvar3cols": )" << getJacobianColsNbr() << R"(, "nparamcols": )" << symbol_table.param_nbr() << R"(, "entries": [)"; for (bool printed_something{false}; @@ -6414,9 +6227,9 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails) auto [eq, var1, var2, var3, param] = vectorToTuple<5>(vidx); - int var1_col = getDynJacobianCol(var1) + 1; - int var2_col = getDynJacobianCol(var2) + 1; - int var3_col = getDynJacobianCol(var3) + 1; + int var1_col = getJacobianCol(var1) + 1; + int var2_col = getJacobianCol(var2) + 1; + int var3_col = getJacobianCol(var3) + 1; int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1; if (writeDetails) diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index a9a31c47..5cdaf360 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -96,10 +96,6 @@ private: //! Nonzero equations in the Hessian set nonzero_hessian_eqs; - //! Number of columns of dynamic jacobian - /*! Set by computeDerivID()s and computeDynJacobianCols() */ - int dynJacobianColsNbr{0}; - //! Creates mapping for variables and equations they are present in map> variableMapping; @@ -124,15 +120,10 @@ private: //! Writes dynamic model file (Matlab version) void writeDynamicMFile(const string &basename) const; //! Writes dynamic model file (Julia version) - void writeDynamicJuliaFile(const string &dynamic_basename) const; + void writeDynamicJuliaFile(const string &basename) const; //! Writes dynamic model file (C version) /*! \todo add third derivatives handling */ void writeDynamicCFile(const string &basename) const; - //! Writes the dynamic model equations and its derivatives - /*! \todo add third derivatives handling in C output */ - void writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const; - void writeDynamicModel(const string &basename, bool use_dll, bool julia) const; - void writeDynamicModel(const string &basename, ostream &DynamicOutput, bool use_dll, bool julia) const; //! Writes the main dynamic function of block decomposed model (MATLAB version) void writeDynamicBlockMFile(const string &basename) const; //! Writes the main dynamic function of block decomposed model (C version) @@ -463,7 +454,22 @@ public: void writeLatexOriginalFile(const string &basename, bool write_equation_tags) const; int getDerivID(int symb_id, int lag) const noexcept(false) override; - int getDynJacobianCol(int deriv_id) const noexcept(false) override; + + int + getJacobianCol(int deriv_id) const override + { + if (auto it = dyn_jacobian_cols_table.find(deriv_id); + it == dyn_jacobian_cols_table.end()) + throw UnknownDerivIDException(); + else + return it->second; + } + int + getJacobianColsNbr() const override + { + return dyn_jacobian_cols_table.size(); + } + void addAllParamDerivId(set &deriv_id_set) override; //! Returns true indicating that this is a dynamic model diff --git a/src/ExprNode.cc b/src/ExprNode.cc index ea461572..1bd72faf 100644 --- a/src/ExprNode.cc +++ b/src/ExprNode.cc @@ -1074,7 +1074,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, case ExprNodeOutputType::juliaDynamicModel: case ExprNodeOutputType::matlabDynamicModel: case ExprNodeOutputType::CDynamicModel: - i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type); + i = datatree.getJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type); output << "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type); break; case ExprNodeOutputType::CStaticModel: diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 8ec265a8..a28f9df0 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -1326,14 +1326,6 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, d output << "]"; } -void -ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) const -{ - temporary_terms_t tt; - temporary_terms_idxs_t ttidxs; - writeModelEquations(output, output_type, tt); -} - void ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const @@ -1599,19 +1591,6 @@ ModelTree::set_cutoff_to_zero() cutoff = 0; } -void -ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const -{ - if (isJuliaOutput(output_type)) - output << " @inbounds "; - output << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type); - if (isMatlabOutput(output_type) || isJuliaOutput(output_type)) - output << eq_nb + 1 << "," << col_nb + 1; - else - output << eq_nb + col_nb *equations.size(); - output << RIGHT_ARRAY_SUBSCRIPT(output_type); -} - void ModelTree::computeParamsDerivatives(int paramsDerivsOrder) { diff --git a/src/ModelTree.hh b/src/ModelTree.hh index de4968dd..c10592a3 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -28,6 +28,7 @@ #include #include #include +#include #include "DataTree.hh" #include "EquationTags.hh" @@ -254,9 +255,13 @@ protected: ostream &output, ExprNodeOutputType output_type, deriv_node_temp_terms_t &tef_terms) const; //! Writes model equations - void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const; void writeModelEquations(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const; + + // Returns outputs for derivatives and temporary terms at each derivation order + template + pair, vector> writeModelFileHelper() const; + //! Writes JSON model equations //! if residuals = true, we are writing the dynamic/static model. //! Otherwise, just the model equations (with line numbers, no tmp terms) @@ -448,9 +453,6 @@ public: void reorderAuxiliaryEquations(); //! Find equations of the form “variable=constant”, excluding equations with “mcp” tag (see dynare#1697) void findConstantEquationsWithoutMcpTag(map &subst_table) const; - //! Helper for writing the Jacobian elements in MATLAB and C - /*! Writes either (i+1,j+1) or [i+j*no_eq] */ - void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const; /* Given an expression, searches for the first equation that has exactly this expression on the LHS, and returns the RHS of that equation. If no such equation can be found, throws an ExprNode::MatchFailureExpression */ @@ -511,4 +513,139 @@ public: } }; +template +pair, vector> +ModelTree::writeModelFileHelper() const +{ + vector d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual) + vector tt_output(derivatives.size()); // Temp terms output (at all orders) + + deriv_node_temp_terms_t tef_terms; + temporary_terms_t temp_term_union; + + writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs, + tt_output[0], output_type, tef_terms); + + writeTemporaryTerms(temporary_terms_derivatives[0], + temp_term_union, + temporary_terms_idxs, + tt_output[0], output_type, tef_terms); + + writeModelEquations(d_output[0], output_type, temp_term_union); + + // Writing Jacobian + if (!derivatives[1].empty()) + { + writeTemporaryTerms(temporary_terms_derivatives[1], + temp_term_union, + temporary_terms_idxs, + tt_output[1], output_type, tef_terms); + + for (const auto &[indices, d1] : derivatives[1]) + { + auto [eq, var] = vectorToTuple<2>(indices); + + if constexpr(isJuliaOutput(output_type)) + d_output[1] << " @inbounds "; + d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type); + if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type)) + d_output[1] << eq + 1 << "," << getJacobianCol(var) + 1; + else + d_output[1] << eq + getJacobianCol(var)*equations.size(); + d_output[1] << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; + d1->writeOutput(d_output[1], output_type, + temp_term_union, temporary_terms_idxs, tef_terms); + d_output[1] << ";" << endl; + } + } + + // Write derivatives for order ≥ 2 + for (size_t i = 2; i < derivatives.size(); i++) + if (!derivatives[i].empty()) + { + writeTemporaryTerms(temporary_terms_derivatives[i], + temp_term_union, + temporary_terms_idxs, + tt_output[i], output_type, tef_terms); + + /* When creating the sparse matrix (in MATLAB or C mode), since storage + is in column-major order, output the first column, then the second, + then the third. This gives a significant performance boost in use_dll + mode (at both compilation and runtime), because it facilitates memory + accesses and expression reusage. */ + ostringstream i_output, j_output, v_output; + + for (int k{0}; // Current line index in the 3-column matrix + const auto &[vidx, d] : derivatives[i]) + { + int eq{vidx[0]}; + + int col_idx{0}; + for (size_t j = 1; j < vidx.size(); j++) + { + col_idx *= getJacobianColsNbr(); + col_idx += getJacobianCol(vidx[j]); + } + + if constexpr(isJuliaOutput(output_type)) + { + d_output[i] << " @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = "; + d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms); + d_output[i] << endl; + } + else + { + i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) + << "=" << eq + 1 << ";" << endl; + j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) + << "=" << col_idx + 1 << ";" << endl; + v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; + d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms); + v_output << ";" << endl; + + k++; + } + + // Output symetric elements at order 2 + if (i == 2 && vidx[1] != vidx[2]) + { + int col_idx_sym{getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1])}; + + if constexpr(isJuliaOutput(output_type)) + d_output[2] << " @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = " + << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl; + else + { + i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) + << "=" << eq + 1 << ";" << endl; + j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) + << "=" << col_idx_sym + 1 << ";" << endl; + v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" + << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) + << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type) + << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl; + + k++; + } + } + } + if constexpr(!isJuliaOutput(output_type)) + d_output[i] << i_output.str() << j_output.str() << v_output.str(); + } + + return { move(d_output), move(tt_output) }; +} + #endif diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 2418cd2d..5777967f 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -939,7 +939,73 @@ StaticModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_co void StaticModel::writeStaticMFile(const string &basename) const { - writeStaticModel(basename, false, false); + auto [d_output, tt_output] = writeModelFileHelper(); + + // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201 + map tmp_paren_vars; + bool message_printed{false}; + for (auto &it : tt_output) + fixNestedParenthesis(it, tmp_paren_vars, message_printed); + for (auto &it : d_output) + fixNestedParenthesis(it, tmp_paren_vars, message_printed); + + ostringstream init_output, end_output; + init_output << "residual = zeros(" << equations.size() << ", 1);"; + end_output << "if ~isreal(residual)" << endl + << " residual = real(residual)+imag(residual).^2;" << endl + << "end"; + writeStaticModelHelper(basename, "static_resid", "residual", "static_resid_tt", + temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(), + "", init_output, end_output, + d_output[0], tt_output[0]); + + init_output.str(""); + end_output.str(""); + init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");"; + end_output << "if ~isreal(g1)" << endl + << " g1 = real(g1)+2*imag(g1);" << endl + << "end"; + writeStaticModelHelper(basename, "static_g1", "g1", "static_g1_tt", + temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(), + "static_resid_tt", + init_output, end_output, + d_output[1], tt_output[1]); + writeWrapperFunctions(basename, "g1"); + + // For order ≥ 2 + int ncols{symbol_table.endo_nbr()}; + int ntt{static_cast(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size())}; + for (size_t i{2}; i < derivatives.size(); i++) + { + ncols *= symbol_table.endo_nbr(); + ntt += temporary_terms_derivatives[i].size(); + string gname{"g" + to_string(i)}; + string gprevname{"g" + to_string(i-1)}; + + init_output.str(""); + end_output.str(""); + if (derivatives[i].size()) + { + init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl + << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl + << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl; + end_output << gname << " = sparse(" + << gname << "_i," << gname << "_j," << gname << "_v," + << equations.size() << "," << ncols << ");"; + } + else + init_output << gname << " = sparse([],[],[]," << equations.size() << "," << ncols << ");"; + writeStaticModelHelper(basename, "static_" + gname, gname, + "static_" + gname + "_tt", + ntt, + "static_" + gprevname + "_tt", + init_output, end_output, + d_output[i], tt_output[i]); + if (i <= 3) + writeWrapperFunctions(basename, gname); + } + + writeStaticMatlabCompatLayer(basename); } void @@ -1098,464 +1164,13 @@ StaticModel::writeStaticMatlabCompatLayer(const string &basename) const output.close(); } -void -StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const -{ - writeStaticModel("", StaticOutput, use_dll, julia); -} - -void -StaticModel::writeStaticModel(const string &basename, bool use_dll, bool julia) const -{ - ofstream StaticOutput; - writeStaticModel(basename, StaticOutput, use_dll, julia); -} - -void -StaticModel::writeStaticModel(const string &basename, - ostream &StaticOutput, bool use_dll, bool julia) const -{ - vector d_output(derivatives.size()); // Derivatives output (at all orders, including 0=residual) - vector tt_output(derivatives.size()); // Temp terms output (at all orders) - - ExprNodeOutputType output_type = (use_dll ? ExprNodeOutputType::CStaticModel : - julia ? ExprNodeOutputType::juliaStaticModel : ExprNodeOutputType::matlabStaticModel); - - deriv_node_temp_terms_t tef_terms; - temporary_terms_t temp_term_union; - - writeModelLocalVariableTemporaryTerms(temp_term_union, temporary_terms_idxs, - tt_output[0], output_type, tef_terms); - - writeTemporaryTerms(temporary_terms_derivatives[0], - temp_term_union, - temporary_terms_idxs, - tt_output[0], output_type, tef_terms); - - writeModelEquations(d_output[0], output_type, temp_term_union); - - int nrows = equations.size(); - int JacobianColsNbr = symbol_table.endo_nbr(); - int hessianColsNbr = JacobianColsNbr*JacobianColsNbr; - - auto getJacobCol = [this](int var) { return symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)); }; - - // Write Jacobian w.r. to endogenous only - if (!derivatives[1].empty()) - { - writeTemporaryTerms(temporary_terms_derivatives[1], - temp_term_union, - temporary_terms_idxs, - tt_output[1], output_type, tef_terms); - - for (const auto & [indices, d1] : derivatives[1]) - { - auto [eq, var] = vectorToTuple<2>(indices); - - jacobianHelper(d_output[1], eq, getJacobCol(var), output_type); - d_output[1] << "="; - d1->writeOutput(d_output[1], output_type, - temp_term_union, temporary_terms_idxs, tef_terms); - d_output[1] << ";" << endl; - } - } - - // Write derivatives for order ≥ 2 - for (size_t i = 2; i < derivatives.size(); i++) - if (!derivatives[i].empty()) - { - writeTemporaryTerms(temporary_terms_derivatives[i], - temp_term_union, - temporary_terms_idxs, - tt_output[i], output_type, tef_terms); - - /* When creating the sparse matrix (in MATLAB or C mode), since storage - is in column-major order, output the first column, then the second, - then the third. This gives a significant performance boost in use_dll - mode (at both compilation and runtime), because it facilitates memory - accesses and expression reusage. */ - ostringstream i_output, j_output, v_output; - - for (int k{0}; // Current line index in the 3-column matrix - const auto &[vidx, d] : derivatives[i]) - { - int eq = vidx[0]; - - int col_idx = 0; - for (size_t j = 1; j < vidx.size(); j++) - { - col_idx *= JacobianColsNbr; - col_idx += getJacobCol(vidx[j]); - } - - if (output_type == ExprNodeOutputType::juliaStaticModel) - { - d_output[i] << " @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = "; - d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms); - d_output[i] << endl; - } - else - { - i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << eq + 1 << ";" << endl; - j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << col_idx + 1 << ";" << endl; - v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << "="; - d->writeOutput(v_output, output_type, temp_term_union, temporary_terms_idxs, tef_terms); - v_output << ";" << endl; - - k++; - } - - // Output symetric elements at order 2 - if (i == 2 && vidx[1] != vidx[2]) - { - int col_idx_sym = getJacobCol(vidx[2]) * JacobianColsNbr + getJacobCol(vidx[1]); - - if (output_type == ExprNodeOutputType::juliaStaticModel) - d_output[2] << " @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = " - << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl; - else - { - i_output << "g" << i << "_i" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << eq + 1 << ";" << endl; - j_output << "g" << i << "_j" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) - << "=" << col_idx_sym + 1 << ";" << endl; - v_output << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" - << "g" << i << "_v" << LEFT_ARRAY_SUBSCRIPT(output_type) - << k-1 + ARRAY_SUBSCRIPT_OFFSET(output_type) - << RIGHT_ARRAY_SUBSCRIPT(output_type) << ";" << endl; - - k++; - } - } - } - if (output_type != ExprNodeOutputType::juliaStaticModel) - d_output[i] << i_output.str() << j_output.str() << v_output.str(); - } - - if (output_type == ExprNodeOutputType::matlabStaticModel) - { - // Check that we don't have more than 32 nested parenthesis because Matlab does not suppor this. See Issue #1201 - map tmp_paren_vars; - bool message_printed = false; - for (auto &it : tt_output) - fixNestedParenthesis(it, tmp_paren_vars, message_printed); - for (auto &it : d_output) - fixNestedParenthesis(it, tmp_paren_vars, message_printed); - - ostringstream init_output, end_output; - init_output << "residual = zeros(" << equations.size() << ", 1);"; - end_output << "if ~isreal(residual)" << endl - << " residual = real(residual)+imag(residual).^2;" << endl - << "end"; - writeStaticModelHelper(basename, "static_resid", "residual", "static_resid_tt", - temporary_terms_mlv.size() + temporary_terms_derivatives[0].size(), - "", init_output, end_output, - d_output[0], tt_output[0]); - - init_output.str(""); - end_output.str(""); - init_output << "g1 = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ");"; - end_output << "if ~isreal(g1)" << endl - << " g1 = real(g1)+2*imag(g1);" << endl - << "end"; - writeStaticModelHelper(basename, "static_g1", "g1", "static_g1_tt", - temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(), - "static_resid_tt", - init_output, end_output, - d_output[1], tt_output[1]); - writeWrapperFunctions(basename, "g1"); - - // For order ≥ 2 - int ncols = JacobianColsNbr; - int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size(); - for (size_t i = 2; i < derivatives.size(); i++) - { - ncols *= JacobianColsNbr; - ntt += temporary_terms_derivatives[i].size(); - string gname = "g" + to_string(i); - string gprevname = "g" + to_string(i-1); - - init_output.str(""); - end_output.str(""); - if (derivatives[i].size()) - { - init_output << gname << "_i = zeros(" << NNZDerivatives[i] << ",1);" << endl - << gname << "_j = zeros(" << NNZDerivatives[i] << ",1);" << endl - << gname << "_v = zeros(" << NNZDerivatives[i] << ",1);" << endl; - end_output << gname << " = sparse(" - << gname << "_i," << gname << "_j," << gname << "_v," - << nrows << "," << ncols << ");"; - } - else - init_output << gname << " = sparse([],[],[]," << nrows << "," << ncols << ");"; - writeStaticModelHelper(basename, "static_" + gname, gname, - "static_" + gname + "_tt", - ntt, - "static_" + gprevname + "_tt", - init_output, end_output, - d_output[i], tt_output[i]); - if (i <= 3) - writeWrapperFunctions(basename, gname); - } - - writeStaticMatlabCompatLayer(basename); - } - else if (output_type == ExprNodeOutputType::CStaticModel) - { - for (size_t i = 0; i < d_output.size(); i++) - { - string funcname = i == 0 ? "resid" : "g" + to_string(i); - StaticOutput << "void static_" << funcname << "_tt(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl - << "{" << endl - << tt_output[i].str() - << "}" << endl - << endl - << "void static_" << funcname << "(const double *restrict y, const double *restrict x, const double *restrict params, const double *restrict T, "; - if (i == 0) - StaticOutput << "double *restrict residual"; - else if (i == 1) - StaticOutput << "double *restrict g1"; - else - StaticOutput << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v"; - StaticOutput << ")" << endl - << "{" << endl; - if (i == 0) - StaticOutput << " double lhs, rhs;" << endl; - StaticOutput << d_output[i].str() - << "}" << endl - << endl; - } - } - else - { - stringstream output; - output << "module " << basename << "Static" << endl - << "#" << endl - << "# NB: this file was automatically generated by Dynare" << endl - << "# from " << basename << ".mod" << endl - << "#" << endl - << "using StatsFuns" << endl << endl - << "export tmp_nbr, static!, staticResid!, staticG1!, staticG2!, staticG3!" << endl << endl - << "#=" << endl - << "# The comments below apply to all functions contained in this module #" << endl - << " NB: The arguments contained on the first line of the function" << endl - << " definition are those that are modified in place" << endl << endl - << "## Exported Functions ##" << endl - << " static! : Wrapper function; computes residuals, Jacobian, Hessian," << endl - << " and third order derivatives matroces depending on the arguments provided" << endl - << " staticResid! : Computes the static model residuals" << endl - << " staticG1! : Computes the static model Jacobian" << endl - << " staticG2! : Computes the static model Hessian" << endl - << " staticG3! : Computes the static model third derivatives" << endl << endl - << "## Exported Variables ##" << endl - << " tmp_nbr : Vector{Int}(4) respectively the number of temporary variables" << endl - << " for the residuals, g1, g2 and g3." << endl << endl - << "## Local Functions ##" << endl - << " staticResidTT! : Computes the static model temporary terms for the residuals" << endl - << " staticG1TT! : Computes the static model temporary terms for the Jacobian" << endl - << " staticG2TT! : Computes the static model temporary terms for the Hessian" << endl - << " staticG3TT! : Computes the static model temporary terms for the third derivatives" << endl << endl - << "## Function Arguments ##" << endl - << " T : Vector{Float64}(num_temp_terms) temporary terms" << endl - << " y : Vector{Float64}(model_.endo_nbr) endogenous variables in declaration order" << endl - << " x : Vector{Float64}(model_.exo_nbr) exogenous variables in declaration order" << endl - << " params : Vector{Float64}(model_.param) parameter values in declaration order" << endl - << " residual : Vector{Float64}(model_.eq_nbr) residuals of the static model equations" << endl - << " in order of declaration of the equations. Dynare may prepend auxiliary equations," << endl - << " see model.aux_vars" << endl - << " g1 : Matrix{Float64}(model.eq_nbr,model_.endo_nbr) Jacobian matrix of the static model equations" << endl - << " The columns and rows respectively correspond to the variables in declaration order and the" << endl - << " equations in order of declaration" << endl - << " g2 : spzeros(model.eq_nbr, model_.endo^2) Hessian matrix of the static model equations" << endl - << " The columns and rows respectively correspond to the variables in declaration order and the" << endl - << " equations in order of declaration" << endl - << " g3 : spzeros(model.eq_nbr, model_.endo^3) Third order derivatives matrix of the static model equations" << endl - << " The columns and rows respectively correspond to the variables in declaration order and the" << endl - << " equations in order of declaration" << endl << endl - << "## Remarks ##" << endl - << " [1] The size of `T`, ie the value of `num_temp_terms`, depends on the version of the static model called. The number of temporary variables" << endl - << " used for the different returned objects (residuals, jacobian, hessian or third order derivatives) is given by the elements in `tmp_nbr`" << endl - << " exported vector. The first element is the number of temporaries used for the computation of the residuals, the second element is the" << endl - << " number of temporaries used for the evaluation of the jacobian matrix, etc. If one calls the version of the static model computing the" << endl - << " residuals, and the jacobian and hessian matrices, then `T` must have at least `sum(tmp_nbr[1:3])` elements." << endl - << "=#" << endl << endl; - - // Write the number of temporary terms - output << "tmp_nbr = zeros(Int,4)" << endl - << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl - << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl - << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl - << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl; - - // staticResidTT! - output << "function staticResidTT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl - << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl - << tt_output[0].str() - << " return nothing" << endl - << "end" << endl << endl; - - // static! - output << "function staticResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl - << " @assert length(y) == " << symbol_table.endo_nbr() << endl - << " @assert length(x) == " << symbol_table.exo_nbr() << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " @assert length(residual) == " << equations.size() << endl - << " if T0_flag" << endl - << " staticResidTT!(T, y, x, params)" << endl - << " end" << endl - << d_output[0].str() - << " if ~isreal(residual)" << endl - << " residual = real(residual)+imag(residual).^2;" << endl - << " end" << endl - << " return nothing" << endl - << "end" << endl << endl; - - // staticG1TT! - output << "function staticG1TT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl - << " if T0_flag" << endl - << " staticResidTT!(T, y, x, params)" << endl - << " end" << endl - << tt_output[1].str() - << " return nothing" << endl - << "end" << endl << endl; - - // staticG1! - output << "function staticG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl - << " @assert length(T) >= " - << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl - << " @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << endl - << " @assert length(y) == " << symbol_table.endo_nbr() << endl - << " @assert length(x) == " << symbol_table.exo_nbr() << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T1_flag" << endl - << " staticG1TT!(T, y, x, params, T0_flag)" << endl - << " end" << endl - << " fill!(g1, 0.0)" << endl - << d_output[1].str() - << " if ~isreal(g1)" << endl - << " g1 = real(g1)+2*imag(g1);" << endl - << " end" << endl - << " return nothing" << endl - << "end" << endl << endl; - - // staticG2TT! - output << "function staticG2TT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl - << " if T1_flag" << endl - << " staticG1TT!(T, y, x, params, TO_flag)" << endl - << " end" << endl - << tt_output[2].str() - << " return nothing" << endl - << "end" << endl << endl; - - // staticG2! - output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl - << " @assert length(T) >= " - << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl - << " @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl - << " @assert length(y) == " << symbol_table.endo_nbr() << endl - << " @assert length(x) == " << symbol_table.exo_nbr() << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T2_flag" << endl - << " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl - << " end" << endl - << " fill!(g2, 0.0)" << endl - << d_output[2].str() - << " return nothing" << endl - << "end" << endl << endl; - - // staticG3TT! - output << "function staticG3TT!(T::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl - << " if T2_flag" << endl - << " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl - << " end" << endl - << tt_output[3].str() - << " return nothing" << endl - << "end" << endl << endl; - - // staticG3! - int ncols = hessianColsNbr * JacobianColsNbr; - output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl - << " @assert length(T) >= " - << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl - << " @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl - << " @assert length(y) == " << symbol_table.endo_nbr() << endl - << " @assert length(x) == " << symbol_table.exo_nbr() << endl - << " @assert length(params) == " << symbol_table.param_nbr() << endl - << " if T3_flag" << endl - << " staticG3TT!(T, y, x, params, T2_flag, T1_flag, T0_flag)" << endl - << " end" << endl - << " fill!(g3, 0.0)" << endl - << d_output[3].str() - << " return nothing" << endl - << "end" << endl << endl; - - // static! - output << "function static!(T::Vector{Float64}, residual::Vector{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl - << " staticResid!(T, residual, y, x, params, true)" << endl - << " return nothing" << endl - << "end" << endl - << endl - << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl - << " staticG1!(T, g1, y, x, params, true, true)" << endl - << " staticResid!(T, residual, y, x, params, false)" << endl - << " return nothing" << endl - << "end" << endl - << endl - << "function static!(T::Vector{Float64}, g1::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl - << " staticG1!(T, g1, y, x, params, true, false)" << endl - << " return nothing" << endl - << "end" << endl - << endl - << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl - << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl - << " staticG2!(T, g2, y, x, params, true, true, true)" << endl - << " staticG1!(T, g1, y, x, params, false, false)" << endl - << " staticResid!(T, residual, y, x, params, false)" << endl - << " return nothing" << endl - << "end" << endl - << endl; - - // Write function definition if BinaryOpcode::powerDeriv is used - writePowerDerivJulia(output); - - output << "end" << endl; - - writeToFileIfModified(output, basename + "Static.jl"); - } -} - void StaticModel::writeStaticCFile(const string &basename) const { // Writing comments and function definition command - string filename = basename + "/model/src/static.c"; + string filename{basename + "/model/src/static.c"}; - int ntt = temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size(); + int ntt{static_cast(temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size())}; ofstream output{filename, ios::out | ios::binary}; if (!output.is_open()) @@ -1581,7 +1196,31 @@ StaticModel::writeStaticCFile(const string &basename) const output << endl; - writeStaticModel(output, true, false); + auto [d_output, tt_output] = writeModelFileHelper(); + + for (size_t i = 0; i < d_output.size(); i++) + { + string funcname{i == 0 ? "resid" : "g" + to_string(i)}; + output << "void static_" << funcname << "_tt(const double *restrict y, const double *restrict x, const double *restrict params, double *restrict T)" << endl + << "{" << endl + << tt_output[i].str() + << "}" << endl + << endl + << "void static_" << funcname << "(const double *restrict y, const double *restrict x, const double *restrict params, const double *restrict T, "; + if (i == 0) + output << "double *restrict residual"; + else if (i == 1) + output << "double *restrict g1"; + else + output << "double *restrict " << funcname << "_i, double *restrict " << funcname << "_j, double *restrict " << funcname << "_v"; + output << ")" << endl + << "{" << endl; + if (i == 0) + output << " double lhs, rhs;" << endl; + output << d_output[i].str() + << "}" << endl + << endl; + } output << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl << "{" << endl @@ -1639,7 +1278,213 @@ StaticModel::writeStaticCFile(const string &basename) const void StaticModel::writeStaticJuliaFile(const string &basename) const { - writeStaticModel(basename, false, true); + auto [d_output, tt_output] = writeModelFileHelper(); + + stringstream output; + output << "module " << basename << "Static" << endl + << "#" << endl + << "# NB: this file was automatically generated by Dynare" << endl + << "# from " << basename << ".mod" << endl + << "#" << endl + << "using StatsFuns" << endl << endl + << "export tmp_nbr, static!, staticResid!, staticG1!, staticG2!, staticG3!" << endl << endl + << "#=" << endl + << "# The comments below apply to all functions contained in this module #" << endl + << " NB: The arguments contained on the first line of the function" << endl + << " definition are those that are modified in place" << endl << endl + << "## Exported Functions ##" << endl + << " static! : Wrapper function; computes residuals, Jacobian, Hessian," << endl + << " and third order derivatives matroces depending on the arguments provided" << endl + << " staticResid! : Computes the static model residuals" << endl + << " staticG1! : Computes the static model Jacobian" << endl + << " staticG2! : Computes the static model Hessian" << endl + << " staticG3! : Computes the static model third derivatives" << endl << endl + << "## Exported Variables ##" << endl + << " tmp_nbr : Vector{Int}(4) respectively the number of temporary variables" << endl + << " for the residuals, g1, g2 and g3." << endl << endl + << "## Local Functions ##" << endl + << " staticResidTT! : Computes the static model temporary terms for the residuals" << endl + << " staticG1TT! : Computes the static model temporary terms for the Jacobian" << endl + << " staticG2TT! : Computes the static model temporary terms for the Hessian" << endl + << " staticG3TT! : Computes the static model temporary terms for the third derivatives" << endl << endl + << "## Function Arguments ##" << endl + << " T : Vector{Float64}(num_temp_terms) temporary terms" << endl + << " y : Vector{Float64}(model_.endo_nbr) endogenous variables in declaration order" << endl + << " x : Vector{Float64}(model_.exo_nbr) exogenous variables in declaration order" << endl + << " params : Vector{Float64}(model_.param) parameter values in declaration order" << endl + << " residual : Vector{Float64}(model_.eq_nbr) residuals of the static model equations" << endl + << " in order of declaration of the equations. Dynare may prepend auxiliary equations," << endl + << " see model.aux_vars" << endl + << " g1 : Matrix{Float64}(model.eq_nbr,model_.endo_nbr) Jacobian matrix of the static model equations" << endl + << " The columns and rows respectively correspond to the variables in declaration order and the" << endl + << " equations in order of declaration" << endl + << " g2 : spzeros(model.eq_nbr, model_.endo^2) Hessian matrix of the static model equations" << endl + << " The columns and rows respectively correspond to the variables in declaration order and the" << endl + << " equations in order of declaration" << endl + << " g3 : spzeros(model.eq_nbr, model_.endo^3) Third order derivatives matrix of the static model equations" << endl + << " The columns and rows respectively correspond to the variables in declaration order and the" << endl + << " equations in order of declaration" << endl << endl + << "## Remarks ##" << endl + << " [1] The size of `T`, ie the value of `num_temp_terms`, depends on the version of the static model called. The number of temporary variables" << endl + << " used for the different returned objects (residuals, jacobian, hessian or third order derivatives) is given by the elements in `tmp_nbr`" << endl + << " exported vector. The first element is the number of temporaries used for the computation of the residuals, the second element is the" << endl + << " number of temporaries used for the evaluation of the jacobian matrix, etc. If one calls the version of the static model computing the" << endl + << " residuals, and the jacobian and hessian matrices, then `T` must have at least `sum(tmp_nbr[1:3])` elements." << endl + << "=#" << endl << endl; + + // Write the number of temporary terms + output << "tmp_nbr = zeros(Int,4)" << endl + << "tmp_nbr[1] = " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << "# Number of temporary terms for the residuals" << endl + << "tmp_nbr[2] = " << temporary_terms_derivatives[1].size() << "# Number of temporary terms for g1 (jacobian)" << endl + << "tmp_nbr[3] = " << temporary_terms_derivatives[2].size() << "# Number of temporary terms for g2 (hessian)" << endl + << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl; + + // staticResidTT! + output << "function staticResidTT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl + << " @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl + << tt_output[0].str() + << " return nothing" << endl + << "end" << endl << endl; + + // static! + output << "function staticResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl + << " @assert length(y) == " << symbol_table.endo_nbr() << endl + << " @assert length(x) == " << symbol_table.exo_nbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " @assert length(residual) == " << equations.size() << endl + << " if T0_flag" << endl + << " staticResidTT!(T, y, x, params)" << endl + << " end" << endl + << d_output[0].str() + << " if ~isreal(residual)" << endl + << " residual = real(residual)+imag(residual).^2;" << endl + << " end" << endl + << " return nothing" << endl + << "end" << endl << endl; + + // staticG1TT! + output << "function staticG1TT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl + << " if T0_flag" << endl + << " staticResidTT!(T, y, x, params)" << endl + << " end" << endl + << tt_output[1].str() + << " return nothing" << endl + << "end" << endl << endl; + + // staticG1! + output << "function staticG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl + << " @assert length(T) >= " + << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl + << " @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << endl + << " @assert length(y) == " << symbol_table.endo_nbr() << endl + << " @assert length(x) == " << symbol_table.exo_nbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T1_flag" << endl + << " staticG1TT!(T, y, x, params, T0_flag)" << endl + << " end" << endl + << " fill!(g1, 0.0)" << endl + << d_output[1].str() + << " if ~isreal(g1)" << endl + << " g1 = real(g1)+2*imag(g1);" << endl + << " end" << endl + << " return nothing" << endl + << "end" << endl << endl; + + // staticG2TT! + output << "function staticG2TT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl + << " if T1_flag" << endl + << " staticG1TT!(T, y, x, params, TO_flag)" << endl + << " end" << endl + << tt_output[2].str() + << " return nothing" << endl + << "end" << endl << endl; + + // staticG2! + int hessianColsNbr{symbol_table.endo_nbr() * symbol_table.endo_nbr()}; + output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl + << " @assert length(T) >= " + << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl + << " @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl + << " @assert length(y) == " << symbol_table.endo_nbr() << endl + << " @assert length(x) == " << symbol_table.exo_nbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T2_flag" << endl + << " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl + << " end" << endl + << " fill!(g2, 0.0)" << endl + << d_output[2].str() + << " return nothing" << endl + << "end" << endl << endl; + + // staticG3TT! + output << "function staticG3TT!(T::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl + << " if T2_flag" << endl + << " staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl + << " end" << endl + << tt_output[3].str() + << " return nothing" << endl + << "end" << endl << endl; + + // staticG3! + int ncols{hessianColsNbr * symbol_table.endo_nbr()}; + output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl + << " @assert length(T) >= " + << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl + << " @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl + << " @assert length(y) == " << symbol_table.endo_nbr() << endl + << " @assert length(x) == " << symbol_table.exo_nbr() << endl + << " @assert length(params) == " << symbol_table.param_nbr() << endl + << " if T3_flag" << endl + << " staticG3TT!(T, y, x, params, T2_flag, T1_flag, T0_flag)" << endl + << " end" << endl + << " fill!(g3, 0.0)" << endl + << d_output[3].str() + << " return nothing" << endl + << "end" << endl << endl; + + // static! + output << "function static!(T::Vector{Float64}, residual::Vector{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl + << " staticResid!(T, residual, y, x, params, true)" << endl + << " return nothing" << endl + << "end" << endl + << endl + << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl + << " staticG1!(T, g1, y, x, params, true, true)" << endl + << " staticResid!(T, residual, y, x, params, false)" << endl + << " return nothing" << endl + << "end" << endl + << endl + << "function static!(T::Vector{Float64}, g1::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl + << " staticG1!(T, g1, y, x, params, true, false)" << endl + << " return nothing" << endl + << "end" << endl + << endl + << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl + << " y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl + << " staticG2!(T, g2, y, x, params, true, true, true)" << endl + << " staticG1!(T, g1, y, x, params, false, false)" << endl + << " staticResid!(T, residual, y, x, params, false)" << endl + << " return nothing" << endl + << "end" << endl + << endl; + + // Write function definition if BinaryOpcode::powerDeriv is used + writePowerDerivJulia(output); + + output << "end" << endl; + + writeToFileIfModified(output, basename + "Static.jl"); } void diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 33c64dbe..28ecc8cc 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -43,9 +43,6 @@ private: //! Writes static model file (Julia version) void writeStaticJuliaFile(const string &basename) const; - //! Writes the static model equations and its derivatives - void writeStaticModel(const string &basename, ostream &StaticOutput, bool use_dll, bool julia) const; - //! Writes the main static function of block decomposed model (MATLAB version) void writeStaticBlockMFile(const string &basename) const; @@ -89,8 +86,18 @@ private: int getLagByDerivID(int deriv_id) const noexcept(false) override; //! Get the symbol ID corresponding to a derivation ID int getSymbIDByDerivID(int deriv_id) const noexcept(false) override; - //! Compute the column indices of the static Jacobian - void computeStatJacobianCols(); + + int + getJacobianCol(int deriv_id) const override + { + return symbol_table.getTypeSpecificID(getSymbIDByDerivID(deriv_id)); + } + int + getJacobianColsNbr() const override + { + return symbol_table.endo_nbr(); + } + //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables void computeChainRuleJacobian(); @@ -106,9 +113,6 @@ private: //! Create a legacy *_static.m file for Matlab/Octave not yet using the temporary terms array interface void writeStaticMatlabCompatLayer(const string &name) const; - void writeStaticModel(ostream &DynamicOutput, bool use_dll, bool julia) const; - void writeStaticModel(const string &dynamic_basename, bool use_dll, bool julia) const; - public: StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants,