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

Bytecode: move computation routines out of the Evaluate class

silicon
Sébastien Villemot 2023-02-24 15:12:07 -05:00
parent 17df930eb9
commit 325e6487d4
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
6 changed files with 233 additions and 232 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -2255,107 +2255,6 @@ Evaluate::compute_block_time(int Per_u_, bool evaluate, bool no_derivative)
#endif
}
void
Evaluate::evaluate_over_periods(bool forward)
{
if (steady_state)
compute_block_time(0, false, false);
else
{
if (forward)
{
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
compute_block_time(0, false, false);
it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
}
else
{
for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
compute_block_time(0, false, false);
it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
}
}
}
void
Evaluate::solve_simple_one_periods()
{
bool cvg = false;
int iter = 0;
double ya;
double slowc = 1;
res1 = 0;
while (!(cvg || iter > maxit_))
{
Per_y_ = it_*y_size;
ya = y[Block_Contain[0].Variable + Per_y_];
compute_block_time(0, false, false);
if (!isfinite(res1))
{
res1 = std::numeric_limits<double>::quiet_NaN();
while ((isinf(res1) || isnan(res1)) && (slowc > 1e-9))
{
compute_block_time(0, false, false);
if (!isfinite(res1))
{
slowc /= 1.5;
mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
y[Block_Contain[0].Variable + Per_y_] = ya - slowc * divide(r[0], g1[0]);
}
}
}
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
if (cvg)
continue;
try
{
y[Block_Contain[0].Variable + Per_y_] += -slowc *divide(rr, g1[0]);
}
catch (FloatingPointException &fpeh)
{
mexPrintf("%s\n \n", fpeh.message.c_str());
mexPrintf(" Singularity in block %d", block_num+1);
}
iter++;
}
if (!cvg)
throw FatalException{"In Solve Forward simple, convergence not achieved in block "
+ to_string(block_num+1) + ", after " + to_string(iter) + " iterations"};
}
void
Evaluate::solve_simple_over_periods(bool forward)
{
g1 = static_cast<double *>(mxMalloc(sizeof(double)));
test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double));
r = static_cast<double *>(mxMalloc(sizeof(double)));
test_mxMalloc(r, __LINE__, __FILE__, __func__, sizeof(double));
if (steady_state)
{
it_ = 0;
solve_simple_one_periods();
}
else
{
if (forward)
{
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
solve_simple_one_periods();
it_= periods+y_kmin-1; // Do not leave it_ in inconsistent state
}
else
{
for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
solve_simple_one_periods();
it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
}
}
mxFree(g1);
mxFree(r);
}
void
Evaluate::gotoBlock(int block)
{
@ -2373,122 +2272,6 @@ Evaluate::gotoBlock(int block)
u_count_int = fb->get_u_count_int();
}
void
Evaluate::compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx)
{
res1 = 0;
*_res1 = 0;
*_res2 = 0;
*_max_res = 0;
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
int shift = (it_-y_kmin) * size;
compute_block_time(Per_u_, false, no_derivatives);
if (!(isnan(res1) || isinf(res1)))
for (int i = 0; i < size; i++)
{
double rr;
rr = r[i];
res[i+shift] = rr;
if (max_res < fabs(rr))
{
*_max_res = fabs(rr);
*_max_res_idx = i;
}
*_res2 += rr*rr;
*_res1 += fabs(rr);
}
else
return;
}
it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
return;
}
bool
Evaluate::compute_complete(bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
{
bool result;
res1 = 0;
compute_block_time(0, false, no_derivatives);
if (!(isnan(res1) || isinf(res1)))
{
_res1 = 0;
_res2 = 0;
_max_res = 0;
for (int i = 0; i < size; i++)
{
double rr;
rr = r[i];
if (max_res < fabs(rr))
{
_max_res = fabs(rr);
_max_res_idx = i;
}
_res2 += rr*rr;
_res1 += fabs(rr);
}
result = true;
}
else
result = false;
return result;
}
bool
Evaluate::compute_complete(double lambda, double *crit)
{
double res1_ = 0, res2_ = 0, max_res_ = 0;
int max_res_idx_ = 0;
if (steady_state)
{
it_ = 0;
for (int i = 0; i < size; i++)
{
int eq = index_vara[i];
y[eq] = ya[eq] + lambda * direction[eq];
}
Per_u_ = 0;
Per_y_ = 0;
if (compute_complete(true, res1, res2, max_res, max_res_idx))
res2_ = res2;
else
return false;
}
else
{
for (int it = y_kmin; it < periods+y_kmin; it++)
for (int i = 0; i < size; i++)
{
int eq = index_vara[i];
y[eq+it*y_size] = ya[eq+it*y_size] + lambda * direction[eq+it*y_size];
}
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
if (compute_complete(true, res1, res2, max_res, max_res_idx))
{
res2_ += res2;
res1_ += res1;
if (max_res > max_res_)
{
max_res = max_res_;
max_res_idx = max_res_idx_;
}
}
else
return false;
}
it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
}
mexPrintf(" lambda=%e, res2=%e\n", lambda, res2_);
*crit = res2_/2;
return true;
}
void
Evaluate::printCurrentBlock()
{

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

@ -93,20 +93,10 @@ protected:
double divide(double a, double b);
double log1(double a);
double log10_1(double a);
void evaluate_over_periods(bool forward);
void solve_simple_one_periods();
void solve_simple_over_periods(bool forward);
void compute_block_time(int Per_u_, bool evaluate, bool no_derivatives);
int Per_u_, Per_y_;
int it_;
int maxit_;
double *direction;
double solve_tolf;
bool print_error;
double res1, res2, max_res;
int max_res_idx;
int *index_vara;
double res1;
int block_num; // Index of the current block
int size; // Size of the current block
@ -161,10 +151,6 @@ public:
Evaluate(int y_size_arg, int y_kmin_arg, int y_kmax_arg, bool steady_state_arg, int periods_arg, BasicSymbolTable &symbol_table_arg);
// TODO: integrate into the constructor
void loadCodeFile(const filesystem::path &codfile);
bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
void compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
bool compute_complete(double lambda, double *crit);
int
get_block_number() const

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

@ -65,6 +65,141 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
print_it = print_it_arg;
}
void
Interpreter::evaluate_over_periods(bool forward)
{
if (steady_state)
compute_block_time(0, false, false);
else
{
if (forward)
{
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
compute_block_time(0, false, false);
it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
}
else
{
for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
compute_block_time(0, false, false);
it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
}
}
}
void
Interpreter::solve_simple_one_periods()
{
bool cvg = false;
int iter = 0;
double ya;
double slowc = 1;
res1 = 0;
while (!(cvg || iter > maxit_))
{
Per_y_ = it_*y_size;
ya = y[Block_Contain[0].Variable + Per_y_];
compute_block_time(0, false, false);
if (!isfinite(res1))
{
res1 = std::numeric_limits<double>::quiet_NaN();
while ((isinf(res1) || isnan(res1)) && (slowc > 1e-9))
{
compute_block_time(0, false, false);
if (!isfinite(res1))
{
slowc /= 1.5;
mexPrintf("Reducing the path length in Newton step slowc=%f\n", slowc);
y[Block_Contain[0].Variable + Per_y_] = ya - slowc * divide(r[0], g1[0]);
}
}
}
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
if (cvg)
continue;
try
{
y[Block_Contain[0].Variable + Per_y_] += -slowc *divide(rr, g1[0]);
}
catch (FloatingPointException &fpeh)
{
mexPrintf("%s\n \n", fpeh.message.c_str());
mexPrintf(" Singularity in block %d", block_num+1);
}
iter++;
}
if (!cvg)
throw FatalException{"In Solve Forward simple, convergence not achieved in block "
+ to_string(block_num+1) + ", after " + to_string(iter) + " iterations"};
}
void
Interpreter::solve_simple_over_periods(bool forward)
{
g1 = static_cast<double *>(mxMalloc(sizeof(double)));
test_mxMalloc(g1, __LINE__, __FILE__, __func__, sizeof(double));
r = static_cast<double *>(mxMalloc(sizeof(double)));
test_mxMalloc(r, __LINE__, __FILE__, __func__, sizeof(double));
if (steady_state)
{
it_ = 0;
solve_simple_one_periods();
}
else
{
if (forward)
{
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
solve_simple_one_periods();
it_= periods+y_kmin-1; // Do not leave it_ in inconsistent state
}
else
{
for (it_ = periods+y_kmin-1; it_ >= y_kmin; it_--)
solve_simple_one_periods();
it_ = y_kmin; // Do not leave it_ in inconsistent state (see #1727)
}
}
mxFree(g1);
mxFree(r);
}
void
Interpreter::compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx)
{
res1 = 0;
*_res1 = 0;
*_res2 = 0;
*_max_res = 0;
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
int shift = (it_-y_kmin) * size;
compute_block_time(Per_u_, false, no_derivatives);
if (!(isnan(res1) || isinf(res1)))
for (int i = 0; i < size; i++)
{
double rr;
rr = r[i];
res[i+shift] = rr;
if (max_res < fabs(rr))
{
*_max_res = fabs(rr);
*_max_res_idx = i;
}
*_res2 += rr*rr;
*_res1 += fabs(rr);
}
else
return;
}
it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
return;
}
void
Interpreter::evaluate_a_block(bool initialization, bool single_block, const string &bin_base_name)
{

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

@ -41,6 +41,10 @@ private:
bool global_temporary_terms;
bool print;
int col_x, col_y;
void evaluate_over_periods(bool forward);
void solve_simple_one_periods();
void solve_simple_over_periods(bool forward);
void compute_complete_2b(bool no_derivatives, double *_res1, double *_res2, double *_max_res, int *_max_res_idx);
protected:
void evaluate_a_block(bool initialization, bool single_block, const string &bin_base_name);
int simulate_a_block(const vector_table_conditional_local_type &vector_table_conditional_local, bool single_block, const string &bin_base_name);

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

@ -1824,6 +1824,88 @@ dynSparseMatrix::Sparse_transpose(const mxArray *A_m)
return C_m;
}
bool
dynSparseMatrix::compute_complete(bool no_derivatives, double &_res1, double &_res2, double &_max_res, int &_max_res_idx)
{
bool result;
res1 = 0;
compute_block_time(0, false, no_derivatives);
if (!(isnan(res1) || isinf(res1)))
{
_res1 = 0;
_res2 = 0;
_max_res = 0;
for (int i = 0; i < size; i++)
{
double rr;
rr = r[i];
if (max_res < fabs(rr))
{
_max_res = fabs(rr);
_max_res_idx = i;
}
_res2 += rr*rr;
_res1 += fabs(rr);
}
result = true;
}
else
result = false;
return result;
}
bool
dynSparseMatrix::compute_complete(double lambda, double *crit)
{
double res1_ = 0, res2_ = 0, max_res_ = 0;
int max_res_idx_ = 0;
if (steady_state)
{
it_ = 0;
for (int i = 0; i < size; i++)
{
int eq = index_vara[i];
y[eq] = ya[eq] + lambda * direction[eq];
}
Per_u_ = 0;
Per_y_ = 0;
if (compute_complete(true, res1, res2, max_res, max_res_idx))
res2_ = res2;
else
return false;
}
else
{
for (int it = y_kmin; it < periods+y_kmin; it++)
for (int i = 0; i < size; i++)
{
int eq = index_vara[i];
y[eq+it*y_size] = ya[eq+it*y_size] + lambda * direction[eq+it*y_size];
}
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
if (compute_complete(true, res1, res2, max_res, max_res_idx))
{
res2_ += res2;
res1_ += res1;
if (max_res > max_res_)
{
max_res = max_res_;
max_res_idx = max_res_idx_;
}
}
else
return false;
}
it_ = periods+y_kmin-1; // Do not leave it_ in inconsistent state
}
mexPrintf(" lambda=%e, res2=%e\n", lambda, res2_);
*crit = res2_/2;
return true;
}
bool
dynSparseMatrix::mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc)
{

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

@ -195,6 +195,17 @@ protected:
int stack_solve_algo, solve_algo;
int minimal_solving_periods;
bool print_it;
int Per_u_, Per_y_;
int maxit_;
double *direction;
double solve_tolf;
double res2, max_res;
int max_res_idx;
int *index_vara;
bool compute_complete(bool no_derivatives, double &res1, double &res2, double &max_res, int &max_res_idx);
bool compute_complete(double lambda, double *crit);
};
#endif // _SPARSEMATRIX_HH