stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Bytecode: the class dynSparseMatrix is no longer derived from Evaluate

mr#2134
Sébastien Villemot 2023-04-18 17:53:27 +02:00
parent e22972849b
commit afe147d88d
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
7 changed files with 52 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
mex/sources/bytecode/Evaluate.cc
View File

@ -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())

23
mex/sources/bytecode/Evaluate.hh
View File

@ -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<string, it_code_type> print_expression(const it_code_type &expr_begin, const optional<it_code_type> &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<int, double> &TEF, map<pair<int, int>, double> &TEFD, map<tuple<int, int, int>, 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<string, it_code_type> print_expression(const it_code_type &expr_begin, const optional<it_code_type> &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<int, double> &TEF, map<pair<int, int>, double> &TEFD, map<tuple<int, int, int>, 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
{

64
mex/sources/bytecode/Interpreter.cc
View File

@ -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<s_plan> &sconstrained_extended_path)
{
vector<int> exogenous {getCurrentBlockExogenous()};
vector<int> endogenous {getCurrentBlockEndogenous()};
vector<int> exogenous {evaluator.getCurrentBlockExogenous()};
vector<int> 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<int> 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<double>(size);
else
residual = vector<double>(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<double>(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<bool, vector<int>>
Interpreter::extended_path(const string &file_name, bool evaluate, int block, int nb_periods, const vector<s_plan> &sextended_path, const vector<s_plan> &sconstrained_extended_path, const vector<string> &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<double *>(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<bool, vector<int>>
Interpreter::compute_blocks(const string &file_name, bool evaluate, int block)
{
ReadCodeFile(file_name);
//The big loop on intructions
vector<s_plan> 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)
{

3
mex/sources/bytecode/Interpreter.hh
View File

@ -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<bool, vector<int>> compute_blocks(const string &file_name, bool evaluate, int block);
void check_for_controlled_exo_validity(int current_block, const vector<s_plan> &sconstrained_extended_path);
pair<bool, vector<int>> MainLoop(const string &bin_basename, bool evaluate, int block, bool constrained, const vector<s_plan> &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

6
mex/sources/bytecode/SparseMatrix.cc
View File

@ -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)
{

6
mex/sources/bytecode/SparseMatrix.hh
View File

@ -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<double> 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;

7
mex/sources/bytecode/bytecode.cc
View File

@ -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,