From a58109d09470f203f2e2f7b85509cb71831616db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 19 Jul 2022 18:24:36 +0200 Subject: [PATCH] Bytecode: refactor methods for writing .cod and .bin files in the block decomposition case --- src/Bytecode.hh | 4 +- src/DynamicModel.cc | 434 ++++++++------------------------------------ src/DynamicModel.hh | 18 +- src/ModelTree.cc | 51 ++++++ src/ModelTree.hh | 238 ++++++++++++++++++++++++ src/StaticModel.cc | 407 ++++------------------------------------- src/StaticModel.hh | 15 +- 7 files changed, 413 insertions(+), 754 deletions(-) diff --git a/src/Bytecode.hh b/src/Bytecode.hh index c7bef351..9d867e55 100644 --- a/src/Bytecode.hh +++ b/src/Bytecode.hh @@ -928,7 +928,7 @@ public: } FBEGINBLOCK_(int size_arg, BlockSimulationType type_arg, int first_element, int block_size, const vector &variable_arg, const vector &equation_arg, - bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg, int nb_col_jacob_arg, + bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg, int nb_col_jacob_arg, int det_exo_size_arg, int nb_col_det_exo_jacob_arg, int exo_size_arg, int nb_col_exo_jacob_arg, int other_endo_size_arg, int nb_col_other_endo_jacob_arg, vector det_exogenous_arg, vector exogenous_arg, vector other_endogenous_arg) : BytecodeInstruction{Tags::FBEGINBLOCK}, @@ -955,7 +955,7 @@ public: } FBEGINBLOCK_(int size_arg, BlockSimulationType type_arg, int first_element, int block_size, const vector &variable_arg, const vector &equation_arg, - bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg, int nb_col_jacob_arg) : + bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int u_count_int_arg, int nb_col_jacob_arg) : BytecodeInstruction{Tags::FBEGINBLOCK}, size{size_arg}, type{type_arg}, diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 8174bc13..a7b140dc 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -165,26 +165,6 @@ DynamicModel::operator=(const DynamicModel &m) return *this; } -void -DynamicModel::writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const -{ - if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), lag) }); - it != derivatives[1].end()) - it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms, temporary_terms_idxs, tef_terms); - else - code_file << FLDZ_{}; -} - -void -DynamicModel::writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const -{ - if (auto it = blocks_derivatives[blk].find({ eq, var, lag }); - it != blocks_derivatives[blk].end()) - it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms, temporary_terms_idxs, tef_terms); - else - code_file << FLDZ_{}; -} - void DynamicModel::additionalBlockTemporaryTerms(int blk, vector> &blocks_temporary_terms, @@ -332,6 +312,39 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const } } +void +DynamicModel::writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block, + const temporary_terms_t &temporary_terms_union, + const deriv_node_temp_terms_t &tef_terms) const +{ + constexpr ExprNodeBytecodeOutputType output_type {ExprNodeBytecodeOutputType::dynamicModel}; + + /* FIXME: there is an inconsistency between endos and the following 3 other + variable types. For the latter, the index of equation within the block is + taken from FNUMEXPR, while it is taken from FSTPG3 for the former. */ + for (const auto &[indices, d] : blocks_derivatives_exo[block]) + { + const auto &[eq, var, lag] {indices}; + code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eq, 0, lag}; + d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag })}; + } + for (const auto &[indices, d] : blocks_derivatives_exo_det[block]) + { + const auto &[eq, var, lag] {indices}; + code_file << FNUMEXPR_{ExpressionType::FirstExodetDerivative, eq, 0, lag}; + d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag })}; + } + for (const auto &[indices, d] : blocks_derivatives_other_endo[block]) + { + const auto &[eq, var, lag] {indices}; + code_file << FNUMEXPR_{ExpressionType::FirstOtherEndoDerivative, eq, 0, lag}; + d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag })}; + } +} + void DynamicModel::writeDynamicPerBlockCFiles(const string &basename) const { @@ -591,301 +604,56 @@ DynamicModel::writeDynamicBytecode(const string &basename) const void DynamicModel::writeDynamicBlockBytecode(const string &basename) const { - struct Uff_l - { - int u, var, lag; - Uff_l *pNext; - }; + BytecodeWriter code_file {basename + "/model/bytecode/dynamic.cod"}; - struct Uff - { - Uff_l *Ufl, *Ufl_First; - }; - - int i, v; - string tmp_s; - ostringstream tmp_output; - expr_t lhs = nullptr, rhs = nullptr; - BinaryOpNode *eq_node; - Uff Uf[symbol_table.endo_nbr()]; - map reference_count; - vector feedback_variables; - bool file_open = false; - BytecodeWriter code_file{basename + "/model/bytecode/dynamic.cod"}; - - //Temporary variables declaration + const string bin_filename {basename + "/model/bytecode/dynamic.bin"}; + ofstream bin_file {bin_filename, ios::out | ios::binary}; + if (!bin_file.is_open()) + { + cerr << R"(Error : Can't open file ")" << bin_filename << R"(" for writing)" << endl; + exit(EXIT_FAILURE); + } + // Temporary variables declaration code_file << FDIMT_{static_cast(blocks_temporary_terms_idxs.size())}; - for (int block = 0; block < static_cast(blocks.size()); block++) + for (int block {0}; block < static_cast(blocks.size()); block++) { - feedback_variables.clear(); - if (block > 0) - code_file << FENDBLOCK_{}; - int count_u; - int u_count_int = 0; - BlockSimulationType simulation_type = blocks[block].simulation_type; - int block_size = blocks[block].size; - int block_mfs = blocks[block].mfs_size; - int block_recursive = blocks[block].getRecursiveSize(); - int block_max_lag = blocks[block].max_lag; - int block_max_lead = blocks[block].max_lead; + const BlockSimulationType simulation_type {blocks[block].simulation_type}; - if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple - || simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveBackwardComplete - || simulation_type == BlockSimulationType::solveForwardComplete) - { - writeBlockBytecodeBinFile(basename, block, u_count_int, file_open, - simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple); - file_open = true; - } + // Write section of .bin file except for evaluate blocks and solve simple blocks + const int u_count {simulation_type == BlockSimulationType::solveTwoBoundariesSimple + || simulation_type == BlockSimulationType::solveTwoBoundariesComplete + || simulation_type == BlockSimulationType::solveBackwardComplete + || simulation_type == BlockSimulationType::solveForwardComplete + ? writeBlockBytecodeBinFile(bin_file, block) + : 0}; - code_file << FBEGINBLOCK_{block_mfs, - simulation_type, - blocks[block].first_equation, - block_size, - endo_idx_block2orig, - eq_idx_block2orig, - blocks[block].linear, - symbol_table.endo_nbr(), - block_max_lag, - block_max_lead, - u_count_int, - static_cast(blocks_jacob_cols_endo[block].size()), - static_cast(blocks_exo_det[block].size()), - static_cast(blocks_jacob_cols_exo_det[block].size()), - static_cast(blocks_exo[block].size()), - static_cast(blocks_jacob_cols_exo[block].size()), - static_cast(blocks_other_endo[block].size()), - static_cast(blocks_jacob_cols_other_endo[block].size()), - vector(blocks_exo_det[block].begin(), blocks_exo_det[block].end()), - vector(blocks_exo[block].begin(), blocks_exo[block].end()), - vector(blocks_other_endo[block].begin(), blocks_other_endo[block].end())}; + code_file << FBEGINBLOCK_{blocks[block].mfs_size, + simulation_type, + blocks[block].first_equation, + blocks[block].size, + endo_idx_block2orig, + eq_idx_block2orig, + blocks[block].linear, + symbol_table.endo_nbr(), + blocks[block].max_lag, + blocks[block].max_lead, + u_count, + static_cast(blocks_jacob_cols_endo[block].size()), + static_cast(blocks_exo_det[block].size()), + static_cast(blocks_jacob_cols_exo_det[block].size()), + static_cast(blocks_exo[block].size()), + static_cast(blocks_jacob_cols_exo[block].size()), + static_cast(blocks_other_endo[block].size()), + static_cast(blocks_jacob_cols_other_endo[block].size()), + { blocks_exo_det[block].begin(), blocks_exo_det[block].end() }, + { blocks_exo[block].begin(), blocks_exo[block].end() }, + { blocks_other_endo[block].begin(), blocks_other_endo[block].end() }}; - temporary_terms_t temporary_terms_union; - - //The Temporary terms - deriv_node_temp_terms_t tef_terms; - - auto write_eq_tt = [&](int eq) - { - for (auto it : blocks_temporary_terms[block][eq]) - { - if (dynamic_cast(it)) - it->writeBytecodeExternalFunctionOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - - code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)}; - it->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPT_{blocks_temporary_terms_idxs.at(it)}; - temporary_terms_union.insert(it); - } - }; - - // The equations - for (i = 0; i < block_size; i++) - { - write_eq_tt(i); - - int variable_ID, equation_ID; - EquationType equ_type; - - switch (simulation_type) - { - evaluation: - case BlockSimulationType::evaluateBackward: - case BlockSimulationType::evaluateForward: - equ_type = getBlockEquationType(block, i); - code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; - if (equ_type == EquationType::evaluate) - { - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - } - else if (equ_type == EquationType::evaluateRenormalized) - { - eq_node = getBlockEquationRenormalizedExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - } - break; - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - case BlockSimulationType::solveTwoBoundariesComplete: - case BlockSimulationType::solveTwoBoundariesSimple: - if (i < block_recursive) - goto evaluation; - variable_ID = getBlockVariableID(block, i); - equation_ID = getBlockEquationID(block, i); - feedback_variables.push_back(variable_ID); - Uf[equation_ID].Ufl = nullptr; - goto end; - default: - end: - code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - - code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive}; - } - } - code_file << FENDEQU_{}; - - /* If the block is not of type “evaluate backward/forward”, then we write - the temporary terms for derivatives at this point, i.e. before the - JMPIFEVAL, because they will be needed in both “simulate” and - “evaluate” modes. */ - if (simulation_type != BlockSimulationType::evaluateBackward - && simulation_type != BlockSimulationType::evaluateForward) - write_eq_tt(blocks[block].size); - - // Get the current code_file position and jump if evaluating - int pos_jmpifeval = code_file.getInstructionCounter(); - code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being - - /* Write the derivatives for the “simulate” mode (not needed if the block - is of type “evaluate backward/forward”) */ - if (simulation_type != BlockSimulationType::evaluateBackward - && simulation_type != BlockSimulationType::evaluateForward) - { - switch (simulation_type) - { - case BlockSimulationType::solveBackwardSimple: - case BlockSimulationType::solveForwardSimple: - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0}; - writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG_{0}; - break; - - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - case BlockSimulationType::solveTwoBoundariesComplete: - case BlockSimulationType::solveTwoBoundariesSimple: - count_u = feedback_variables.size(); - for (const auto &[indices, ignore] : blocks_derivatives[block]) - { - auto [eq, var, lag] = indices; - int eqr = getBlockEquationID(block, eq); - int varr = getBlockVariableID(block, var); - if (eq >= block_recursive and var >= block_recursive) - { - if (lag != 0 - && (simulation_type == BlockSimulationType::solveForwardComplete - || simulation_type == BlockSimulationType::solveBackwardComplete)) - continue; - if (!Uf[eqr].Ufl) - { - Uf[eqr].Ufl = static_cast(malloc(sizeof(Uff_l))); - Uf[eqr].Ufl_First = Uf[eqr].Ufl; - } - else - { - Uf[eqr].Ufl->pNext = static_cast(malloc(sizeof(Uff_l))); - Uf[eqr].Ufl = Uf[eqr].Ufl->pNext; - } - Uf[eqr].Ufl->pNext = nullptr; - Uf[eqr].Ufl->u = count_u; - Uf[eqr].Ufl->var = varr; - Uf[eqr].Ufl->lag = lag; - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag}; - writeBytecodeChainRuleDerivative(code_file, block, eq, var, lag, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPU_{count_u}; - count_u++; - } - } - for (i = 0; i < block_size; i++) - { - if (i >= block_recursive) - { - code_file << FLDR_{i-block_recursive} << FLDZ_{}; - - v = getBlockEquationID(block, i); - for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext) - code_file << FLDU_{Uf[v].Ufl->u} - << FLDV_{SymbolType::endogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag} - << FBINARY_{BinaryOpcode::times} - << FCUML_{}; - Uf[v].Ufl = Uf[v].Ufl_First; - while (Uf[v].Ufl) - { - Uf[v].Ufl_First = Uf[v].Ufl->pNext; - free(Uf[v].Ufl); - Uf[v].Ufl = Uf[v].Ufl_First; - } - code_file << FBINARY_{BinaryOpcode::minus} - << FSTPU_{i - block_recursive}; - } - } - break; - default: - break; - } - } - - // Jump unconditionally after the block - int pos_jmp = code_file.getInstructionCounter(); - code_file << FJMP_{0}; // Use 0 as jump offset for the time being - // Update jump offset for previous JMPIFEVAL - code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval}); - - /* If the block is of type “evaluate backward/forward”, then write the - temporary terms for derivatives at this point, because they have not - been written before the JMPIFEVAL. */ - if (simulation_type == BlockSimulationType::evaluateBackward - || simulation_type == BlockSimulationType::evaluateForward) - write_eq_tt(blocks[block].size); - - // Write the derivatives for the “evaluate” mode - for (const auto &[indices, d] : blocks_derivatives[block]) - { - auto [eq, var, lag] = indices; - int eqr = getBlockEquationID(block, eq); - int varr = getBlockVariableID(block, var); - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag}; - writeBytecodeDerivative(code_file, eqr, varr, lag, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_endo[block].at({ var, lag })}; - } - for (const auto &[indices, d] : blocks_derivatives_exo[block]) - { - auto [eqr, var, lag] = indices; - int eq = getBlockEquationID(block, eqr); - int varr = 0; // Dummy value, actually unused by the bytecode MEX - code_file << FNUMEXPR_{ExpressionType::FirstExoDerivative, eqr, varr, lag}; - d->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo[block].at({ var, lag })}; - } - for (const auto &[indices, d] : blocks_derivatives_exo_det[block]) - { - auto [eqr, var, lag] = indices; - int eq = getBlockEquationID(block, eqr); - int varr = 0; // Dummy value, actually unused by the bytecode MEX - code_file << FNUMEXPR_{ExpressionType::FirstExodetDerivative, eqr, varr, lag}; - d->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_exo_det[block].at({ var, lag })}; - } - for (const auto &[indices, d] : blocks_derivatives_other_endo[block]) - { - auto [eqr, var, lag] = indices; - int eq = getBlockEquationID(block, eqr); - int varr = 0; // Dummy value, actually unused by the bytecode MEX - code_file << FNUMEXPR_{ExpressionType::FirstOtherEndoDerivative, eqr, varr, lag}; - d->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::dynamicModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG3_{eq, var, lag, blocks_jacob_cols_other_endo[block].at({ var, lag })}; - } - - // Update jump offset for previous JMP - int pos_end_block = code_file.getInstructionCounter(); - code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1}); + writeBlockBytecodeHelper(code_file, block); } - code_file << FENDBLOCK_{} << FEND_{}; + code_file << FEND_{}; } void @@ -1333,60 +1101,6 @@ DynamicModel::printNonZeroHessianEquations(ostream &output) const output << "]"; } -void -DynamicModel::writeBlockBytecodeBinFile(const string &basename, int num, int &u_count_int, - bool &file_open, bool is_two_boundaries) const -{ - int j; - std::ofstream SaveCode; - string filename = basename + "/model/bytecode/dynamic.bin"; - - if (file_open) - SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate); - else - SaveCode.open(filename, ios::out | ios::binary); - if (!SaveCode.is_open()) - { - cerr << R"(Error : Can't open file ")" << filename << R"(" for writing)" << endl; - exit(EXIT_FAILURE); - } - u_count_int = 0; - int block_size = blocks[num].size; - int block_mfs = blocks[num].mfs_size; - int block_recursive = blocks[num].getRecursiveSize(); - for (const auto &[indices, ignore] : blocks_derivatives[num]) - { - auto [eq, var, lag] = indices; - if (lag != 0 && !is_two_boundaries) - continue; - if (eq >= block_recursive && var >= block_recursive) - { - int v = eq - block_recursive; - SaveCode.write(reinterpret_cast(&v), sizeof(v)); - int varr = var - block_recursive + lag * block_mfs; - SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); - SaveCode.write(reinterpret_cast(&lag), sizeof(lag)); - int u = u_count_int + block_mfs; - SaveCode.write(reinterpret_cast(&u), sizeof(u)); - u_count_int++; - } - } - - if (is_two_boundaries) - u_count_int += block_mfs; - for (j = block_recursive; j < block_size; j++) - { - int varr = getBlockVariableID(num, j); - SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); - } - for (j = block_recursive; j < block_size; j++) - { - int eqr = getBlockEquationID(num, j); - SaveCode.write(reinterpret_cast(&eqr), sizeof(eqr)); - } - SaveCode.close(); -} - void DynamicModel::writeDynamicBlockMFile(const string &basename) const { diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 4d44f883..491cbeef 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -140,11 +140,12 @@ private: void writeDynamicPerBlockCFiles(const string &basename) const; //! Writes the code of the block-decomposed model in virtual machine bytecode void writeDynamicBlockBytecode(const string &basename) const; + // Writes derivatives w.r.t. exo, exo det and other endogenous + void writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block, + const temporary_terms_t &temporary_terms_union, + const deriv_node_temp_terms_t &tef_terms) const override; //! Writes the code of the model in virtual machine bytecode void writeDynamicBytecode(const string &basename) const; - //! Adds per-block information for bytecode simulation in a separate .bin file - void writeBlockBytecodeBinFile(const string &basename, int num, int &u_count_int, bool &file_open, - bool is_two_boundaries) const; void writeSetAuxiliaryVariables(const string &basename, bool julia) const; void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const; @@ -175,11 +176,6 @@ private: vector> &blocks_temporary_terms, map> &reference_count) const override; - //! Write derivative bytecode of an equation w.r. to a variable - void writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const; - //! Write chain rule derivative bytecode of an equation w.r. to a variable - void writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const; - SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override; int getLagByDerivID(int deriv_id) const noexcept(false) override; int getSymbIDByDerivID(int deriv_id) const noexcept(false) override; @@ -295,6 +291,12 @@ private: (with a lag), and also as a contemporaneous diff lag auxvar). */ vector getVARDerivIDs(int lhs_symb_id, int lead_lag) const; + int + getBlockJacobianEndoCol(int blk, int var, int lag) const override + { + return blocks_jacob_cols_endo[blk].at({ var, lag }); + } + public: DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 0e70d8da..b03581ca 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -1262,6 +1262,49 @@ ModelTree::writeBytecodeBinFile(const string &filename, bool is_two_boundaries) return u_count; } +int +ModelTree::writeBlockBytecodeBinFile(ofstream &bin_file, int block) const +{ + int u_count {0}; + int block_size {blocks[block].size}; + int block_mfs {blocks[block].mfs_size}; + int block_recursive {blocks[block].getRecursiveSize()}; + BlockSimulationType simulation_type {blocks[block].simulation_type}; + bool is_two_boundaries {simulation_type == BlockSimulationType::solveTwoBoundariesComplete + || simulation_type == BlockSimulationType::solveTwoBoundariesSimple}; + for (const auto &[indices, ignore] : blocks_derivatives[block]) + { + const auto &[eq, var, lag] {indices}; + if (lag != 0 && !is_two_boundaries) + continue; + if (eq >= block_recursive && var >= block_recursive) + { + int v {eq - block_recursive}; + bin_file.write(reinterpret_cast(&v), sizeof v); + int varr {var - block_recursive + lag * block_mfs}; + bin_file.write(reinterpret_cast(&varr), sizeof varr); + bin_file.write(reinterpret_cast(&lag), sizeof lag); + int u {u_count + block_mfs}; + bin_file.write(reinterpret_cast(&u), sizeof u); + u_count++; + } + } + + if (is_two_boundaries) + u_count += block_mfs; + for (int j {block_recursive}; j < block_size; j++) + { + int varr {getBlockVariableID(block, j)}; + bin_file.write(reinterpret_cast(&varr), sizeof varr); + } + for (int j {block_recursive}; j < block_size; j++) + { + int eqr {getBlockEquationID(block, j)}; + bin_file.write(reinterpret_cast(&eqr), sizeof eqr); + } + return u_count; +} + void ModelTree::writeLatexModelFile(const string &mod_basename, const string &latex_basename, ExprNodeOutputType output_type, bool write_equation_tags) const { @@ -1823,3 +1866,11 @@ ModelTree::getRHSFromLHS(expr_t lhs) const return eq->arg2; throw ExprNode::MatchFailureException{"Cannot find an equation with the requested LHS"}; } + +void +ModelTree::writeBlockBytecodeAdditionalDerivatives([[maybe_unused]] BytecodeWriter &code_file, + [[maybe_unused]] int block, + [[maybe_unused]] const temporary_terms_t &temporary_terms_union, + [[maybe_unused]] const deriv_node_temp_terms_t &tef_terms) const +{ +} diff --git a/src/ModelTree.hh b/src/ModelTree.hh index a195cf79..f11ac7c9 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -251,6 +251,9 @@ protected: file. Returns the number of first derivatives w.r.t. endogenous variables */ int writeBytecodeBinFile(const string &filename, bool is_two_boundaries) const; + //! Adds per-block information for bytecode simulation in a separate .bin file + int writeBlockBytecodeBinFile(ofstream &bin_file, int block) const; + //! Fixes output when there are more than 32 nested parens, Issue #1201 void fixNestedParenthesis(ostringstream &output, map &tmp_paren_vars, bool &message_printed) const; //! Tests if string contains more than 32 nested parens, Issue #1201 @@ -279,6 +282,17 @@ protected: template void writeBytecodeHelper(BytecodeWriter &code_file) const; + // Helper for writing blocks in bytecode + template + void writeBlockBytecodeHelper(BytecodeWriter &code_file, int block) const; + + /* Write additional derivatives w.r.t. to exogenous, exogenous det and other endo + in block+bytecode mode. Does nothing by default, but overriden by + DynamicModel which needs those. */ + virtual void writeBlockBytecodeAdditionalDerivatives(BytecodeWriter &code_file, int block, + const temporary_terms_t &temporary_terms_union, + const deriv_node_temp_terms_t &tef_terms) const; + /* Helper for writing JSON output for residuals and derivatives. Returns mlv and derivatives output at each derivation order. */ template @@ -430,6 +444,10 @@ protected: Assumes that derivatives have already been computed. */ map, expr_t> collectFirstOrderDerivativesEndogenous(); + /* Get column number within Jacobian of a given block. + “var” is the block-specific endogenous variable index. */ + virtual int getBlockJacobianEndoCol(int blk, int var, int lag) const = 0; + private: //! Internal helper for the copy constructor and assignment operator /*! Copies all the structures that contain ExprNode*, by the converting the @@ -1205,6 +1223,226 @@ ModelTree::writeBytecodeHelper(BytecodeWriter &code_file) const code_file << FENDBLOCK_{} << FEND_{}; } +template +void +ModelTree::writeBlockBytecodeHelper(BytecodeWriter &code_file, int block) const +{ + constexpr ExprNodeBytecodeOutputType output_type + { dynamic ? ExprNodeBytecodeOutputType::dynamicModel : ExprNodeBytecodeOutputType::staticModel }; + constexpr ExprNodeBytecodeOutputType assignment_lhs_output_type + { dynamic ? ExprNodeBytecodeOutputType::dynamicAssignmentLHS : ExprNodeBytecodeOutputType::staticAssignmentLHS }; + + const BlockSimulationType simulation_type {blocks[block].simulation_type}; + const int block_size {blocks[block].size}; + const int block_mfs {blocks[block].mfs_size}; + const int block_recursive {blocks[block].getRecursiveSize()}; + + temporary_terms_t temporary_terms_union; + deriv_node_temp_terms_t tef_terms; + + auto write_eq_tt = [&](int eq) + { + for (auto it : blocks_temporary_terms[block][eq]) + { + if (dynamic_cast(it)) + it->writeBytecodeExternalFunctionOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + + code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)}; + it->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + if constexpr(dynamic) + code_file << FSTPT_{blocks_temporary_terms_idxs.at(it)}; + else + code_file << FSTPST_{blocks_temporary_terms_idxs.at(it)}; + temporary_terms_union.insert(it); + } + }; + + // The equations + for (int i {0}; i < block_size; i++) + { + write_eq_tt(i); + + switch (simulation_type) + { + evaluation: + case BlockSimulationType::evaluateBackward: + case BlockSimulationType::evaluateForward: + code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; + if (EquationType equ_type {getBlockEquationType(block, i)}; + equ_type == EquationType::evaluate) + { + BinaryOpNode *eq_node {getBlockEquationExpr(block, i)}; + expr_t lhs {eq_node->arg1}; + expr_t rhs {eq_node->arg2}; + rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + lhs->writeBytecodeOutput(code_file, assignment_lhs_output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + } + else if (equ_type == EquationType::evaluateRenormalized) + { + BinaryOpNode *eq_node {getBlockEquationRenormalizedExpr(block, i)}; + expr_t lhs {eq_node->arg1}; + expr_t rhs {eq_node->arg2}; + rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + lhs->writeBytecodeOutput(code_file, assignment_lhs_output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + } + break; + case BlockSimulationType::solveBackwardComplete: + case BlockSimulationType::solveForwardComplete: + case BlockSimulationType::solveTwoBoundariesComplete: + case BlockSimulationType::solveTwoBoundariesSimple: + if (i < block_recursive) + goto evaluation; + [[fallthrough]]; + default: + code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; + BinaryOpNode *eq_node {getBlockEquationExpr(block, i)}; + expr_t lhs {eq_node->arg1}; + expr_t rhs {eq_node->arg2}; + lhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + rhs->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive}; + } + } + code_file << FENDEQU_{}; + + /* If the block is not of type “evaluate backward/forward”, then we write + the temporary terms for derivatives at this point, i.e. before the + JMPIFEVAL, because they will be needed in both “simulate” and + “evaluate” modes. */ + if (simulation_type != BlockSimulationType::evaluateBackward + && simulation_type != BlockSimulationType::evaluateForward) + write_eq_tt(blocks[block].size); + + // Get the current code_file position and jump if evaluating + int pos_jmpifeval {code_file.getInstructionCounter()}; + code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being + + /* Write the derivatives for the “simulate” mode (not needed if the block + is of type “evaluate backward/forward”) */ + if (simulation_type != BlockSimulationType::evaluateBackward + && simulation_type != BlockSimulationType::evaluateForward) + { + switch (simulation_type) + { + case BlockSimulationType::solveBackwardSimple: + case BlockSimulationType::solveForwardSimple: + { + int eqr {getBlockEquationID(block, 0)}; + int varr {getBlockVariableID(block, 0)}; + code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, 0}; + // Get contemporaneous derivative of the single variable in the block + if (auto it { blocks_derivatives[block].find({ 0, 0, 0 }) }; + it != blocks_derivatives[block].end()) + it->second->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + else + code_file << FLDZ_{}; + code_file << FSTPG_{0}; + } + break; + + case BlockSimulationType::solveBackwardComplete: + case BlockSimulationType::solveForwardComplete: + case BlockSimulationType::solveTwoBoundariesComplete: + case BlockSimulationType::solveTwoBoundariesSimple: + { + // For each equation, stores a list of tuples (index_u, var, lag) + vector>> Uf(symbol_table.endo_nbr()); + + for (int count_u {block_mfs}; + const auto &[indices, ignore] : blocks_derivatives[block]) + { + const auto &[eq, var, lag] {indices}; + int eqr {getBlockEquationID(block, eq)}; + int varr {getBlockVariableID(block, var)}; + if (eq >= block_recursive && var >= block_recursive) + { + if constexpr(dynamic) + if (lag != 0 + && (simulation_type == BlockSimulationType::solveForwardComplete + || simulation_type == BlockSimulationType::solveBackwardComplete)) + continue; + code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag}; + if (auto it { blocks_derivatives[block].find({ eq, var, lag }) }; + it != blocks_derivatives[block].end()) + it->second->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + else + code_file << FLDZ_{}; + if constexpr(dynamic) + code_file << FSTPU_{count_u}; + else + code_file << FSTPSU_{count_u}; + Uf[eqr].emplace_back(count_u, varr, lag); + count_u++; + } + } + for (int i {0}; i < block_size; i++) + if (i >= block_recursive) + { + code_file << FLDR_{i-block_recursive} << FLDZ_{}; + + int eqr {getBlockEquationID(block, i)}; + for (const auto &[index_u, var, lag] : Uf[eqr]) + { + if constexpr(dynamic) + code_file << FLDU_{index_u} + << FLDV_{SymbolType::endogenous, var, lag}; + else + code_file << FLDSU_{index_u} + << FLDSV_{SymbolType::endogenous, var}; + code_file << FBINARY_{BinaryOpcode::times} + << FCUML_{}; + } + code_file << FBINARY_{BinaryOpcode::minus}; + if constexpr(dynamic) + code_file << FSTPU_{i - block_recursive}; + else + code_file << FSTPSU_{i - block_recursive}; + } + } + break; + default: + break; + } + } + + // Jump unconditionally after the block + int pos_jmp {code_file.getInstructionCounter()}; + code_file << FJMP_{0}; // Use 0 as jump offset for the time being + // Update jump offset for previous JMPIFEVAL + code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval}); + + /* If the block is of type “evaluate backward/forward”, then write the + temporary terms for derivatives at this point, because they have not + been written before the JMPIFEVAL. */ + if (simulation_type == BlockSimulationType::evaluateBackward + || simulation_type == BlockSimulationType::evaluateForward) + write_eq_tt(blocks[block].size); + + // Write the derivatives for the “evaluate” mode + for (const auto &[indices, d] : blocks_derivatives[block]) + { + const auto &[eq, var, lag] {indices}; + int eqr {getBlockEquationID(block, eq)}; + int varr {getBlockVariableID(block, var)}; + code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, lag}; + d->writeBytecodeOutput(code_file, output_type, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); + if constexpr(dynamic) + code_file << FSTPG3_{eq, var, lag, getBlockJacobianEndoCol(block, var, lag)}; + else + code_file << FSTPG2_{eq, getBlockJacobianEndoCol(block, var, lag)}; + } + + /* Write derivatives w.r.t. exo, exo det and other endogenous, but only in + dynamic mode */ + writeBlockBytecodeAdditionalDerivatives(code_file, block, temporary_terms_union, tef_terms); + + // Update jump offset for previous JMP + int pos_end_block {code_file.getInstructionCounter()}; + code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1}); + + code_file << FENDBLOCK_{}; +} + template pair> ModelTree::writeJsonComputingPassOutputHelper(bool writeDetails) const diff --git a/src/StaticModel.cc b/src/StaticModel.cc index ad32e007..45d8723c 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -94,26 +94,6 @@ StaticModel::StaticModel(const DynamicModel &m) : user_set_compiler = m.user_set_compiler; } -void -StaticModel::writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const -{ - if (auto it = derivatives[1].find({ eq, getDerivID(symbol_table.getID(SymbolType::endogenous, symb_id), 0) }); - it != derivatives[1].end()) - it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms, temporary_terms_idxs, tef_terms); - else - code_file << FLDZ_{}; -} - -void -StaticModel::writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const -{ - if (auto it = blocks_derivatives[blk].find({ eq, var, lag }); - it != blocks_derivatives[blk].end()) - it->second->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms, temporary_terms_idxs, tef_terms); - else - code_file << FLDZ_{}; -} - void StaticModel::writeStaticPerBlockMFiles(const string &basename) const { @@ -288,368 +268,45 @@ StaticModel::writeStaticBytecode(const string &basename) const void StaticModel::writeStaticBlockBytecode(const string &basename) const { - struct Uff_l - { - int u, var, lag; - Uff_l *pNext; - }; + BytecodeWriter code_file {basename + "/model/bytecode/static.cod"}; - struct Uff - { - Uff_l *Ufl, *Ufl_First; - }; - - int i, v; - string tmp_s; - ostringstream tmp_output; - expr_t lhs = nullptr, rhs = nullptr; - BinaryOpNode *eq_node; - Uff Uf[symbol_table.endo_nbr()]; - map reference_count; - vector feedback_variables; - bool file_open = false; - - BytecodeWriter code_file{basename + "/model/bytecode/static.cod"}; - - //Temporary variables declaration - code_file << FDIMST_{static_cast(blocks_temporary_terms_idxs.size())}; - - temporary_terms_t temporary_terms_union; - - for (int block = 0; block < static_cast(blocks.size()); block++) + const string bin_filename {basename + "/model/bytecode/static.bin"}; + ofstream bin_file {bin_filename, ios::out | ios::binary}; + if (!bin_file.is_open()) { - feedback_variables.clear(); - if (block > 0) - code_file << FENDBLOCK_{}; - int count_u; - int u_count_int = 0; - BlockSimulationType simulation_type = blocks[block].simulation_type; - int block_size = blocks[block].size; - int block_mfs = blocks[block].mfs_size; - int block_recursive = blocks[block].getRecursiveSize(); - - if (simulation_type == BlockSimulationType::solveTwoBoundariesSimple - || simulation_type == BlockSimulationType::solveTwoBoundariesComplete - || simulation_type == BlockSimulationType::solveBackwardComplete - || simulation_type == BlockSimulationType::solveForwardComplete) - { - writeBlockBytecodeBinFile(basename, block, u_count_int, file_open); - file_open = true; - } - - code_file << FBEGINBLOCK_{block_mfs, - simulation_type, - blocks[block].first_equation, - block_size, - endo_idx_block2orig, - eq_idx_block2orig, - blocks[block].linear, - symbol_table.endo_nbr(), - 0, - 0, - u_count_int, - block_size}; - - // Get the current code_file position and jump if evaluating - int pos_jmpifeval = code_file.getInstructionCounter(); - code_file << FJMPIFEVAL_{0}; // Use 0 as jump offset for the time being - - //The Temporary terms - deriv_node_temp_terms_t tef_terms; - /* Keep a backup of temporary_terms_union here, since temp. terms are - written a second time below. This is probably unwanted… */ - temporary_terms_t ttu_old = temporary_terms_union; - - auto write_eq_tt = [&](int eq) - { - for (auto it : blocks_temporary_terms[block][eq]) - { - if (dynamic_cast(it)) - it->writeBytecodeExternalFunctionOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - - code_file << FNUMEXPR_{ExpressionType::TemporaryTerm, blocks_temporary_terms_idxs.at(it)}; - it->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPST_{blocks_temporary_terms_idxs.at(it)}; - temporary_terms_union.insert(it); - } - }; - - for (i = 0; i < block_size; i++) - { - write_eq_tt(i); - - // The equations - int variable_ID, equation_ID; - EquationType equ_type; - switch (simulation_type) - { - evaluation: - case BlockSimulationType::evaluateBackward: - case BlockSimulationType::evaluateForward: - equ_type = getBlockEquationType(block, i); - code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; - if (equ_type == EquationType::evaluate) - { - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - } - else if (equ_type == EquationType::evaluateRenormalized) - { - eq_node = getBlockEquationRenormalizedExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - } - break; - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - if (i < block_recursive) - goto evaluation; - variable_ID = getBlockVariableID(block, i); - equation_ID = getBlockEquationID(block, i); - feedback_variables.push_back(variable_ID); - Uf[equation_ID].Ufl = nullptr; - goto end; - default: - end: - code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - - code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive}; - } - } - code_file << FENDEQU_{}; - - // The Jacobian if we have to solve the block - if (simulation_type != BlockSimulationType::evaluateBackward - && simulation_type != BlockSimulationType::evaluateForward) - { - // Write temporary terms for derivatives - write_eq_tt(blocks[block].size); - - switch (simulation_type) - { - case BlockSimulationType::solveBackwardSimple: - case BlockSimulationType::solveForwardSimple: - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, 0, 0}; - writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG_{0}; - break; - - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - count_u = feedback_variables.size(); - for (const auto &[indices, ignore2] : blocks_derivatives[block]) - { - auto [eq, var, ignore] = indices; - int eqr = getBlockEquationID(block, eq); - int varr = getBlockVariableID(block, var); - if (eq >= block_recursive && var >= block_recursive) - { - if (!Uf[eqr].Ufl) - { - Uf[eqr].Ufl = static_cast(malloc(sizeof(Uff_l))); - Uf[eqr].Ufl_First = Uf[eqr].Ufl; - } - else - { - Uf[eqr].Ufl->pNext = static_cast(malloc(sizeof(Uff_l))); - Uf[eqr].Ufl = Uf[eqr].Ufl->pNext; - } - Uf[eqr].Ufl->pNext = nullptr; - Uf[eqr].Ufl->u = count_u; - Uf[eqr].Ufl->var = varr; - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr}; - writeBytecodeChainRuleDerivative(code_file, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPSU_{count_u}; - count_u++; - } - } - for (i = 0; i < block_size; i++) - if (i >= block_recursive) - { - code_file << FLDR_{i-block_recursive} << FLDZ_{}; - - v = getBlockEquationID(block, i); - for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext) - code_file << FLDSU_{Uf[v].Ufl->u} - << FLDSV_{SymbolType::endogenous, Uf[v].Ufl->var} - << FBINARY_{BinaryOpcode::times} - << FCUML_{}; - Uf[v].Ufl = Uf[v].Ufl_First; - while (Uf[v].Ufl) - { - Uf[v].Ufl_First = Uf[v].Ufl->pNext; - free(Uf[v].Ufl); - Uf[v].Ufl = Uf[v].Ufl_First; - } - code_file << FBINARY_{BinaryOpcode::minus} - << FSTPSU_{i - block_recursive}; - } - break; - default: - break; - } - } - - // Jump unconditionally after the block - int pos_jmp = code_file.getInstructionCounter(); - code_file << FJMP_{0}; // Use 0 as jump offset for the time being - // Update jump offset for previous JMPIFEVAL - code_file.overwriteInstruction(pos_jmpifeval, FJMPIFEVAL_{pos_jmp-pos_jmpifeval}); - - tef_terms.clear(); - temporary_terms_union = ttu_old; - - for (i = 0; i < block_size; i++) - { - write_eq_tt(i); - - // The equations - int variable_ID, equation_ID; - EquationType equ_type; - switch (simulation_type) - { - evaluation_l: - case BlockSimulationType::evaluateBackward: - case BlockSimulationType::evaluateForward: - equ_type = getBlockEquationType(block, i); - code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; - if (equ_type == EquationType::evaluate) - { - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - } - else if (equ_type == EquationType::evaluateRenormalized) - { - eq_node = getBlockEquationRenormalizedExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticAssignmentLHS, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - } - break; - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - if (i < block_recursive) - goto evaluation_l; - variable_ID = getBlockVariableID(block, i); - equation_ID = getBlockEquationID(block, i); - feedback_variables.push_back(variable_ID); - Uf[equation_ID].Ufl = nullptr; - goto end_l; - default: - end_l: - code_file << FNUMEXPR_{ExpressionType::ModelEquation, getBlockEquationID(block, i)}; - eq_node = getBlockEquationExpr(block, i); - lhs = eq_node->arg1; - rhs = eq_node->arg2; - lhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - rhs->writeBytecodeOutput(code_file, ExprNodeBytecodeOutputType::staticModel, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - - code_file << FBINARY_{BinaryOpcode::minus} << FSTPR_{i - block_recursive}; - } - } - code_file << FENDEQU_{}; - - // The Jacobian if we have to solve the block determinsitic bloc - - // Write temporary terms for derivatives - write_eq_tt(blocks[block].size); - - switch (simulation_type) - { - case BlockSimulationType::solveBackwardSimple: - case BlockSimulationType::solveForwardSimple: - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, 0, 0}; - writeBytecodeDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG2_{0, 0}; - break; - case BlockSimulationType::evaluateBackward: - case BlockSimulationType::evaluateForward: - case BlockSimulationType::solveBackwardComplete: - case BlockSimulationType::solveForwardComplete: - count_u = feedback_variables.size(); - for (const auto &[indices, ignore2] : blocks_derivatives[block]) - { - auto &[eq, var, ignore] = indices; - int eqr = getBlockEquationID(block, eq); - int varr = getBlockVariableID(block, var); - code_file << FNUMEXPR_{ExpressionType::FirstEndoDerivative, eqr, varr, 0}; - writeBytecodeChainRuleDerivative(code_file, block, eq, var, 0, temporary_terms_union, blocks_temporary_terms_idxs, tef_terms); - code_file << FSTPG2_{eq, var}; - } - break; - default: - break; - } - // Set codefile position to previous JMP_ and set the number of instructions to jump - // Update jump offset for previous JMP - int pos_end_block = code_file.getInstructionCounter(); - code_file.overwriteInstruction(pos_jmp, FJMP_{pos_end_block-pos_jmp-1}); - } - code_file << FENDBLOCK_{} << FEND_{}; -} - -void -StaticModel::writeBlockBytecodeBinFile(const string &basename, int num, - int &u_count_int, bool &file_open) const -{ - int j; - std::ofstream SaveCode; - string filename = basename + "/model/bytecode/static.bin"; - if (file_open) - SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate); - else - SaveCode.open(filename, ios::out | ios::binary); - if (!SaveCode.is_open()) - { - cerr << "Error : Can't open file " << filename << " for writing" << endl; + cerr << R"(Error : Can't open file ")" << bin_filename << R"(" for writing)" << endl; exit(EXIT_FAILURE); } - u_count_int = 0; - int block_size = blocks[num].size; - int block_mfs = blocks[num].mfs_size; - int block_recursive = blocks[num].getRecursiveSize(); - for (const auto &[indices, ignore2] : blocks_derivatives[num]) - { - auto [eq, var, ignore] = indices; - int lag = 0; - if (eq >= block_recursive && var >= block_recursive) - { - int v = eq - block_recursive; - SaveCode.write(reinterpret_cast(&v), sizeof(v)); - int varr = var - block_recursive; - SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); - SaveCode.write(reinterpret_cast(&lag), sizeof(lag)); - int u = u_count_int + block_mfs; - SaveCode.write(reinterpret_cast(&u), sizeof(u)); - u_count_int++; - } - } - for (j = block_recursive; j < block_size; j++) + // Temporary variables declaration + code_file << FDIMST_{static_cast(blocks_temporary_terms_idxs.size())}; + + for (int block {0}; block < static_cast(blocks.size()); block++) { - int varr = getBlockVariableID(num, j); - SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); + const BlockSimulationType simulation_type {blocks[block].simulation_type}; + const int block_size {blocks[block].size}; + + const int u_count {simulation_type == BlockSimulationType::solveBackwardComplete + || simulation_type == BlockSimulationType::solveForwardComplete + ? writeBlockBytecodeBinFile(bin_file, block) + : 0}; + + code_file << FBEGINBLOCK_{blocks[block].mfs_size, + simulation_type, + blocks[block].first_equation, + block_size, + endo_idx_block2orig, + eq_idx_block2orig, + blocks[block].linear, + symbol_table.endo_nbr(), + 0, + 0, + u_count, + block_size}; + + writeBlockBytecodeHelper(code_file, block); } - for (j = block_recursive; j < block_size; j++) - { - int eqr = getBlockEquationID(num, j); - SaveCode.write(reinterpret_cast(&eqr), sizeof(eqr)); - } - SaveCode.close(); + code_file << FEND_{}; } void diff --git a/src/StaticModel.hh b/src/StaticModel.hh index 8948d846..facb4e09 100644 --- a/src/StaticModel.hh +++ b/src/StaticModel.hh @@ -65,10 +65,6 @@ private: //! Writes the code of the model in virtual machine bytecode void writeStaticBytecode(const string &basename) const; - //! Adds per-block information for bytecode simulation in a separate .bin file - void writeBlockBytecodeBinFile(const string &basename, int num, - int &u_count_int, bool &file_open) const; - //! Computes jacobian and prepares for equation normalization /*! Using values from initval/endval blocks and parameter initializations: - computes the jacobian for the model w.r. to contemporaneous variables @@ -76,11 +72,6 @@ private: */ void evaluateJacobian(const eval_context_t &eval_context, jacob_map_t *j_m, bool dynamic); - //! Write derivative bytecode of an equation w.r. to a variable - void writeBytecodeDerivative(BytecodeWriter &code_file, int eq, int symb_id, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const; - //! Write chain rule derivative bytecode of an equation w.r. to a variable - void writeBytecodeChainRuleDerivative(BytecodeWriter &code_file, int blk, int eq, int var, int lag, const temporary_terms_t &temporary_terms, const temporary_terms_idxs_t &temporary_terms_idxs, const deriv_node_temp_terms_t &tef_terms) const; - SymbolType getTypeByDerivID(int deriv_id) const noexcept(false) override; int getLagByDerivID(int deriv_id) const noexcept(false) override; int getSymbIDByDerivID(int deriv_id) const noexcept(false) override; @@ -113,6 +104,12 @@ private: //! Create a legacy *_static.m file for MATLAB/Octave not yet using the temporary terms array interface void writeStaticMCompatFile(const string &name) const; + int + getBlockJacobianEndoCol([[maybe_unused]] int blk, int var, [[maybe_unused]] int lag) const override + { + return var; + } + public: StaticModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants,