diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc index b1b3ec3fd..216c22c62 100644 --- a/mex/sources/bytecode/Evaluate.cc +++ b/mex/sources/bytecode/Evaluate.cc @@ -29,14 +29,9 @@ #include "CommonEnums.hh" #include "ErrorHandling.hh" -Evaluate::Evaluate(bool steady_state_arg, const BasicSymbolTable &symbol_table_arg) : +Evaluate::Evaluate(const filesystem::path &codfile, bool steady_state_arg, const BasicSymbolTable &symbol_table_arg) : symbol_table {symbol_table_arg}, steady_state {steady_state_arg} -{ -} - -void -Evaluate::loadCodeFile(const filesystem::path &codfile) { ifstream CompiledCode {codfile, ios::in | ios::binary | ios::ate}; if (!CompiledCode.is_open()) diff --git a/mex/sources/bytecode/Evaluate.hh b/mex/sources/bytecode/Evaluate.hh index 2d6a0f869..dcb004d7c 100644 --- a/mex/sources/bytecode/Evaluate.hh +++ b/mex/sources/bytecode/Evaluate.hh @@ -65,6 +65,13 @@ private: string error_location(it_code_type expr_begin, it_code_type faulty_op, int it_) const; + /* Prints a bytecode expression in human readable form. + If faulty_op is not default constructed, it should point to a tag within + the expression that created a floating point exception, in which case the + corresponding mathematical operator will be printed within braces. + The second output argument points to the tag past the expression. */ + pair print_expression(const it_code_type &expr_begin, const optional &faulty_op = nullopt) const; + FBEGINBLOCK_ * currentBlockTag() const { @@ -78,15 +85,10 @@ private: return instructions_list.begin() + begin_block[block_num] + 1; } -protected: - void evaluateBlock(int it_, double *__restrict__ y, const double *__restrict__ ya, int y_size, double *__restrict__ x, int nb_row_x, double *__restrict__ params, const double *__restrict__ steady_y, double *__restrict__ u, int Per_u_, double *__restrict__ T, int T_nrows, map &TEF, map, double> &TEFD, map, double> &TEFDD, double *__restrict__ r, double *__restrict__ g1, double *__restrict__ jacob, double *__restrict__ jacob_exo, double *__restrict__ jacob_exo_det, bool evaluate, bool no_derivatives); +public: + Evaluate(const filesystem::path &codfile, bool steady_state_arg, const BasicSymbolTable &symbol_table_arg); - /* Prints a bytecode expression in human readable form. - If faulty_op is not default constructed, it should point to a tag within - the expression that created a floating point exception, in which case the - corresponding mathematical operator will be printed within braces. - The second output argument points to the tag past the expression. */ - pair print_expression(const it_code_type &expr_begin, const optional &faulty_op = nullopt) const; + void evaluateBlock(int it_, double *__restrict__ y, const double *__restrict__ ya, int y_size, double *__restrict__ x, int nb_row_x, double *__restrict__ params, const double *__restrict__ steady_y, double *__restrict__ u, int Per_u_, double *__restrict__ T, int T_nrows, map &TEF, map, double> &TEFD, map, double> &TEFDD, double *__restrict__ r, double *__restrict__ g1, double *__restrict__ jacob, double *__restrict__ jacob_exo, double *__restrict__ jacob_exo_det, bool evaluate, bool no_derivatives); // Prints current block void printCurrentBlock(); @@ -146,11 +148,6 @@ protected: return currentBlockTag()->get_det_exo_size(); } -public: - Evaluate(bool steady_state_arg, const BasicSymbolTable &symbol_table_arg); - // TODO: integrate into the constructor - void loadCodeFile(const filesystem::path &codfile); - int get_block_number() const { diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index 60476d725..3e0c2bc78 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -28,14 +28,14 @@ constexpr double BIG = 1.0e+8, SMALL = 1.0e-5; -Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, +Interpreter::Interpreter(Evaluate &evaluator_arg, double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *direction_arg, size_t y_size_arg, size_t nb_row_x_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, size_t size_of_direction_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg, bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg, bool steady_state_arg, bool block_decomposed_arg, bool print_it_arg, int col_x_arg, int col_y_arg, const BasicSymbolTable &symbol_table_arg) -: dynSparseMatrix {y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, block_decomposed_arg, periods_arg, minimal_solving_periods_arg, symbol_table_arg, print_error_arg} +: dynSparseMatrix {evaluator_arg, y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, block_decomposed_arg, periods_arg, minimal_solving_periods_arg, symbol_table_arg, print_error_arg} { params = params_arg; y = y_arg; @@ -615,20 +615,11 @@ Interpreter::simulate_a_block(const vector_table_conditional_local_type &vector_ return NO_ERROR_ON_EXIT; } -void -Interpreter::ReadCodeFile(const string &file_name) -{ - filesystem::path codfile {file_name + "/model/bytecode/" + (block_decomposed ? "block/" : "") - + (steady_state ? "static" : "dynamic") + ".cod"}; - - loadCodeFile(codfile); -} - void Interpreter::check_for_controlled_exo_validity(int current_block, const vector &sconstrained_extended_path) { - vector exogenous {getCurrentBlockExogenous()}; - vector endogenous {getCurrentBlockEndogenous()}; + vector exogenous {evaluator.getCurrentBlockExogenous()}; + vector endogenous {evaluator.getCurrentBlockEndogenous()}; for (auto & it : sconstrained_extended_path) { if (find(endogenous.begin(), endogenous.end(), it.exo_num) != endogenous.end() @@ -660,23 +651,25 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool { initializeTemporaryTerms(global_temporary_terms); - if (block >= get_block_number()) + int nb_blocks {evaluator.get_block_number()}; + + if (block >= nb_blocks) throw FatalException {"Interpreter::MainLoop: Input argument block = " + to_string(block+1) + " is greater than the number of blocks in the model (" - + to_string(get_block_number()) + " see M_.block_structure" + (steady_state ? "_stat" : "") + ".block)"}; + + to_string(nb_blocks) + " see M_.block_structure" + (steady_state ? "_stat" : "") + ".block)"}; vector blocks; if (block < 0) { - blocks.resize(get_block_number()); + blocks.resize(nb_blocks); iota(blocks.begin(), blocks.end(), 0); } else blocks.push_back(block); - jacobian_block.resize(get_block_number()); - jacobian_exo_block.resize(get_block_number()); - jacobian_det_exo_block.resize(get_block_number()); + jacobian_block.resize(nb_blocks); + jacobian_exo_block.resize(nb_blocks); + jacobian_det_exo_block.resize(nb_blocks); double max_res_local = 0; int max_res_idx_local = 0; @@ -691,13 +684,13 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool for (int current_block : blocks) { - gotoBlock(current_block); + evaluator.gotoBlock(current_block); block_num = current_block; - size = getCurrentBlockSize(); - type = getCurrentBlockType(); - is_linear = isCurrentBlockLinear(); - Block_Contain = getCurrentBlockEquationsAndVariables(); - u_count_int = getCurrentBlockUCount(); + size = evaluator.getCurrentBlockSize(); + type = evaluator.getCurrentBlockType(); + is_linear = evaluator.isCurrentBlockLinear(); + Block_Contain = evaluator.getCurrentBlockEquationsAndVariables(); + u_count_int = evaluator.getCurrentBlockUCount(); if (constrained) check_for_controlled_exo_validity(current_block, sconstrained_extended_path); @@ -707,23 +700,23 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool residual = vector(size); else residual = vector(size*(periods+y_kmin)); - printCurrentBlock(); + evaluator.printCurrentBlock(); } else if (evaluate) { #ifdef DEBUG mexPrintf("jacobian_block=mxCreateDoubleMatrix(%d, %d, mxREAL)\n", size, getCurrentBlockNbColJacob()); #endif - jacobian_block[current_block] = mxCreateDoubleMatrix(size, getCurrentBlockNbColJacob(), mxREAL); + jacobian_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockNbColJacob(), mxREAL); if (!steady_state) { #ifdef DEBUG - mexPrintf("allocates jacobian_exo_block( %d, %d, mxREAL)\n", size, getCurrentBlockExoSize()); + mexPrintf("allocates jacobian_exo_block( %d, %d, mxREAL)\n", size, evaluator.getCurrentBlockExoSize()); mexPrintf("(0) Allocating Jacobian\n"); #endif - jacobian_exo_block[current_block] = mxCreateDoubleMatrix(size, getCurrentBlockExoSize(), mxREAL); - jacobian_det_exo_block[current_block] = mxCreateDoubleMatrix(size, getCurrentBlockExoDetSize(), mxREAL); + jacobian_exo_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockExoSize(), mxREAL); + jacobian_det_exo_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockExoDetSize(), mxREAL); } if (block >= 0) { @@ -743,9 +736,9 @@ Interpreter::MainLoop(const string &bin_basename, bool evaluate, int block, bool bool result; if (sconstrained_extended_path.size()) { - jacobian_block[current_block] = mxCreateDoubleMatrix(size, getCurrentBlockNbColJacob(), mxREAL); - jacobian_exo_block[current_block] = mxCreateDoubleMatrix(size, getCurrentBlockExoSize(), mxREAL); - jacobian_det_exo_block[current_block] = mxCreateDoubleMatrix(size, getCurrentBlockExoDetSize(), mxREAL); + jacobian_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockNbColJacob(), mxREAL); + jacobian_exo_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockExoSize(), mxREAL); + jacobian_det_exo_block[current_block] = mxCreateDoubleMatrix(size, evaluator.getCurrentBlockExoDetSize(), mxREAL); residual = vector(size*(periods+y_kmin)); result = simulate_a_block(vector_table_conditional_local, block >= 0, bin_basename); } @@ -810,7 +803,6 @@ Interpreter::elastic(string str, unsigned int len, bool left) pair> Interpreter::extended_path(const string &file_name, bool evaluate, int block, int nb_periods, const vector &sextended_path, const vector &sconstrained_extended_path, const vector &dates, const table_conditional_global_type &table_conditional_global) { - ReadCodeFile(file_name); size_t size_of_direction = y_size*col_y*sizeof(double); auto *y_save = static_cast(mxMalloc(size_of_direction)); test_mxMalloc(y_save, __LINE__, __FILE__, __func__, size_of_direction); @@ -907,8 +899,6 @@ Interpreter::extended_path(const string &file_name, bool evaluate, int block, in pair> Interpreter::compute_blocks(const string &file_name, bool evaluate, int block) { - ReadCodeFile(file_name); - //The big loop on intructions vector s_plan_junk; vector_table_conditional_local_type vector_table_conditional_local_junk; @@ -923,7 +913,7 @@ Interpreter::compute_blocks(const string &file_name, bool evaluate, int block) void Interpreter::initializeTemporaryTerms(bool global_temporary_terms) { - int ntt { getNumberOfTemporaryTerms() }; + int ntt { evaluator.getNumberOfTemporaryTerms() }; if (steady_state) { diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh index efe90d342..528d4160e 100644 --- a/mex/sources/bytecode/Interpreter.hh +++ b/mex/sources/bytecode/Interpreter.hh @@ -51,7 +51,7 @@ protected: int simulate_a_block(const vector_table_conditional_local_type &vector_table_conditional_local, bool single_block, const string &bin_base_name); string elastic(string str, unsigned int len, bool left); public: - Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, + Interpreter(Evaluate &evaluator_arg, double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *direction_arg, size_t y_size_arg, size_t nb_row_x_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, size_t size_of_direction_arg, int y_decal_arg, double markowitz_c_arg, @@ -62,7 +62,6 @@ public: pair> compute_blocks(const string &file_name, bool evaluate, int block); void check_for_controlled_exo_validity(int current_block, const vector &sconstrained_extended_path); pair> MainLoop(const string &bin_basename, bool evaluate, int block, bool constrained, const vector &sconstrained_extended_path, const vector_table_conditional_local_type &vector_table_conditional_local); - void ReadCodeFile(const string &file_name); inline mxArray * get_jacob(int block_num) const diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index f9f2ccfa3..56828f97c 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -24,13 +24,13 @@ #include "SparseMatrix.hh" -dynSparseMatrix::dynSparseMatrix(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, +dynSparseMatrix::dynSparseMatrix(Evaluate &evaluator_arg, int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg, bool print_error_arg) : - Evaluate {steady_state_arg, symbol_table_arg}, symbol_table {symbol_table_arg}, steady_state {steady_state_arg}, block_decomposed {block_decomposed_arg}, + evaluator {evaluator_arg}, minimal_solving_periods {minimal_solving_periods_arg}, print_it {print_it_arg}, y_size {y_size_arg}, @@ -1845,7 +1845,7 @@ dynSparseMatrix::compute_block_time(int Per_u_, bool evaluate, bool no_derivativ try { - evaluateBlock(it_, y, ya, y_size, x, nb_row_x, params, steady_y, u, Per_u_, T, periods+y_kmin+y_kmax, TEF, TEFD, TEFDD, r, g1, jacob, jacob_exo, jacob_exo_det, evaluate, no_derivatives); + evaluator.evaluateBlock(it_, y, ya, y_size, x, nb_row_x, params, steady_y, u, Per_u_, T, periods+y_kmin+y_kmax, TEF, TEFD, TEFDD, r, g1, jacob, jacob_exo, jacob_exo_det, evaluate, no_derivatives); } catch (FloatingPointException &e) { diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index 6d26f6d7e..e896f1d4c 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -65,10 +65,10 @@ constexpr double eps = 1e-15, very_big = 1e24; constexpr int alt_symbolic_count_max = 1; constexpr double mem_increasing_factor = 1.1; -class dynSparseMatrix : public Evaluate +class dynSparseMatrix { public: - dynSparseMatrix(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg, bool print_error_arg); + dynSparseMatrix(Evaluate &evaluator_arg, int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool print_it_arg, bool steady_state_arg, bool block_decomposed_arg, int periods_arg, int minimal_solving_periods_arg, const BasicSymbolTable &symbol_table_arg, bool print_error_arg); void Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin, int y_kmax, int Size, int periods, bool cvg, int minimal_solving_periods, int stack_solve_algo, const vector_table_conditional_local_type &vector_table_conditional_local); void Simulate_Newton_One_Boundary(bool forward); void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1); @@ -147,6 +147,8 @@ protected: // Whether to use the block-decomposed version of the bytecode file bool block_decomposed; + Evaluate &evaluator; + stack Stack; int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u; int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y; diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 71e24f223..97156848a 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -716,7 +716,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) size_t nb_row_x = row_x; clock_t t0 = clock(); - Interpreter interprete {params, y, ya, x, steady_yd, direction, y_size, nb_row_x, + + const filesystem::path codfile {file_name + "/model/bytecode/" + (block_decomposed ? "block/" : "") + + (steady_state ? "static" : "dynamic") + ".cod"}; + Evaluate evaluator {codfile, steady_state, symbol_table}; + + Interpreter interprete {evaluator, params, y, ya, x, steady_yd, direction, y_size, nb_row_x, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms,