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

- compute the jacobian matrix for exogenous, det_exogenous and previous blocks endogenous 

- preprocessor provides informations about the dynamic structure of each block
- extends the algorithms available with bytecode: stack_solve_algo = 1, 2, 3 and 4 is compatible the bytecode. Speed and memory requirement are improved with stack_solve_algo = 1 or 4 for large scale models.
- bytecode can be used to evaluate a model
time-shift
Ferhat Mihoubi 2010-07-23 11:20:24 +02:00 committed by Sébastien Villemot
parent f268513ffb
commit 0a3c8a4b0c
21 changed files with 2813 additions and 1164 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/model_info.m
View File

@ -34,7 +34,7 @@ if(isfield(M_,'block_structure'))
end;
for j=1:size_block
if(j==1)
fprintf('| %3d (%4d) | %10d | %30s | %14d | %-6d %24s |\n',i,M_.block_structure.block(i).num,size_block,Sym_type(M_.block_structure.block(i).Simulation_Type),M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
fprintf('| %10d | %10d | %30s | %14d | %-6d %24s |\n',i,size_block,Sym_type(M_.block_structure.block(i).Simulation_Type),M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
else
fprintf('| %10s | %10s | %30s | %14d | %-6d %24s |\n','','','',M_.block_structure.block(i).equation(j),M_.block_structure.block(i).variable(j),M_.endo_names(M_.block_structure.block(i).variable(j),:));
end;

2
matlab/read_data_.m
View File

@ -43,7 +43,7 @@ if size(oo_.endo_simul,2) < M_.maximum_lag+M_.maximum_lead+options_.periods
positions = [positions ; strmatch(chopped,M_.endo_names,'exact')];
end
Values=fscanf(fid,'%f',inf);
Values=reshape(Values,M_.endo_nbr,size(Values,1)/M_.endo_nbr);
Values=reshape(Values,M_.orig_endo_nbr,size(Values,1)/M_.orig_endo_nbr);
oo_.endo_simul=Values(positions,:);
fclose(fid);
end

4
matlab/simul.m
View File

@ -63,8 +63,8 @@ end
if options_.block && ~options_.bytecode && (options_.stack_solve_algo == 0 || options_.stack_solve_algo == 5)
error('SIMUL: for the moment, you must use stack_solve_algo={1,2,3,4} when using block without bytecode option')
end
if options_.bytecode && options_.stack_solve_algo ~= 5
error('SIMUL: for the moment, you must use stack_solve_algo=5 with bytecode option')
if options_.bytecode && (options_.stack_solve_algo ~= 1 && options_.stack_solve_algo ~= 2 && options_.stack_solve_algo ~= 3 && options_.stack_solve_algo ~= 4 && options_.stack_solve_algo ~= 5)
error('SIMUL: for the moment, you must use stack_solve_algo= 1, 2, 3, 4 or 5 with bytecode option')
end
if exist('OCTAVE_VERSION') && options_.stack_solve_algo == 2

7
matlab/steady_.m
View File

@ -29,7 +29,8 @@ function steady_()
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ oo_ it_ options_
if options_.bytecode && options_.solve_algo ~= 5
if options_.bytecode && ...
(options_.solve_algo ~= 1 && options_.solve_algo ~= 2 && options_.solve_algo ~= 3 && options_.solve_algo ~= 4 && options_.solve_algo ~= 5)
error('STEADY: for the moment, you must use solve_algo=5 with bytecode option')
end
if ~options_.bytecode && options_.solve_algo == 5
@ -75,7 +76,9 @@ if options_.steadystate_flag
check1 = 1;
end
elseif options_.block && options_.bytecode
[residuals, check1] = bytecode('evaluate','static');
[residuals, check1] = bytecode('evaluate','static',oo_.steady_state,...
[oo_.exo_steady_state; ...
oo_.exo_det_steady_state], M_.params, 1);
else
check1 = 0;
check1 = max(abs(feval([M_.fname '_static'],...

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

@ -38,7 +38,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
double *direction_arg, int y_size_arg,
int nb_row_x_arg, int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg,
int maxit_arg_, double solve_tolf_arg, int size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
string &filename_arg, int minimal_solving_periods_arg)
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg)
{
params = params_arg;
y = y_arg;
@ -64,6 +64,8 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
T = NULL;
error_not_printed = true;
minimal_solving_periods = minimal_solving_periods_arg;
stack_solve_algo = stack_solve_algo_arg;
solve_algo = solve_algo_arg;
}
string
@ -105,7 +107,7 @@ Interpreter::add_underscore_to_fpe(string str)
string
Interpreter::get_variable(SymbolType variable_type, int variable_num)
Interpreter::get_variable(SymbolType variable_type, unsigned int variable_num)
{
ostringstream res;
#ifndef DEBUG_EX
@ -364,6 +366,9 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2();
dvar3 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2();
break;
default:
mexPrintf("Dérivatives %d not implemented yet\n", it_code->first);
mexErrMsgTxt("end of bytecode\n");
}
break;
case FLDV:
@ -1098,15 +1103,31 @@ Interpreter::print_expression(it_code_type it_code, bool evaluate)
}
void
Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state)
{
int var, lag = 0, op;
int var = 0, lag = 0, op;
unsigned int eq, pos_col;
ostringstream tmp_out;
double v1, v2, v3;
bool go_on = true;
double ll;
double rr;
double *jacob = NULL, *jacob_other_endo = NULL, *jacob_exo = NULL, *jacob_exo_det = NULL;
EQN_block = block_num;
#ifdef DEBUG
mexPrintf("compute_block_time\n");
#endif
#ifndef DEBUG_EX
if (evaluate && !steady_state)
{
jacob = mxGetPr(jacobian_block[block_num]);
mexPrintf("jacobian_block[%d]=%x\n",block_num, jacobian_block[block_num]);
jacob_other_endo = mxGetPr(jacobian_other_endo_block[block_num]);
jacob_exo = mxGetPr(jacobian_exo_block[block_num]);
jacob_exo_det = mxGetPr(jacobian_det_exo_block[block_num]);
}
#endif
//feclearexcept (FE_ALL_EXCEPT);
while (go_on)
{
@ -1114,41 +1135,74 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
switch (it_code->first)
{
case FNUMEXPR:
#ifdef DEBUG
mexPrintf("FNUMEXPR\n");
#endif
it_code_expr = it_code;
switch (((FNUMEXPR_ *) it_code->second)->get_expression_type())
{
case TemporaryTerm:
#ifdef DEBUG
mexPrintf("TemporaryTerm\n");
#endif
EQN_type = TemporaryTerm;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
break;
case ModelEquation:
#ifdef DEBUG
mexPrintf("ModelEquation\n");
#endif
EQN_type = ModelEquation;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
break;
case FirstEndoDerivative:
#ifdef DEBUG
mexPrintf("FirstEndoDerivative\n");
#endif
EQN_type = FirstEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
case FirstOtherEndoDerivative:
#ifdef DEBUG
mexPrintf("FirstOtherEndoDerivative\n");
#endif
EQN_type = FirstOtherEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
case FirstExoDerivative:
#ifdef DEBUG
mexPrintf("FirstExoDerivative\n");
#endif
EQN_type = FirstExoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
case FirstExodetDerivative:
#ifdef DEBUG
mexPrintf("FirstExodetDerivative\n");
#endif
EQN_type = FirstExodetDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_lag1 = ((FNUMEXPR_ *) it_code->second)->get_lag1();
break;
case FirstParamDerivative:
#ifdef DEBUG
mexPrintf("FirstParamDerivative\n");
#endif
EQN_type = FirstParamDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
break;
case SecondEndoDerivative:
#ifdef DEBUG
mexPrintf("SecondEndoDerivative\n");
#endif
EQN_type = SecondEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1157,6 +1211,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2();
break;
case SecondExoDerivative:
#ifdef DEBUG
mexPrintf("SecondExoDerivative\n");
#endif
EQN_type = SecondExoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1165,6 +1222,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2();
break;
case SecondExodetDerivative:
#ifdef DEBUG
mexPrintf("SecondExodetDerivative\n");
#endif
EQN_type = SecondExodetDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1173,12 +1233,18 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag2 = ((FNUMEXPR_ *) it_code->second)->get_lag2();
break;
case SecondParamDerivative:
#ifdef DEBUG
mexPrintf("SecondParamDerivative\n");
#endif
EQN_type = SecondParamDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
EQN_dvar2 = ((FNUMEXPR_ *) it_code->second)->get_dvariable2();
break;
case ThirdEndoDerivative:
#ifdef DEBUG
mexPrintf("ThirdEndoDerivative\n");
#endif
EQN_type = ThirdEndoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1189,6 +1255,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3();
break;
case ThirdExoDerivative:
#ifdef DEBUG
mexPrintf("ThirdExoDerivative\n");
#endif
EQN_type = ThirdExoDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1199,6 +1268,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3();
break;
case ThirdExodetDerivative:
#ifdef DEBUG
mexPrintf("ThirdExodetDerivative\n");
#endif
EQN_type = ThirdExodetDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1209,6 +1281,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
EQN_lag3 = ((FNUMEXPR_ *) it_code->second)->get_lag3();
break;
case ThirdParamDerivative:
#ifdef DEBUG
mexPrintf("ThirdParamDerivative\n");
#endif
EQN_type = ThirdParamDerivative;
EQN_equation = ((FNUMEXPR_ *) it_code->second)->get_equation();
EQN_dvar1 = ((FNUMEXPR_ *) it_code->second)->get_dvariable1();
@ -1223,14 +1298,18 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
case eParameter:
var = ((FLDV_ *) it_code->second)->get_pos();
Stack.push(params[var]);
#ifdef DEBUG
mexPrintf("FLDV Param[var=%d]\n",var);
tmp_out << " params[" << var << "](" << params[var] << ")";
#endif
Stack.push(params[var]);
break;
case eEndogenous:
var = ((FLDV_ *) it_code->second)->get_pos();
lag = ((FLDV_ *) it_code->second)->get_lead_lag();
#ifdef DEBUG
mexPrintf("FLDV y[var=%d, lag=%d, it_=%d], y_size=%d evaluate=%d\n",var, lag, it_, y_size, evaluate);
#endif
if (evaluate)
Stack.push(ya[(it_+lag)*y_size+var]);
else
@ -1242,10 +1321,11 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
case eExogenous:
var = ((FLDV_ *) it_code->second)->get_pos();
lag = ((FLDV_ *) it_code->second)->get_lead_lag();
Stack.push(x[it_+lag+var*nb_row_x]);
#ifdef DEBUG
mexPrintf("FLDV x[var=%d, lag=%d, it_=%d], nb_row_x=%d evaluate=%d\n",var, lag, it_, nb_row_x, evaluate);
tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")";
#endif
Stack.push(x[it_+lag+var*nb_row_x]);
break;
case eExogenousDet:
var = ((FLDV_ *) it_code->second)->get_pos();
@ -1268,18 +1348,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
case eParameter:
var = ((FLDSV_ *) it_code->second)->get_pos();
Stack.push(params[var]);
#ifdef DEBUG
mexPrintf("FLDSV Param[var=%d]\n",var);
tmp_out << " params[" << var << "](" << params[var] << ")";
#endif
Stack.push(params[var]);
break;
case eEndogenous:
var = ((FLDSV_ *) it_code->second)->get_pos();
if (evaluate)
Stack.push(ya[var]);
else
Stack.push(y[var]);
#ifdef DEBUG
mexPrintf("FLDSV y[var=%d]\n",var);
tmp_out << " y[" << var << "](" << y[var] << ")";
if(var<0 || var>= y_size)
{
@ -1287,16 +1365,24 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
mexErrMsgTxt("End of bytecode");
}
#endif
if (evaluate)
Stack.push(ya[var]);
else
Stack.push(y[var]);
break;
case eExogenous:
var = ((FLDSV_ *) it_code->second)->get_pos();
Stack.push(x[var]);
#ifdef DEBUG
mexPrintf("FLDSV x[var=%d]\n",var);
tmp_out << " x[" << var << "](" << x[var] << ")";
#endif
Stack.push(x[var]);
break;
case eExogenousDet:
var = ((FLDSV_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("FLDSV xd[var=%d]\n",var);
#endif
Stack.push(x[var]);
break;
case eModelLocalVariable:
@ -1316,28 +1402,28 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
case eParameter:
var = ((FLDVS_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("params[%d]=%f\n", var, params[var]);
mexPrintf("params[%d]\n", var);
#endif
Stack.push(params[var]);
break;
case eEndogenous:
var = ((FLDVS_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf(" steady_y[%d]=%f\n", var, steady_y[var]);
mexPrintf("FLDVS steady_y[%d]\n", var);
#endif
Stack.push(steady_y[var]);
break;
case eExogenous:
var = ((FLDVS_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf(" x[%d] = %f\n", var, x[var]);
mexPrintf("FLDVS x[%d] \n", var);
#endif
Stack.push(x[var]);
break;
case eExogenousDet:
var = ((FLDVS_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf(" xd[%d] = %f\n", var, x[var]);
mexPrintf("FLDVS xd[%d]\n", var);
#endif
Stack.push(x[var]);
break;
@ -1364,6 +1450,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
//load a temporary variable in the processor
var = ((FLDST_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("FLDST T[%d]\n",var);
tmp_out << " T[" << var << "](" << T[var] << ")";
#endif
Stack.push(T[var]);
@ -1373,6 +1460,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
var = ((FLDU_ *) it_code->second)->get_pos();
var += Per_u_;
#ifdef DEBUG
mexPrintf("FLDU u[%d]\n",var);
tmp_out << " u[" << var << "](" << u[var] << ")";
#endif
Stack.push(u[var]);
@ -1381,6 +1469,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
//load u variable in the processor
var = ((FLDSU_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("FLDSU u[%d]\n",var);
tmp_out << " u[" << var << "](" << u[var] << ")";
#endif
Stack.push(u[var]);
@ -1388,10 +1477,16 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
case FLDR:
//load u variable in the processor
var = ((FLDR_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("FLDR r[%d]\n",var);
#endif
Stack.push(r[var]);
break;
case FLDZ:
//load 0 in the processor
#ifdef DEBUG
mexPrintf("FLDZ\n");
#endif
Stack.push(0.0);
#ifdef DEBUG
tmp_out << " 0";
@ -1401,6 +1496,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
//load a numerical constant in the processor
ll = ((FLDC_ *) it_code->second)->get_value();
#ifdef DEBUG
mexPrintf("FLDC = %f\n",ll);
tmp_out << " " << ll;
#endif
@ -1412,6 +1508,9 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
{
case eParameter:
var = ((FSTPV_ *) it_code->second)->get_pos();
#ifdef DEBUG
mexPrintf("FSTPV params[%d]\n",var);
#endif
params[var] = Stack.top();
Stack.pop();
break;
@ -1559,6 +1658,63 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
#endif
Stack.pop();
break;
case FSTPG3:
//store in derivative (g) variable from the processor
#ifdef DEBUG
mexPrintf("FSTPG3\n");
mexEvalString("drawnow;");
#endif
rr = Stack.top();
switch(EQN_type)
{
case FirstEndoDerivative:
eq = ((FSTPG3_ *) it_code->second)->get_row();
var = ((FSTPG3_ *) it_code->second)->get_col();
lag = ((FSTPG3_ *) it_code->second)->get_lag();
pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
mexPrintf("jacob[%d(size=%d*pos_col=%d + eq=%d )]=%f\n",eq + size*pos_col, size, pos_col, eq, rr);
jacob[eq + size*pos_col] = rr;
break;
case FirstOtherEndoDerivative:
//eq = ((FSTPG3_ *) it_code->second)->get_row();
eq = EQN_equation;
var = ((FSTPG3_ *) it_code->second)->get_col();
lag = ((FSTPG3_ *) it_code->second)->get_lag();
pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
mexPrintf("jacob_other_endo[%d(size=%d*pos_col=%d + eq=%d)]=%f\n",size*pos_col + eq, size, pos_col, eq, rr);
jacob_other_endo[eq + size*pos_col] = rr;
break;
case FirstExoDerivative:
//eq = ((FSTPG3_ *) it_code->second)->get_row();
eq = EQN_equation;
var = ((FSTPG3_ *) it_code->second)->get_col();
lag = ((FSTPG3_ *) it_code->second)->get_lag();
pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
mexPrintf("jacob_exo[%d(size=%d*pos_col=%dr + eq=%d)]=%f\n",size*pos_col+eq, size, pos_col, eq, rr);
jacob_exo[eq + size*pos_col] = rr;
break;
case FirstExodetDerivative:
//eq = ((FSTPG3_ *) it_code->second)->get_row();
eq = EQN_equation;
var = ((FSTPG3_ *) it_code->second)->get_col();
lag = ((FSTPG3_ *) it_code->second)->get_lag();
pos_col = ((FSTPG3_ *) it_code->second)->get_col_pos();
mexPrintf("jacob_exo_det[%d(size=%d*pos_col=%dr + eq=%d)]=%f\n",size*pos_col+eq, size, pos_col, eq, rr);
jacob_exo_det[eq + size*pos_col] = rr;
break;
default:
mexPrintf("Variable %d not used yet\n", EQN_type);
mexErrMsgTxt("end of bytecode\n");
}
#ifdef DEBUG
tmp_out << "=>";
mexPrintf(" g1[%d](%f)=%s\n", var, g1[var], tmp_out.str().c_str());
tmp_out.str("");
#endif
Stack.pop();
break;
case FBINARY:
op = ((FBINARY_ *) it_code->second)->get_op_type();
v2 = Stack.top();
@ -1691,7 +1847,6 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
Stack.push(log1(v1, evaluate));
if (isnan(res1))
go_on = false;
#ifdef DEBUG
tmp_out << " |log(" << v1 << ")|";
#endif
@ -1808,10 +1963,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
switch (op)
{
case oNormcdf:
#ifndef _MSC_VER
//mexPrintf("normcdf(v1=%f, v2=%f, v3=%f)=%f\n", v1, v2, v3, 0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
Stack.push(0.5*(1+erf((v1-v2)/v3/M_SQRT2)));
# ifdef DEBUG
#ifndef _MSC_VER Stack.push(0.5*(1+erf((v1-v2)/v3/M_SQRT2)));# ifdef DEBUG
tmp_out << " |normcdf(" << v1 << ", " << v2 << ", " << v3 << ")|";
# endif
#else
@ -1837,10 +1989,30 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
break;
case FENDBLOCK:
//it's the block end
#ifdef DEBUG
mexPrintf("FENDBLOCK\n");
#endif
go_on = false;
break;
case FENDEQU:
break;
case FJMPIFEVAL:
if (evaluate)
{
#ifdef DEBUG
mexPrintf("FJMPIFEVAL length=%d\n",((FJMPIFEVAL_ *) it_code->second)->get_pos());
mexEvalString("drawnow;");
#endif
it_code += ((FJMPIFEVAL_ *) it_code->second)->get_pos()/* - 1*/;
}
break;
case FJMP:
#ifdef DEBUG
mexPrintf("FJMP length=%d\n",((FJMP_ *) it_code->second)->get_pos());
mexEvalString("drawnow;");
#endif
it_code += ((FJMP_ *) it_code->second)->get_pos() /*- 1 */;
break;
case FOK:
op = ((FOK_ *) it_code->second)->get_arg();
if (Stack.size() > 0)
@ -1856,18 +2028,23 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num)
}
it_code++;
}
#ifdef DEBUG
mexPrintf("==> end of compute_block_time Block = %d\n",block_num);
mexEvalString("drawnow;");
#endif
}
void
Interpreter::evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag, const int Block_List_Max_Lead, const int u_count_int)
const bool is_linear, const int symbol_table_endo_nbr, const int Block_List_Max_Lag,
const int Block_List_Max_Lead, const int u_count_int)
{
it_code_type begining;
switch (type)
{
case EVALUATE_FORWARD:
if (steady_state)
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
else
{
begining = it_code;
@ -1875,7 +2052,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
}
}
break;
@ -1884,7 +2061,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@ -1893,10 +2070,9 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
begining = it_code;
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
it_code = begining;
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@ -1906,11 +2082,11 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
break;
case SOLVE_FORWARD_COMPLETE:
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@ -1921,7 +2097,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@ -1930,7 +2106,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
break;
case EVALUATE_BACKWARD:
if (steady_state)
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
else
{
begining = it_code;
@ -1938,7 +2114,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
}
}
break;
@ -1947,7 +2123,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@ -1958,7 +2134,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@ -1968,11 +2144,11 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
break;
case SOLVE_BACKWARD_COMPLETE:
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
r = (double *) mxMalloc(size*sizeof(double));
if (steady_state)
{
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[Block_Contain[j].Variable] += r[j];
}
@ -1983,7 +2159,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, true, block_num);
compute_block_time(0, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@ -1993,7 +2169,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true);
Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true, stack_solve_algo, solve_algo);
u_count = u_count_int*(periods+y_kmax+y_kmin);
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@ -2002,7 +2178,7 @@ Interpreter::evaluate_a_block(const int size, const int type, string bin_basenam
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
it_code = begining;
compute_block_time(Per_u_, true, block_num);
compute_block_time(Per_u_, true, block_num, size, steady_state);
for (int j = 0; j < size; j++)
y[it_*y_size+Block_Contain[j].Variable] += r[j];
}
@ -2021,11 +2197,18 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
int giter;
bool result = true;
double *y_save;
res1 = 0;
#ifdef DEBUG
mexPrintf("simulate_a_block\n");
#endif
switch (type)
{
case EVALUATE_FORWARD:
#ifdef DEBUG
mexPrintf("EVALUATE_FORWARD\n");
#endif
if (steady_state)
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
else
{
begining = it_code;
@ -2033,13 +2216,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
}
}
break;
case EVALUATE_BACKWARD:
#ifdef DEBUG
mexPrintf("EVALUATE_BACKWARD\n");
#endif
if (steady_state)
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
else
{
begining = it_code;
@ -2047,11 +2233,14 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
}
}
break;
case SOLVE_FORWARD_SIMPLE:
#ifdef DEBUG
mexPrintf("SOLVE_FORWARD_SIMPLE size=%d\n",size);
#endif
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@ -2063,7 +2252,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
@ -2099,7 +2288,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
double rr;
if (fabs(1+y[Per_y_+Block_Contain[0].Variable]) > eps)
rr = r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
@ -2131,6 +2320,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
mxFree(r);
break;
case SOLVE_BACKWARD_SIMPLE:
#ifdef DEBUG
mexPrintf("SOLVE_BACKWARD_SIMPLE\n");
#endif
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@ -2142,7 +2334,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
double rr;
rr = r[0];
cvg = (fabs(rr) < solve_tolf);
@ -2177,7 +2369,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
{
it_code = begining;
Per_y_ = it_*y_size;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
double rr;
if (fabs(1+y[Per_y_+Block_Contain[0].Variable]) > eps)
rr = r[0]/(1+y[Per_y_+Block_Contain[0].Variable]);
@ -2209,8 +2401,11 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
mxFree(r);
break;
case SOLVE_FORWARD_COMPLETE:
#ifdef DEBUG
mexPrintf("SOLVE_FORWARD_COMPLETE\n");
#endif
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@ -2231,7 +2426,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@ -2253,7 +2448,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2273,11 +2468,14 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
iter = 0;
res1 = res2 = max_res = 0; max_res_idx = 0;
res1 = 0;
res2 = 0;
max_res = 0;
max_res_idx = 0;
error_not_printed = true;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
if (!result)
{
mexPrintf("Convergence not achieved in block %d\n", Block_Count);
@ -2304,7 +2502,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@ -2329,7 +2527,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2352,11 +2550,13 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
iter = 0;
res1 = res2 = max_res = 0; max_res_idx = 0;
res1 = 0;
res2 = 0;
max_res = 0; max_res_idx = 0;
error_not_printed = true;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
}
}
}
@ -2368,8 +2568,11 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
mxFree(u);
break;
case SOLVE_BACKWARD_COMPLETE:
#ifdef DEBUG
mexPrintf("SOLVE_BACKWARD_COMPLETE\n");
#endif
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false);
Read_SparseMatrix(bin_basename, size, 1, 0, 0, steady_state, false, stack_solve_algo, solve_algo);
g1 = (double *) mxMalloc(size*size*sizeof(double));
r = (double *) mxMalloc(size*sizeof(double));
begining = it_code;
@ -2389,7 +2592,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@ -2411,7 +2614,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2431,11 +2634,13 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
iter = 0;
res1 = res2 = max_res = 0; max_res_idx = 0;
res1 = 0;
res2 = 0;
max_res = 0; max_res_idx = 0;
error_not_printed = true;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, 0, 0, 0, size, false, cvg, iter, true, EQN_block_number, solve_algo);
if (!result)
{
mexPrintf("Convergence not achieved in block %d\n", Block_Count);
@ -2462,7 +2667,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
res2 = 0;
res1 = 0;
max_res = 0;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
if (!(isnan(res1) || isinf(res1)))
{
for (i = 0; i < size; i++)
@ -2487,7 +2692,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
if (cvg)
continue;
int prev_iter = iter;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2510,9 +2715,9 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
it_code = begining;
Per_y_ = it_*y_size;
error_not_printed = true;
compute_block_time(0, false, block_num);
compute_block_time(0, false, block_num, size, steady_state);
cvg = false;
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number);
result = simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter, false, EQN_block_number, solve_algo);
}
}
}
@ -2525,13 +2730,16 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
break;
case SOLVE_TWO_BOUNDARIES_SIMPLE:
case SOLVE_TWO_BOUNDARIES_COMPLETE:
#ifdef DEBUG
mexPrintf("SOLVE_TWO_BOUNDARIES\n");
#endif
if (steady_state)
{
mexPrintf("SOLVE_TWO_BOUNDARIES in a steady state model: impossible case\n");
return false;
}
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true);
Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax, steady_state, true, stack_solve_algo, solve_algo);
u_count = u_count_int*(periods+y_kmax+y_kmin);
r = (double *) mxMalloc(size*sizeof(double));
y_save = (double *) mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@ -2560,7 +2768,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
Per_y_ = it_*y_size;
it_code = begining;
error_not_printed = true;
compute_block_time(Per_u_, false, block_num);
compute_block_time(Per_u_, false, block_num, size, steady_state);
if (isnan(res1) || isinf(res1))
{
memcpy(y, y_save, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
@ -2588,7 +2796,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
cvg = (max_res < solve_tolf);
u_count = u_count_saved;
int prev_iter = iter;
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number);
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number, stack_solve_algo);
iter++;
if (iter > prev_iter)
{
@ -2605,13 +2813,15 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
else
{
res1 = res2 = max_res = 0; max_res_idx = 0;
res1 = 0;
res2 = 0;
max_res = 0; max_res_idx = 0;
for (it_ = y_kmin; it_ < periods+y_kmin; it_++)
{
Per_u_ = (it_-y_kmin)*u_count_int;
Per_y_ = it_*y_size;
it_code = begining;
compute_block_time(Per_u_, false, block_num);
compute_block_time(Per_u_, false, block_num, size, steady_state);
for (i = 0; i < size; i++)
{
double rr;
@ -2626,7 +2836,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
}
cvg = false;
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number);
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter, minimal_solving_periods, EQN_block_number, stack_solve_algo);
}
mxFree(r);
mxFree(y_save);
@ -2644,7 +2854,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name,
}
bool
Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate)
Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks)
{
bool result = true;
@ -2686,8 +2896,23 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
Block_Contain = fb->get_Block_Contain();
it_code++;
if (evaluate)
{
#ifndef DEBUG_EX
if (!steady_state)
{
jacobian_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_jacob(),mxREAL));
mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_nb_col_jacob(), sizeof(jacobian_block[Block_Count]));
jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_exo_size(),mxREAL));
mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_exo_block.size()=%d fb->get_exo_size()=%d\n",Block_Count, fb->get_size(), fb->get_exo_size(), sizeof(jacobian_exo_block[Block_Count]), fb->get_exo_size());
jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_det_exo_size(),mxREAL));
mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_det_exo_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_det_exo_size(), sizeof(jacobian_det_exo_block[Block_Count]));
jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(),mxREAL));
mexPrintf("at block = %d, mxCreateDoubleMatrix(%d, %d, mxREAL) jacobian_other_endo_block.size()=%d\n",Block_Count, fb->get_size(), fb->get_nb_col_other_endo_jacob(), sizeof(jacobian_other_endo_block[Block_Count]));
}
#endif
evaluate_a_block(fb->get_size(), fb->get_type(), bin_basename, steady_state, Block_Count,
fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int());
}
else
result = simulate_a_block(fb->get_size(), fb->get_type(), file_name, bin_basename, true, steady_state, Block_Count,
fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int());
@ -2724,13 +2949,14 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool steady_s
it_code++;
break;
default:
mexPrintf("Unknown command \n");
mexPrintf("Unknown command %d\n",it_code->first);
mexEvalString("drawnow;");
mexErrMsgTxt("End of bytecode");
break;
}
}
mxFree(Init_Code->second);
nb_blocks = Block_Count+1;
if (T)
mxFree(T);
return result;

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

@ -41,12 +41,15 @@
using namespace std;
#define pow_ pow
typedef vector<pair<Tags, void * > >::const_iterator it_code_type;
typedef vector<pair<Tags, void * > > code_liste_type;
typedef code_liste_type::const_iterator it_code_type;
class Interpreter : SparseMatrix
{
private:
#ifndef DEBUG_EX
vector<mxArray*> jacobian_block, jacobian_other_endo_block, jacobian_exo_block, jacobian_det_exo_block;
#endif
ExpressionType EQN_type;
unsigned int EQN_equation, EQN_block, EQN_block_number;
unsigned int EQN_dvar1, EQN_dvar2, EQN_dvar3;
@ -58,9 +61,9 @@ private:
double log10_1(double a, bool evaluate);
string remove_white(string str);
string add_underscore_to_fpe(string str);
string get_variable(SymbolType variable_type, int variable_num);
string get_variable(SymbolType variable_type, unsigned int variable_num);
string error_location(bool evaluate);
void compute_block_time(int Per_u_, bool evaluate, int block_num);
void compute_block_time(int Per_u_, bool evaluate, int block_num, int size, bool steady_state);
string print_expression(it_code_type it_code, bool evaluate);
void evaluate_a_block(const int size, const int type, string bin_basename, bool steady_state, int block_num,
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
@ -68,7 +71,7 @@ private:
const bool is_linear = false, const int symbol_table_endo_nbr = 0, const int Block_List_Max_Lag = 0, const int Block_List_Max_Lead = 0, const int u_count_int = 0);
double *T;
vector<Block_contain_type> Block_Contain;
vector<pair<Tags, void * > > code_liste;
code_liste_type code_liste;
it_code_type it_code;
stack<double> Stack;
int Block_Count, Per_u_, Per_y_;
@ -82,13 +85,20 @@ private:
int equation, derivative_equation, derivative_variable;
string filename;
int minimal_solving_periods;
int stack_solve_algo, solve_algo;
public:
~Interpreter();
Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg,
double *direction_arg, int y_size_arg, int nb_row_x_arg,
int nb_row_xd_arg, int periods_arg, int y_kmin_arg, int y_kmax_arg, int maxit_arg_, double solve_tolf_arg, int size_o_direction_arg,
double slowc_arg, int y_decal_arg, double markowitz_c_arg, string &filename_arg, int minimal_solving_periods_arg);
bool compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate);
double slowc_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 compute_blocks(string file_name, string bin_basename, bool steady_state, bool evaluate, bool block, int &nb_blocks);
#ifndef DEBUG_EX
inline mxArray* get_jacob(int block_num) {return jacobian_block[block_num];};
inline mxArray* get_jacob_exo(int block_num) {return jacobian_exo_block[block_num];};
inline mxArray* get_jacob_exo_det(int block_num) {return jacobian_det_exo_block[block_num];};
inline mxArray* get_jacob_other_endo(int block_num) {return jacobian_other_endo_block[block_num];};
#endif
};
#endif

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

@ -27,7 +27,7 @@ Mem_Mngr::Mem_Mngr()
void
Mem_Mngr::Print_heap()
{
int i;
unsigned int i;
mexPrintf("i :");
for (i = 0; i < CHUNK_SIZE; i++)
mexPrintf("%3d ", i);
@ -61,7 +61,7 @@ Mem_Mngr::init_CHUNK_BLCK_SIZE(int u_count)
NonZeroElem *
Mem_Mngr::mxMalloc_NZE()
{
long int i;
long unsigned int i;
if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/
{
NonZeroElem *p1 = Chunk_Stack.back();
@ -102,7 +102,7 @@ Mem_Mngr::mxMalloc_NZE()
void
Mem_Mngr::mxFree_NZE(void *pos)
{
int i;
unsigned int i;
size_t gap;
for (i = 0; i < Nb_CHUNK; i++)
{

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

@ -52,8 +52,8 @@ public:
bool swp_f;
private:
v_NonZeroElem Chunk_Stack;
int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK;
int CHUNK_heap_pos;
unsigned int CHUNK_SIZE, CHUNK_BLCK_SIZE, Nb_CHUNK;
unsigned int CHUNK_heap_pos;
NonZeroElem **NZE_Mem_add;
NonZeroElem *NZE_Mem;
vector<NonZeroElem *> NZE_Mem_Allocated;

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

@ -46,6 +46,8 @@ SparseMatrix::SparseMatrix()
start_compare = 0;
restart = 0;
IM_i.clear();
lu_inc_tol = 1e-10;
reduced = true;
}
int
@ -278,7 +280,7 @@ SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_
}
void
SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries)
SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo)
{
unsigned int eq, var;
int i, j, lag;
@ -302,6 +304,8 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
}
IM_i.clear();
if (two_boundaries)
{
if (stack_solve_algo == 5)
{
for (i = 0; i < u_count_init-Size; i++)
{
@ -314,7 +318,24 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
for (j = 0; j < Size; j++)
IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
}
else if (stack_solve_algo >= 1 || stack_solve_algo <= 4)
{
for (i = 0; i < u_count_init-Size; i++)
{
SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
IM_i[make_pair(make_pair(var - lag*Size, -lag), eq)] = j;
}
for (j = 0; j < Size; j++)
IM_i[make_pair(make_pair(Size*(periods+y_kmax), 0), j)] = j;
}
}
else
{
if ((stack_solve_algo == 5 && !steady_state) || (solve_algo == 5 && steady_state))
{
for (i = 0; i < u_count_init; i++)
{
@ -325,6 +346,18 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i
IM_i[make_pair(make_pair(eq, var), lag)] = j;
}
}
else if ( ((stack_solve_algo >= 1 || stack_solve_algo <= 4) && !steady_state) || ((solve_algo >= 1 || solve_algo <= 4) && steady_state) )
{
for (i = 0; i < u_count_init; i++)
{
SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
IM_i[make_pair(make_pair(var - lag*Size, -lag), eq)] = j;
}
}
}
index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int));
for (j = 0; j < Size; j++)
SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara));
@ -415,7 +448,230 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pa
}
void
SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
SparseMatrix::Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m)
{
int i, eq, var;
double *b = mxGetPr(b_m);
if (!b)
{
mexPrintf("Can't retrieve b vector in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
mwIndex *Ai = mxGetIr(A_m);
if (!Ai)
{
mexPrintf("Can't allocate Ai index vector in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
mwIndex *Aj = mxGetJc(A_m);
if (!Aj)
{
mexPrintf("Can't allocate Aj index vector in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
double *A = mxGetPr(A_m);
if (!A)
{
mexPrintf("Can't retrieve A matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
map<pair<pair<int, int>, int>, int>::iterator it4;
for (i = 0; i < y_size*(periods+y_kmin); i++)
ya[i] = y[i];
#if DEBUG
unsigned int max_nze = mxGetNzmax(A_m);
#endif
unsigned int NZE = 0;
int last_var = 0;
for (i = 0; i < Size; i++)
b[i] = u[i];
Aj[0] = 0;
last_var = 0;
it4 = IM.begin();
while (it4 != IM.end())
{
var = it4->first.first.first;
if (var != last_var)
{
Aj[1+last_var ] = NZE;
last_var = var;
}
eq = it4->first.second;
int index = it4->second;
#if DEBUG
if (index<0 || index >= u_count_alloc)
{
mexPrintf("index (%d) out of range for u vector (0)\n",index);
mexErrMsgTxt("end of bytecode\n");
}
if (NZE >= max_nze)
{
mexPrintf("Exceeds the capacity of A_m sparse matrix in LU (0)\n");
mexErrMsgTxt("end of bytecode\n");
}
#endif
A[NZE] = u[index];
Ai[NZE] = eq;
NZE++;
#if DEBUG
if ((index+lag*u_count_init) < 0 || (index+lag*u_count_init) >= u_count_alloc)
{
mexPrintf("index (%d) out of range for u vector (1)\n",index+lag*u_count_init);
mexErrMsgTxt("end of bytecode\n");
}
if (eq < 0 || eq >= (Size*periods))
{
mexPrintf("index (%d) out of range for b vector (0)\n",eq);
mexErrMsgTxt("end of bytecode\n");
}
if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
{
mexPrintf("index (%d) out of range for index_vara vector (0)\n",var+Size*(y_kmin+t+lag));
mexErrMsgTxt("end of bytecode\n");
}
if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
{
mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var+Size*(y_kmin+t+lag)], y_size*(periods+y_kmin+y_kmax));
mexErrMsgTxt("end of bytecode\n");
}
#endif
it4++;
}
Aj[Size] = NZE;
}
void
SparseMatrix::Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m)
{
int t, i, eq, var, lag, ti_y_kmin, ti_y_kmax;
double *b = mxGetPr(b_m);
if (!b)
{
mexPrintf("Can't retrieve b vector in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
mwIndex *Ai = mxGetIr(A_m);
if (!Ai)
{
mexPrintf("Can't allocate Ai index vector in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
mwIndex *Aj = mxGetJc(A_m);
if (!Aj)
{
mexPrintf("Can't allocate Aj index vector in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
double *A = mxGetPr(A_m);
if (!A)
{
mexPrintf("Can't retrieve A matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
map<pair<pair<int, int>, int>, int>::iterator it4;
for (i = 0; i < y_size*(periods+y_kmin); i++)
ya[i] = y[i];
#if DEBUG
unsigned int max_nze = mxGetNzmax(A_m);
#endif
unsigned int NZE = 0;
int last_var = 0;
for (i = 0; i < periods*Size; i++)
b[i] = 0;
Aj[0] = 0;
for (t = 0; t < periods; t++)
{
last_var = 0;
it4 = IM.begin();
while (it4 != IM.end())
{
var = it4->first.first.first;
if (var != last_var)
{
Aj[1+last_var + t * Size] = NZE;
last_var = var;
}
eq = it4->first.second+Size*t;
lag = -it4->first.first.second;
int index = it4->second+ (t-lag) * u_count_init;
if (var < (periods+y_kmax)*Size)
{
ti_y_kmin = -min(t, y_kmin);
ti_y_kmax = min(periods-(t +1 ), y_kmax);
int ti_new_y_kmax = min(t, y_kmax);
int ti_new_y_kmin = -min(periods-(t+1), y_kmin);
if (lag <= ti_new_y_kmax && lag >= ti_new_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/
{
#if DEBUG
if (index<0 || index >= u_count_alloc)
{
mexPrintf("index (%d) out of range for u vector (0)\n",index);
mexErrMsgTxt("end of bytecode\n");
}
if (NZE >= max_nze)
{
mexPrintf("Exceeds the capacity of A_m sparse matrix in LU (0)\n");
mexErrMsgTxt("end of bytecode\n");
}
#endif
A[NZE] = u[index];
Ai[NZE] = eq - lag * Size;
NZE++;
}
if (lag > ti_y_kmax || lag < ti_y_kmin)
{
#if DEBUG
if ((index+lag*u_count_init) < 0 || (index+lag*u_count_init) >= u_count_alloc)
{
mexPrintf("index (%d) out of range for u vector (1)\n",index+lag*u_count_init);
mexErrMsgTxt("end of bytecode\n");
}
if (eq < 0 || eq >= (Size*periods))
{
mexPrintf("index (%d) out of range for b vector (0)\n",eq);
mexErrMsgTxt("end of bytecode\n");
}
if (var+Size*(y_kmin+t+lag) < 0 || var+Size*(y_kmin+t+lag) >= Size*(periods+y_kmin+y_kmax))
{
mexPrintf("index (%d) out of range for index_vara vector (0)\n",var+Size*(y_kmin+t+lag));
mexErrMsgTxt("end of bytecode\n");
}
if (index_vara[var+Size*(y_kmin+t+lag)] < 0 || index_vara[var+Size*(y_kmin+t+lag)] >= y_size*(periods+y_kmin+y_kmax))
{
mexPrintf("index (%d) out of range for y vector max=%d (0)\n",index_vara[var+Size*(y_kmin+t+lag)], y_size*(periods+y_kmin+y_kmax));
mexErrMsgTxt("end of bytecode\n");
}
#endif
b[eq] += u[index+lag*u_count_init]*y[index_vara[var+Size*(y_kmin+t+lag)]];
}
}
else /* ...and store it in the u vector*/
{
#if DEBUG
if (index < 0 || index >= u_count_alloc)
{
mexPrintf("index (%d) out of range for u vector (2)\n",index);
mexErrMsgTxt("end of bytecode\n");
}
if (eq < 0 || eq >= (Size*periods))
{
mexPrintf("index (%d) out of range for b vector (1)\n",eq);
mexErrMsgTxt("end of bytecode\n");
}
#endif
b[eq] += u[index];
}
it4++;
}
}
Aj[Size*periods] = NZE;
}
void
SparseMatrix::Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
{
int t, i, eq, var, lag, ti_y_kmin, ti_y_kmax;
double tmp_b = 0.0;
@ -524,47 +780,6 @@ SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<
mxFree(temp_NZE_C);
}
void
SparseMatrix::ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
{
int t, eq, var, lag, ti_y_kmin, ti_y_kmax;
double tmp_b = 0.0;
map<pair<pair<int, int>, int>, int>::iterator it4;
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(it4, ti_y_kmin, ti_y_kmax, eq, var, lag, tmp_b) schedule(dynamic)
for (t = 0; t < periods; t++)
{
ti_y_kmin = -min(t, y_kmin);
ti_y_kmax = min(periods-(t+1), y_kmax);
it4 = IM.begin();
eq = -1;
while (it4 != IM.end())
{
var = it4->first.first.second;
if (eq != it4->first.first.first+Size*t)
tmp_b = 0;
eq = it4->first.first.first+Size*t;
if (var < (periods+y_kmax)*Size)
{
lag = it4->first.second;
if (lag <= ti_y_kmax && lag >= ti_y_kmin)
{
var += Size*t;
}
else
{
tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
}
}
else
{
b[eq] = it4->second+u_count_init*t;
u[b[eq]] += tmp_b;
}
it4++;
}
}
}
int
SparseMatrix::Get_u()
{
@ -619,7 +834,7 @@ SparseMatrix::Print_u()
}
void
SparseMatrix::End(int Size)
SparseMatrix::End_GE(int Size)
{
mem_mngr.Free_All();
mxFree(FNZE_R);
@ -995,7 +1210,7 @@ SparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
}
bool
SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number)
SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number, int solve_algo)
{
int i, j, k;
int pivj = 0, pivk = 0;
@ -1005,6 +1220,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
int *pivj_v, *pivk_v, *NR;
int l, N_max;
bool one;
mxArray *b_m, *A_m;
Clear_u();
piv_v = (double *) mxMalloc(Size*sizeof(double));
pivj_v = (int *) mxMalloc(Size*sizeof(int));
@ -1131,7 +1347,27 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
}
if (solve_algo == 5)
Simple_Init(it_, y_kmin, y_kmax, Size, IM_i);
else
{
b_m = mxCreateDoubleMatrix(Size,1,mxREAL);
if (!b_m)
{
mexPrintf("Can't allocate b_m matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
A_m = mxCreateSparse(Size, Size, IM_i.size()*2, mxREAL);
if (!A_m)
{
mexPrintf("Can't allocate A_m matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
Init_Matlab_Sparse_Simple(Size, IM_i, A_m, b_m);
}
if (solve_algo == 5)
{
NonZeroElem **bc;
bc = (NonZeroElem **) mxMalloc(Size*sizeof(*bc));
for (i = 0; i < Size; i++)
@ -1261,10 +1497,6 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
}
}
}
/*if (fabs(piv) < eps)
mexPrintf("==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f\n",NR_max, N_max, piv, piv_abs, markovitz_max);
if (NR_max == 0)
mexPrintf("==> Error NR_max=0 and piv=%f, markovitz_max=%f\n",piv, markovitz_max);*/
pivot[i] = pivj;
pivotk[i] = pivk;
pivotv[i] = piv;
@ -1372,12 +1604,19 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
ya[i+it_*y_size] = y[i+it_*y_size];
slowc_save = slowc;
res1bx = simple_bksub(it_, Size, slowc_lbx);
End(Size);
End_GE(Size);
mxFree(piv_v);
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
mxFree(bc);
}
else if (solve_algo == 1 || solve_algo == 4)
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc);
else if (solve_algo == 2)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck);
else if (solve_algo == 3)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck);
return true;
}
@ -1456,7 +1695,7 @@ void
SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b)
{
const double epsilon = 1e-10;
Init(periods, y_kmin, y_kmax, Size, IM_i);
Init_GE(periods, y_kmin, y_kmax, Size, IM_i);
NonZeroElem *first;
int cal_y = y_kmin*Size;
mexPrintf(" ");
@ -1492,8 +1731,170 @@ SparseMatrix::Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size,
mexErrMsgTxt(filename.c_str());
}
void
SparseMatrix::Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l)
{
int n = mxGetM(A_m);
mxArray *z;
mxArray *rhs[2];
rhs[0] = A_m;
rhs[1] = b_m;
mexCallMATLAB(1,&z,2, rhs, "mldivide");
double *res = mxGetPr(z);
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc_l * yy;
}
mxDestroyArray(A_m);
mxDestroyArray(b_m);
mxDestroyArray(z);
}
void
SparseMatrix::Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block)
{
int n = mxGetM(A_m);
/*[L1, U1]=luinc(g1a,luinc_tol);*/
mxArray *lhs0[2];
mxArray *rhs0[2];
rhs0[0] = A_m;
rhs0[1] = mxCreateDoubleScalar(lu_inc_tol);
mexCallMATLAB(2, lhs0, 2, rhs0, "luinc");
mxArray *L1 = lhs0[0];
mxArray *U1 = lhs0[1];
/*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
mxArray *rhs[7];
rhs[0] = A_m;
rhs[1] = b_m;
rhs[2] = mxCreateDoubleScalar(Size);
rhs[3] = mxCreateDoubleScalar(1e-6);
rhs[4] = mxCreateDoubleScalar(n);
rhs[5] = L1;
rhs[6] = U1;
mxArray *lhs[2];
mexCallMATLAB(2,lhs, 7, rhs, "gmres");
mxArray *z = lhs[0];
mxArray *flag = lhs[1];
double *flag1 = mxGetPr(flag);
mxDestroyArray(rhs0[1]);
mxDestroyArray(rhs[2]);
mxDestroyArray(rhs[3]);
mxDestroyArray(rhs[4]);
mxDestroyArray(rhs[5]);
mxDestroyArray(rhs[6]);
if (*flag1 > 0 || reduced)
{
ostringstream tmp;
if (*flag1 == 1)
{
tmp << "Error in simul: No convergence inside GMRES, in block " << block;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 2)
{
tmp << "Error in simul: Preconditioner is ill-conditioned, in block " << block;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 3)
{
tmp << "Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block " << block;
mexWarnMsgTxt(tmp.str().c_str());
}
lu_inc_tol /= 10;
reduced = false;
}
else
{
double *res = mxGetPr(z);
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
}
mxDestroyArray(A_m);
mxDestroyArray(b_m);
mxDestroyArray(z);
mxDestroyArray(flag);
}
void
SparseMatrix::Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block)
{
int n = mxGetM(A_m);
/*[L1, U1]=luinc(g1a,luinc_tol);*/
mxArray *lhs0[2];
mxArray *rhs0[2];
rhs0[0] = A_m;
rhs0[1] = mxCreateDoubleScalar(lu_inc_tol);
mexCallMATLAB(2, lhs0, 2, rhs0, "luinc");
mxArray *L1 = lhs0[0];
mxArray *U1 = lhs0[1];
/*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
mxArray *rhs[6];
rhs[0] = A_m;
rhs[1] = b_m;
rhs[2] = mxCreateDoubleScalar(1e-6);
rhs[3] = mxCreateDoubleScalar(n);
rhs[4] = L1;
rhs[5] = U1;
mxArray *lhs[2];
mexCallMATLAB(2,lhs, 6, rhs, "bicgstab");
mxArray *z = lhs[0];
mxArray *flag = lhs[1];
double *flag1 = mxGetPr(flag);
mxDestroyArray(rhs0[1]);
mxDestroyArray(rhs[2]);
mxDestroyArray(rhs[3]);
mxDestroyArray(rhs[4]);
mxDestroyArray(rhs[5]);
if (*flag1 > 0 || reduced)
{
ostringstream tmp;
if (*flag1 == 1)
{
tmp << "Error in simul: No convergence inside BiCGStab, in block " << block;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 2)
{
tmp << "Error in simul: Preconditioner is ill-conditioned, in block " << block;
mexWarnMsgTxt(tmp.str().c_str());
}
else if (*flag1 == 3)
{
tmp << "Error in simul: BiCGStab stagnated (Two consecutive iterates were the same.), in block " << block;
mexWarnMsgTxt(tmp.str().c_str());
}
lu_inc_tol /= 10;
reduced = false;
}
else
{
double *res = mxGetPr(z);
for (int i = 0; i < n; i++)
{
int eq = index_vara[i+Size*y_kmin];
double yy = - (res[i] + y[eq]);
direction[eq] = yy;
y[eq] += slowc * yy;
}
}
mxDestroyArray(A_m);
mxDestroyArray(b_m);
mxDestroyArray(z);
mxDestroyArray(flag);
}
int
SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number)
SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo)
{
/*Triangularisation at each period of a block using a simple gaussian Elimination*/
t_save_op_s *save_op_s;
@ -1513,6 +1914,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
clock_t t1 = clock();
nop1 = 0;
error_not_printed = true;
mxArray *b_m, *A_m;
if (iter > 0)
{
mexPrintf("Sim : %f ms\n", (1000.0*(double (clock())-double (time00)))/double (CLOCKS_PER_SEC));
@ -1565,7 +1967,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
filename += " stopped";
mexErrMsgTxt(filename.c_str());
}
if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)))
if(!(isnan(res1) || isinf(res1)) && !(isnan(g0) || isinf(g0)) && (stack_solve_algo == 1 || stack_solve_algo == 5))
{
if (try_at_iteration == 0)
{
@ -1583,6 +1985,16 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
glambda2 = res2;
try_at_iteration ++;
if (slowc_save<=0.1)
{
for (i = 0; i < y_size*(periods+y_kmin); i++)
y[i] = ya[i]+direction[i];
g0 = res2;
gp0 = -res2;
try_at_iteration = 0;
iter--;
return 0;
}
}
else
{
@ -1596,7 +2008,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
iter--;
return 0;
}
/*if (isnan(res1) || isinf(res1))
if (isnan(res1) || isinf(res1))
{
if (iter == 0)
{
@ -1629,8 +2041,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
y[i] = ya[i]+slowc_save*direction[i];
iter--;
return (0);
}*/
}
u_count += u_count_init;
if (stack_solve_algo == 5)
{
if (alt_symbolic && alt_symbolic_count < alt_symbolic_count_max)
{
mexPrintf("Pivoting method will be applied only to the first periods.\n");
@ -1661,16 +2075,41 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
start_compare = max(y_kmin, minimal_solving_periods);
restart = 0;
}
}
res1a = res1;
if (print_it)
{
if (iter == 0)
{
switch(stack_solve_algo)
{
case 1:
mexPrintf("MODEL SIMULATION: (method=Sparse LU)\n");
break;
case 2:
mexPrintf("MODEL SIMULATION: (method=GMRES)\n");
break;
case 3:
mexPrintf("MODEL SIMULATION: (method=BiCGStab)\n");
break;
case 4:
mexPrintf("MODEL SIMULATION: (method=Sparse LU & optimal path length)\n");
break;
case 5:
mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n");
break;
default:
mexPrintf("MODEL SIMULATION: (method=Unknown - %d - )\n", stack_solve_algo);
}
}
mexPrintf("-----------------------------------\n");
mexPrintf(" Simulate iteration no %d \n", iter+1);
mexPrintf(" max. error=%.10e \n", double (max_res));
mexPrintf(" sqr. error=%.10e \n", double (res2));
mexPrintf(" abs. error=%.10e \n", double (res1));
mexPrintf("-----------------------------------\n");
mexEvalString("drawnow;");
}
if (cvg)
{
@ -1678,7 +2117,26 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
else
{
Init(periods, y_kmin, y_kmax, Size, IM_i);
if (stack_solve_algo == 5)
Init_GE(periods, y_kmin, y_kmax, Size, IM_i);
else
{
b_m = mxCreateDoubleMatrix(periods*Size,1,mxREAL);
if (!b_m)
{
mexPrintf("Can't allocate b_m matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
A_m = mxCreateSparse(periods*Size, periods*Size, IM_i.size()* periods*2, mxREAL);
if (!A_m)
{
mexPrintf("Can't allocate A_m matrix in LU solver\n");
mexErrMsgTxt("end of bytecode\n");
}
Init_Matlab_Sparse(periods, y_kmin, y_kmax, Size, IM_i, A_m, b_m);
}
if (stack_solve_algo == 5)
{
double *piv_v;
int *pivj_v, *pivk_v, *NR;
piv_v = (double *) mxMalloc(Size*sizeof(double));
@ -1997,19 +2455,12 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
else /*first_sub->c_index==first_piv->c_index*/
{
if (i == sub_c_index)
{
//#pragma omp barrier
//#pragma omp single
//#pragma omp critical
{
NonZeroElem *firsta = first;
NonZeroElem *first_suba = first_sub->NZE_R_N;
Delete(first_sub->r_index, first_sub->c_index);
first = firsta->NZE_C_N;
first_sub = first_suba;
}
if (first_sub)
sub_c_index = first_sub->c_index;
else
@ -2137,7 +2588,6 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mxFree(pivj_v);
mxFree(pivk_v);
mxFree(NR);
}
nop_all += nop;
if (symbolic)
{
@ -2156,7 +2606,15 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
slowc_save = slowc;
res1bx = bksub(tbreak, last_period, Size, slowc_lbx);
t01 = clock();
End(Size);
End_GE(Size);
}
else if (stack_solve_algo == 1 || stack_solve_algo == 4)
Solve_Matlab_LU_UMFPack(A_m, b_m, Size, slowc);
else if (stack_solve_algo == 2)
Solve_Matlab_GMRES(A_m, b_m, Size, slowc, blck);
else if (stack_solve_algo == 3)
Solve_Matlab_BiCGStab(A_m, b_m, Size, slowc, blck);
}
if (print_it)
{
clock_t t2 = clock();

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

@ -33,7 +33,10 @@ using namespace std;
#ifdef _MSC_VER
# include <limits>
#include <boost/math/special_functions/erf.hpp>
#define M_SQRT2 1.4142135623730950488016887242097
#define M_PI 3.1415926535897932384626433832795
#define erf(x) boost::math::erf(x)
extern unsigned long _nan[2];
extern double NAN;
@ -100,18 +103,23 @@ class SparseMatrix
{
public:
SparseMatrix();
int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number);
bool simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number);
int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo);
bool simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number, int solve_algo);
void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter);
void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries);
void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo);
void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u);
double g0, gp0, glambda2, try_at_iteration;
private:
void Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM);
void End(int Size);
void End_GE(int Size);
void Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l);
void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER
, long int *ndiv, long int *nsub
@ -173,6 +181,8 @@ private:
protected:
int u_count_alloc, u_count_alloc_save;
double *u, *y, *ya;
vector<double*> jac;
double *jcb;
double res1, res2, max_res, max_res_idx;
double slowc, slowc_save, prev_slowc_save, markowitz_c;
int y_kmin, y_kmax, y_size, periods, y_decal;
@ -184,6 +194,8 @@ protected:
int restart;
bool error_not_printed;
double g_lambda1, g_lambda2, gp_0;
double lu_inc_tol;
bool reduced;
};
#endif

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

@ -54,15 +54,8 @@ main(int argc, const char *argv[])
FILE *fid;
bool steady_state = false;
bool evaluate = false;
printf("argc=%d\n", argc);
/*fexcept_t *flagp;
flagp = (fexcept_t*) mxMalloc(sizeof(fexcept_t));
if (fegetexceptflag(flagp, FE_ALL_EXCEPT))
mexPrintf("fegetexceptflag failed\n");
if (fesetexceptflag(flagp,FE_INVALID | FE_DIVBYZERO))
mexPrintf("fesetexceptflag failed\n");
mxFree(flagp);
feclearexcept (FE_ALL_EXCEPT);*/
bool block = false;
if (argc < 2)
{
mexPrintf("model filename expected\n");
@ -72,7 +65,7 @@ main(int argc, const char *argv[])
float f_tmp;
ostringstream tmp_out("");
tmp_out << argv[1] << "_options.txt";
cout << tmp_out.str().c_str() << "\n";
int nb_params;
int i, row_y, col_y, row_x, col_x;
double *yd, *xd;
@ -81,6 +74,7 @@ main(int argc, const char *argv[])
string file_name(argv[1]);
int count_array_argument = 0;
for (i = 2; i < argc; i++)
{
if (Get_Argument(argv[i]) == "static")
@ -89,6 +83,8 @@ main(int argc, const char *argv[])
steady_state = false;
else if (Get_Argument(argv[i]) == "evaluate")
evaluate = true;
else if (Get_Argument(argv[i]) == "block")
block = true;
else
{
mexPrintf("Unknown argument : ");
@ -99,6 +95,7 @@ main(int argc, const char *argv[])
mexErrMsgTxt(f.c_str());
}
}
fid = fopen(tmp_out.str().c_str(), "r");
int periods = 1;
if (!steady_state)
@ -143,10 +140,12 @@ main(int argc, const char *argv[])
params[i] = f_tmp;
}
fclose(fid);
yd = (double *) malloc(row_y*col_y*sizeof(yd[0]));
xd = (double *) malloc(row_x*col_x*sizeof(xd[0]));
tmp_out.str("");
tmp_out << argv[1] << "_oo.txt";
fid = fopen(tmp_out.str().c_str(), "r");
for (i = 0; i < col_y*row_y; i++)
{
@ -193,7 +192,8 @@ main(int argc, const char *argv[])
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods);
string f(file_name);
interprete.compute_blocks(f, f, steady_state, evaluate);
int nb_blocks = 0;
interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks);
clock_t t1 = clock();
if (!evaluate)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
@ -235,29 +235,56 @@ void
mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mxArray *M_, *oo_, *options_;
int i, row_y, col_y, row_x, col_x, nb_row_xd;
mxArray *block_structur = NULL;
int i, row_y = 0, col_y = 0, row_x = 0, col_x = 0, nb_row_xd = 0;
int steady_row_y, steady_col_y, steady_row_x, steady_col_x, steady_nb_row_xd;
int y_kmin = 0, y_kmax = 0, y_decal = 0, periods = 1;
double *pind;
double *direction;
bool steady_state = false;
bool evaluate = false;
/*fexcept_t *flagp;
flagp = (fexcept_t*) mxMalloc(sizeof(fexcept_t));
if(fegetexceptflag(flagp, FE_ALL_EXCEPT))
mexPrintf("fegetexceptflag failed\n");
if(fesetexceptflag(flagp,FE_INVALID | FE_DIVBYZERO))
mexPrintf("fesetexceptflag failed\n");
mxFree(flagp);
feclearexcept (FE_ALL_EXCEPT);*/
bool block = false;
double *params = NULL;
double *yd = NULL, *xd = NULL;
int count_array_argument = 0;
for (i = 0; i < nrhs; i++)
{
if (Get_Argument(prhs[i]) == "static")
if (!mxIsChar(prhs[i]))
{
switch (count_array_argument)
{
case 0:
yd = mxGetPr(prhs[i]);
row_y = mxGetM(prhs[i]);
col_y = mxGetN(prhs[i]);
break;
case 1:
xd = mxGetPr(prhs[i]);
row_x = mxGetM(prhs[i]);
col_x = mxGetN(prhs[i]);
break;
case 2:
params = mxGetPr(prhs[i]);
break;
case 3:
periods = mxGetScalar(prhs[i]);
break;
case 4:
block_structur = mxDuplicateArray(prhs[i]);
break;
default:
mexPrintf("Unknown argument\n");
}
count_array_argument++;
}
else if (Get_Argument(prhs[i]) == "static")
steady_state = true;
else if (Get_Argument(prhs[i]) == "dynamic")
steady_state = false;
else if (Get_Argument(prhs[i]) == "evaluate")
evaluate = true;
else if (Get_Argument(prhs[i]) == "block")
block = true;
else
{
mexPrintf("Unknown argument : ");
@ -267,6 +294,12 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexErrMsgTxt(f.c_str());
}
}
if (count_array_argument > 0 && count_array_argument < 4)
{
mexPrintf("Missing arguments. All the following arguments have to be indicated y, x, params, it_\n");
mexErrMsgTxt("end of bytecode\n");
}
M_ = mexGetVariable("global", "M_");
if (M_ == NULL)
{
@ -286,21 +319,29 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mexPrintf("Global variable not found : ");
mexErrMsgTxt("options_ \n");
}
double *params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
double *yd, *xd, *steady_yd = NULL, *steady_xd = NULL;
if (!count_array_argument)
params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "params")));
double *steady_yd = NULL, *steady_xd = NULL;
if (!steady_state)
{
if (!count_array_argument)
{
yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));
col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "endo_simul")));;
}
if (!count_array_argument)
{
xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_simul")));
}
nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
y_kmin = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_lag"))))));
y_kmax = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_lead"))))));
y_decal = max(0, y_kmin-int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "maximum_endo_lag")))))));
if (!count_array_argument)
periods = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "periods"))))));
steady_yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
@ -312,19 +353,27 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
steady_nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
}
else
{
if (!count_array_argument)
{
yd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
row_y = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));
col_y = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "steady_state")));;
}
if (!count_array_argument)
{
xd = mxGetPr(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
row_x = mxGetM(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
col_x = mxGetN(mxGetFieldByNumber(oo_, 0, mxGetFieldNumber(oo_, "exo_steady_state")));
}
nb_row_xd = int (floor(*(mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_, "exo_det_nbr"))))));
}
int maxit_ = int (floor(*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "maxit_"))))));
double slowc = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "slowc")))));
double markowitz_c = double (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "markowitz")))));
int minimal_solving_periods = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "minimal_solving_periods")))));
int stack_solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "stack_solve_algo")))));
int solve_algo = int (*(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_algo")))));
double solve_tolf;
if (steady_state)
solve_tolf = *(mxGetPr(mxGetFieldByNumber(options_, 0, mxGetFieldNumber(options_, "solve_tolf"))));
@ -356,9 +405,10 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
int nb_row_x = row_x;
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods);
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal, markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo);
string f(fname);
bool result = interprete.compute_blocks(f, f, steady_state, evaluate);
int nb_blocks = 0;
bool result = interprete.compute_blocks(f, f, steady_state, evaluate, block, nb_blocks);
clock_t t1 = clock();
if (!steady_state && !evaluate)
mexPrintf("Simulation Time=%f milliseconds\n", 1000.0*(double (t1)-double (t0))/double (CLOCKS_PER_SEC));
@ -380,6 +430,62 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
pind[0] = 0;
else
pind[0] = 1;
if (nlhs > 2)
{
int jacob_field_number = 0, jacob_exo_field_number = 0, jacob_exo_det_field_number = 0, jacob_other_endo_field_number = 0;
if (!block_structur)
{
const char *field_names[] = {"jacob","jacob_exo","jacob_exo_det","jacob_other_endo"};
jacob_field_number=0;
jacob_exo_field_number=1;
jacob_exo_det_field_number=2;
jacob_other_endo_field_number=2;
mwSize dims[1] = {nb_blocks };
plhs[2] = mxCreateStructArray(1, dims, 4, field_names);
mexPrintf("the structure has been created\n");
}
else if (!mxIsStruct(block_structur))
{
mexPrintf("The third output argument must be a structure\n");
mexErrMsgTxt("end of bytecode\n");
}
else
{
mexPrintf("Adding Fields\n");
jacob_field_number = mxAddField(block_structur, "jacob");
if (jacob_field_number == -1)
{
mexPrintf("Cannot add extra field to the structArray\n");
mexErrMsgTxt("end of bytecode\n");
}
jacob_exo_field_number = mxAddField(block_structur, "jacob_exo");
if (jacob_exo_field_number == -1)
{
mexPrintf("Cannot add extra field to the structArray\n");
mexErrMsgTxt("end of bytecode\n");
}
jacob_exo_det_field_number = mxAddField(block_structur, "jacob_exo_det");
if (jacob_exo_det_field_number == -1)
{
mexPrintf("Cannot add extra field to the structArray\n");
mexErrMsgTxt("end of bytecode\n");
}
jacob_other_endo_field_number = mxAddField(block_structur, "jacob_other_endo");
if (jacob_other_endo_field_number == -1)
{
mexPrintf("Cannot add extra field to the structArray\n");
mexErrMsgTxt("end of bytecode\n");
}
}
for (int i = 0; i < nb_blocks; i++)
{
mxSetFieldByNumber(block_structur,i,jacob_field_number,interprete.get_jacob(i));
mxSetFieldByNumber(block_structur,i,jacob_exo_field_number,interprete.get_jacob_exo(i));
mxSetFieldByNumber(block_structur,i,jacob_exo_det_field_number,interprete.get_jacob_exo_det(i));
mxSetFieldByNumber(block_structur,i,jacob_other_endo_field_number,interprete.get_jacob_other_endo(i));
}
plhs[2] = block_structur;
}
}
}
if (x)

298
preprocessor/CodeInterpreter.hh
View File

@ -57,47 +57,53 @@ using namespace std;
*/
enum Tags
{
FLDZ, //!< Stores zero in the stack - 0
FLDC, //!< Stores a constant term in the stack - 1
FLDZ, //!< Stores zero in the stack - 0 (0)
FLDC, //!< Stores a constant term in the stack - 1 (1)
FDIMT, //!< Defines the number of temporary terms - dynamic context (the period has to be indicated) - 2
FDIMST, //!< Defines the number of temporary terms - static context (the period hasn't to be indicated) - 3
FLDT, //!< Stores a temporary term in the stack - dynamic context (the period has to be indicated) - 4
FLDST, //!< Stores a temporary term in the stack - static context (the period hasn't to be indicated) - 5
FSTPT, //!< Loads a temporary term from the stack - dynamic context (the period has to be indicated) - 6
FSTPST, //!< Loads a temporary term from the stack - static context (the period hasn't to be indicated) - 7
FDIMT, //!< Defines the number of temporary terms - dynamic context (the period has to be indicated) - 2 (2)
FDIMST, //!< Defines the number of temporary terms - static context (the period hasn't to be indicated) - 3 (3)
FLDT, //!< Stores a temporary term in the stack - dynamic context (the period has to be indicated) - 4 (4)
FLDST, //!< Stores a temporary term in the stack - static context (the period hasn't to be indicated) - 5 (5)
FSTPT, //!< Loads a temporary term from the stack - dynamic context (the period has to be indicated) - 6 (6)
FSTPST, //!< Loads a temporary term from the stack - static context (the period hasn't to be indicated) - 7 (7)
FLDU, //!< Stores an element of the vector U in the stack - dynamic context (the period has to be indicated) - 8
FLDSU, //!< Stores an element of the vector U in the stack - static context (the period hasn't to be indicated) - 9
FSTPU, //!< Loads an element of the vector U from the stack - dynamic context (the period has to be indicated) - A
FSTPSU, //!< Loads an element of the vector U from the stack - static context (the period hasn't to be indicated) - B
FLDU, //!< Stores an element of the vector U in the stack - dynamic context (the period has to be indicated) - 8 (8)
FLDSU, //!< Stores an element of the vector U in the stack - static context (the period hasn't to be indicated) - 9 (9)
FSTPU, //!< Loads an element of the vector U from the stack - dynamic context (the period has to be indicated) - A (10)
FSTPSU, //!< Loads an element of the vector U from the stack - static context (the period hasn't to be indicated) - B (11)
FLDV, //!< Stores a variable (described in SymbolType) in the stack - dynamic context (the period has to be indicated) - C
FLDSV, //!< Stores a variable (described in SymbolType) in the stack - static context (the period hasn't to be indicated) - D
FLDVS, //!< Stores a variable (described in SymbolType) in the stack - dynamic context but inside the STEADYSTATE function (the period hasn't to be indicated) - E
FSTPV, //!< Loads a variable (described in SymbolType) from the stack - dynamic context (the period has to be indicated) - F
FSTPSV, //!< Loads a variable (described in SymbolType) from the stack - static context (the period hasn't to be indicated) - 10
FLDV, //!< Stores a variable (described in SymbolType) in the stack - dynamic context (the period has to be indicated) - C (12)
FLDSV, //!< Stores a variable (described in SymbolType) in the stack - static context (the period hasn't to be indicated) - D (13)
FLDVS, //!< Stores a variable (described in SymbolType) in the stack - dynamic context but inside the STEADYSTATE function (the period hasn't to be indicated) - E (14)
FSTPV, //!< Loads a variable (described in SymbolType) from the stack - dynamic context (the period has to be indicated) - F (15)
FSTPSV, //!< Loads a variable (described in SymbolType) from the stack - static context (the period hasn't to be indicated) - 10 (16)
FLDR, //!< Stores a residual in the stack - 11
FSTPR, //!< Loads a residual from the stack - 12
FLDR, //!< Stores a residual in the stack - 11 (17)
FSTPR, //!< Loads a residual from the stack - 12 (18)
FSTPG, //!< Loads a derivative from the stack - 13
FSTPG2, //!< Loads a derivative matrix from the stack - 13
FSTPG, //!< Loads a derivative from the stack - 13 (19)
FSTPG2, //!< Loads a derivative matrix for static model from the stack - 14 (20)
FSTPG3, //!< Loads a derivative matrix for a dynamic model from the stack - 15 (21)
FSTPG4, //!< Loads a second order derivative matrix for a dynamic model from the stack - 16 (22)
FUNARY, //!< A Unary operator - 14
FBINARY, //!< A binary operator - 15
FTRINARY, //!< A trinary operator - 15'
FCUML, //!< Cumulates the result - 16
FUNARY, //!< A Unary operator - 17 (23)
FBINARY, //!< A binary operator - 18 (24)
FTRINARY, //!< A trinary operator - 19 (25)
FBEGINBLOCK, //!< Defines the begining of a model block - 17
FENDBLOCK, //!< Defines the end of a model block - 18
FENDEQU, //!< Defines the last equation of the block. For block that has to be solved, the derivatives appear just after this flag - 19
FEND, //!< Defines the end of the model code - 1A
FCUML, //!< Cumulates the result - 1A (26)
FOK, //!< Used for debugging purpose - 1B
FJMPIFEVAL, //!< Jump if evaluate = true - 1B (27)
FJMP, //!< Jump - 1C (28)
FNUMEXPR //!< Store the expression type and references
FBEGINBLOCK, //!< Defines the begining of a model block - 1D (29)
FENDBLOCK, //!< Defines the end of a model block - 1E (30)
FENDEQU, //!< Defines the last equation of the block. For block that has to be solved, the derivatives appear just after this flag - 1F (31)
FEND, //!< Defines the end of the model code - 20 (32)
FOK, //!< Used for debugging purpose - 21 (33)
FNUMEXPR //!< Store the expression type and references - 22 (34)
};
@ -148,6 +154,7 @@ enum ExpressionType
TemporaryTerm,
ModelEquation,
FirstEndoDerivative,
FirstOtherEndoDerivative,
FirstExoDerivative,
FirstExodetDerivative,
FirstParamDerivative,
@ -224,9 +231,10 @@ public:
{
};
inline void
write(ostream &CompileCode)
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
instruction_number++;
};
};
@ -244,9 +252,10 @@ public:
{
};
inline void
write(ostream &CompileCode)
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(this), sizeof(TagWithOneArgument));
instruction_number++;
};
};
@ -265,9 +274,10 @@ public:
{
};
inline void
write(ostream &CompileCode)
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
instruction_number++;
};
};
@ -287,12 +297,40 @@ public:
{
};
inline void
write(ostream &CompileCode)
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
instruction_number++;
};
};
template < class T1, class T2, class T3, class T4 >
class TagWithFourArguments
{
protected:
uint8_t op_code;
T1 arg1;
T2 arg2;
T3 arg3;
T4 arg4;
public:
inline TagWithFourArguments(uint8_t op_code_arg) : op_code(op_code_arg)
{
};
inline TagWithFourArguments(uint8_t op_code_arg, T1 arg_arg1, T2 arg_arg2, T3 arg_arg3, T4 arg_arg4) : op_code(op_code_arg), arg1(arg_arg1), arg2(arg_arg2), arg3(arg_arg3), arg4(arg_arg4)
{
};
inline void
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(this), sizeof(*this));
instruction_number++;
};
};
class FLDZ_ : public TagWithoutArgument
{
public:
@ -578,6 +616,36 @@ public:
};
};
class FSTPG3_ : public TagWithFourArguments<unsigned int, unsigned int, int, unsigned int>
{
public:
inline FSTPG3_() : TagWithFourArguments<unsigned int, unsigned int, int, unsigned int>::TagWithFourArguments(FSTPG3, 0, 0, 0, 0)
{
};
inline FSTPG3_(const unsigned int pos_arg1, const unsigned int pos_arg2, const int pos_arg3, const unsigned int pos_arg4) : TagWithFourArguments<unsigned int, unsigned int, int, unsigned int>::TagWithFourArguments(FSTPG3, pos_arg1, pos_arg2, pos_arg3, pos_arg4)
{
};
inline unsigned int
get_row()
{
return arg1;
};
inline unsigned int
get_col()
{
return arg2;
};
inline int
get_lag()
{
return arg2;
};
inline unsigned int
get_col_pos()
{
return arg4;
};
};
class FUNARY_ : public TagWithOneArgument<uint8_t>
{
@ -643,6 +711,38 @@ public:
};
};
class FJMPIFEVAL_ : public TagWithOneArgument<unsigned int>
{
public:
inline FJMPIFEVAL_() : TagWithOneArgument<unsigned int>::TagWithOneArgument(FJMPIFEVAL)
{
};
inline FJMPIFEVAL_(unsigned int arg_pos) : TagWithOneArgument<unsigned int>::TagWithOneArgument(FJMPIFEVAL, arg_pos)
{
};
inline unsigned int
get_pos()
{
return arg1;
}
};
class FJMP_ : public TagWithOneArgument<unsigned int>
{
public:
inline FJMP_() : TagWithOneArgument<unsigned int>::TagWithOneArgument(FJMP)
{
};
inline FJMP_(unsigned int arg_pos) : TagWithOneArgument<unsigned int>::TagWithOneArgument(FJMP, arg_pos)
{
};
inline unsigned int
get_pos()
{
return arg1;
}
};
class FLDVS_ : public TagWithTwoArguments<uint8_t, unsigned int>
{
public:
@ -873,9 +973,10 @@ public:
return lag3;
};
inline void
write(ostream &CompileCode)
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(this), sizeof(FNUMEXPR_));
instruction_number++;
};
};
@ -887,27 +988,52 @@ private:
uint8_t type;
vector<int> variable;
vector<int> equation;
vector<unsigned int> other_endogenous;
vector<unsigned int> exogenous;
vector<unsigned int> det_exogenous;
bool is_linear;
vector<Block_contain_type> Block_Contain_;
int endo_nbr;
int Max_Lag;
int Max_Lead;
int u_count_int;
int nb_col_jacob;
unsigned int det_exo_size, exo_size, other_endo_size;
unsigned int nb_col_other_endo_jacob;
public:
inline FBEGINBLOCK_()
{
op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; /*variable = NULL; equation = NULL;*/
is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0;
is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0; nb_col_jacob = 0;
};
inline FBEGINBLOCK_(unsigned int size_arg, BlockSimulationType type_arg, int unsigned first_element, int unsigned block_size,
const vector<int> &variable_arg, const vector<int> &equation_arg,
bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_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,
unsigned int det_exo_size_arg, unsigned int exo_size_arg, unsigned int other_endo_size_arg, unsigned int nb_col_other_endo_jacob_arg,
const vector<unsigned int> &det_exogenous_arg, const vector<unsigned int> &exogenous_arg, const vector<unsigned int> &other_endogenous_arg)
{
op_code = FBEGINBLOCK; size = size_arg; type = type_arg;
variable = vector<int>(variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size));
equation = vector<int>(equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size));
is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg; /*Block_Contain.clear();*/
det_exogenous = vector<unsigned int>(det_exogenous_arg);
exogenous = vector<unsigned int>(exogenous_arg);
other_endogenous = vector<unsigned int>(other_endogenous_arg);
is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg;
nb_col_jacob = nb_col_jacob_arg; det_exo_size = det_exo_size_arg; exo_size = exo_size_arg; other_endo_size = other_endo_size_arg;
nb_col_other_endo_jacob = nb_col_other_endo_jacob_arg;
};
inline FBEGINBLOCK_(unsigned int size_arg, BlockSimulationType type_arg, int unsigned first_element, int unsigned block_size,
const vector<int> &variable_arg, const vector<int> &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)
{
op_code = FBEGINBLOCK; size = size_arg; type = type_arg;
variable = vector<int>(variable_arg.begin()+first_element, variable_arg.begin()+(first_element+block_size));
equation = vector<int>(equation_arg.begin()+first_element, equation_arg.begin()+(first_element+block_size));
is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg;
nb_col_jacob = nb_col_jacob_arg;
det_exo_size = 0; exo_size = 0; other_endo_size = 0;
nb_col_other_endo_jacob = 0;
}
inline unsigned int
get_size()
{
@ -948,8 +1074,33 @@ public:
{
return Block_Contain_;
};
inline int
get_nb_col_jacob()
{
return nb_col_jacob;
};
inline unsigned int
get_exo_size()
{
return exo_size;
};
inline unsigned int
get_det_exo_size()
{
return det_exo_size;
};
inline unsigned int
get_other_endo_size()
{
return other_endo_size;
};
inline unsigned int
get_nb_col_other_endo_jacob()
{
return nb_col_other_endo_jacob;
};
inline void
write(ostream &CompileCode)
write(ostream &CompileCode, unsigned int &instruction_number)
{
CompileCode.write(reinterpret_cast<char *>(&op_code), sizeof(op_code));
CompileCode.write(reinterpret_cast<char *>(&size), sizeof(size));
@ -968,6 +1119,19 @@ public:
CompileCode.write(reinterpret_cast<char *>(&Max_Lead), sizeof(Max_Lead));
CompileCode.write(reinterpret_cast<char *>(&u_count_int), sizeof(u_count_int));
}
CompileCode.write(reinterpret_cast<char *>(&nb_col_jacob), sizeof(nb_col_jacob));
CompileCode.write(reinterpret_cast<char *>(&det_exo_size), sizeof(det_exo_size));
CompileCode.write(reinterpret_cast<char *>(&exo_size), sizeof(exo_size));
CompileCode.write(reinterpret_cast<char *>(&other_endo_size), sizeof(other_endo_size));
CompileCode.write(reinterpret_cast<char *>(&nb_col_other_endo_jacob), sizeof(nb_col_other_endo_jacob));
for (unsigned int i = 0; i < det_exo_size; i++)
CompileCode.write(reinterpret_cast<char *>(&det_exogenous[i]), sizeof(det_exogenous[0]));
for (unsigned int i = 0; i < exo_size; i++)
CompileCode.write(reinterpret_cast<char *>(&exogenous[i]), sizeof(exogenous[0]));
for (unsigned int i = 0; i < other_endo_size; i++)
CompileCode.write(reinterpret_cast<char *>(&other_endogenous[i]), sizeof(other_endogenous[0]));
instruction_number++;
};
#ifdef BYTE_CODE
@ -993,6 +1157,30 @@ public:
memcpy(&Max_Lead, code, sizeof(Max_Lead)); code += sizeof(Max_Lead);
memcpy(&u_count_int, code, sizeof(u_count_int)); code += sizeof(u_count_int);
}
memcpy(&nb_col_jacob, code, sizeof(nb_col_jacob)); code += sizeof(nb_col_jacob);
memcpy(&det_exo_size, code, sizeof(det_exo_size)); code += sizeof(det_exo_size);
memcpy(&exo_size, code, sizeof(exo_size)); code += sizeof(exo_size);
memcpy(&other_endo_size, code, sizeof(other_endo_size)); code += sizeof(other_endo_size);
memcpy(&nb_col_other_endo_jacob, code, sizeof(nb_col_other_endo_jacob)); code += sizeof(nb_col_other_endo_jacob);
for (unsigned int i = 0; i < det_exo_size; i++)
{
unsigned int tmp_i;
memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
det_exogenous.push_back(tmp_i);
}
for (unsigned int i = 0; i < exo_size; i++)
{
unsigned int tmp_i;
memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
exogenous.push_back(tmp_i);
}
for (unsigned int i = 0; i < other_endo_size; i++)
{
unsigned int tmp_i;
memcpy(&tmp_i, code, sizeof(tmp_i)); code += sizeof(tmp_i);
other_endogenous.push_back(tmp_i);
}
return code;
};
#endif
@ -1036,6 +1224,7 @@ public:
CompiledCode.close();
nb_blocks = 0;
bool done = false;
int instruction = 0;
while (!done)
{
switch (*code)
@ -1181,6 +1370,20 @@ public:
tags_liste.push_back(make_pair(FSTPG, code));
code += sizeof(FSTPG_);
break;
case FSTPG2:
# ifdef DEBUGL
mexPrintf("FSTPG2\n");
# endif
tags_liste.push_back(make_pair(FSTPG2, code));
code += sizeof(FSTPG2_);
break;
case FSTPG3:
# ifdef DEBUGL
mexPrintf("FSTPG3\n");
# endif
tags_liste.push_back(make_pair(FSTPG3, code));
code += sizeof(FSTPG3_);
break;
case FUNARY:
# ifdef DEBUGL
mexPrintf("FUNARY\n");
@ -1257,10 +1460,25 @@ public:
nb_blocks++;
}
break;
case FJMPIFEVAL:
# ifdef DEBUGL
mexPrintf("FJMPIFEVAL\n");
# endif
tags_liste.push_back(make_pair(FJMPIFEVAL, code));
code += sizeof(FJMPIFEVAL_);
break;
case FJMP:
# ifdef DEBUGL
mexPrintf("FJMP\n");
# endif
tags_liste.push_back(make_pair(FJMP, code));
code += sizeof(FJMP_);
break;
default:
mexPrintf("Unknown Tag value=%d code=%x\n", *code, code);
done = true;
}
instruction++;
}
return tags_liste;
};

601
preprocessor/DynamicModel.cc
View File

@ -57,28 +57,28 @@ DynamicModel::AddVariable(int symb_id, int lag)
}
void
DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
DynamicModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const
{
first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, symb_id), lag)));
if (it != first_derivatives.end())
(it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
else
{
FLDZ_ fldz;
fldz.write(code_file);
fldz.write(code_file, instruction_number);
}
}
void
DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, const map_idx_t &map_idx) const
DynamicModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, const map_idx_t &map_idx) const
{
map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
if (it != first_chain_rule_derivatives.end())
(it->second)->compile(code_file, false, temporary_terms, map_idx, true, false);
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
else
{
FLDZ_ fldz;
fldz.write(code_file);
fldz.write(code_file, instruction_number);
}
}
@ -140,7 +140,6 @@ DynamicModel::computeTemporaryTermsOrdered()
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
set<int> temporary_terms_in_use;
temporary_terms_in_use.clear();
v_temporary_terms_inuse[block] = temporary_terms_in_use;
@ -471,7 +470,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
goto evaluation;
feedback_variables.push_back(variable_ID);
output << " % equation " << equation_ID+1 << " variable : " << sModel
<< " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
<< " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
output << " " << "residual(" << i+1-block_recursive << ") = (";
goto end;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
@ -480,7 +479,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
goto evaluation;
feedback_variables.push_back(variable_ID);
output << " % equation " << equation_ID+1 << " variable : " << sModel
<< " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << endl;
<< " (" << variable_ID+1 << ") " << c_Equation_Type(equ_type) << " symb_id=" << symbol_table.getID(eEndogenous, variable_ID) << endl;
Uf[equation_ID] << " b(" << i+1-block_recursive << "+Per_J_) = -residual(" << i+1-block_recursive << ", it_)";
output << " residual(" << i+1-block_recursive << ", it_) = (";
goto end;
@ -526,6 +525,34 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.first;
int var = it->first.second.second;
int eqr = getBlockInitialEquationID(block, eq);
expr_t id = it->second;
output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
<< "(" << lag
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.first;
int var = it->first.second.second;
int eqr = getBlockInitialEquationID(block, eq);
expr_t id = it->second;
output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
<< "(" << lag
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
int lag = it->first.first;
@ -564,7 +591,36 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.first;
int var = it->first.second.second;
int eqr = getBlockInitialEquationID(block, eq);
expr_t id = it->second;
output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
<< "(" << lag
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.first;
int var = it->first.second.second;
int eqr = getBlockInitialEquationID(block, eq);
expr_t id = it->second;
output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
<< "(" << lag
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
int lag = it->first.first;
@ -695,6 +751,36 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != derivative_exo[block].end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.first;
int var = it->first.second.second;
int eqr = getBlockInitialEquationID(block, eq);
expr_t id = it->second;
output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
<< "(" << lag
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != derivative_exo_det[block].end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.first;
int var = it->first.second.second;
int eqr = getBlockInitialEquationID(block, eq);
expr_t id = it->second;
output << " g1_x(" << eqr+1 << ", " << var+1+(lag+block_max_lag)*block_size << ") = ";
id->writeOutput(output, local_output_type, local_temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eExogenous, var))
<< "(" << lag
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
}
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != derivative_other_endo[block].end(); it++)
{
int lag = it->first.first;
@ -726,6 +812,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
{
ostringstream tmp_output;
ofstream code_file;
unsigned int instruction_number = 0;
bool file_open = false;
string main_name = file_name;
@ -737,7 +824,6 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
exit(EXIT_FAILURE);
}
int count_u;
int u_count_int = 0;
BlockSimulationType simulation_type;
@ -753,8 +839,15 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
fdimt.write(code_file);
fdimt.write(code_file, instruction_number);
int other_endo_size = 0;
vector<unsigned int> exo, exo_det, other_endo;
for(int i = 0; i < symbol_table.exo_det_nbr(); i++)
exo_det.push_back(i);
for(int i = 0; i < symbol_table.exo_nbr(); i++)
exo.push_back(i);
FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
simulation_type,
0,
@ -765,16 +858,24 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
symbol_table.endo_nbr(),
0,
0,
u_count_int
u_count_int,
0,
symbol_table.exo_det_nbr(),
symbol_table.exo_nbr(),
other_endo_size,
0,
exo_det,
exo,
other_endo
);
fbeginblock.write(code_file);
fbeginblock.write(code_file, instruction_number);
compileTemporaryTerms(code_file, temporary_terms, map_idx, true, false);
compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false);
compileModelEquations(code_file, temporary_terms, map_idx, true, false);
compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, true, false);
FENDEQU_ fendequ;
fendequ.write(code_file);
fendequ.write(code_file, instruction_number);
vector<vector<pair<pair<int, int>, int > > > derivatives;
derivatives.resize(symbol_table.endo_nbr());
count_u = symbol_table.endo_nbr();
@ -790,45 +891,45 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
unsigned int var = symbol_table.getTypeSpecificID(symb);
int lag = getLagByDerivID(deriv_id);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var, lag);
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
if (!derivatives[eq].size())
derivatives[eq].clear();
derivatives[eq].push_back(make_pair(make_pair(var, lag), count_u));
d1->compile(code_file, false, temporary_terms, map_idx, true, false);
d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FSTPU_ fstpu(count_u);
fstpu.write(code_file);
fstpu.write(code_file, instruction_number);
count_u++;
}
}
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
FLDR_ fldr(i);
fldr.write(code_file);
fldr.write(code_file, instruction_number);
for(vector<pair<pair<int, int>, int> >::const_iterator it = derivatives[i].begin();
it != derivatives[i].end(); it++)
{
FLDU_ fldu(it->second);
fldu.write(code_file);
fldu.write(code_file, instruction_number);
FLDV_ fldv(eEndogenous, it->first.first, it->first.second);
fldv.write(code_file);
fldv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
if (it != derivatives[i].begin())
{
FBINARY_ fbinary(oPlus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
}
}
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPU_ fstpu(i);
fstpu.write(code_file);
fstpu.write(code_file, instruction_number);
}
FENDBLOCK_ fendblock;
fendblock.write(code_file);
fendblock.write(code_file, instruction_number);
FEND_ fend;
fend.write(code_file);
fend.write(code_file, instruction_number);
code_file.close();
}
@ -852,6 +953,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
unsigned int instruction_number = 0;
expr_t lhs = NULL, rhs = NULL;
BinaryOpNode *eq_node;
Uff Uf[symbol_table.endo_nbr()];
@ -870,7 +972,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
fdimt.write(code_file);
fdimt.write(code_file, instruction_number);
for (unsigned int block = 0; block < getNbBlocks(); block++)
{
@ -878,7 +980,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
if (block > 0)
{
FENDBLOCK_ fendblock;
fendblock.write(code_file);
fendblock.write(code_file, instruction_number);
}
int count_u;
int u_count_int = 0;
@ -886,6 +988,8 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
unsigned int block_size = getBlockSize(block);
unsigned int block_mfs = getBlockMfs(block);
unsigned int block_recursive = block_size - block_mfs;
unsigned int block_exo_det_size = exo_det_block[block].size();
unsigned int block_other_endo_size = other_endo_block[block].size();
int block_max_lag = max_leadlag_block[block].first;
if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE
@ -895,6 +999,45 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE);
file_open = true;
}
map<pair<int, pair<int, int> >, expr_t> tmp_block_endo_derivative;
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
tmp_block_endo_derivative[make_pair(it->second.first, make_pair(it->first.second, it->first.first) )] = it->second.second ;
map<pair<int, pair<int, int> >, expr_t> tmp_exo_derivative;
for (derivative_t::const_iterator it = derivative_exo[block].begin(); it != (derivative_exo[block]).end(); it++)
tmp_exo_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first) )] = it->second ;
map<pair<int, pair<int, int> >, expr_t> tmp_exo_det_derivative;
for (derivative_t::const_iterator it = derivative_exo_det[block].begin(); it != (derivative_exo_det[block]).end(); it++)
tmp_exo_det_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first) )] = it->second;
map<pair<int, pair<int, int> >, expr_t> tmp_other_endo_derivative;
for (derivative_t::const_iterator it = derivative_other_endo[block].begin(); it != (derivative_other_endo[block]).end(); it++)
tmp_other_endo_derivative[make_pair(it->first.first, make_pair(it->first.second.second, it->first.second.first) )] = it->second;
int prev_var = -1;
int prev_lag = -999999999;
int count_col_endo = 0;
for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
{
int lag = it->first.first;
int var = it->first.second.first;
if(prev_var != var || prev_lag != lag)
{
prev_var = var;
prev_lag = lag;
count_col_endo++;
}
}
vector<unsigned int> exo_det;
for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++)
exo_det.push_back(it->first);
vector<unsigned int> exo;
for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
exo.push_back(it->first);
vector<unsigned int> other_endo;
unsigned int count_col_other_endo = 0;
for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
{
other_endo.push_back(it->first);
count_col_other_endo += it->second.size();
}
FBEGINBLOCK_ fbeginblock(block_mfs,
simulation_type,
getBlockFirstEquation(block),
@ -905,9 +1048,17 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
symbol_table.endo_nbr(),
block_max_lag,
block_max_lag,
u_count_int
u_count_int,
count_col_endo,
block_exo_det_size,
getBlockExoColSize(block),
block_other_endo_size,
count_col_other_endo,
exo_det,
exo,
other_endo
);
fbeginblock.write(code_file);
fbeginblock.write(code_file, instruction_number);
// The equations
for (i = 0; i < (int) block_size; i++)
@ -921,14 +1072,15 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
it != v_temporary_terms[block][i].end(); it++)
{
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
fnumexpr.write(code_file);
(*it)->compile(code_file, false, tt2, map_idx, true, false);
fnumexpr.write(code_file, instruction_number);
(*it)->compile(code_file, instruction_number, false, tt2, map_idx, true, false);
FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
fstpt.write(code_file);
fstpt.write(code_file, instruction_number);
// Insert current node into tt2
tt2.insert(*it);
#ifdef DEBUGC
cout << "FSTPT " << v << "\n";
instruction_number++;
code_file.write(&FOK, sizeof(FOK));
code_file.write(reinterpret_cast<char *>(&k), sizeof(k));
ki++;
@ -956,23 +1108,23 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
equ_type = getBlockEquationType(block, i);
{
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
}
if (equ_type == E_EVALUATE)
{
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
}
else if (equ_type == E_EVALUATE_S)
{
eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
lhs->compile(code_file, true, temporary_terms, map_idx, true, false);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, true, false);
}
break;
case SOLVE_BACKWARD_COMPLETE:
@ -989,22 +1141,28 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
default:
end:
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
lhs->compile(code_file, false, temporary_terms, map_idx, true, false);
rhs->compile(code_file, false, temporary_terms, map_idx, true, false);
lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(i - block_recursive);
fstpr.write(code_file);
fstpr.write(code_file, instruction_number);
}
}
FENDEQU_ fendequ;
fendequ.write(code_file);
// The Jacobian if we have to solve the block
fendequ.write(code_file, instruction_number);
// Get the current code_file position and jump if eval = true
streampos pos1 = code_file.tellp();
FJMPIFEVAL_ fjmp_if_eval(0);
fjmp_if_eval.write(code_file, instruction_number);
int prev_instruction_number = instruction_number;
// The Jacobian if we have to solve the block determinsitic bloc
if (simulation_type != EVALUATE_BACKWARD
&& simulation_type != EVALUATE_FORWARD)
{
@ -1013,13 +1171,13 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
case SOLVE_BACKWARD_SIMPLE:
case SOLVE_FORWARD_SIMPLE:
{
FNUMEXPR_ fnumexpr(FirstEndoDerivative, getBlockEquationID(block, i), getBlockVariableID(block, 0), 0);
fnumexpr.write(code_file);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0);
fnumexpr.write(code_file, instruction_number);
}
compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
{
FSTPG_ fstpg(0);
fstpg.write(code_file);
fstpg.write(code_file, instruction_number);
}
break;
@ -1030,12 +1188,12 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
count_u = feedback_variables.size();
for (block_derivatives_equation_variable_laglead_nodeid_t::const_iterator it = blocks_derivatives[block].begin(); it != (blocks_derivatives[block]).end(); it++)
{
int lag = it->second.first;
unsigned int eq = it->first.first;
unsigned int var = it->first.second;
unsigned int eqr = getBlockEquationID(block, eq);
unsigned int varr = getBlockVariableID(block, var);
int lag = it->second.first;
if (eq >= block_recursive && var >= block_recursive)
if (eq >= block_recursive and var >= block_recursive)
{
if (!Uf[eqr].Ufl)
{
@ -1052,10 +1210,10 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
Uf[eqr].Ufl->var = varr;
Uf[eqr].Ufl->lag = lag;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
fnumexpr.write(code_file);
compileChainRuleDerivative(code_file, eqr, varr, lag, map_idx);
fnumexpr.write(code_file, instruction_number);
compileChainRuleDerivative(code_file, instruction_number, eqr, varr, lag, map_idx);
FSTPU_ fstpu(count_u);
fstpu.write(code_file);
fstpu.write(code_file, instruction_number);
count_u++;
}
}
@ -1064,24 +1222,24 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
if (i >= (int) block_recursive)
{
FLDR_ fldr(i-block_recursive);
fldr.write(code_file);
fldr.write(code_file, instruction_number);
FLDZ_ fldz;
fldz.write(code_file);
fldz.write(code_file, instruction_number);
v = getBlockEquationID(block, i);
for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
{
FLDU_ fldu(Uf[v].Ufl->u);
fldu.write(code_file);
fldu.write(code_file, instruction_number);
FLDV_ fldv(eEndogenous, Uf[v].Ufl->var, Uf[v].Ufl->lag);
fldv.write(code_file);
fldv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FCUML_ fcuml;
fcuml.write(code_file);
fcuml.write(code_file, instruction_number);
}
Uf[v].Ufl = Uf[v].Ufl_First;
while (Uf[v].Ufl)
@ -1091,10 +1249,10 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
Uf[v].Ufl = Uf[v].Ufl_First;
}
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPU_ fstpu(i - block_recursive);
fstpu.write(code_file);
fstpu.write(code_file, instruction_number);
}
}
break;
@ -1102,11 +1260,125 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
break;
}
}
// Get the current code_file position and jump = true
streampos pos2 = code_file.tellp();
FJMP_ fjmp(0);
fjmp.write(code_file, instruction_number);
// Set code_file position to previous JMPIFEVAL_ and set the number of instructions to jump
streampos pos3 = code_file.tellp();
code_file.seekp(pos1);
FJMPIFEVAL_ fjmp_if_eval1(instruction_number - prev_instruction_number);
fjmp_if_eval1.write(code_file, instruction_number);
code_file.seekp(pos3);
prev_instruction_number = instruction_number ;
// The Jacobian if we have to solve the block determinsitic bloc
prev_var = -1;
prev_lag = -999999999;
count_col_endo = 0;
for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_block_endo_derivative.begin(); it != tmp_block_endo_derivative.end(); it++)
{
int lag = it->first.first;
unsigned int eq = it->first.second.second;
int var = it->first.second.first;
unsigned int eqr = getBlockEquationID(block, eq);
unsigned int varr = getBlockVariableID(block, var);
if(prev_var != var || prev_lag != lag)
{
prev_var = var;
prev_lag = lag;
count_col_endo++;
}
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr, lag);
fnumexpr.write(code_file, instruction_number);
compileDerivative(code_file, instruction_number, eqr, varr, lag, map_idx);
FSTPG3_ fstpg3(eq, var, lag, count_col_endo-1);
fstpg3.write(code_file, instruction_number);
}
prev_var = -1;
prev_lag = -999999999;
int count_col_exo = 0;
for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_exo_derivative.begin(); it != tmp_exo_derivative.end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.second;
int var = it->first.second.first;
int eqr = getBlockInitialEquationID(block, eq);
int varr = getBlockInitialExogenousID(block, var);
if(prev_var != var || prev_lag != lag)
{
prev_var = var;
prev_lag = lag;
count_col_exo++;
}
expr_t id = it->second;
FNUMEXPR_ fnumexpr(FirstExoDerivative, eqr, varr, lag);
fnumexpr.write(code_file, instruction_number);
id->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FSTPG3_ fstpg3(eq, var, lag, /*var*/count_col_exo-1);
fstpg3.write(code_file, instruction_number);
}
prev_var = -1;
prev_lag = -999999999;
int count_col_exo_det = 0;
for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_exo_det_derivative.begin(); it != tmp_exo_det_derivative.end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.second;
int var = it->first.second.first;
int eqr = getBlockInitialEquationID(block, eq);
int varr = getBlockInitialDetExogenousID(block, var);
if(prev_var != var || prev_lag != lag)
{
prev_var = var;
prev_lag = lag;
count_col_exo_det++;
}
expr_t id = it->second;
FNUMEXPR_ fnumexpr(FirstExodetDerivative, eqr, varr, lag);
fnumexpr.write(code_file, instruction_number);
id->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FSTPG3_ fstpg3(eq, var, lag, count_col_exo_det-1);
fstpg3.write(code_file, instruction_number);
}
prev_var = -1;
prev_lag = -999999999;
count_col_other_endo = 0;
for (map<pair<int, pair<int, int> >, expr_t>::const_iterator it = tmp_other_endo_derivative.begin(); it != tmp_other_endo_derivative.end(); it++)
{
int lag = it->first.first;
int eq = it->first.second.second;
int var = it->first.second.first;
int eqr = getBlockInitialEquationID(block, eq);
int varr = getBlockInitialOtherEndogenousID(block, var);;
if(prev_var != var || prev_lag != lag)
{
prev_var = var;
prev_lag = lag;
count_col_other_endo++;
}
expr_t id = it->second;
FNUMEXPR_ fnumexpr(FirstOtherEndoDerivative, eqr, varr, lag);
fnumexpr.write(code_file, instruction_number);
id->compile(code_file, instruction_number, false, temporary_terms, map_idx, true, false);
FSTPG3_ fstpg3(eq, var, lag, count_col_other_endo-1);
fstpg3.write(code_file, instruction_number);
}
// Set codefile position to previous JMP_ and set the number of instructions to jump
pos1 = code_file.tellp();
code_file.seekp(pos2);
FJMP_ fjmp1(instruction_number - prev_instruction_number);
fjmp1.write(code_file, instruction_number);
code_file.seekp(pos1);
}
FENDBLOCK_ fendblock;
fendblock.write(code_file);
fendblock.write(code_file, instruction_number);
FEND_ fend;
fend.write(code_file);
fend.write(code_file, instruction_number);
code_file.close();
}
@ -1909,7 +2181,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
if (block_decomposition)
{
int count_lead_lag_incidence = 0;
int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo;
int max_lead, max_lag, max_lag_endo, max_lead_endo, max_lag_exo, max_lead_exo, max_lag_exo_det, max_lead_exo_det;
unsigned int nb_blocks = getNbBlocks();
for (unsigned int block = 0; block < nb_blocks; block++)
{
@ -1921,11 +2193,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
max_lead = max_leadlag_block[block].second;
max_lag_endo = endo_max_leadlag_block[block].first;
max_lead_endo = endo_max_leadlag_block[block].second;
max_lag_exo = max(exo_max_leadlag_block[block].first, exo_det_max_leadlag_block[block].first);
max_lead_exo = max(exo_max_leadlag_block[block].second, exo_det_max_leadlag_block[block].second);
vector<int> exogenous(symbol_table.exo_nbr(), -1);
vector<int>::iterator it_exogenous;
exogenous.clear();
max_lag_exo = exo_max_leadlag_block[block].first;
max_lead_exo = exo_max_leadlag_block[block].second;
max_lag_exo_det = exo_det_max_leadlag_block[block].first;
max_lead_exo_det = exo_det_max_leadlag_block[block].second;
ostringstream tmp_s, tmp_s_eq;
tmp_s.str("");
tmp_s_eq.str("");
@ -1934,10 +2205,21 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
tmp_s << " " << getBlockVariableID(block, i)+1;
tmp_s_eq << " " << getBlockEquationID(block, i)+1;
}
it_exogenous = exogenous.begin();
set<int> exogenous;
exogenous.clear();
for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
it_exogenous = set_union(it->second.begin(), it->second.end(), exogenous.begin(), exogenous.end(), it_exogenous);
output << "M_.block_structure.block(" << block+1 << ").num = " << block+1 << ";\n";
for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
exogenous.insert(*it1);
set<int> exogenous_det;
exogenous_det.clear();
for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++)
for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
exogenous_det.insert(*it1);
set<int> other_endogenous;
other_endogenous.clear();
for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++)
for(var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
other_endogenous.insert(*it1);
output << "M_.block_structure.block(" << block+1 << ").Simulation_Type = " << simulation_type << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_lag = " << max_lag << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_lead = " << max_lead << ";\n";
@ -1945,21 +2227,43 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
output << "M_.block_structure.block(" << block+1 << ").maximum_endo_lead = " << max_lead_endo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lag = " << max_lag_exo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_lead = " << max_lead_exo << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_det_lag = " << max_lag_exo_det << ";\n";
output << "M_.block_structure.block(" << block+1 << ").maximum_exo_det_lead = " << max_lead_exo_det << ";\n";
output << "M_.block_structure.block(" << block+1 << ").endo_nbr = " << block_size << ";\n";
output << "M_.block_structure.block(" << block+1 << ").mfs = " << getBlockMfs(block) << ";\n";
output << "M_.block_structure.block(" << block+1 << ").equation = [" << tmp_s_eq.str() << "];\n";
output << "M_.block_structure.block(" << block+1 << ").variable = [" << tmp_s.str() << "];\n";
output << "M_.block_structure.block(" << block+1 << ").exo_nbr = " << getBlockExoSize(block) << ";\n";
output << "M_.block_structure.block(" << block+1 << ").exogenous = [";
int i = 0;
for (it_exogenous = exogenous.begin(); it_exogenous != exogenous.end(); it_exogenous++)
for (set<int>::iterator it_exogenous = exogenous.begin(); it_exogenous != exogenous.end(); it_exogenous++)
if (*it_exogenous >= 0)
{
output << " " << *it_exogenous+1;
i++;
}
output << "];\n";
output << "M_.block_structure.block(" << block+1 << ").exo_nbr = " << i << ";\n";
output << "M_.block_structure.block(" << block+1 << ").exo_det_nbr = " << i << ";\n";
output << "M_.block_structure.block(" << block+1 << ").exogenous_det = [";
i = 0;
for (set<int>::iterator it_exogenous_det = exogenous_det.begin(); it_exogenous_det != exogenous_det.end(); it_exogenous_det++)
if (*it_exogenous_det >= 0)
{
output << " " << *it_exogenous_det+1;
i++;
}
output << "];\n";
output << "M_.block_structure.block(" << block+1 << ").exo_det_nbr = " << exo_det_block.size() << ";\n";
output << "M_.block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";\n";
output << "M_.block_structure.block(" << block+1 << ").other_endogenous = [";
i = 0;
for (set<int>::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++)
if (*it_other_endogenous >= 0)
{
output << " " << *it_other_endogenous+1;
i++;
}
output << "];\n";
tmp_s.str("");
count_lead_lag_incidence = 0;
@ -1970,7 +2274,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
int last_var = -1;
for (int lag = -max_lag_endo; lag < max_lead_endo+1; lag++)
{
last_var = 0;
last_var = -1;
for (dynamic_jacob_map_t::const_iterator it = reordered_dynamic_jacobian.begin(); it != reordered_dynamic_jacobian.end(); it++)
{
if (lag == it->first.first && last_var != it->first.second.first)
@ -1989,6 +2293,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
output << "M_.block_structure.block(" << block+1 << ").lead_lag_incidence = [ M_.block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << "\n";
tmp_s.str("");
}
output << "M_.block_structure.block(" << block+1 << ").n_static = " << block_col_type[block].first.first << ";\n";
output << "M_.block_structure.block(" << block+1 << ").n_forward = " << block_col_type[block].first.second << ";\n";
output << "M_.block_structure.block(" << block+1 << ").n_backward = " << block_col_type[block].second.first << ";\n";
output << "M_.block_structure.block(" << block+1 << ").n_mixed = " << block_col_type[block].second.second << ";\n";
}
string cst_s;
int nb_endo = symbol_table.endo_nbr();
@ -2027,7 +2335,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").lead_lag = " << prev_lag << ";\n";
output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").sparse_IM = [";
}
output << it->first.second.first << " " << it->first.second.second << ";\n";
output << it->first.second.first+1 << " " << it->first.second.second+1 << ";\n";
}
output << "];\n";
}
@ -2142,6 +2450,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
if (block)
{
vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
jacob_map_t contemporaneous_jacobian, static_jacobian;
// for each block contains pair<Size, Feddback_variable>
@ -2159,10 +2468,11 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
cout << "Finding the optimal block decomposition of the model ...\n";
if (prologue+epilogue < (unsigned int) equation_number())
computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered);
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type);
printBlockDecomposition(blocks);
@ -2172,10 +2482,24 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
collect_block_first_order_derivatives();
collectBlockVariables();
global_temporary_terms = true;
if (!no_tmp_terms)
computeTemporaryTermsOrdered();
int k = 0;
equation_block = vector<int>(equation_number());
variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equation_number());
for (unsigned int i = 0; i < getNbBlocks(); i++)
{
for(unsigned int j = 0; j< getBlockSize(i); j++)
{
equation_block[equation_reordered[k]] = i;
int l = variable_reordered[k];
variable_block_lead_lag[l] = make_pair(i, make_pair(variable_lag_lead[l].first, variable_lag_lead[l].second));
k++;
}
}
}
else
if (!no_tmp_terms)
@ -2241,22 +2565,19 @@ DynamicModel::get_Derivatives(int block)
}
void
DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_endo_derivatives)
{
map<int, expr_t> recursive_variables;
unsigned int nb_blocks = getNbBlocks();
blocks_derivatives = blocks_derivatives_t(nb_blocks);
blocks_endo_derivatives = blocks_derivatives_t(nb_blocks);
for (unsigned int block = 0; block < nb_blocks; block++)
{
block_derivatives_equation_variable_laglead_nodeid_t tmp_derivatives;
recursive_variables.clear();
BlockSimulationType simulation_type = getBlockSimulationType(block);
int block_size = getBlockSize(block);
int block_nb_mfs = getBlockMfs(block);
int block_nb_recursives = block_size - block_nb_mfs;
if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
{
blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
blocks_endo_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == E_EVALUATE_S)
@ -2289,34 +2610,7 @@ DynamicModel::computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives)
}
tmp_derivatives.push_back(make_pair(make_pair(eq, var), make_pair(lag, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))])));
}
}
else if (simulation_type == SOLVE_BACKWARD_SIMPLE || simulation_type == SOLVE_FORWARD_SIMPLE
|| simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE)
{
blocks_derivatives.push_back(block_derivatives_equation_variable_laglead_nodeid_t(0));
for (int i = 0; i < block_nb_recursives; i++)
{
if (getBlockEquationType(block, i) == E_EVALUATE_S)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationRenormalizedExpr(block, i);
else
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, getBlockVariableID(block, i)), 0)] = getBlockEquationExpr(block, i);
}
for (int eq = block_nb_recursives; eq < block_size; eq++)
{
int eqr = getBlockEquationID(block, eq);
for (int var = block_nb_recursives; var < block_size; var++)
{
int varr = getBlockVariableID(block, var);
expr_t d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
if (d1 == Zero)
continue;
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
tmp_derivatives.push_back(
make_pair(make_pair(eq, var), make_pair(0, first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))])));
}
}
}
blocks_derivatives[block] = tmp_derivatives;
blocks_endo_derivatives[block] = tmp_derivatives;
}
}
@ -2379,6 +2673,20 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > other_endo_max_leadlag_block[block_eq].second)
other_endo_max_leadlag_block[block_eq] = make_pair(other_endo_max_leadlag_block[block_eq].first, lag);
tmp_derivative = derivative_other_endo[block_eq];
{
map< int, map<int, int> >::const_iterator it = block_other_endo_index.find(block_eq);
if (it == block_other_endo_index.end())
block_other_endo_index[block_eq][var] = 0;
else
{
map<int, int>::const_iterator it1 = it->second.find(var);
if (it1 == it->second.end())
{
int size = block_other_endo_index[block_eq].size();
block_other_endo_index[block_eq][var] = size;
}
}
}
tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eEndogenous, var), lag))];
derivative_other_endo[block_eq] = tmp_derivative;
lag_var = other_endo_block[block_eq];
@ -2394,6 +2702,20 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > exo_max_leadlag_block[block_eq].second)
exo_max_leadlag_block[block_eq] = make_pair(exo_max_leadlag_block[block_eq].first, lag);
tmp_derivative = derivative_exo[block_eq];
{
map< int, map<int, int> >::const_iterator it = block_exo_index.find(block_eq);
if (it == block_exo_index.end())
block_exo_index[block_eq][var] = 0;
else
{
map<int, int>::const_iterator it1 = it->second.find(var);
if (it1 == it->second.end())
{
int size = block_exo_index[block_eq].size();
block_exo_index[block_eq][var] = size;
}
}
}
tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eExogenous, var), lag))];
derivative_exo[block_eq] = tmp_derivative;
lag_var = exo_block[block_eq];
@ -2408,6 +2730,20 @@ DynamicModel::collect_block_first_order_derivatives()
if (lag > 0 && lag > exo_det_max_leadlag_block[block_eq].second)
exo_det_max_leadlag_block[block_eq] = make_pair(exo_det_max_leadlag_block[block_eq].first, lag);
tmp_derivative = derivative_exo_det[block_eq];
{
map< int, map<int, int> >::const_iterator it = block_det_exo_index.find(block_eq);
if (it == block_det_exo_index.end())
block_det_exo_index[block_eq][var] = 0;
else
{
map<int, int>::const_iterator it1 = it->second.find(var);
if (it1 == it->second.end())
{
int size = block_det_exo_index[block_eq].size();
block_det_exo_index[block_eq][var] = size;
}
}
}
tmp_derivative[make_pair(lag, make_pair(eq, var))] = first_derivatives[make_pair(eq, getDerivID(symbol_table.getID(eExogenous, var), lag))];
derivative_exo_det[block_eq] = tmp_derivative;
lag_var = exo_det_block[block_eq];
@ -2427,6 +2763,35 @@ DynamicModel::collect_block_first_order_derivatives()
}
void DynamicModel::collectBlockVariables()
{
for (unsigned int block = 0; block < getNbBlocks(); block++)
{
int prev_var = -1;
int prev_lag = -999999999;
int count_col_exo = 0;
var_t tmp_var_exo;
for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++)
{
int lag = it->first;
for(var_t::const_iterator it2 = it->second.begin(); it2 != it->second.end(); it2++)
{
int var = *it2;
tmp_var_exo.insert(var);
if(prev_var != var || prev_lag != lag)
{
prev_var = var;
prev_lag = lag;
count_col_exo++;
}
}
}
block_var_exo.push_back(make_pair(tmp_var_exo, count_col_exo));
}
}
void
DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const
{

98
preprocessor/DynamicModel.hh
View File

@ -142,9 +142,9 @@ private:
//! creates a mapping from the index of temporary terms to a natural index
void computeTemporaryTermsMapping();
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, const map_idx_t &map_idx) const;
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, const map_idx_t &map_idx) const;
//! Get the type corresponding to a derivation ID
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
@ -180,6 +180,9 @@ private:
//! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
void collect_block_first_order_derivatives();
//! Collecte the informations about exogenous, deterministic exogenous and endogenous from the previous block for each block
void collectBlockVariables();
//! Factorized code for substitutions of leads/lags
/*! \param[in] type determines which type of variables is concerned */
void substituteLeadLagInternal(aux_var_t type);
@ -211,11 +214,26 @@ private:
//! Vector of derivative for each blocks
vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
//!List for each block and for each lag-leag all the other endogenous variables and exogenous variables
//!List for each block and for each lag-lead all the other endogenous variables and exogenous variables
typedef set<int> var_t;
typedef map<int, var_t> lag_var_t;
vector<lag_var_t> other_endo_block, exo_block, exo_det_block;
//!List for each block the exogenous variables
vector<pair<var_t, int> > block_var_exo;
map< int, map<int , int> > block_exo_index, block_det_exo_index, block_other_endo_index;
//! for each block described the number of static, forward, backward and mixed variables in the block
/*! pair< pair<static, forward>, pair<backward,mixed> > */
vector<pair< pair<int, int>, pair<int,int> > > block_col_type;
//! List for each variable its block number and its maximum lag and lead inside the block
vector<pair<int, pair<int, int> > > variable_block_lead_lag;
//! List for each equation its block number
vector<int> equation_block;
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
@ -323,6 +341,18 @@ public:
{
return (block_type_firstequation_size_mfs[block_number].second.first);
};
//! Return the number of exogenous variable in the block block_number
virtual unsigned int
getBlockExoSize(int block_number) const
{
return (block_var_exo[block_number].first.size());
};
//! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
virtual unsigned int
getBlockExoColSize(int block_number) const
{
return (block_var_exo[block_number].second);
};
//! Return the number of feedback variable of the block block_number
virtual unsigned int
getBlockMfs(int block_number) const
@ -377,6 +407,13 @@ public:
{
return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
};
//! Return the original number of the exogenous variable varexo_number belonging to the block block_number
virtual int
getBlockVariableExoID(int block_number, int variable_number) const
{
map<int, var_t>::const_iterator it = exo_block[block_number].find(variable_number);
return (it->first);
};
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int
getBlockInitialEquationID(int block_number, int equation_number) const
@ -389,7 +426,60 @@ public:
{
return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
};
//! Return the block number containing the endogenous variable variable_number
int
getBlockVariableID(int variable_number) const
{
return(variable_block_lead_lag[variable_number].first);
};
//! Return the position of the exogenous variable_number in the block number belonging to the block block_number
virtual int
getBlockInitialExogenousID(int block_number, int variable_number) const
{
map< int, map<int, int> >::const_iterator it = block_exo_index.find(block_number);
if (it != block_exo_index.end())
{
map<int, int>::const_iterator it1 = it->second.find(variable_number);
if( it1 != it->second.end())
return it1->second;
else
return -1;
}
else
return (-1);
};
//! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
virtual int
getBlockInitialDetExogenousID(int block_number, int variable_number) const
{
map< int, map<int, int> >::const_iterator it = block_det_exo_index.find(block_number);
if (it != block_det_exo_index.end())
{
map<int, int>::const_iterator it1 = it->second.find(variable_number);
if( it1 != it->second.end())
return it1->second;
else
return -1;
}
else
return (-1);
};
//! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
virtual int
getBlockInitialOtherEndogenousID(int block_number, int variable_number) const
{
map< int, map<int, int> >::const_iterator it = block_other_endo_index.find(block_number);
if (it != block_other_endo_index.end())
{
map<int, int>::const_iterator it1 = it->second.find(variable_number);
if( it1 != it->second.end())
return it1->second;
else
return -1;
}
else
return (-1);
};
};
#endif

62
preprocessor/ExprNode.cc
View File

@ -300,10 +300,10 @@ NumConstNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
}
void
NumConstNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
NumConstNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
FLDC_ fldc(datatree.num_constants.getDouble(id));
fldc.write(CompileCode);
fldc.write(CompileCode, instruction_number);
}
void
@ -689,10 +689,10 @@ VariableNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
}
void
VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
VariableNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
if (type == eModelLocalVariable || type == eModFileLocalVariable)
datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
datatree.local_variables_table[symb_id]->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
else
{
int tsid = datatree.symbol_table.getTypeSpecificID(symb_id);
@ -705,26 +705,26 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
if (steady_dynamic) // steady state values in a dynamic model
{
FLDVS_ fldvs(type, tsid);
fldvs.write(CompileCode);
fldvs.write(CompileCode, instruction_number);
}
else
{
if (type == eParameter)
{
FLDV_ fldv(type, tsid);
fldv.write(CompileCode);
fldv.write(CompileCode, instruction_number);
}
else
{
FLDV_ fldv(type, tsid, lag);
fldv.write(CompileCode);
fldv.write(CompileCode, instruction_number);
}
}
}
else
{
FLDSV_ fldsv(type, tsid);
fldsv.write(CompileCode);
fldsv.write(CompileCode, instruction_number);
}
}
else
@ -741,19 +741,19 @@ VariableNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
if (type == eParameter)
{
FSTPV_ fstpv(type, tsid);
fstpv.write(CompileCode);
fstpv.write(CompileCode, instruction_number);
}
else
{
FSTPV_ fstpv(type, tsid, lag);
fstpv.write(CompileCode);
fstpv.write(CompileCode, instruction_number);
}
}
}
else
{
FSTPSV_ fstpsv(type, tsid);
fstpsv.write(CompileCode);
fstpsv.write(CompileCode, instruction_number);
}
}
}
@ -1567,7 +1567,7 @@ UnaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalException
}
void
UnaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
UnaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<UnaryOpNode *>(this));
if (it != temporary_terms.end())
@ -1576,23 +1576,23 @@ UnaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDT_ fldt(ii->second);
fldt.write(CompileCode);
fldt.write(CompileCode, instruction_number);
}
else
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDST_ fldst(ii->second);
fldst.write(CompileCode);
fldst.write(CompileCode, instruction_number);
}
return;
}
if (op_code == oSteadyState)
arg->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, true);
arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, true);
else
{
arg->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FUNARY_ funary(op_code);
funary.write(CompileCode);
funary.write(CompileCode, instruction_number);
}
}
@ -2271,7 +2271,7 @@ BinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExceptio
}
void
BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
BinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
// If current node is a temporary term
temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<BinaryOpNode *>(this));
@ -2281,20 +2281,20 @@ BinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDT_ fldt(ii->second);
fldt.write(CompileCode);
fldt.write(CompileCode, instruction_number);
}
else
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDST_ fldst(ii->second);
fldst.write(CompileCode);
fldst.write(CompileCode, instruction_number);
}
return;
}
arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FBINARY_ fbinary(op_code);
fbinary.write(CompileCode);
fbinary.write(CompileCode, instruction_number);
}
void
@ -3205,7 +3205,7 @@ TrinaryOpNode::eval(const eval_context_t &eval_context) const throw (EvalExcepti
}
void
TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
TrinaryOpNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
// If current node is a temporary term
temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<TrinaryOpNode *>(this));
@ -3215,21 +3215,21 @@ TrinaryOpNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDT_ fldt(ii->second);
fldt.write(CompileCode);
fldt.write(CompileCode, instruction_number);
}
else
{
map_idx_t::const_iterator ii = map_idx.find(idx);
FLDST_ fldst(ii->second);
fldst.write(CompileCode);
fldst.write(CompileCode, instruction_number);
}
return;
}
arg1->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg2->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg3->compile(CompileCode, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg1->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg2->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
arg3->compile(CompileCode, instruction_number, lhs_rhs, temporary_terms, map_idx, dynamic, steady_dynamic);
FTRINARY_ ftrinary(op_code);
ftrinary.write(CompileCode);
ftrinary.write(CompileCode, instruction_number);
}
void
@ -3665,7 +3665,7 @@ ExternalFunctionNode::eval(const eval_context_t &eval_context) const throw (Eval
}
void
ExternalFunctionNode::compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
ExternalFunctionNode::compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
cerr << "ExternalFunctionNode::compile: operation impossible!" << endl;
exit(EXIT_FAILURE);

19
preprocessor/ExprNode.hh
View File

@ -241,7 +241,7 @@ public:
};
virtual double eval(const eval_context_t &eval_context) const throw (EvalException) = 0;
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const = 0;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const = 0;
//! Creates a static version of this node
/*!
This method duplicates the current node by creating a similar node from which all leads/lags have been stripped,
@ -387,7 +387,7 @@ public:
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
@ -429,8 +429,13 @@ public:
int equation) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
SymbolType
get_type() const
{
return type;
};
int
get_symb_id() const
{
@ -486,7 +491,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(UnaryOpcode op_code, double v) throw (EvalException);
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
//! Returns operand
expr_t
get_arg() const
@ -549,7 +554,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (EvalException);
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const;
//! Returns first operand
expr_t
@ -620,7 +625,7 @@ public:
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
static double eval_opcode(double v1, TrinaryOpcode op_code, double v2, double v3) throw (EvalException);
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
@ -678,7 +683,7 @@ public:
virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
virtual double eval(const eval_context_t &eval_context) const throw (EvalException);
virtual void compile(ostream &CompileCode, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual void compile(ostream &CompileCode, unsigned int &instruction_number, bool lhs_rhs, const temporary_terms_t &temporary_terms, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
virtual expr_t toStatic(DataTree &static_datatree) const;
virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const;
virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);

2
preprocessor/ModFile.cc
View File

@ -302,7 +302,7 @@ ModFile::computingPass(bool no_tmp_terms)
if (mod_file_struct.simul_present)
{
dynamic_model.initializeVariablesAndEquations();
dynamic_model.computingPass(false, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
dynamic_model.computingPass(true, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
}
else
{

136
preprocessor/ModelTree.cc
View File

@ -207,6 +207,8 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
if (first_derivatives.find(make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))) == first_derivatives.end())
first_derivatives[make_pair(it->first.first, getDerivID(symbol_table.getID(eEndogenous, it->first.second), 0))] = Zero;
}
}
}
@ -486,7 +488,7 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
for (dynamic_jacob_map_t::const_iterator it = dynamic_jacobian.begin(); it != dynamic_jacobian.end(); it++)
{
int lag = it->first.first;
int j_1 = it->first.second.second;
int j_1 = it->first.second.first;
int i_1 = it->first.second.second;
if (variable_blck[i_1] == equation_blck[j_1])
{
@ -503,7 +505,7 @@ ModelTree::getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian
}
void
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const
ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const
{
int nb_var = variable_reordered.size();
int n = nb_var - prologue - epilogue;
@ -570,8 +572,6 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
components_set[endo2block[i]].first.insert(i);
}
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
getVariableLeadLagByBlock(dynamic_jacobian, endo2block, num, equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered);
vector<int> tmp_equation_reordered(equation_reordered), tmp_variable_reordered(variable_reordered);
@ -594,6 +594,23 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
if (Equation_Type[equation_reordered[i+prologue]].first == E_SOLVE || mfs == 0)
add_edge(i, i, G2);
}
//Determines the dynamic structure of each equation
n_static = vector<unsigned int>(prologue+num+epilogue, 0);
n_forward = vector<unsigned int>(prologue+num+epilogue, 0);
n_backward = vector<unsigned int>(prologue+num+epilogue, 0);
n_mixed = vector<unsigned int>(prologue+num+epilogue, 0);
for (int i = 0; i < prologue; i++)
{
if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
n_mixed[i]++;
else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second != 0)
n_forward[i]++;
else if (variable_lag_lead[tmp_variable_reordered[i]].first != 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
n_backward[i]++;
else if (variable_lag_lead[tmp_variable_reordered[i]].first == 0 && variable_lag_lead[tmp_variable_reordered[i]].second == 0)
n_static[i]++;
}
//For each block, the minimum set of feedback variable is computed
// and the non-feedback variables are reordered to get
// a sub-recursive block without feedback variables
@ -611,21 +628,88 @@ ModelTree::computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob
Reorder_the_recursive_variables(G, feed_back_vertices, Reordered_Vertice);
//First we have the recursive equations conditional on feedback variables
for (int j = 0; j < 4; j++)
{
for (vector<int>::iterator its = Reordered_Vertice.begin(); its != Reordered_Vertice.end(); its++)
{
bool something_done = false;
if (j == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
{
n_mixed[prologue+i]++;
something_done = true;
}
else if (j == 1 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second != 0)
{
n_forward[prologue+i]++;
something_done = true;
}
else if (j == 2 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
{
n_backward[prologue+i]++;
something_done = true;
}
else if (j == 3 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[*its +prologue]].second == 0)
{
n_static[prologue+i]++;
something_done = true;
}
if (something_done)
{
equation_reordered[order] = tmp_equation_reordered[*its+prologue];
variable_reordered[order] = tmp_variable_reordered[*its+prologue];
order++;
}
}
}
components_set[i].second.second = Reordered_Vertice;
//Second we have the equations related to the feedback variables
for (int j = 0; j < 4; j++)
{
for (set<int>::iterator its = feed_back_vertices.begin(); its != feed_back_vertices.end(); its++)
{
bool something_done = false;
if (j == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
{
n_mixed[prologue+i]++;
something_done = true;
}
else if (j == 1 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second != 0)
{
n_forward[prologue+i]++;
something_done = true;
}
else if (j == 2 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first != 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
{
n_backward[prologue+i]++;
something_done = true;
}
else if (j == 3 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].first == 0 && variable_lag_lead[tmp_variable_reordered[v_index[vertex(*its, G)]+prologue]].second == 0)
{
n_static[prologue+i]++;
something_done = true;
}
if (something_done)
{
equation_reordered[order] = tmp_equation_reordered[v_index[vertex(*its, G)]+prologue];
variable_reordered[order] = tmp_variable_reordered[v_index[vertex(*its, G)]+prologue];
order++;
}
}
}
}
for (int i = 0; i < epilogue; i++)
{
if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second != 0)
n_mixed[prologue+num+i]++;
else if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second != 0)
n_forward[prologue+num+i]++;
else if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first != 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second == 0)
n_backward[prologue+num+i]++;
else if (variable_lag_lead[tmp_variable_reordered[prologue+num+i]].first == 0 && variable_lag_lead[tmp_variable_reordered[prologue+num+i]].second == 0)
n_static[prologue+num+i]++;
}
inv_equation_reordered = vector<int>(nb_var);
inv_variable_reordered = vector<int>(nb_var);
for (int i = 0; i < nb_var; i++)
@ -665,7 +749,7 @@ ModelTree::printBlockDecomposition(const vector<pair<int, int> > &blocks) const
}
block_type_firstequation_size_mfs_t
ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, const vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered)
ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int,int> > > &block_col_type)
{
int i = 0;
int count_equ = 0, blck_count_simult = 0;
@ -674,6 +758,10 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
block_type_firstequation_size_mfs_t block_type_size_mfs;
BlockSimulationType Simulation_Type, prev_Type = UNKNOWN;
int eq = 0;
unsigned int l_n_static = 0;
unsigned int l_n_forward = 0;
unsigned int l_n_backward = 0;
unsigned int l_n_mixed = 0;
for (i = 0; i < prologue+(int) blocks.size()+epilogue; i++)
{
int first_count_equ = count_equ;
@ -736,6 +824,10 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
else
Simulation_Type = SOLVE_FORWARD_SIMPLE;
}
l_n_static = n_static[i];
l_n_forward = n_forward[i];
l_n_backward = n_backward[i];
l_n_mixed = n_mixed[i];
if (Blck_Size == 1)
{
if (Equation_Type[equation_reordered[eq]].first == E_EVALUATE || Equation_Type[equation_reordered[eq]].first == E_EVALUATE_S)
@ -754,29 +846,35 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j
BlockSimulationType c_Type = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.first;
int c_Size = (block_type_size_mfs[block_type_size_mfs.size()-1]).second.first;
int first_equation = (block_type_size_mfs[block_type_size_mfs.size()-1]).first.second;
block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(++c_Size, block_type_size_mfs[block_type_size_mfs.size()-1].second.second));
c_Size++;
block_type_size_mfs[block_type_size_mfs.size()-1] = make_pair(make_pair(c_Type, first_equation), make_pair(c_Size, c_Size));
if (block_lag_lead[block_type_size_mfs.size()-1].first > Lag)
Lag = block_lag_lead[block_type_size_mfs.size()-1].first;
if (block_lag_lead[block_type_size_mfs.size()-1].second > Lead)
Lead = block_lag_lead[block_type_size_mfs.size()-1].second;
block_lag_lead[block_type_size_mfs.size()-1] = make_pair(Lag, Lead);
pair< pair< unsigned int, unsigned int>, pair<unsigned int, unsigned int> > tmp = block_col_type[block_col_type.size()-1];
block_col_type[block_col_type.size()-1] = make_pair( make_pair(tmp.first.first+l_n_static, tmp.first.second+l_n_forward), make_pair(tmp.second.first+l_n_backward, tmp.second.second+l_n_mixed) );
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
}
}
else
{
block_type_size_mfs.push_back(make_pair(make_pair(Simulation_Type, eq), make_pair(Blck_Size, MFS_Size)));
block_lag_lead.push_back(make_pair(Lag, Lead));
block_col_type.push_back(make_pair( make_pair(l_n_static, l_n_forward), make_pair(l_n_backward, l_n_mixed) ));
}
prev_Type = Simulation_Type;
eq += Blck_Size;
@ -1015,7 +1113,7 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, ostream &output,
}
void
ModelTree::compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const
{
// Local var used to keep track of temp nodes already written
temporary_terms_t tt2;
@ -1023,17 +1121,17 @@ ModelTree::compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt
it != tt.end(); it++)
{
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
fnumexpr.write(code_file);
(*it)->compile(code_file, false, tt2, map_idx, dynamic, steady_dynamic);
fnumexpr.write(code_file, instruction_number);
(*it)->compile(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic);
if (dynamic)
{
FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second));
fstpt.write(code_file);
fstpt.write(code_file, instruction_number);
}
else
{
FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
fstpst.write(code_file);
fstpst.write(code_file, instruction_number);
}
// Insert current node into tt2
tt2.insert(*it);
@ -1108,7 +1206,7 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
}
void
ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
ModelTree::compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const
{
for (int eq = 0; eq < (int) equations.size(); eq++)
{
@ -1116,7 +1214,7 @@ ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt
expr_t lhs = eq_node->get_arg1();
expr_t rhs = eq_node->get_arg2();
FNUMEXPR_ fnumexpr(ModelEquation, eq);
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
// Test if the right hand side of the equation is empty.
double vrhs = 1.0;
try
@ -1129,20 +1227,20 @@ ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_t &tt
if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
{
lhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic);
rhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic);
lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(eq);
fstpr.write(code_file);
fstpr.write(code_file, instruction_number);
}
else // The right hand side of the equation is empty ==> residual=lhs;
{
lhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic);
lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, dynamic, steady_dynamic);
FSTPR_ fstpr(eq);
fstpr.write(code_file);
fstpr.write(code_file, instruction_number);
}
}
}

21
preprocessor/ModelTree.hh
View File

@ -110,7 +110,7 @@ protected:
//! Writes temporary terms
void writeTemporaryTerms(const temporary_terms_t &tt, ostream &output, ExprNodeOutputType output_type) const;
//! Compiles temporary terms
void compileTemporaryTerms(ostream &code_file, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
void compileTemporaryTerms(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, map_idx_t map_idx, bool dynamic, bool steady_dynamic) const;
//! Adds informations for simulation in a binary file
void Write_Inf_To_Bin_File(const string &basename, int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const;
@ -120,7 +120,7 @@ protected:
//! Writes model equations
void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const;
//! Compiles model equations
void compileModelEquations(ostream &code_file, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
void compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
//! Writes LaTeX model file
void writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const;
@ -167,9 +167,9 @@ protected:
//! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations
equation_type_and_normalized_equation_t equationTypeDetermination(const map<pair<int, pair<int, int> >, expr_t> &first_order_endo_derivatives, const vector<int> &Index_Var_IM, const vector<int> &Index_Equ_IM, int mfs) const;
//! Compute the block decomposition and for a non-recusive block find the minimum feedback set
void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered) const;
void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector<int> &equation_reordered, vector<int> &variable_reordered, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector<int> &inv_equation_reordered, vector<int> &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed) const;
//! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, const vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered);
block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector<pair<int, int> > &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector<int> &variable_reordered, const vector<int> &equation_reordered, vector<unsigned int> &n_static, vector<unsigned int> &n_forward, vector<unsigned int> &n_backward, vector<unsigned int> &n_mixed, vector<pair< pair<int, int>, pair<int,int> > > &block_col_type);
//! Determine the maximum number of lead and lag for the endogenous variable in a bloc
void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector<int> &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector<int> &equation_reordered, const vector<int> &variable_reordered) const;
//! Print an abstract of the block structure of the model
@ -189,6 +189,10 @@ protected:
virtual unsigned int getBlockFirstEquation(int block_number) const = 0;
//! Return the size of the block block_number
virtual unsigned int getBlockSize(int block_number) const = 0;
//! Return the number of exogenous variable in the block block_number
virtual unsigned int getBlockExoSize(int block_number) const = 0;
//! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
virtual unsigned int getBlockExoColSize(int block_number) const = 0;
//! Return the number of feedback variable of the block block_number
virtual unsigned int getBlockMfs(int block_number) const = 0;
//! Return the maximum lag in a block
@ -207,11 +211,18 @@ protected:
virtual int getBlockEquationID(int block_number, int equation_number) const = 0;
//! Return the original number of variable variable_number belonging to the block block_number
virtual int getBlockVariableID(int block_number, int variable_number) const = 0;
//! Return the original number of the exogenous variable varexo_number belonging to the block block_number
virtual int getBlockVariableExoID(int block_number, int variable_number) const = 0;
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int getBlockInitialEquationID(int block_number, int equation_number) const = 0;
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialVariableID(int block_number, int variable_number) const = 0;
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int getBlockInitialExogenousID(int block_number, int variable_number) const = 0;
//! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
virtual int getBlockInitialDetExogenousID(int block_number, int variable_number) const = 0;
//! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
virtual int getBlockInitialOtherEndogenousID(int block_number, int variable_number) const = 0;
public:
ModelTree(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_arg);
//! Declare a node as an equation of the model

130
preprocessor/StaticModel.cc
View File

@ -46,28 +46,28 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
}
void
StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_t &map_idx) const
StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx) const
{
first_derivatives_t::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id)));
if (it != first_derivatives.end())
(it->second)->compile(code_file, false, temporary_terms, map_idx, false, false);
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
{
FLDZ_ fldz;
fldz.write(code_file);
fldz.write(code_file, instruction_number);
}
}
void
StaticModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, map_idx_t &map_idx) const
StaticModel::compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eqr, int varr, int lag, map_idx_t &map_idx) const
{
map<pair<int, pair<int, int> >, expr_t>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
if (it != first_chain_rule_derivatives.end())
(it->second)->compile(code_file, false, temporary_terms, map_idx, false, false);
(it->second)->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
else
{
FLDZ_ fldz;
fldz.write(code_file);
fldz.write(code_file, instruction_number);
}
}
@ -406,6 +406,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
ostringstream tmp_output;
ofstream code_file;
unsigned int instruction_number = 0;
bool file_open = false;
string main_name = file_name;
@ -424,8 +425,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
fdimt.write(code_file);
fdimt.write(code_file, instruction_number);
FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(),
SOLVE_FORWARD_COMPLETE,
0,
@ -436,22 +436,22 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
symbol_table.endo_nbr(),
0,
0,
u_count_int
u_count_int,
symbol_table.endo_nbr()
);
fbeginblock.write(code_file);
fbeginblock.write(code_file, instruction_number);
// Add a mapping form node ID to temporary terms order
int j = 0;
for (temporary_terms_t::const_iterator it = temporary_terms.begin();
it != temporary_terms.end(); it++)
map_idx[(*it)->idx] = j++;
compileTemporaryTerms(code_file, temporary_terms, map_idx, false, false);
compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, false, false);
compileModelEquations(code_file, temporary_terms, map_idx, false, false);
compileModelEquations(code_file, instruction_number, temporary_terms, map_idx, false, false);
FENDEQU_ fendequ;
fendequ.write(code_file);
fendequ.write(code_file, instruction_number);
vector<vector<pair<int, int> > > derivatives;
derivatives.resize(symbol_table.endo_nbr());
@ -467,46 +467,46 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
int symb = getSymbIDByDerivID(deriv_id);
unsigned int var = symbol_table.getTypeSpecificID(symb);
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eq, var);
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
if (!derivatives[eq].size())
derivatives[eq].clear();
derivatives[eq].push_back(make_pair(var, count_u));
d1->compile(code_file, false, temporary_terms, map_idx, false, false);
d1->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
FSTPSU_ fstpsu(count_u);
fstpsu.write(code_file);
fstpsu.write(code_file, instruction_number);
count_u++;
}
}
for (int i = 0; i < symbol_table.endo_nbr(); i++)
{
FLDR_ fldr(i);
fldr.write(code_file);
fldr.write(code_file, instruction_number);
for(vector<pair<int, int> >::const_iterator it = derivatives[i].begin();
it != derivatives[i].end(); it++)
{
FLDSU_ fldsu(it->second);
fldsu.write(code_file);
fldsu.write(code_file, instruction_number);
FLDSV_ fldsv(eEndogenous, it->first);
fldsv.write(code_file);
fldsv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
if (it != derivatives[i].begin())
{
FBINARY_ fbinary(oPlus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
}
}
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPSU_ fstpsu(i);
fstpsu.write(code_file);
fstpsu.write(code_file, instruction_number);
}
FENDBLOCK_ fendblock;
fendblock.write(code_file);
fendblock.write(code_file, instruction_number);
FEND_ fend;
fend.write(code_file);
fend.write(code_file, instruction_number);
code_file.close();
}
@ -528,6 +528,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
unsigned int instruction_number = 0;
expr_t lhs = NULL, rhs = NULL;
BinaryOpNode *eq_node;
Uff Uf[symbol_table.endo_nbr()];
@ -546,7 +547,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
//Temporary variables declaration
FDIMT_ fdimt(temporary_terms.size());
fdimt.write(code_file);
fdimt.write(code_file, instruction_number);
for (unsigned int block = 0; block < getNbBlocks(); block++)
{
@ -554,7 +555,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
if (block > 0)
{
FENDBLOCK_ fendblock;
fendblock.write(code_file);
fendblock.write(code_file, instruction_number);
}
int count_u;
int u_count_int = 0;
@ -580,9 +581,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
symbol_table.endo_nbr(),
0,
0,
u_count_int
u_count_int,
symbol_table.endo_nbr()
);
fbeginblock.write(code_file);
fbeginblock.write(code_file, instruction_number);
// The equations
for (i = 0; i < (int) block_size; i++)
@ -596,10 +598,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
it != v_temporary_terms[block][i].end(); it++)
{
FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find((*it)->idx)->second));
fnumexpr.write(code_file);
(*it)->compile(code_file, false, tt2, map_idx, false, false);
fnumexpr.write(code_file, instruction_number);
(*it)->compile(code_file, instruction_number, false, tt2, map_idx, false, false);
FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second));
fstpst.write(code_file);
fstpst.write(code_file, instruction_number);
// Insert current node into tt2
tt2.insert(*it);
}
@ -615,23 +617,23 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
equ_type = getBlockEquationType(block, i);
{
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
}
if (equ_type == E_EVALUATE)
{
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
lhs->compile(code_file, true, temporary_terms, map_idx, false, false);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
}
else if (equ_type == E_EVALUATE_S)
{
eq_node = (BinaryOpNode *) getBlockEquationRenormalizedExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
lhs->compile(code_file, true, temporary_terms, map_idx, false, false);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
lhs->compile(code_file, instruction_number, true, temporary_terms, map_idx, false, false);
}
break;
case SOLVE_BACKWARD_COMPLETE:
@ -646,22 +648,22 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
default:
end:
FNUMEXPR_ fnumexpr(ModelEquation, getBlockEquationID(block, i));
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
eq_node = (BinaryOpNode *) getBlockEquationExpr(block, i);
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
lhs->compile(code_file, false, temporary_terms, map_idx, false, false);
rhs->compile(code_file, false, temporary_terms, map_idx, false, false);
lhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
rhs->compile(code_file, instruction_number, false, temporary_terms, map_idx, false, false);
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPR_ fstpr(i - block_recursive);
fstpr.write(code_file);
fstpr.write(code_file, instruction_number);
}
}
FENDEQU_ fendequ;
fendequ.write(code_file);
fendequ.write(code_file, instruction_number);
// The Jacobian if we have to solve the block
if (simulation_type != EVALUATE_BACKWARD
&& simulation_type != EVALUATE_FORWARD)
@ -672,12 +674,12 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
case SOLVE_FORWARD_SIMPLE:
{
FNUMEXPR_ fnumexpr(FirstEndoDerivative, 0, 0);
fnumexpr.write(code_file);
fnumexpr.write(code_file, instruction_number);
}
compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
compileDerivative(code_file, instruction_number, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
{
FSTPG_ fstpg(0);
fstpg.write(code_file);
fstpg.write(code_file, instruction_number);
}
break;
@ -706,10 +708,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
Uf[eqr].Ufl->u = count_u;
Uf[eqr].Ufl->var = varr;
FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
fnumexpr.write(code_file);
compileChainRuleDerivative(code_file, eqr, varr, 0, map_idx);
fnumexpr.write(code_file, instruction_number);
compileChainRuleDerivative(code_file, instruction_number, eqr, varr, 0, map_idx);
FSTPSU_ fstpsu(count_u);
fstpsu.write(code_file);
fstpsu.write(code_file, instruction_number);
count_u++;
}
}
@ -718,24 +720,24 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
if (i >= (int) block_recursive)
{
FLDR_ fldr(i-block_recursive);
fldr.write(code_file);
fldr.write(code_file, instruction_number);
FLDZ_ fldz;
fldz.write(code_file);
fldz.write(code_file, instruction_number);
v = getBlockEquationID(block, i);
for (Uf[v].Ufl = Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl = Uf[v].Ufl->pNext)
{
FLDSU_ fldsu(Uf[v].Ufl->u);
fldsu.write(code_file);
fldsu.write(code_file, instruction_number);
FLDSV_ fldsv(eEndogenous, Uf[v].Ufl->var);
fldsv.write(code_file);
fldsv.write(code_file, instruction_number);
FBINARY_ fbinary(oTimes);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FCUML_ fcuml;
fcuml.write(code_file);
fcuml.write(code_file, instruction_number);
}
Uf[v].Ufl = Uf[v].Ufl_First;
while (Uf[v].Ufl)
@ -745,10 +747,10 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
Uf[v].Ufl = Uf[v].Ufl_First;
}
FBINARY_ fbinary(oMinus);
fbinary.write(code_file);
fbinary.write(code_file, instruction_number);
FSTPSU_ fstpsu(i - block_recursive);
fstpsu.write(code_file);
fstpsu.write(code_file, instruction_number);
}
}
@ -759,9 +761,9 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
}
}
FENDBLOCK_ fendblock;
fendblock.write(code_file);
fendblock.write(code_file, instruction_number);
FEND_ fend;
fend.write(code_file);
fend.write(code_file, instruction_number);
code_file.close();
}
@ -858,6 +860,7 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
if (block)
{
jacob_map_t contemporaneous_jacobian, static_jacobian;
vector<unsigned int> n_static, n_forward, n_backward, n_mixed;
// for each block contains pair<Size, Feddback_variable>
vector<pair<int, int> > blocks;
@ -874,10 +877,11 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
cout << "Finding the optimal block decomposition of the model ...\n";
if (prologue+epilogue < (unsigned int) equation_number())
computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered);
lag_lead_vector_t equation_lag_lead, variable_lag_lead;
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered);
computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed);
block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type);
printBlockDecomposition(blocks);

49
preprocessor/StaticModel.hh
View File

@ -82,9 +82,9 @@ private:
void computeTemporaryTermsMapping();
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_t &map_idx) const;
void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx) const;
//! Write chain rule derivative code of an equation w.r. to a variable
void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_t &map_idx) const;
void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, map_idx_t &map_idx) const;
//! Get the type corresponding to a derivation ID
virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
@ -147,6 +147,16 @@ protected:
typedef map<int, var_t> lag_var_t;
vector<lag_var_t> other_endo_block, exo_block, exo_det_block;
//! for each block described the number of static, forward, backward and mixed variables in the block
/*! pair< pair<static, forward>, pair<backward,mixed> > */
vector<pair< pair<int, int>, pair<int,int> > > block_col_type;
//! List for each variable its block number and its maximum lag and lead inside the block
vector<pair<int, pair<int, int> > > variable_block_lead_lag;
//! List for each equation its block number
vector<int> equation_block;
//!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
@ -219,6 +229,16 @@ public:
{
return (block_type_firstequation_size_mfs[block_number].second.first);
};
//! Return the number of exogenous variable in the block block_number
virtual unsigned int getBlockExoSize(int block_number) const
{
return 0;
};
//! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
virtual unsigned int getBlockExoColSize(int block_number) const
{
return 0;
}
//! Return the number of feedback variable of the block block_number
virtual unsigned int
getBlockMfs(int block_number) const
@ -273,6 +293,12 @@ public:
{
return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
};
//! Return the original number of the exogenous variable varexo_number belonging to the block block_number
virtual int
getBlockVariableExoID(int block_number, int variable_number) const
{
return 0;
};
//! Return the position of equation_number in the block number belonging to the block block_number
virtual int
getBlockInitialEquationID(int block_number, int equation_number) const
@ -285,7 +311,24 @@ public:
{
return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
};
//! Return the position of variable_number in the block number belonging to the block block_number
virtual int
getBlockInitialExogenousID(int block_number, int variable_number) const
{
return -1;
};
//! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
virtual int
getBlockInitialDetExogenousID(int block_number, int variable_number) const
{
return -1;
};
//! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
virtual int
getBlockInitialOtherEndogenousID(int block_number, int variable_number) const
{
return -1;
};
};
#endif