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

- sparse_dll option works fine with feedback variables 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@2851 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2009-07-21 15:50:12 +00:00
parent f5920bbc34
commit 1d9260251d
20 changed files with 971 additions and 681 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
matlab/solve_two_boundaries.m
View File

@ -58,7 +58,7 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global oo_ M_;
global oo_ M_ T9025 T1149 T11905;
cvg=0;
iter=0;
Per_u_=0;
@ -73,8 +73,33 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_
reduced = 0;
while ~(cvg==1 | iter>maxit_),
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, y_kmin, Blck_size);
% for i=1:periods;
% disp([sprintf('%5.14f ',[T9025(i) T1149(i) T11905(i)])]);
% end;
% return;
residual = r(:,y_kmin+1:y_kmin+1+y_kmax_l);
%num2str(residual,' %1.6f')
%jac_ = g1(1:(y_kmin)*Blck_size, 1:(y_kmin+1+y_kmax_l)*Blck_size);
%jac_
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
b = b -g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';
term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
term2 = g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';
b = b - term1 - term2;
% fid = fopen(['result' num2str(iter)],'w');
% fg1a = full(g1a);
% fprintf(fid,'%d\n',size(fg1a,1));
% fprintf(fid,'%d\n',size(fg1a,2));
% fprintf(fid,'%5.14f\n',fg1a);
% fprintf(fid,'%d\n',size(b,1));
% fprintf(fid,'%5.14f\n',b);
% fclose(fid);
% return;
%ipconfigb_ = b(1:(1+y_kmin)*Blck_size);
%b_
[max_res, max_indx]=max(max(abs(r')));
if(~isreal(r))
max_res = (-max_res^2)^0.5;
@ -151,6 +176,14 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_
dx = g1a\b- ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
% v = '';
% for i=1:(size(y_index,2))
% v = [v ' %1.6f'];
% end;
% v = [v '\n'];
% v
% sprintf(v,y(:,y_index)')
% return;
elseif(simulation_method==2),
flag1=1;
while(flag1>0)

209
mex/sources/simulate/Interpreter.cc
View File

@ -144,7 +144,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
//load a temporary variable in the processor
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FLDT %d\n",var);
mexPrintf("FLDT %d [%d]=%f Stack.size()=%d\n",var,var*(periods+y_kmin+y_kmax)+it_,T[var*(periods+y_kmin+y_kmax)+it_], Stack.size());
mexEvalString("drawnow;");
#endif
Stack.push(T[var*(periods+y_kmin+y_kmax)+it_]);
@ -161,11 +161,11 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
break;
case FLDR :
//load u variable in the processor
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FLDR\n");
mexPrintf("FLDR r[%d]=%f\n",var, r[var] );
mexEvalString("drawnow;");
#endif
var=get_code_int;
Stack.push(r[var]);
break;
case FLDZ :
@ -212,7 +212,7 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
lag=get_code_int;
#ifdef DEBUGC
//mexPrintf("y[%d(it_=%d, lag=%d, y_size=%d, var=%d)](%d)=",(it_+lag)*y_size+var,it_, lag, y_size, var, Stack.size());
mexPrintf("y[it_=%d, lag=%d, y_size=%d, var=%d, block=%d)]=",it_, lag, y_size, var+1, Block_Count+1);
mexPrintf("FSTP y[it_=%d, lag=%d, y_size=%d, var=%d, block=%d)]=",it_, lag, y_size, var+1, Block_Count+1);
mexEvalString("drawnow;");
#endif
y[(it_+lag)*y_size+var] = Stack.top();
@ -226,11 +226,15 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
#ifdef DEBUGC
mexPrintf("Exogenous\n");
#endif
var=get_code_int;
//var=get_code_int;
var=get_code_int;
lag=get_code_int;
/*mexPrintf("FSTP x[it_=%d, lag=%d, y_size=%d, var=%d, block=%d)]=",it_, lag, y_size, var+1, Block_Count+1);
mexEvalString("drawnow;");*/
x[it_+lag+var*nb_row_x] = Stack.top();
Stack.pop();
/*mexPrintf("%f\n",x[it_+lag+var*nb_row_x]);
mexEvalString("drawnow;");*/
break;
case eExogenousDet :
#ifdef DEBUGC
@ -278,9 +282,13 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
case FSTPR :
//load u variable in the processor
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FSTPR residual[%d]=",var);
mexEvalString("drawnow;");
#endif
r[var] = Stack.top();
#ifdef DEBUGC
mexPrintf("FSTPR residual[%d]=%f\n",var,r[var]);
mexPrintf("%f\n",r[var]);
mexEvalString("drawnow;");
#endif
Stack.pop();
@ -288,9 +296,13 @@ Interpreter::compute_block_time(int Per_u_) /*throw(EvalException)*/
case FSTPG :
//load u variable in the processor
var=get_code_int;
#ifdef DEBUGC
mexPrintf("FSTPG g1[%d]=",var);
mexEvalString("drawnow;");
#endif
g1[var] = Stack.top();
#ifdef DEBUGC
mexPrintf("FSTPG g1[%d)=%f\n",var,g1[var]);
mexPrintf("%f\n",g1[var]);
mexEvalString("drawnow;");
#endif
Stack.pop();
@ -628,10 +640,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
{
case EVALUATE_FORWARD :
//case EVALUATE_FORWARD_R :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("EVALUATE_FORWARD\n");
mexEvalString("drawnow;");
//#endif
#endif
begining=get_code_pointer;
//mexPrintf("y_kmin=%d periods=%d",y_kmin, periods);
for (it_=y_kmin;it_<periods+y_kmin;it_++)
@ -641,14 +653,19 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Per_y_=it_*y_size;
//mexPrintf("bef compute_block_time()\n");
compute_block_time(0);
//mexPrintf("x(%d, 638)=%1.14f/x(%d, 116)=%1.14f=%1.14f\n",it_,x[it_+638*nb_row_x],it_,x[it_+116*nb_row_x],x[it_+638*nb_row_x]/x[it_+116*nb_row_x]);
#ifdef PRINT_OUT
for (int j = 0; j<size; j++)
mexPrintf("y[%d, %d] = %1.14f\n", Block_Contain[j].Variable, it_, y[Per_y_ + Block_Contain[j].Variable]);
#endif
}
break;
case EVALUATE_BACKWARD :
//case EVALUATE_BACKWARD_R :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("EVALUATE_BACKWARD\n");
mexEvalString("drawnow;");
//#endif
#endif
begining=get_code_pointer;
for (it_=periods+y_kmin-1;it_>=y_kmin;it_--)
{
@ -656,16 +673,16 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
Per_y_=it_*y_size;
compute_block_time(0);
#ifdef PRINT_OUT
for (j = 0; j<size; j++)
mexPrintf("y[%d, %d] = %f\n", Block_Contain[j].Variable, it_, y[Per_y_ + Block_Contain[j].Variable]);
for (int j = 0; j<size; j++)
mexPrintf("y[%d, %d] = %1.14f\n", Block_Contain[j].Variable, it_, y[Per_y_ + Block_Contain[j].Variable]);
#endif
}
break;
case SOLVE_FORWARD_SIMPLE :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("SOLVE_FORWARD_SIMPLE\n");
mexEvalString("drawnow;");
//#endif
#endif
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
@ -700,10 +717,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mxFree(r);
break;
case SOLVE_BACKWARD_SIMPLE :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("SOLVE_BACKWARD_SIMPLE\n");
mexEvalString("drawnow;");
//#endif
#endif
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
@ -739,28 +756,28 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mxFree(r);
break;
case SOLVE_FORWARD_COMPLETE :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("SOLVE FORWARD_COMPLETE\n");
mexEvalString("drawnow;");
//#endif
#endif
is_linear=get_code_bool;
//mexPrintf("is_linear=%d\n",is_linear);
//mexEvalString("drawnow;");
/*mexPrintf("is_linear=%d\n",is_linear);
mexEvalString("drawnow;");*/
max_lag_plus_max_lead_plus_1=get_code_int;
//mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1);
//mexEvalString("drawnow;");
/*mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1);
mexEvalString("drawnow;");*/
symbol_table_endo_nbr=get_code_int;
//mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr);
//mexEvalString("drawnow;");
/*mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr);
mexEvalString("drawnow;");*/
Block_List_Max_Lag=get_code_int;
//mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag);
//mexEvalString("drawnow;");
/*mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag);
mexEvalString("drawnow;");*/
Block_List_Max_Lead=get_code_int;
//mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead);
//mexEvalString("drawnow;");
/*mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead);
mexEvalString("drawnow;");*/
u_count_int=get_code_int;
//mexPrintf("u_count_int=%d\n",u_count_int);
//mexEvalString("drawnow;");
/*mexPrintf("u_count_int=%d\n",u_count_int);
mexEvalString("drawnow;");*/
fixe_u(&u, u_count_int, u_count_int);
//Read_file(file_name, periods, 0, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, nb_endo, u_count, u_count_init, u);
//sparse_matrix.initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);
@ -850,28 +867,39 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
/*mexPrintf("end of forward solve\n");
mexEvalString("drawnow;");*/
mxFree(index_equa);
mxFree(index_vara);
memset(direction,0,size_of_direction);
mxFree(g1);
mxFree(r);
mxFree(u);
break;
case SOLVE_BACKWARD_COMPLETE :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("SOLVE_BACKWARD_COMPLETE\n");
mexEvalString("drawnow;");
//#endif
#endif
is_linear=get_code_bool;
max_lag_plus_max_lead_plus_1=get_code_int;
symbol_table_endo_nbr=get_code_int;
Block_List_Max_Lag=get_code_int;
Block_List_Max_Lead=get_code_int;
//Read_file(file_name, periods, 0, symbol_table_endo_nbr, Block_List_Max_Lag, Block_List_Max_Lead, nb_endo, u_count, u_count_init, u);
//sparse_matrix.initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);
u_count_int=get_code_int;
//mexPrintf("u_count_int=%d\n",u_count_int);
//mexPrintf("is_linear=%d\n",is_linear);
//mexEvalString("drawnow;");
max_lag_plus_max_lead_plus_1=get_code_int;
//mexPrintf("max_lag_plus_max_lead_plus_1=%d\n",max_lag_plus_max_lead_plus_1);
//mexEvalString("drawnow;");
symbol_table_endo_nbr=get_code_int;
//mexPrintf("symbol_table_endo_nbr=%d\n",symbol_table_endo_nbr);
//mexEvalString("drawnow;");
Block_List_Max_Lag=get_code_int;
//mexPrintf("Block_List_Max_Lag=%d\n",Block_List_Max_Lag);
//mexEvalString("drawnow;");
Block_List_Max_Lead=get_code_int;
//mexPrintf("Block_List_Max_Lead=%d\n",Block_List_Max_Lead);
//mexEvalString("drawnow;");
u_count_int=get_code_int;
/*mexPrintf("u_count_int=%d\n",u_count_int);
mexEvalString("drawnow;");*/
fixe_u(&u, u_count_int, u_count_int);
Read_SparseMatrix(bin_basename, size, 1, 0, 0);
//mexPrintf("size=%d\n",size);
//mexEvalString("drawnow;");
g1=(double*)mxMalloc(size*size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
begining=get_code_pointer;
@ -929,6 +957,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
simulate_NG(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, false, cvg, iter);
}
}
mxFree(index_equa);
mxFree(index_vara);
memset(direction,0,size_of_direction);
mxFree(g1);
mxFree(r);
@ -939,10 +969,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
#if GNUVER >= 432
//mexPrintf("omp_get_max_threads=%d\n",omp_get_max_threads());
#endif
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n");
mexEvalString("drawnow;");
//#endif
#endif
is_linear=get_code_bool;
#ifdef DEBUGC
mexPrintf("is_linear=%d\n",is_linear);
@ -977,16 +1007,24 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
//sparse_matrix.initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);
fixe_u(&u, u_count_int, max_lag_plus_max_lead_plus_1);
//fixe_u(&u, u_count_int, max_lag_plus_max_lead_plus_1);
fixe_u(&u, u_count_int, u_count_int);
#ifdef DEBUGC
mexPrintf("u=%x\n",u);
mexPrintf("size=%d\n",size);
#endif
Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax);
#ifdef DEBUGC
mexPrintf("size=%d\n",size);
mexEvalString("drawnow;");
#endif
//mexPrintf("aft reading_sparse_matrix\n");
//mexEvalString("drawnow;");
u_count=u_count_int*(periods+y_kmax+y_kmin);
//g1=(double*)mxMalloc(size*size*sizeof(double));
//mexPrintf("r=(double*)mxMalloc(%d)=",size*sizeof(double));
r=(double*)mxMalloc(size*sizeof(double));
//mexPrintf("%x\n",r);
y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin));
#ifdef DEBUGC
mexPrintf("u_count=%d\n",u_count);
@ -1023,9 +1061,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
max_res=0;
max_res_idx=0;
memcpy(y_save, y, y_size*sizeof(double)*(periods+y_kmax+y_kmin));
//double res[size][13];
for (it_=y_kmin;it_<periods+y_kmin;it_++)
{
Per_u_=(it_-y_kmin)*max_lag_plus_max_lead_plus_1;
Per_u_=(it_-y_kmin)*/*max_lag_plus_max_lead_plus_1*/u_count_int;
//mexPrintf("Per_u_=%d\n",Per_u_);
Per_y_=it_*y_size;
//mexPrintf("ok\n");
@ -1033,6 +1072,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mexEvalString("drawnow;");*/
set_code_pointer(begining);
compute_block_time(Per_u_);
/*mexPrintf("periods it_=%d ",it_);
for (i=0; i< size; i++)
mexPrintf(" %f ",r[i]);
mexPrintf("\n");*/
/*mexPrintf("end of compute_block_time it_=%d\n",it_);*/
/*if(Gaussian_Elimination)
initialize(periods, nb_endo, y_kmin, y_kmax, y_size, u_count, u_count_init, u, y, ya, slowc, y_decal, markowitz_c, res1, res2, max_res);*/
@ -1045,6 +1088,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}
for (i=0; i< size; i++)
{
/*if(it_<13)
res[i][it_-y_kmin]=r[i];*/
double rr;
/*if(fabs(y[Per_y_+Block_Contain[i].Variable])>solve_tolf)*/
//mexPrintf("res[%d]=%f\n",i,r[i]);
@ -1069,9 +1114,19 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}*/
/*mexPrintf("r[%d] (i=%d)",i+size*(it_-y_kmin),i);
mexPrintf("=%f\n",r[i]);*/
//mexPrintf("u[b[%d]=%d]=%f\n",i+(it_-y_kmin)*u_count_int,b[i+(it_-y_kmin)*u_count_int],u[b[i+(it_-y_kmin)*u_count_int]]);
}
//mexPrintf("---------------------------------------------------------------------------------\n");
//mexPrintf("(log(y(%d, 151))) - (x(%d, 338)=%f\n",it_, it_,y[Per_y_+ 150] - exp(x[it_+337*nb_row_x]));
}
/*for(i=0;i<size;i++)
{
for(int j=0;j<5;j++)
mexPrintf(" % 1.6f ",res[i][j]);
mexPrintf("\n");
}*/
cvg=(max_res<solve_tolf);
//mexPrintf("it_=%d\n",it_);
if(Gaussian_Elimination)
{
/*mexPrintf("bef simulate_NG1\n");
@ -1079,6 +1134,13 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
u_count=u_count_saved;
/*mexPrintf("u_count=%d &u_count=%x\n",u_count,&u_count);
mexEvalString("drawnow;");*/
/*for(int t=0;t<periods;t++)
{
mexPrintf("%5.14f %5.14f %5.14f %5.14f \n",T[437*(periods+y_kmin+y_kmax)+t], T[72*(periods+y_kmin+y_kmax)+t], T[473*(periods+y_kmin+y_kmax)+t],
T[437*(periods+y_kmin+y_kmax)+t] * T[72*(periods+y_kmin+y_kmax)+t] * T[473*(periods+y_kmin+y_kmax)+t]);
}
filename=" stopped";
mexErrMsgTxt(filename.c_str());*/
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
/*mexPrintf("after simulate_NG1\n");
mexEvalString("drawnow;");*/
@ -1103,12 +1165,17 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}
else
{
res1=res2=max_res=0;max_res_idx=0;
for (it_=y_kmin;it_<periods+y_kmin;it_++)
{
Per_u_=(it_-y_kmin)*max_lag_plus_max_lead_plus_1;
Per_u_=(it_-y_kmin)*/*max_lag_plus_max_lead_plus_1*/u_count_int;
Per_y_=it_*y_size;
set_code_pointer(begining);
compute_block_time(Per_u_);
/*mexPrintf("periods it_=%d ",it_);
for (i=0; i< size; i++)
mexPrintf(" %f ",r[i]);
mexPrintf("\n");*/
#ifdef PRINT_OUT
for (j=0; j<max_lag_plus_max_lead_plus_1; j++)
{
@ -1120,9 +1187,30 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
for(i=0; i<y_size; i++)
mexPrintf(" %f",y[i]);
mexPrintf("\n");*/
for (i=0; i< size; i++)
{
double rr;
/*if(fabs(y[Per_y_+Block_Contain[i].Variable])>solve_tolf)*/
//mexPrintf("res[%d]=%f\n",i,r[i]);
/*if(fabs(1+y[Per_y_+Block_Contain[i].Variable])>eps)
rr=r[i]/(1+y[Per_y_+Block_Contain[i].Variable]);
else*/
rr=r[i];
/*else
rr=r[i];*/
if (max_res<fabs(rr))
{
max_res=fabs(rr);
max_res_idx=i;
}
res2+=rr*rr;
res1+=fabs(rr);
}
}
res1=res2=max_res=0;max_res_idx=0;
//res1=res2=max_res=0;max_res_idx=0;
cvg = false;
//mexPrintf("it_=%d\n",it_);
if(Gaussian_Elimination)
{
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
@ -1238,46 +1326,46 @@ Interpreter::compute_blocks(string file_name, string bin_basename)
mexEvalString("drawnow;");
#endif
lBlock.size=get_code_int;
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("Block[Block_Count].size=%d\n",lBlock.size);
mexEvalString("drawnow;");
//#endif
#endif
lBlock.type=get_code_int;
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("Block[Block_Count].type=%d\n",lBlock.type);
mexEvalString("drawnow;");
//#endif
#endif
Block.push_back(lBlock);
for (int i=0;i</*Block[Block_Count].size*/lBlock.size;i++)
{
lBlock_Contain.Variable=get_code_int;
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("Block_Contain[%d].Variable=%d\n",i,lBlock_Contain.Variable);
mexEvalString("drawnow;");
//#endif
#endif
lBlock_Contain.Equation=get_code_int;
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("Block_Contain[%d].Equation=%d\n",i,lBlock_Contain.Equation);
mexEvalString("drawnow;");
//#endif
#endif
lBlock_Contain.Own_Derivative=get_code_int;
//mexPrintf("Block_Contain[%d].Own_Derivative=%d\n",i,lBlock_Contain.Own_Derivative);
Block_Contain.push_back(lBlock_Contain);
}
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("Block Completed\n");
mexEvalString("drawnow;");
//#endif
#endif
simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,true);
/*mexPrintf("after simulate_a_block\n");
mexEvalString("drawnow;");*/
//simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,false);
break;
case FEND :
//#ifdef DEBUGC
#ifdef DEBUGC
mexPrintf("FEND\n");
mexEvalString("drawnow;");
//#endif
#endif
go_on=false;
break;
case FDIMT :
@ -1304,6 +1392,3 @@ Interpreter::compute_blocks(string file_name, string bin_basename)
/*mexPrintf("compute_blocks\n");
mexEvalString("drawnow;");*/
}

6
mex/sources/simulate/Interpreter.hh
View File

@ -32,7 +32,11 @@
#ifdef LINBCG
# include "linbcg.hh"
#endif
#include "mex.h"
#ifndef DEBUG_EX
#include "mex.h"
#else
#include "mex_interface.hh"
#endif
//#define DEBUGC

78
mex/sources/simulate/Mem_Mngr.cc
View File

@ -23,7 +23,6 @@ Mem_Mngr::Mem_Mngr()
{
swp_f=false;
swp_f_b=0;
//verbose=false;
}
void
Mem_Mngr::Print_heap()
@ -44,6 +43,7 @@ Mem_Mngr::init_Mem()
NZE_Mem=NULL;
NZE_Mem_add=NULL;
CHUNK_heap_pos=0;
NZE_Mem_Allocated.clear();
}
void Mem_Mngr::fixe_file_name(string filename_arg)
@ -60,7 +60,7 @@ Mem_Mngr::init_CHUNK_BLCK_SIZE(int u_count)
NonZeroElem*
Mem_Mngr::mxMalloc_NZE()
{
int i;
long int i;
if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/
{
NonZeroElem* p1 = Chunk_Stack.back();
@ -69,32 +69,27 @@ Mem_Mngr::mxMalloc_NZE()
}
else if (CHUNK_heap_pos<CHUNK_SIZE) /*there is enough allocated memory space available we keep it at the top of the heap*/
{
int i=CHUNK_heap_pos++;
i=CHUNK_heap_pos++;
return(NZE_Mem_add[i]);
}
else /*We have to allocate extra memory space*/
{
//mexPrintf("CHUNK_SIZE=%d CHUNK_BLCK_SIZE=%d Nb_CHUNK=%d\n",CHUNK_SIZE,CHUNK_BLCK_SIZE,Nb_CHUNK);
CHUNK_SIZE+=CHUNK_BLCK_SIZE;
/*mexPrintf("Allocate %f Ko\n",double(CHUNK_BLCK_SIZE)*double(sizeof(NonZeroElem))/double(1024));
mexEvalString("drawnow;");*/
Nb_CHUNK++;
#ifdef MEM_ALLOC_CHK
mexPrintf("CHUNK_BLCK_SIZE=%d\n",CHUNK_BLCK_SIZE);
#endif
NZE_Mem=(NonZeroElem*)mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem));
//mexPrintf("in mxMalloc NZE_Mem=%x CHUNK_heap_pos=%d CHUNK_BLCK_SIZE=%d Nb_CHUNK=%d\n",NZE_Mem, CHUNK_heap_pos, CHUNK_BLCK_SIZE, Nb_CHUNK);
NZE_Mem=(NonZeroElem*)mxMalloc(CHUNK_BLCK_SIZE*sizeof(NonZeroElem)); /*The block of memory allocated*/
NZE_Mem_Allocated.push_back(NZE_Mem);
if(!NZE_Mem)
{
mexPrintf("Not enough memory available\n");
mexEvalString("drawnow;");
}
#ifdef MEM_ALLOC_CHK
mexPrintf("CHUNK_SIZE=%d\n",CHUNK_SIZE);
#endif
NZE_Mem_add=(NonZeroElem**)mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem*));
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
#endif
NZE_Mem_add=(NonZeroElem**)mxRealloc(NZE_Mem_add, CHUNK_SIZE*sizeof(NonZeroElem*)); /*We have to redefine the size of pointer on the memory*/
if(!NZE_Mem_add)
{
mexPrintf("Not enough memory available\n");
mexEvalString("drawnow;");
}
for (i=CHUNK_heap_pos;i<CHUNK_SIZE;i++)
{
NZE_Mem_add[i]=(NonZeroElem*)(NZE_Mem+(i-CHUNK_heap_pos));
@ -109,36 +104,13 @@ void
Mem_Mngr::mxFree_NZE(void* pos)
{
int i, gap;
/*if(verbose)
{
mexPrintf("pos=%x Nb_CHUNK=%d CHUNK_BLCK_SIZE=%d\n",pos,Nb_CHUNK, CHUNK_BLCK_SIZE);
mexEvalString("drawnow;");
}
*/
for (i=0;i<Nb_CHUNK;i++)
{
/*if(verbose)
{
mexPrintf("i=%d\n",i);
mexEvalString("drawnow;");
mexPrintf("NZE_Mem_add[i*CHUNK_BLCK_SIZE]=%d\n",NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
mexEvalString("drawnow;");
}*/
gap=((uint64_t)(pos)-(uint64_t)(NZE_Mem_add[i*CHUNK_BLCK_SIZE]))/sizeof(NonZeroElem);
if ((gap<CHUNK_BLCK_SIZE) && (gap>=0))
break;
}
/*if(verbose)
{
mexPrintf("push_back()\n");
mexEvalString("drawnow;");
}*/
Chunk_Stack.push_back((NonZeroElem*)pos);
/*if(verbose)
{
mexPrintf("End\n");
mexEvalString("drawnow;");
}*/
}
@ -151,9 +123,6 @@ Mem_Mngr::write_swp_f(int *save_op_all,long int *nop_all)
if (!SaveCode_swp.is_open())
{
mexPrintf("open the swp file for writing\n");
#ifdef PRINT_OUT
mexPrintf("file opened\n");
#endif
SaveCode_swp.open((filename + ".swp").c_str(), std::ios::out | std::ios::binary);
if (!SaveCode_swp.is_open())
{
@ -161,9 +130,6 @@ Mem_Mngr::write_swp_f(int *save_op_all,long int *nop_all)
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
#ifdef PRINT_OUT
mexPrintf("done\n");
#endif
}
SaveCode_swp.write(reinterpret_cast<char *>(nop_all), sizeof(*nop_all));
SaveCode_swp.write(reinterpret_cast<char *>(save_op_all), (*nop_all)*sizeof(int));
@ -177,9 +143,6 @@ Mem_Mngr::read_swp_f(int **save_op_all,long int *nop_all)
swp_f=true;
if (!SaveCode_swp.is_open())
{
#ifdef PRINT_OUT
mexPrintf("file opened\n");
#endif
mexPrintf("open the file %s\n",(filename + ".swp").c_str());
SaveCode_swp.open((filename + ".swp").c_str(), std::ios::in | std::ios::binary);
j=SaveCode_swp.is_open();
@ -191,9 +154,6 @@ Mem_Mngr::read_swp_f(int **save_op_all,long int *nop_all)
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
#ifdef PRINT_OUT
mexPrintf("done\n");
#endif
SaveCode_swp.seekg(0);
}
@ -277,15 +237,11 @@ Mem_Mngr::chk_avail_mem(int **save_op_all,long int *nop_all,long int *nopa_all,i
void
Mem_Mngr::Free_All()
{
int i;
/*mexPrintf("Nb_CHUNK=%d\n",Nb_CHUNK);
mexEvalString("drawnow;");*/
for (i=0;i<Nb_CHUNK;i++)
{
/*mexPrintf("NZE_Mem_add[%d]=%x\n",i*CHUNK_BLCK_SIZE,NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
mexEvalString("drawnow;");*/
mxFree(NZE_Mem_add[i*CHUNK_BLCK_SIZE]);
}
while(NZE_Mem_Allocated.size())
{
mxFree(NZE_Mem_Allocated.back());
NZE_Mem_Allocated.pop_back();
}
mxFree(NZE_Mem_add);
init_Mem();
}

8
mex/sources/simulate/Mem_Mngr.hh
View File

@ -20,9 +20,14 @@
#ifndef MEM_MNGR_HH_INCLUDED
#define MEM_MNGR_HH_INCLUDED
//#include <stack>
#include <vector>
#include <fstream>
#include "mex.h"
#ifndef DEBUG_EX
#include "mex.h"
#else
#include "mex_interface.hh"
#endif
using namespace std;
struct NonZeroElem
@ -59,6 +64,7 @@ private:
int CHUNK_heap_pos/*, CHUNK_heap_max_size*/;
NonZeroElem** NZE_Mem_add;
NonZeroElem* NZE_Mem;
vector<NonZeroElem*> NZE_Mem_Allocated;
int swp_f_b;
fstream SaveCode_swp;
string filename;

509
mex/sources/simulate/SparseMatrix.cc
View File

@ -18,7 +18,7 @@
*/
#include <cstring>
#include <sstream>
#include "SparseMatrix.hh"
SparseMatrix::SparseMatrix()
@ -33,10 +33,13 @@ SparseMatrix::SparseMatrix()
max_u=0;
min_u=0x7FFFFFFF;
res1a=9.0e60;
tbreak_g=0;
start_compare=0;
restart = 0;
}
//
int SparseMatrix::NRow(int r)
@ -156,6 +159,7 @@ double tdelete1=0, tdelete2=0, tdelete21=0, tdelete22=0, tdelete221=0, tdelete22
void SparseMatrix::Delete(const int r,const int c, const int Size)
{
//mexPrintf("Delete r=%d c=%d\n",r,c);
NonZeroElem *first=FNZE_R[r], *firsta=NULL;
#ifdef PROFILER
clock_t td0, td1, td2;
@ -220,6 +224,23 @@ void SparseMatrix::Delete(const int r,const int c, const int Size)
tdelete22+=clock()-td1;
tdelete2+=clock()-td0;
#endif
/*Check the deletition*/
/*int nb_var=NbNZRow[r];
first=FNZE_R[r];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in Delete (Row) r=%d and c=%d \n",r,c);
first=first->NZE_R_N;
}
nb_var=NbNZCol[c];
first=FNZE_C[c];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in Delete (Col) r=%d and c=%d \n",r,c);
first=first->NZE_C_N;
}*/
}
@ -274,6 +295,7 @@ void SparseMatrix::Print(int Size, int *b)
void SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_index)
{
//mexPrintf("Insert r=%d c=%d\n",r,c);
#ifdef PRINT_OUT
mexPrintf("In Insert r=%d, c=%d, u=%d, lag=%d \n",r,c,u_index,lag_index);
#endif
@ -315,6 +337,8 @@ void SparseMatrix::Insert(const int r, const int c, const int u_index, const int
}
else /*first.c_index<c*/
{
/*if(first->c_index==c)
mexPrintf("Error in Insert (r=%d, c=%d -Row-) already exist!!\n");*/
first->NZE_R_N=firstn;
firstn->NZE_R_N=NULL;
}
@ -336,10 +360,29 @@ void SparseMatrix::Insert(const int r, const int c, const int u_index, const int
}
else /*first.r_index<r*/
{
/*if(first->r_index==r)
mexPrintf("Error in Insert (r=%d, c=%d -Col-) already exist!!\n");*/
first->NZE_C_N=firstn;
firstn->NZE_C_N=NULL;
}
NbNZCol[c]++;
/*Check the insertion*/
/*int nb_var=NbNZRow[r];
first=FNZE_R[r];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in insert (Row) r=%d and c=%d \n",r,c);
first=first->NZE_R_N;
}
nb_var=NbNZCol[c];
first=FNZE_C[c];
for(int j=0;j<nb_var;j++)
{
if(!first)
mexPrintf("Error in insert (Col) r=%d and c=%d \n",r,c);
first=first->NZE_C_N;
}*/
}
void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, int y_kmin, int y_kmax)
@ -364,12 +407,14 @@ void SparseMatrix::Read_SparseMatrix(string file_name, int Size, int periods, in
#endif
}
IM_i.clear();
//mexPrintf("u_count_init=%d\n",u_count_init);
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));
//mexPrintf("eq=%d var=%d lag=%d j=%d\n",eq, var, lag, j);
IM_i[std::make_pair(std::make_pair(eq, var), lag)] = j;
}
#ifdef MEM_ALLOC_CHK
@ -416,6 +461,7 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m
NonZeroElem* first;
//mexPrintf("periods=%d, y_kmin=%d, y_kmax=%d, SizeInit=%d, IM.size()=%d\n",periods, y_kmin, y_kmax, Size, IM.size());
pivot=(int*)mxMalloc(Size*sizeof(int));
pivot_save=(int*)mxMalloc(Size*sizeof(int));
pivotk=(int*)mxMalloc(Size*sizeof(int));
pivotv=(double*)mxMalloc(Size*sizeof(double));
pivotva=(double*)mxMalloc(Size*sizeof(double));
@ -445,7 +491,7 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m
it4=IM.begin();
eq=-1;
double tmp_b[Size];
///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for(i=0; i< Size;i++)
{
tmp_b[i]=0;//u[i];
@ -501,12 +547,11 @@ void SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::m
}
it4++;
}
///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for(i=0;i<Size;i++)
{
b[i]=u_count1+i;
u[b[i]]=-tmp_b[i];
//mexPrintf("b[%d]=%f\n",i,u[b[i]]);
}
//mexEvalString("Init");
mxFree(temp_NZE_R);
@ -527,6 +572,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
mexPrintf("pivot=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods);
#endif
pivot=(int*)mxMalloc(Size*periods*sizeof(int));
pivot_save=(int*)mxMalloc(Size*periods*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("pivota=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods);
#endif
@ -592,13 +638,13 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
//i=periods*Size*sizeof(*b);
//memset(b,0,i);
///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for(i=0; i< periods*Size;i++)
{
b[i]=0;
line_done[i]=0;
}
///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for(i=0; i< (periods+y_kmax+1)*Size;i++)
{
FNZE_C[i]=0;
@ -630,7 +676,9 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
if (eq!=it4->first.first.first+Size*t)
tmp_b=0;
eq=it4->first.first.first+Size*t;
lag=it4->first.second;
#ifdef PRINT_OUT
mexPrintf("=) eq=%d var=%d lag=%d t=%d\n",eq,var, lag, t);
mexPrintf("eq=%d, var=%d",eq,var);
mexEvalString("drawnow;");
#endif
@ -650,6 +698,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
#else
first=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem));
#endif
//mexPrintf("=> eq=%d var=%d lag=%d u=%d\n",eq,var, lag, it4->second+u_count_init*t);
first->NZE_C_N=NULL;
first->NZE_R_N=NULL;
first->u_index=it4->second+u_count_init*t;
@ -659,9 +708,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
/*if(eq==0 && var==0)
mexPrintf("alloc FNZE_R[0]=%x\n",first);*/
if (FNZE_R[eq]==NULL)
{
FNZE_R[eq]=first;
}
FNZE_R[eq]=first;
if (FNZE_C[var]==NULL)
FNZE_C[var]=first;
if (temp_NZE_R[eq]!=NULL)
@ -674,15 +721,30 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
mexPrintf("=> ");
#endif
}
else /*Build the additive terms ooutside the simulation periods related to the first lags and the las leads...*/
else /*Build the additive terms ooutside the simulation periods related to the first lags and the last leads...*/
{
if(lag<ti_y_kmin)
{
#ifdef PRINT_OUT
mexPrintf("nn ");
mexPrintf("tmp_b+=u[%d]*y[index_var[%d]]\n",it4->second+u_count_init*t,var+Size*(y_kmin+t));
mexPrintf("tmp_b+=u[%d](%f)*y[%d(%d)](%f)",it4->second+u_count_init*t,u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t),y[index_vara[var+Size*(y_kmin+t)]]);
mexEvalString("drawnow;");
mexPrintf("nn var=%d, Size=%d, t=%d, y_kmin=%d, y_kmax=%d\n", var, Size, t, y_kmin, y_kmax);
mexPrintf(" tmp_b+=u[%d]*y[index_var[%d]]\n", it4->second+u_count_init*t, var+Size*(y_kmin+t));
mexPrintf(" tmp_b+=u[%d](%f)*y[%d(%d)](f)\n", it4->second+u_count_init*t, u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t)/*,y[index_vara[var+Size*(y_kmin+t)]]*/);
mexEvalString("drawnow;");
#endif
tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
}
else
{
#ifdef PRINT_OUT
var -= Size;
mexPrintf("nn var=%d, Size=%d, t=%d, y_kmin=%d, y_kmax=%d\n", var, Size, t, y_kmin, y_kmax);
mexPrintf(" tmp_b+=u[%d]*y[index_var[%d]]\n", it4->second+u_count_init*t, var+Size*(y_kmin+t));
mexPrintf(" tmp_b+=u[%d](%f)*y[%d(%d)](f)\n", it4->second+u_count_init*t, u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t)/*,y[index_vara[var+Size*(y_kmin+t)]]*/);
mexEvalString("drawnow;");
#endif
tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
}
}
}
else /* ...and store it in the u vector*/
@ -692,9 +754,10 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<
#endif
b[eq]=it4->second+u_count_init*t;
u[b[eq]]+=tmp_b;
tmp_b = 0;
//mexPrintf("u[%d]=%f corr=%f\n",b[eq],u[b[eq]],tmp_b);
#ifdef PRINT_OUT
mexPrintf("=> b[%d]=%f\n", eq, u[b[eq]]);
mexPrintf("=> u[b[%d]=%d]=%f\n", eq, b[eq], u[b[eq]]);
mexEvalString("drawnow;");
#endif
}
@ -867,6 +930,7 @@ void SparseMatrix::End(int Size)
mxFree(b);
mxFree(line_done);
mxFree(pivot);
mxFree(pivot_save);
mxFree(pivotk);
mxFree(pivotv);
mxFree(pivotva);
@ -879,22 +943,13 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
#endif
)
{
long int i,j,/*nop=nop4/4*/nop=nop4/2, t, index_d, k;
long int i,j,nop=nop4/2, t, index_d, k;
double r=0.0;
bool OK=true;
t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s;
int *diff1, *diff2;
#ifdef MEM_ALLOC_CHK
mexPrintf("diff1=(int*)mxMalloc(%f)\n",double(nop)*double(sizeof(int))/1024);
#endif
diff1=(int*)mxMalloc(nop*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("diff2=(int*)mxMalloc(%f)\n",double(nop)*double(sizeof(int))/1024);
#endif
diff2=(int*)mxMalloc(nop*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
#endif
int max_save_ops_first=-1;
j=k=i=0;
while (i<nop4 && OK)
@ -932,16 +987,9 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
}
j++;
}
#ifdef PROFILER
if(OK)
mexPrintf("at %d same construction\n",beg_t);
else
mexPrintf("at %d different construction\n",beg_t);
mexEvalString("drawnow;");
#endif
// the same pivot for all remaining periods
if (OK)
///#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(j) schedule(dynamic)
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(j) schedule(dynamic)
for (i=beg_t;i<periods;i++)
{
for (j=0;j<Size;j++)
@ -952,17 +1000,9 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
}
if (OK)
{
/*#ifdef WRITE_u
long int i_toto=0;
fstream toto;
toto.open("compare_s.txt", std::ios::out);
#endif*/
if (max_save_ops_first>=u_count_alloc)
{
u_count_alloc+=5*u_count_alloc_save;
/*#ifdef MEM_ALLOC_CHK
mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))=%d, t=%d, omp_get_thread_num()=%d\n",u_count_alloc,t,omp_get_thread_num());
#endif*/
u=(double*)mxRealloc(u,u_count_alloc*sizeof(double));
if (!u)
{
@ -971,13 +1011,24 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
mexErrMsgTxt("Exit from Dynare");
}
}
for (t=1;t<periods-beg_t-max(y_kmax,y_kmin);t++)
for (t=1;t<periods-beg_t-y_kmax/*max(y_kmax,y_kmin)*/;t++)
{
i=j=0;
while (i<nop4)
{
save_op_s=(t_save_op_s*)(&(save_op[i]));
index_d=save_op_s->first+t*diff1[j];
if (index_d>u_count_alloc)
{
u_count_alloc+=2*u_count_alloc_save;
u=(double*)mxRealloc(u,u_count_alloc*sizeof(double));
if (!u)
{
mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n",u_count_alloc*sizeof(double));
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
}
switch (save_op_s->operat)
{
case IFLD :
@ -998,28 +1049,12 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
break;
}
j++;
if (index_d+3>=u_count_alloc)
{
u_count_alloc+=2*u_count_alloc_save;
u=(double*)mxRealloc(u,u_count_alloc*sizeof(double));
if (!u)
{
mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n",u_count_alloc*sizeof(double));
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
}
}
}
int t1=periods-beg_t-max(y_kmax,y_kmin);
#ifdef PROFILER
mexPrintf("first step done\n");
mexEvalString("drawnow;");
#endif
int t1=max(1,periods-beg_t-y_kmax);
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(t, i,j, save_op_s, index_d, r) schedule(dynamic)
for (t=t1;t<periods-beg_t;t++)
{
//mexPrintf("omp_in_parallel=%hd, omp_get_thread_num=%d, t=%d\n",omp_in_parallel(), omp_get_thread_num(), t);
i=j=0;
//#pragma omp ordered
while (i<nop4)
@ -1028,47 +1063,10 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
if (save_op_s->lag<((periods-beg_t)-t))
{
index_d=save_op_s->first+t*diff1[j];
switch (save_op_s->operat)
{
case IFLD :
r=u[index_d];
#ifdef PRINT_u
mexPrintf("FLD u[%d] (%f)\n",index_d,u[index_d]);
#endif
i+=2;
break;
case IFDIV :
u[index_d]/=r;
#ifdef PRINT_u
mexPrintf("FDIV u[%d](%f)/=r(%f)=(%f)\n",index_d,u[index_d],r,u[index_d]);
#endif
i+=2;
break;
case IFSUB :
u[index_d]-=u[save_op_s->second+t*diff2[j]]*r;
#ifdef PRINT_u
mexPrintf("FSUB u[%d]-=u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] );
#endif
i+=3;
break;
case IFLESS:
u[index_d]=-u[save_op_s->second+t*diff2[j]]*r;
#ifdef PRINT_u
mexPrintf("FLESS u[%d]=-u[%d](%f)*r(%f)=(%f)\n",index_d,save_op_s->second+t*diff2[j],u[save_op_s->second+t*diff2[j]],r,u[index_d] );
#endif
i+=3;
break;
}
if (index_d+3>=u_count_alloc)
if (index_d>u_count_alloc)
{
u_count_alloc+=2*u_count_alloc_save;
#ifdef MEM_ALLOC_CHK
mexPrintf("u=(double*)mxRealloc(u,u_count_alloc*sizeof(double))=%d, t=%d, omp_get_thread_num()=%d\n",u_count_alloc,t,omp_get_thread_num());
#endif
u=(double*)mxRealloc(u,u_count_alloc*sizeof(double));
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
#endif
if (!u)
{
mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n",u_count_alloc*sizeof(double));
@ -1076,6 +1074,25 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
mexErrMsgTxt("Exit from Dynare");
}
}
switch (save_op_s->operat)
{
case IFLD :
r=u[index_d];
i+=2;
break;
case IFDIV :
u[index_d]/=r;
i+=2;
break;
case IFSUB :
u[index_d]-=u[save_op_s->second+t*diff2[j]]*r;
i+=3;
break;
case IFLESS:
u[index_d]=-u[save_op_s->second+t*diff2[j]]*r;
i+=3;
break;
}
}
else
{
@ -1094,21 +1111,9 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i
j++;
}
}
#ifdef WRITE_u
toto.close();
filename+=" stopped";
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt(filename.c_str());
#endif
}
//mexPrintf("mxFree(diff1)\n");
mxFree(diff1);
//mexPrintf("mxFree(diff2)\n");
mxFree(diff2);
#ifdef PROFILER
mexPrintf("end of compare\n");
mexEvalString("drawnow;");
#endif
return OK;
}
@ -1238,7 +1243,7 @@ void
SparseMatrix::run_triangular(int nop_all,int *op_all)
{
int j=0;
mexPrintf("begining of run_triangular nop_all=%d\n",nop_all);
//mexPrintf("begining of run_triangular nop_all=%d\n",nop_all);
if (mem_mngr.swp_f)
{
bool OK=true;
@ -1592,7 +1597,8 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
//mexEvalString("drawnow;");
/*finding the max-pivot*/
double piv=piv_abs=0;
int nb_eq=At_Col(i, 0, &first);
//int nb_eq=At_Col(i, 0, &first);
int nb_eq=At_Col(i, &first);
//mexPrintf("nb_eq=%d\n",nb_eq);
//mexEvalString("drawnow;");
#ifdef MARKOVITZ
@ -1727,8 +1733,11 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
nb_eq=At_Col(i,&first);
NonZeroElem *first_piva;
int nb_var_piva=At_Row(pivj,&first_piva);
for (j=0;j<Size and first;j++)
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for (j=0;j<Size /*and first*/;j++)
{
if(first)
{
int row=first->r_index;
//mexPrintf("j=%d row=%d line_done[row]=%d\n",j, row, line_done[row]);
if (!line_done[row])
@ -1810,6 +1819,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
first=first->NZE_C_N;
/*first=first->NZE_R_N;*/
}
}
}
/*mexPrintf("before bcksub\n");
mexEvalString("drawnow;");*/
@ -1821,12 +1831,156 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax,
//res1bx=bksub( NULL, 1, y_size, slowc_lbx);
res1bx=simple_bksub(it_,Size,slowc_lbx);
//mexPrintf("End of simulate_NG\n");
End(Size);
return(0);
}
void
SparseMatrix::CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter)
{
const double epsilon=1e-7;
fstream SaveResult;
ostringstream out;
out << "Result" << iter;
SaveResult.open(out.str().c_str(), std::ios::in );
if (!SaveResult.is_open())
{
mexPrintf("Error : Can't open file \"%s\" for reading\n", "Result");
mexEvalString("st=fclose('all');clear all;");
mexErrMsgTxt("Exit from Dynare");
}
mexPrintf("Reading Result...");
int row, col;
SaveResult >> row;
mexPrintf("row=%d\n",row);
SaveResult >> col;
mexPrintf("col=%d\n",col);
//double G1a[row][col];
double G1a;
mexPrintf("Allocated\n");
NonZeroElem *first;
for(int j=0; j< col; j++)
{
mexPrintf("j=%d ",j);
int nb_equ=At_Col(j,&first);
mexPrintf("nb_equ=%d\n",nb_equ);
int line;
if(first)
line = first->r_index;
else
line = -9999999;
for(int i=0; i< row; i++)
{
SaveResult >> G1a;
if(line == i)
{
if(abs(u[first->u_index]/G1a-1)>epsilon)
mexPrintf("Problem at r=%d c=%d u[first->u_index]=%5.14f G1a[i][j]=%5.14f %f\n",i,j,u[first->u_index],G1a, u[first->u_index]/G1a-1);
first=first->NZE_C_N;
if(first)
line = first->r_index;
else
line = -9999999;
}
else
{
if(G1a!=0.0)
mexPrintf("Problem at r=%d c=%d G1a[i][j]=%f\n",i,j,G1a);
}
}
}
mexPrintf("G1a red done\n");
SaveResult >> row;
mexPrintf("row(2)=%d\n",row);
double B[row];
for(int i=0; i< row; i++)
SaveResult >> B[i];
SaveResult.close();
mexPrintf("done\n");
mexPrintf("Comparing...");
/*NonZeroElem *first;
for(int i=0;i<row;i++)
{
mexPrintf("i=%d ",i);
int nb_var=At_Row(i,&first);
mexPrintf("nb_var=%d\n",nb_var);
int column;
if(first)
column = first->c_index;
else
column = -9999999;
for(int j=0;j<col;j++)
{
if(column == j)
{
if(abs(u[first->u_index]-G1a[i][j])>epsilon)
mexPrintf("Problem at r=%d c=%d u[first->u_index]=%f G1a[i][j]=%f\n",i,j,u[first->u_index],G1a[i][j]);
first=first->NZE_R_N;
if(first)
column = first->c_index;
else
column = -9999999;
}
else
{
if(G1a[i][j]!=0)
mexPrintf("Problem at r=%d c=%d G1a[i][j]=%f\n",i,j,G1a[i][j]);
}
}
}*/
for(int i=0; i<row; i++)
{
if(abs(u[b[i]]+B[i])>epsilon)
mexPrintf("Problem at i=%d u[b[i]]=%f B[i]=%f\n",i,u[b[i]],B[i]);
}
}
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;
//std::map<std::pair<std::pair<int, int> ,int>, int> IM_i;
Init(periods, y_kmin, y_kmax, Size, IM_i);
NonZeroElem *first;
int cal_y = y_kmin*Size;
mexPrintf(" ");
for(int i=0; i<Size; i++)
mexPrintf(" %8d",i);
mexPrintf("\n");
for(int t=y_kmin; t<periods+y_kmin; t++)
{
mexPrintf("t=%5d",t);
for(int i=0; i<Size; i++)
mexPrintf(" %d %1.6f",t*y_size+index_vara[i], y[t*y_size+index_vara[i]]);
mexPrintf("\n");
}
for(int i=0;i<Size*periods;i++)
{
double res=0;
int pos = pivot[i];
mexPrintf("pos[%d]=%d",i,pos);
int nb_var = At_Row(pos, &first);
mexPrintf(" nb_var=%d\n",nb_var);
for(int j=0;j<nb_var; j++)
{
mexPrintf("(y[%d]=%f)*(u[%d]=%f)(r=%d, c=%d)\n",index_vara[first->c_index]+cal_y, y[index_vara[first->c_index]+cal_y], first->u_index, u[first->u_index], first->r_index, first->c_index);
res += y[index_vara[first->c_index]+cal_y]*u[first->u_index];
first=first->NZE_R_N;
}
double tmp_=res;
res += u[b[pos]];
if(abs(res)>epsilon)
mexPrintf("Error for equation %d => res=%f y[%d]=%f u[b[%d]]=%f somme(y*u)=%f\n",pos,res,pos,y[index_vara[pos]], pos, u[b[pos]], tmp_);
}
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
}
@ -1835,7 +1989,6 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{
/*Triangularisation at each period of a block using a simple gaussian Elimination*/
t_save_op_s *save_op_s;
int start_compare=y_kmin;
bool record=false;
int *save_op=NULL, *save_opa=NULL, *save_opaa=NULL;
long int nop=0, nopa=0;
@ -1845,9 +1998,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
int row, nb_var_piv, nb_var_sub, l_sub, sub_c_index, tmp_lag, l_piv, piv_c_index, tmp_u_count, lag;
NonZeroElem *first, *firsta, *first_sub, *first_piv, *first_suba;
double piv_abs, first_elem;
//SparseMatrix sparse_matrix;
//mexPrintf("->u_count=%d &u_count=%x\n",u_count,&u_count);
//mexPrintf("GNU version=%d\n",GNUVER);
if(start_compare==0)
start_compare=y_kmin;;
#ifdef RECORD_ALL
int save_u_count=u_count;
#endif
@ -1876,6 +2028,17 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
#endif
if (isnan(res1) || isinf(res1))
{
if(iter==0)
{
for(j=0;j<y_size; j++)
mexPrintf("variable %d at time %d = %f, %f\n",j+1, it_, y[j+it_*y_size], y[j+(it_+1)*y_size]);
mexPrintf("The initial values of endogenous variables are too far from the solution.\n");
mexPrintf("Change them!\n");
mexEvalString("drawnow;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
}
if (slowc_save<1e-8)
{
mexPrintf("slowc_save=%g\n", slowc_save);
@ -1883,7 +2046,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mexPrintf("variable %d at time %d = %f, %f\n",j+1, it_, y[j+it_*y_size], y[j+(it_+1)*y_size]);
mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx);
mexEvalString("drawnow;");
//mexEvalString("st=fclose('all');clear all;");
mexEvalString("st=fclose('all');clear all;");
filename+=" stopped";
mexErrMsgTxt(filename.c_str());
#ifdef DEBUG_EX
@ -1906,23 +2069,32 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
markowitz_c=markowitz_c_s;
alt_symbolic_count++;
}
if (((res1/res1a-1)>-0.3) && symbolic)
if (((res1/res1a-1)>-0.3) && symbolic && iter>0)
{
if (start_compare==y_kmin)
if(restart>2)
{
mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n");
symbolic=false;
alt_symbolic=true;
markowitz_c_s=markowitz_c;
markowitz_c=0;
}
else
{
mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied for a longer period.\n");
start_compare=min(tbreak_g,periods);
}
else
{
mexPrintf("Divergence or slowdown occured during simulation.\nIn the next iteration, pivoting method will be applied to all periods.\n");
symbolic=false;
alt_symbolic=true;
markowitz_c_s=markowitz_c;
markowitz_c=0;
restart++;
}
}
else
{
start_compare=y_kmin;
restart = 0;
}
res1a=res1;
if(print_it)
{
mexPrintf("-----------------------------------\n");
@ -1932,6 +2104,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mexPrintf(" abs. error=%.10e \n",double(res1));
mexPrintf("-----------------------------------\n");
}
//Print(Size, b);
if (cvg)
{
/*mexPrintf("End of simulate_NG1\n");
@ -1987,20 +2160,20 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
#ifdef RECORD_ALL
if (record_all && nop_all)
{
#ifdef PRINT_OUT
//#ifdef PRINT_OUT
mexPrintf("ShortInit\n");
mexEvalString("drawnow;");
#endif
//#endif
ShortInit(periods, y_kmin, y_kmax, Size, IM_i);
#ifdef PRINT_OUT
//#ifdef PRINT_OUT
mexPrintf("run_triangular\n");
mexEvalString("drawnow;");
#endif
//#endif
run_triangular(nop_all,save_op_all);
#ifdef PRINT_OUT
//#ifdef PRINT_OUT
mexPrintf("OK\n");
mexEvalString("drawnow;");
#endif
//#endif
}
else
#endif
@ -2023,12 +2196,26 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mexEvalString("drawnow;");
#endif
Init(periods, y_kmin, y_kmax, Size, IM_i);
/*ua = (double*)mxMalloc(u_count_init * periods*sizeof(double));
for(i=0; i< u_count_init * periods;i++)
ua[i] = u[i];*/
#ifdef PRINT_OUT
mexPrintf("done\n");
mexEvalString("drawnow;");
#endif
//Print(Size, b);
//CheckIt(y_size, y_kmin, y_kmax, Size, periods, iter);
//mexErrMsgTxt("Exit from Dynare");
/*for(int i=0; i<row; i++)
{
u[b[i]] = - u[b[i]];
}*/
for (int t=0;t<periods;t++)
{
//mexPrintf("t=%d periods=%d\n",t,periods);
#ifdef WRITE_u
if (!symbolic && ((periods-t)<=y_kmax))
{
@ -2084,7 +2271,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
#ifdef PRINT_OUT
mexPrintf("nb_eq=%d\n",nb_eq);
#endif
if /*(t<=y_kmin) */((symbolic && t<=start_compare) || !symbolic)
//mexPrintf("symbolic=%d t=%d start_compare=%d\n",symbolic, t, start_compare);
if ((symbolic && t<=start_compare) || !symbolic)
{
#ifdef MARKOVITZ
double piv_v[Size];
@ -2197,8 +2385,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
mexPrintf("pivot[%d]=%d\n",i,pivj);
mexEvalString("drawnow;");
#endif
if(iter>0 && t>start_compare)
{
if(pivot_save[i-Size]+Size!=pivj)
mexPrintf("At t=%d at line i=%d pivj=%d and pivot_save[i-Size]+Size=%d\n",t,i,pivj, pivot_save[i-Size]+Size);
}
pivot[i]=pivj;
pivot_save[i]=pivj;
pivotk[i]=pivk;
pivotv[i]=piv;
}
@ -2256,6 +2449,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nop_all+=2;
}
#endif
/*mexPrintf("piv_abs=%f\n",piv_abs);
mexEvalString("drawnow;");*/
if (piv_abs<eps)
{
mexPrintf("Error: singular system in Simulate_NG1\n");
@ -2267,9 +2462,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
/*divide all the non zeros elements of the line pivj by the max_pivot*/
int nb_var=At_Row(pivj,&first);
NonZeroElem* bb[nb_var];
/*mexPrintf("nb_var=%d\n",nb_var);
mexEvalString("drawnow;");*/
for(j=0;j<nb_var;j++)
{
bb[j]=first;
/*mexPrintf("j=%d",j);
mexPrintf(" first->NZE_R_N=%x\n",first->NZE_R_N);*/
first=first->NZE_R_N;
}
@ -2327,7 +2526,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
save_op_s->first=first->u_index;
save_op_s->lag=first->lag_index;
}
//nopi+=2;
//nop+=2; ///!!
}
#ifdef RECORD_ALL
else if (record_all)
@ -2400,6 +2599,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nb_eq=At_Col(i,&first);
NonZeroElem *first_piva;
int nb_var_piva=At_Row(pivj,&first_piva);
//mexPrintf("pivj=%d\n",pivj);
#ifdef PRINT_OUT
if(iter>0)
{
@ -2414,10 +2614,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
bc[j]=first;
first=first->NZE_C_N;
}
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) private(first, row, first_elem, nopa, save_op_s, nb_var_piv, nb_var_piva, first_piv, first_piva, first_sub, nb_var_sub, l_sub, l_piv, sub_c_index, piv_c_index, tmp_lag)
for (j=0;j<nb_eq;j++)
{
first=bc[j];
row=first->r_index;
/*mexPrintf("-------------------\n");
mexPrintf("j=%d line_done[row=%d]=%d\n",j,row, line_done[row]);*/
#ifdef PRINT_OUT
mexPrintf("t=%d, j=%d, line_done[%d]=%hd process=%d\n", t, j, row, line_done[row],omp_get_thread_num());
#endif
@ -2470,6 +2673,22 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
nop_all+=2;
}
#endif
/*mexPrintf("For equ=9\n");
int nb_var__=At_Row(9,&first_piv);
for(int uu=0; uu<nb_var__; uu++)
{
mexPrintf("-> first_piv->c_index=%d\n",first_piv->c_index);
first_piv=first_piv->NZE_R_N;
}
first_piv = first_piva;
mexPrintf("OK\n");
for(int uu=0; uu<nb_var_piva; uu++)
{
mexPrintf("-> first_piv->c_index=%d\n",first_piv->c_index);
first_piv=first_piv->NZE_R_N;
}*/
nb_var_piv=nb_var_piva;
first_piv=first_piva;
nb_var_sub=At_Row(row,&first_sub);
@ -2490,8 +2709,11 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
if (l_sub<nb_var_sub)
mexPrintf(" sub eq=%d lag=%d var=%d l0=%d",first_sub->r_index, first_sub->lag_index, first_sub->c_index,l_sub);
#endif
//mexPrintf("sub_c_index=%d piv_c_index=%d, l_sub=%d nb_var_sub=%d, l_piv=%d nb_var_piv=%d\n",sub_c_index, piv_c_index, l_sub, nb_var_sub, l_piv, nb_var_piv);
if (l_sub<nb_var_sub && (sub_c_index<piv_c_index || l_piv>=nb_var_piv))
{
//There is no nonzero element at line pivot for this column=> Nothing to do for the current element got to next column
//mexPrintf("Nothing\n");
first_sub=first_sub->NZE_R_N;
if (first_sub)
sub_c_index=first_sub->c_index;
@ -2501,6 +2723,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
else if (sub_c_index>piv_c_index || l_sub>=nb_var_sub)
{
// There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row
//mexPrintf("Insert\n");
tmp_u_count=Get_u();
#ifdef PROFILER
clock_t td0=clock();
@ -2567,6 +2791,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
{
if (i==sub_c_index)
{
//mexPrintf("Delete\n");
#ifdef PRINT_OUT
/*if(iter>0)
{
@ -2602,6 +2827,7 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
}
else
{
//mexPrintf("Substract\n");
#ifdef PRINT_OUT
mexPrintf(" u[%d]-=u[%d]*%f\n",first_sub->u_index,first_piv->u_index,double(first_elem));
#endif
@ -2955,6 +3181,13 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax
#endif
close_swp_file();
time00=clock();
if(tbreak_g==0)
tbreak_g=periods;
/*Check the solution*/
/*Check_the_Solution(periods, y_kmin, y_kmax, Size, ua, pivot, b);
mxFree(ua);*/
return(0);
}

35
mex/sources/simulate/SparseMatrix.hh
View File

@ -41,7 +41,7 @@ using namespace std;
struct t_save_op_s
{
short int lag, operat;
int first, second;
long int first, second;
};
const int IFLD =0;
@ -52,24 +52,10 @@ const int IFLDZ =4;
const int IFMUL =5;
const int IFSTP =6;
const int IFADD =7;
const double eps=1e-7;
const double eps=1e-10;
const double very_big=1e24;
const int alt_symbolic_count_max=1;
struct t_table_y
{
int index, nb;
int *u_index, *y_index;
};
struct t_table_u
{
t_table_u* pNext;
unsigned char type;
int index;
int op1, op2;
};
class SparseMatrix
{
@ -107,6 +93,8 @@ class SparseMatrix
void Delete_u(int pos);
void Clear_u();
void Print_u();
void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter);
void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int* pivot, int* b);
int complete(int beg_t, int Size, int periods, int *b);
double bksub( int tbreak, int last_period, int Size, double slowc_l
#ifdef PROFILER
@ -118,19 +106,11 @@ class SparseMatrix
void run_it(int nop_all,int *op_all);
void run_u_period1(int periods);
void close_swp_file();
void read_file_table_u(t_table_u **table_u, t_table_u **F_table_u, t_table_u **i_table_u, t_table_u **F_i_table_u, int *nb_table_u, bool i_to_do, bool shifting, int *nb_add_u_count, int y_kmin, int y_kmax, int u_size);
void read_file_table_y(t_table_y **table_y, t_table_y **i_table_y, int *nb_table_y, bool i_to_do, bool shifting, int y_kmin, int y_kmax, int u_size, int y_size);
void close_SaveCode();
stack<double> Stack;
int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
int middle_count_loop;
char type;
t_table_u *prologue_table_u, *first_table_u, *first_i_table_u, *middle_table_u, *middle_i_table_u, *last_table_u;
t_table_y *prologue_table_y, *first_table_y, *middle_table_y, *middle_i_table_y, *last_table_y;
t_table_u *F_prologue_table_u, *F_first_table_u, *F_first_i_table_u, *F_middle_table_u, *F_middle_i_table_u, *F_last_table_u;
fstream SaveCode;
string filename;
int max_u, min_u;
@ -142,7 +122,7 @@ class SparseMatrix
NonZeroElem **FNZE_R, **FNZE_C;
int nb_endo, u_count_init;
int *pivot, *pivotk;
int *pivot, *pivotk, *pivot_save;
double *pivotv, *pivotva;
int *b;
bool *line_done;
@ -154,7 +134,7 @@ class SparseMatrix
int u_count_alloc, u_count_alloc_save;
double markowitz_c_s;
double res1a;
long int nop_all, /*nopa_all,*/ nop1, nop2;
long int nop_all, nop1, nop2;
map<pair<pair<int, int> ,int>, int> IM_i;
protected:
double *u, *y, *ya;
@ -166,7 +146,8 @@ protected:
int u_count, tbreak_g;
int iter;
double *direction;
int start_compare;
int restart;
};

3
mex/sources/simulate/simulate.cc
View File

@ -178,10 +178,7 @@ main( int argc, const char* argv[] )
mxFree(ya);
if(direction)
mxFree(direction);
free(params);
}
#else

81
mex/sources/simulate/simulate.hh
View File

@ -20,87 +20,20 @@
#ifndef SIMULATE_HH_INCLUDED
#define SIMULATE_HH_INCLUDED
/*#include <stack>
#include <set>
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>*/
//#include <ctime>
//#include <string>
/*#include <map>
#include <algorithm>
#include "CodeInterpreter.hh"
#include "SymbolTableTypes.hh"*/
#include "Interpreter.hh"
//#include "Mem_Mngr.hh"
/*#include "LinBCG.hh"*/
#include "mex.h"
/*#include "ExprNode.hh"*/
//#define pow pow1
/*typedef struct Variable_l
{
int* Index;
};
typedef struct tBlock
{
int Size, Sized, Type, Max_Lead, Max_Lag, Simulation_Type, Nb_Lead_Lag_Endo;
int *Variable, *dVariable, *Equation;
int *variable_dyn_index, *variable_dyn_leadlag;
IM_compact *IM_lead_lag;
};
typedef struct tModel_Block
{
int Size;
tBlock * List;
};
#define MARKOVITZ
#define PRINT_OUT_p
//#define RECORD_ALL
//#define DEBUGC
//#define PRINT_OUT
//#define PRINT_OUT_y1
//#define PRINT_u
//#define PRINT_OUT_b
//#define PRINT_OUT_y
//#define DEBUG
//#define EXTENDED
//#define FLOAT
//#define WRITE_u
//#define MEMORY_LEAKS
//#define N_MX_ALLOC
//#define MEM_ALLOC_CHK
#define NEW_ALLOC
//#define PROFILER
//#ifdef EXTENDED
//typedef long double double;
//#else
//typedef double double;
//#endif
*/
#ifndef DEBUG_EX
#include "mex.h"
#else
#include "mex_interface.hh"
#endif
using namespace std;
/*std::multimap<std::pair<int,int>, int> var_in_equ_and_lag_i, equ_in_var_and_lag_i;
int *pivota=NULL, *save_op_all=NULL;
#ifdef RECORD_ALL
bool record_all=false;
#endif
long int nopa_all;
*/
int /*Per_y_, Per_u_, it_, */nb_row_x, nb_row_xd, u_size, y_size, x_size, y_kmin, y_kmax, y_decal;
int nb_row_x, nb_row_xd, u_size, y_size, x_size, y_kmin, y_kmax, y_decal;
int periods, maxit_;
double *params, markowitz_c, slowc, slowc_save;
double *u, *y, *x, *r, *g1, *g2, *ya;

15
mex/sources/simulate/testing/mex_interface.cc
View File

@ -5,9 +5,9 @@
using namespace std;
int
mexPrintf(string str, ...)
mexPrintf(/*string */const char *str, ...)
{
va_list vl;
/*va_list vl;
size_t found, found_p=0;
found=str.find_first_of("%");
va_start(vl,str);
@ -49,7 +49,16 @@ mexPrintf(string str, ...)
found=str.find_first_of("%",found_p);
}
printf(str.substr(found_p, str.size()-found_p+1).c_str());
return 0;
return 0;*/
va_list args;
int retval;
va_start (args, str);
retval = vprintf (str, args);
va_end (args);
return retval;
}
void

2
mex/sources/simulate/testing/mex_interface.hh
View File

@ -5,7 +5,7 @@
#include <stdarg.h>
using namespace std;
int mexPrintf(const string str, ...);
int mexPrintf(/*const string*/const char* str, ...);
void mexErrMsgTxt(const string str);
void* mxMalloc(int amount);
void* mxRealloc(void* to_extend, int amount);

36
preprocessor/BlockTriangular.cc
View File

@ -328,7 +328,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
ModelBlock->Block_List[count_Block].Type = type;
ModelBlock->Block_List[count_Block].Nb_Recursives = recurs_Size;
ModelBlock->Block_List[count_Block].Temporary_InUse = new temporary_terms_inuse_type();
ModelBlock->Block_List[count_Block].Chaine_Rule_Derivatives = new chaine_rule_derivatives_type();
ModelBlock->Block_List[count_Block].Chain_Rule_Derivatives = new chain_rule_derivatives_type();
ModelBlock->Block_List[count_Block].Temporary_InUse->clear();
ModelBlock->Block_List[count_Block].Simulation_Type = SimType;
ModelBlock->Block_List[count_Block].Equation = (int *) malloc(ModelBlock->Block_List[count_Block].Size * sizeof(int));
@ -341,7 +341,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
first_count_equ = *count_Equ;
tmp_var = (int *) malloc(size * sizeof(int));
tmp_endo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_other_endo = (int *) malloc(symbol_table.endo_nbr() * sizeof(int));
tmp_other_endo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size_other_endo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
tmp_size_exo = (int *) malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int));
@ -349,7 +349,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
memset(tmp_size_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_size, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
memset(tmp_other_endo, 0, symbol_table.endo_nbr()*sizeof(int));
memset(tmp_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int));
nb_lead_lag_endo = 0;
Lag_Endo = Lead_Endo = Lag_Other_Endo = Lead_Other_Endo = Lag_Exo = Lead_Exo = 0;
@ -438,7 +438,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int count_Block, Block
{
if (!tmp_variable_evaluated[j])
{
tmp_other_endo[j] = 1;
tmp_other_endo[incidencematrix.Model_Max_Lag + k]++;
tmp_nb_other_endo++;
}
if (k > 0 && k > Lead_Other_Endo)
@ -679,7 +679,7 @@ BlockTriangular::Free_Block(Model_Block *ModelBlock) const
delete ModelBlock->Block_List[blk].Temporary_Terms_in_Equation[i];
free(ModelBlock->Block_List[blk].Temporary_Terms_in_Equation);
delete (ModelBlock->Block_List[blk].Temporary_InUse);
delete ModelBlock->Block_List[blk].Chaine_Rule_Derivatives;
delete ModelBlock->Block_List[blk].Chain_Rule_Derivatives;
}
free(ModelBlock->Block_List);
free(ModelBlock);
@ -728,12 +728,13 @@ BlockTriangular::Equation_Type_determination(vector<BinaryOpNode *> &equations,
}
else
{
//vector<pair<int, NodeID> > List_of_Op_RHS;
//the equation could be normalized by a permutation of the rhs and the lhs
if (d_endo_variable == result.end()) //the equation is linear in the endogenous and could be normalized using the derivative
{
Equation_Simulation_Type = E_EVALUATE_S;
//cout << " gone normalized : ";
res = equations[eq]->normalizeLinearInEndoEquation(var, derivative->second);
res = equations[eq]->normalizeLinearInEndoEquation(var, derivative->second/*, List_of_Op_RHS*/);
/*res.second->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);
cout << " done\n";*/
}
@ -854,10 +855,11 @@ BlockTriangular::Reduce_Blocks_and_type_determination(int prologue, int epilogue
return (Type);
}
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>
BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
{
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int> Derivatives;
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> Derivatives;
Derivatives.clear();
int nb_endo = symbol_table.endo_nbr();
/*ModelBlock.Block_List[Blck].first_order_determinstic_simulation_derivatives = new*/
@ -876,11 +878,8 @@ BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
cout << "varr=" << varr << " eqr=" << eqr << " lag=" << lag << "\n";*/
if(IM[varr+eqr*nb_endo])
{
/*if(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)
{*/
bool OK = true;
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag)));
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr)));
if(its!=Derivatives.end())
{
if(its->second == 2)
@ -890,20 +889,21 @@ BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck)
if(OK)
{
if (ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
Derivatives[make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag))] = 1;
//It's a normalized equation, we have to recompute the derivative using chain rule derivative function*/
Derivatives[make_pair(make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr))] = 1;
else
Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varr, var), lag))] = 0;
//It's a feedback equation we can use the derivatives
Derivatives[make_pair(make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr))] = 0;
}
/*}
else if(eq<ModelBlock->Block_List[blck].Nb_Recursives and var<ModelBlock->Block_List[blck].Nb_Recursives)*/
if(var<ModelBlock->Block_List[blck].Nb_Recursives)
{
int eqs = ModelBlock->Block_List[blck].Equation[var];
for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; vars<ModelBlock->Block_List[blck].Size; vars++)
{
int varrs = ModelBlock->Block_List[blck].Variable[vars];
if(Derivatives.find(make_pair(make_pair(eqs, var), make_pair(make_pair(varrs, vars), lag)))!=Derivatives.end())
Derivatives[make_pair(make_pair(eqr, eq),make_pair(make_pair(varrs, vars), lag))] = 2;
//A new derivative need to be computed using the chain rule derivative function (a feedback variable appear in a recursive equation)
if(Derivatives.find(make_pair(make_pair(lag, make_pair(var, vars)), make_pair(eqs, varrs)))!=Derivatives.end())
Derivatives[make_pair(make_pair(lag, make_pair(eq, vars)), make_pair(eqr, varrs))] = 2;
}
}
}

6
preprocessor/BlockTriangular.hh
View File

@ -44,7 +44,7 @@ typedef vector<pair< int, int> > t_vtype;
typedef set<int> temporary_terms_inuse_type;
typedef vector<pair< pair<int, int>, pair<int, pair<int, int> > > > chaine_rule_derivatives_type;
typedef vector<pair< pair<int, pair<int, int> >, pair<int, int> > > chain_rule_derivatives_type;
//! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model
struct IM_compact
@ -73,7 +73,7 @@ struct Block
//temporary_terms_type *Temporary_terms;
temporary_terms_inuse_type *Temporary_InUse;
IM_compact *IM_lead_lag;
chaine_rule_derivatives_type *Chaine_Rule_Derivatives;
chain_rule_derivatives_type *Chain_Rule_Derivatives;
int Code_Start, Code_Length;
};
@ -118,7 +118,7 @@ public:
//! Frees the Model structure describing the content of each block
void Free_Block(Model_Block* ModelBlock) const;
map<pair<pair<int, int>, pair<pair<int, int>, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck);
void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector<BinaryOpNode *> &equations, t_etype &V_Equation_Type, map<pair<int, pair<int, int> >, NodeID> &first_order_endo_derivatives);

552
preprocessor/DynamicModel.cc
View File

@ -67,6 +67,18 @@ DynamicModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int la
code_file.write(&FLDZ, sizeof(FLDZ));
}
void
DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, int lag, map_idx_type &map_idx) const
{
map<pair<int, pair<int, int> >, NodeID>::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);
else
code_file.write(&FLDZ, sizeof(FLDZ));
}
void
DynamicModel::BuildIncidenceMatrix()
{
@ -105,7 +117,7 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
ostringstream tmp_output;
BinaryOpNode *eq_node;
first_derivatives_type::const_iterator it;
first_chaine_rule_derivatives_type::const_iterator it_chr;
first_chain_rule_derivatives_type::const_iterator it_chr;
ostringstream tmp_s;
temporary_terms.clear();
@ -134,14 +146,14 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
it->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
}
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
eq=it.first.second;
var=it.second.second.first;
lag=it.second.second.second;
it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
//it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
lag=it.first.first;
eq=it.first.second.first;
var=it.first.second.second;
it_chr=first_chain_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
//it_chr->second->writeChainRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
it_chr->second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, j, ModelBlock, ModelBlock->Block_List[j].Size-1, map_idx);
}
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
@ -199,14 +211,14 @@ DynamicModel::computeTemporaryTermsOrdered(Model_Block *ModelBlock)
it->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
}
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
eq=it.first.second;
var=it.second.second.first;
lag=it.second.second.second;
it_chr=first_chaine_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
//it_chr->second->writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
lag=it.first.first;
eq=it.first.second.first;
var=it.first.second.second;
it_chr=first_chain_rule_derivatives.find(make_pair(eq, make_pair( var, lag)));
//it_chr->second->writeChainRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
it_chr->second->collectTemporary_terms(temporary_terms, ModelBlock, j);
}
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
@ -337,12 +349,12 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
output << " g1 = spalloc(" << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)*ModelBlock->Periods
<< ", " << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)*(ModelBlock->Periods+ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1)
<< ", " << nze*ModelBlock->Periods << ");\n";
output << " g1_tmp_r = spalloc(" << (ModelBlock->Block_List[j].Nb_Recursives)
/*output << " g1_tmp_r = spalloc(" << (ModelBlock->Block_List[j].Nb_Recursives)
<< ", " << (ModelBlock->Block_List[j].Size)*(ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1)
<< ", " << nze << ");\n";
output << " g1_tmp_b = spalloc(" << (ModelBlock->Block_List[j].Size-ModelBlock->Block_List[j].Nb_Recursives)
<< ", " << (ModelBlock->Block_List[j].Size)*(ModelBlock->Block_List[j].Max_Lag+ModelBlock->Block_List[j].Max_Lead+1)
<< ", " << nze << ");\n";
<< ", " << nze << ");\n";*/
}
else
{
@ -409,10 +421,10 @@ DynamicModel::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const strin
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
tmp_output.str("");
if((ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) and (i<ModelBlock->Block_List[j].Nb_Recursives))
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
else
/*if((ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_BACKWARD or ModelBlock->Block_List[j].Simulation_Type!=EVALUATE_FORWARD) and (i<ModelBlock->Block_List[j].Nb_Recursives))
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
else*/
lhs->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms);
switch (ModelBlock->Block_List[j].Simulation_Type)
{
case EVALUATE_BACKWARD:
@ -598,16 +610,16 @@ end:
k=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
output << " g1(" << eqr+1 << ", "
<< varr+1 + m*(ModelBlock->Block_List[j].Size) << ") = ";
writeDerivative(output, eq, symbol_table.getID(eEndogenous, var), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k << ") " << var+1
<< ", equation=" << eq+1 << endl;
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];
output << " g1(" << eq+1 << ", "
<< var+1 + m*(ModelBlock->Block_List[j].Size) << ") = ";
writeDerivative(output, eqr, symbol_table.getID(eEndogenous, varr), k, oMatlabDynamicModelSparse, temporary_terms);
output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
<< "(" << k << ") " << varr+1
<< ", equation=" << eqr+1 << endl;
}
}
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
@ -652,15 +664,16 @@ end:
m=ModelBlock->Block_List[j].Max_Lag;
//cout << "\nDerivatives in Block " << j << "\n";
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
//Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
int eqr=it.first.first;
int eq=it.first.second;
int varr=it.second.first;
int var=it.second.second.first;
k=it.second.second.second;
//Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
k=it.first.first;
int eq=it.first.second.first;
int var=it.first.second.second;
int eqr=it.second.first;
int varr=it.second.second;
/*for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
@ -672,14 +685,14 @@ end:
{
if (varr>=ModelBlock->Block_List[j].Nb_Recursives)
{*/
output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
output << " g1(" << eq+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< var+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
writeChainRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;
<< ") " << varr+1
<< ", equation=" << eqr+1 << endl;
}
/*}
}
@ -698,17 +711,15 @@ end:
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i];*/
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
//Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
int eqr=it.first.first;
int eq=it.first.second;
int varr=it.second.first;
int var=it.second.second.first;
k=it.second.second.second;
//Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
k=it.first.first;
int eq=it.first.second.first;
int var=it.first.second.second;
int eqr=it.second.first;
int varr=it.second.second;
//bool derivative_exist;
ostringstream tmp_output;
@ -761,12 +772,12 @@ end:
output << " " << tmp_output.str();
writeChaineRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
writeChainRuleDerivative(output, eqr, varr, k, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr))
<< "(" << k << ") " << varr+1
<< ", equation=" << eqr+1 << endl;
<< ", equation=" << eqr+1 << " (" << eq+1 << ")" << endl;
}
//cout << " done\n";
/* }
@ -888,7 +899,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
int eqr;
};
int i,j,k,m, v, ModelBlock_Aggregated_Count, k0, k1;
int i,j,k,m, v;
string tmp_s;
ostringstream tmp_output;
ofstream code_file;
@ -897,9 +908,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
bool lhs_rhs_done;
Uff Uf[symbol_table.endo_nbr()];
map<NodeID, int> reference_count;
//map<int,int> ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number;
vector<int> feedback_variables;
int prev_Simulation_Type=-1;
bool file_open=false;
string main_name=file_name;
main_name+=".cod";
@ -914,20 +923,19 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const Model
k=temporary_terms.size();
code_file.write(reinterpret_cast<char *>(&k),sizeof(k));
ModelBlock_Aggregated_Count = ModelBlock->Size;
//For each block
for (j = 0; j < ModelBlock->Size ;j++)
{
feedback_variables.clear();
if (j>0)
code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
code_file.write(&FBEGINBLOCK, sizeof(FBEGINBLOCK));
v=ModelBlock->Block_List[j].Size;
v=ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives;
//cout << "v (Size) = " << v << "\n";
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=ModelBlock->Block_List[j].Simulation_Type;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
for (i=0; i < ModelBlock->Block_List[j].Size;i++)
int count_u;
for (i=ModelBlock->Block_List[j].Nb_Recursives; i < ModelBlock->Block_List[j].Size;i++)
{
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i]));
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].Equation[i]),sizeof(ModelBlock->Block_List[j].Equation[i]));
@ -936,8 +944,14 @@ ModelBlock_Aggregated_Count = ModelBlock->Size;
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE ||
ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_COMPLETE)
{
int u_count_int=0;
//cout << "ModelBlock->Block_List[j].Nb_Recursives = " << ModelBlock->Block_List[j].Nb_Recursives << "\n";
Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
//cout << "u_count_int=" << u_count_int << "\n";
code_file.write(reinterpret_cast<char *>(&ModelBlock->Block_List[j].is_linear),sizeof(ModelBlock->Block_List[j].is_linear));
v=block_triangular.ModelBlock->Block_List[j].IM_lead_lag[block_triangular.ModelBlock->Block_List[j].Max_Lag + block_triangular.ModelBlock->Block_List[j].Max_Lead].u_finish + 1;
//v=block_triangular.ModelBlock->Block_List[j].IM_lead_lag[block_triangular.ModelBlock->Block_List[j].Max_Lag + block_triangular.ModelBlock->Block_List[j].Max_Lead].u_finish + 1;
v = u_count_int ;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=symbol_table.endo_nbr();
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
@ -945,16 +959,12 @@ ModelBlock_Aggregated_Count = ModelBlock->Size;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=block_triangular.ModelBlock->Block_List[j].Max_Lead;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
int u_count_int=0;
Write_Inf_To_Bin_File(file_name, bin_basename, j, u_count_int,file_open,
ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE);
v=u_count_int;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
file_open=true;
//}
}
//For a block composed of a single equation determines whether we have to evaluate or to solve the equation
if (ModelBlock->Block_List[j].Size==1)
/*if (ModelBlock->Block_List[j].Size==1)
{
lhs_rhs_done=true;
eq_node = equations[ModelBlock->Block_List[j].Equation[0]];
@ -962,12 +972,15 @@ ModelBlock_Aggregated_Count = ModelBlock->Size;
rhs = eq_node->get_arg2();
}
else
lhs_rhs_done=false;
lhs_rhs_done=false;*/
// The equations
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
//The Temporary terms
temporary_terms_type tt2;
tt2.clear();
/*if (ModelBlock->Block_List[j].Temporary_Terms_in_Equation[i]->size())
output << " " << sps << "% //Temporary variables" << endl;*/
#ifdef DEBUGC
k=0;
#endif
@ -997,12 +1010,12 @@ ModelBlock_Aggregated_Count = ModelBlock->Size;
cout << "map_idx[" << (*it)->idx <<"]=" << ii->second << "\n";
}
#endif
if (!lhs_rhs_done)
{
/*if (!lhs_rhs_done)
{*/
eq_node = equations[ModelBlock->Block_List[j].Equation[i]];
lhs = eq_node->get_arg1();
rhs = eq_node->get_arg2();
}
/*}*/
switch (ModelBlock->Block_List[j].Simulation_Type)
{
evaluation:
@ -1042,11 +1055,9 @@ end:
int v=oMinus;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
code_file.write(&FSTPR, sizeof(FSTPR));
code_file.write(reinterpret_cast<char *>(&i), sizeof(i));
#ifdef CONDITION
if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE)
output << " condition[" << i << "]=0;\n";
#endif
v = i - ModelBlock->Block_List[j].Nb_Recursives;
//cout << "residual[" << v << "]\n";
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
}
}
code_file.write(&FENDEQU, sizeof(FENDEQU));
@ -1068,112 +1079,121 @@ end:
break;
case SOLVE_BACKWARD_COMPLETE:
case SOLVE_FORWARD_COMPLETE:
if(feedback_variable)
count_u = feedback_variables.size();
//cout << "todo SOLVE_COMPLETE\n";
for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
int u = feedback_variables.size();
for(i=0; i<ModelBlock->Block_List[j].Chaine_Rule_Derivatives->size();i++)
//Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
k=it.first.first;
int eq=it.first.second.first;;
int var=it.first.second.second;
int eqr=it.second.first;
int varr=it.second.second;
int v=ModelBlock->Block_List[j].Equation[eq];
if (!Uf[v].Ufl)
{
//Chaine_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, int>, pair<int, pair<int, int> > > it = ModelBlock->Block_List[j].Chaine_Rule_Derivatives->at(i);
int eqr=it.first.first;
int eq=it.first.second;
int varr=it.second.first;
int var=it.second.second.first;
int v=ModelBlock->Block_List[j].Equation[eqr];
k=it.second.second.second;
if (!Uf[v].Ufl)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
}
else
{
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=u;
Uf[v].Ufl->var=var;
compileDerivative(code_file, eq, var, 0, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
u++;
/*output << " g1(" << eqr+1-ModelBlock->Block_List[j].Nb_Recursives << ", "
<< varr+1-ModelBlock->Block_List[j].Nb_Recursives << ") = ";
writeChaineRuleDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms);
output << ";";
output << " %2 variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var))
<< "(" << k
<< ") " << var+1
<< ", equation=" << eq+1 << endl;*/
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
}
}
else
{
m=ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
else
{
int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index[i];
int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index[i];
int u=ModelBlock->Block_List[j].IM_lead_lag[m].us[i];
int eqr=ModelBlock->Block_List[j].IM_lead_lag[m].Equ[i];
int v=ModelBlock->Block_List[j].Equation[eqr];
if (!Uf[v].Ufl)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
}
else
{
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=u;
Uf[v].Ufl->var=var;
compileDerivative(code_file, eq, var, 0, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=count_u;
Uf[v].Ufl->var=varr;
compileChainRuleDerivative(code_file, eqr, varr, 0, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&count_u), sizeof(count_u));
count_u++;
}
//cout << "done SOLVE_COMPLETE\n";
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
code_file.write(&FLDR, sizeof(FLDR));
code_file.write(reinterpret_cast<char *>(&i), sizeof(i));
code_file.write(&FLDZ, sizeof(FLDZ));
int v=ModelBlock->Block_List[j].Equation[i];
for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
{
code_file.write(&FLDU, sizeof(FLDU));
code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
code_file.write(&FLDV, sizeof(FLDV));
char vc=eEndogenous;
code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc));
code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var));
int v1=0;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FBINARY, sizeof(FBINARY));
v1=oTimes;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FCUML, sizeof(FCUML));
}
Uf[v].Ufl=Uf[v].Ufl_First;
while (Uf[v].Ufl)
{
Uf[v].Ufl_First=Uf[v].Ufl->pNext;
free(Uf[v].Ufl);
if(i>=ModelBlock->Block_List[j].Nb_Recursives)
{
code_file.write(&FLDR, sizeof(FLDR));
v = i-ModelBlock->Block_List[j].Nb_Recursives;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FLDZ, sizeof(FLDZ));
int v=ModelBlock->Block_List[j].Equation[i];
for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
{
code_file.write(&FLDU, sizeof(FLDU));
code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
code_file.write(&FLDV, sizeof(FLDV));
char vc=eEndogenous;
code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc));
code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var));
int v1=0;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FBINARY, sizeof(FBINARY));
v1=oTimes;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FCUML, sizeof(FCUML));
}
Uf[v].Ufl=Uf[v].Ufl_First;
}
code_file.write(&FBINARY, sizeof(FBINARY));
v=oMinus;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&i), sizeof(i));
while (Uf[v].Ufl)
{
Uf[v].Ufl_First=Uf[v].Ufl->pNext;
free(Uf[v].Ufl);
Uf[v].Ufl=Uf[v].Ufl_First;
}
code_file.write(&FBINARY, sizeof(FBINARY));
v=oMinus;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FSTPU, sizeof(FSTPU));
v = i - ModelBlock->Block_List[j].Nb_Recursives;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
}
}
break;
case SOLVE_TWO_BOUNDARIES_COMPLETE:
case SOLVE_TWO_BOUNDARIES_SIMPLE:
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
count_u=ModelBlock->Block_List[j].Size - ModelBlock->Block_List[j].Nb_Recursives;
//cout << "todo SOLVE_TWO_BOUNDARIES\n";
//cout << "ModelBlock->Block_List[j].Nb_Recursives=" << ModelBlock->Block_List[j].Nb_Recursives << "\n";
for(i=0; i<ModelBlock->Block_List[j].Chain_Rule_Derivatives->size();i++)
{
//Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, pair<int, int> >, pair<int, int> > it = ModelBlock->Block_List[j].Chain_Rule_Derivatives->at(i);
k=it.first.first;
int eq=it.first.second.first;
int var=it.first.second.second;
int eqr=it.second.first;
int varr=it.second.second;
//cout << "k=" << k << " eq=" << eq << " (" << eq-ModelBlock->Block_List[j].Nb_Recursives << ") var=" << var << " (" << var-ModelBlock->Block_List[j].Nb_Recursives << ") eqr=" << eqr << " varr=" << varr << " count_u=" << count_u << "\n";
int v=ModelBlock->Block_List[j].Equation[eq];
/*m = ModelBlock->Block_List[j].Max_Lag + k;
int u=ModelBlock->Block_List[j].IM_lead_lag[m].u[i];*/
if(eq>=ModelBlock->Block_List[j].Nb_Recursives and var>=ModelBlock->Block_List[j].Nb_Recursives)
{
if (!Uf[v].Ufl)
{
Uf[v].Ufl=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl_First=Uf[v].Ufl;
}
else
{
Uf[v].Ufl->pNext=(Uff_l*)malloc(sizeof(Uff_l));
Uf[v].Ufl=Uf[v].Ufl->pNext;
}
Uf[v].Ufl->pNext=NULL;
Uf[v].Ufl->u=count_u;
Uf[v].Ufl->var=varr;
Uf[v].Ufl->lag=k;
//writeChainRuleDerivative(cout, eqr, varr, k, oMatlabDynamicModelSparse, /*map_idx*/temporary_terms);
//cout <<"\n";
compileChainRuleDerivative(code_file, eqr, varr, k, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&count_u), sizeof(count_u));
count_u++;
}
}
//cout << "done it SOLVE_TWO_BOUNDARIES\n";
/*for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
{
k=m-ModelBlock->Block_List[j].Max_Lag;
for (i=0;i<ModelBlock->Block_List[j].IM_lead_lag[m].size;i++)
@ -1200,50 +1220,47 @@ end:
compileDerivative(code_file, eq, var, k, map_idx);
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&u), sizeof(u));
#ifdef CONDITION
output << " if (fabs(condition[" << eqr << "])<fabs(u[" << u << "+Per_u_]))\n";
output << " condition[" << eqr << "]=u[" << u << "+Per_u_];\n";
#endif
}
}
}*/
for (i = 0;i < ModelBlock->Block_List[j].Size;i++)
{
code_file.write(&FLDR, sizeof(FLDR));
code_file.write(reinterpret_cast<char *>(&i), sizeof(i));
code_file.write(&FLDZ, sizeof(FLDZ));
int v=ModelBlock->Block_List[j].Equation[i];
for (Uf[v].Ufl=Uf[v].Ufl_First;Uf[v].Ufl;Uf[v].Ufl=Uf[v].Ufl->pNext)
{
code_file.write(&FLDU, sizeof(FLDU));
code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
code_file.write(&FLDV, sizeof(FLDV));
char vc=eEndogenous;
code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc));
int v1=Uf[v].Ufl->var;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
v1=Uf[v].Ufl->lag;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FBINARY, sizeof(FBINARY));
v1=oTimes;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FCUML, sizeof(FCUML));
}
Uf[v].Ufl=Uf[v].Ufl_First;
while (Uf[v].Ufl)
{
Uf[v].Ufl_First=Uf[v].Ufl->pNext;
free(Uf[v].Ufl);
if(i>=ModelBlock->Block_List[j].Nb_Recursives)
{
code_file.write(&FLDR, sizeof(FLDR));
v = i-ModelBlock->Block_List[j].Nb_Recursives;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FLDZ, sizeof(FLDZ));
int v=ModelBlock->Block_List[j].Equation[i];
for (Uf[v].Ufl=Uf[v].Ufl_First; Uf[v].Ufl; Uf[v].Ufl=Uf[v].Ufl->pNext)
{
code_file.write(&FLDU, sizeof(FLDU));
code_file.write(reinterpret_cast<char *>(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u));
code_file.write(&FLDV, sizeof(FLDV));
char vc=eEndogenous;
code_file.write(reinterpret_cast<char *>(&vc), sizeof(vc));
int v1=Uf[v].Ufl->var;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
v1=Uf[v].Ufl->lag;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FBINARY, sizeof(FBINARY));
v1=oTimes;
code_file.write(reinterpret_cast<char *>(&v1), sizeof(v1));
code_file.write(&FCUML, sizeof(FCUML));
}
Uf[v].Ufl=Uf[v].Ufl_First;
}
code_file.write(&FBINARY, sizeof(FBINARY));
v=oMinus;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FSTPU, sizeof(FSTPU));
code_file.write(reinterpret_cast<char *>(&i), sizeof(i));
#ifdef CONDITION
output << " if (fabs(condition[" << i << "])<fabs(u[" << i << "+Per_u_]))\n";
output << " condition[" << i << "]=u[" << i << "+Per_u_];\n";
#endif
while (Uf[v].Ufl)
{
Uf[v].Ufl_First=Uf[v].Ufl->pNext;
free(Uf[v].Ufl);
Uf[v].Ufl=Uf[v].Ufl_First;
}
code_file.write(&FBINARY, sizeof(FBINARY));
v=oMinus;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
code_file.write(&FSTPU, sizeof(FSTPU));
v = i - ModelBlock->Block_List[j].Nb_Recursives;
code_file.write(reinterpret_cast<char *>(&v), sizeof(v));
}
}
#ifdef CONDITION
for (m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++)
@ -1265,8 +1282,6 @@ end:
default:
break;
}
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
}
}
code_file.write(&FENDBLOCK, sizeof(FENDBLOCK));
@ -1415,7 +1430,34 @@ DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string
exit(EXIT_FAILURE);
}
u_count_int=0;
for (int m=0;m<=block_triangular.ModelBlock->Block_List[num].Max_Lead+block_triangular.ModelBlock->Block_List[num].Max_Lag;m++)
int Size = block_triangular.ModelBlock->Block_List[num].Size - block_triangular.ModelBlock->Block_List[num].Nb_Recursives;
for(int i=0; i<block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->size();i++)
{
//Chain_Rule_Derivatives.insert(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, lag))));
pair< pair<int, pair<int, int> >, pair<int, int> > it = block_triangular.ModelBlock->Block_List[num].Chain_Rule_Derivatives->at(i);
int k=it.first.first;
int eq=it.first.second.first;
int var_init=it.first.second.second;
/*int eqr=it.second.first;
int varr=it.second.second;*/
if(eq>=block_triangular.ModelBlock->Block_List[num].Nb_Recursives and var_init>=block_triangular.ModelBlock->Block_List[num].Nb_Recursives)
{
int v=eq-block_triangular.ModelBlock->Block_List[num].Nb_Recursives;
SaveCode.write(reinterpret_cast<char *>(&v), sizeof(v));
int var=it.first.second.second-block_triangular.ModelBlock->Block_List[num].Nb_Recursives + k * Size;
SaveCode.write(reinterpret_cast<char *>(&var), sizeof(var));
SaveCode.write(reinterpret_cast<char *>(&k), sizeof(k));
int u = u_count_int + Size;
SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
//cout << "eq=" << eq << " var=" << var << " k=" << k << " u=" << u << "\n";
u_count_int++;
}
}
/*for (int m=0;m<=block_triangular.ModelBlock->Block_List[num].Max_Lead+block_triangular.ModelBlock->Block_List[num].Max_Lag;m++)
{
int k1=m-block_triangular.ModelBlock->Block_List[num].Max_Lag;
for (j=0;j<block_triangular.ModelBlock->Block_List[num].IM_lead_lag[m].size;j++)
@ -1429,13 +1471,13 @@ DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string
SaveCode.write(reinterpret_cast<char *>(&u), sizeof(u));
u_count_int++;
}
}
}*/
if (is_two_boundaries)
{
for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
for (j=0;j<Size;j++)
{
int eqr1=j;
int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods
int varr=/*block_triangular.ModelBlock->Block_List[num].Size*/Size*(block_triangular.periods
+block_triangular.incidencematrix.Model_Max_Lead_Endo);
int k1=0;
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
@ -1446,12 +1488,12 @@ DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string
}
}
//cout << "u_count_int=" << u_count_int << "\n";
for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
int varr=block_triangular.ModelBlock->Block_List[num].Variable[j];
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
}
for (j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j];
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
@ -1548,7 +1590,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << tmp1.str();
tmp1.str("");
}
for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1;
@ -1735,7 +1777,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " g1=0;\n";
mDynamicModelFile << " r=0;\n";
tmp.str("");
for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
@ -1766,10 +1808,9 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
mDynamicModelFile << " g1=0;\n";
mDynamicModelFile << " r=0;\n";
tmp.str("");
for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
if (ik>=block_triangular.ModelBlock->Block_List[i].Nb_Recursives)
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
mDynamicModelFile << " y_index = [" << tmp.str() << "];\n";
int nze, m;
@ -1800,10 +1841,9 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
for (nze=0,m=0;m<=block_triangular.ModelBlock->Block_List[i].Max_Lead+block_triangular.ModelBlock->Block_List[i].Max_Lag;m++)
nze+=block_triangular.ModelBlock->Block_List[i].IM_lead_lag[m].size;
mDynamicModelFile << " y_index=[";
for (int ik=0;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ik<block_triangular.ModelBlock->Block_List[i].Size;ik++)
{
if (ik>=block_triangular.ModelBlock->Block_List[i].Nb_Recursives)
mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
mDynamicModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1;
}
mDynamicModelFile << " ];\n";
mDynamicModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n";
@ -2524,7 +2564,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
if (!no_tmp_terms)
computeTemporaryTermsOrdered(block_triangular.ModelBlock);
computeChaineRuleJacobian(block_triangular.ModelBlock);
computeChainRuleJacobian(block_triangular.ModelBlock);
}
else
if (!no_tmp_terms)
@ -2747,12 +2787,12 @@ DynamicModel::getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDExcepti
void
DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
DynamicModel::computeChainRuleJacobian(Model_Block *ModelBlock)
{
//cout << "computeChaineRuleJacobian\n";
//cout << "computeChainRuleJacobian\n";
//clock_t t1 = clock();
map<int, NodeID> recursive_variables;
first_chaine_rule_derivatives.clear();
first_chain_rule_derivatives.clear();
for(int blck = 0; blck<ModelBlock->Size; blck++)
{
//cout << "blck=" << blck << "\n";
@ -2760,7 +2800,7 @@ DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
if (ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_TWO_BOUNDARIES_SIMPLE)
{
//cout << "SOLVE_TWO_BOUNDARIES_COMPLETE \n";
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->clear();
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->clear();
for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
{
if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S)
@ -2769,27 +2809,27 @@ DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
recursive_variables[getDerivID(symbol_table.getID(eEndogenous, ModelBlock->Block_List[blck].Variable[i]), 0)] = equations[ModelBlock->Block_List[blck].Equation[i]];
}
//cout << "After recursive_alloc\n";
map<pair<pair<int, int >, pair<pair<int, int>,int> > , int > Derivatives = block_triangular.get_Derivatives(ModelBlock, blck);
for(map<pair<pair<int, int >, pair<pair<int, int>,int> > , int >::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++)
map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> Derivatives = block_triangular.get_Derivatives(ModelBlock, blck);
for(map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++)
{
int eqr = it->first.first.first;
int eq = it->first.first.second;
int varr = it->first.second.first.first;
int var = it->first.second.first.second;
int lag = it->first.second.second;
int lag = it->first.first.first;
int eq = it->first.first.second.first;
int var = it->first.first.second.second;
int eqr = it->first.second.first;
int varr = it->first.second.second;
int Deriv_type = it->second;
if(Deriv_type == 0)
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = first_derivatives[make_pair(eqr, getDerivID(symbol_table.getID(eEndogenous, varr), lag))];
else if (Deriv_type == 1)
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
else if (Deriv_type == 2)
{
if(ModelBlock->Block_List[blck].Equation_Type[eq] == E_EVALUATE_S and eq<ModelBlock->Block_List[blck].Nb_Recursives)
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = ModelBlock->Block_List[blck].Equation_Normalized[eq]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
else
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), lag), recursive_variables);
}
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(lag, make_pair(eq, var)), make_pair(eqr, varr)));
}
@ -2803,11 +2843,11 @@ DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
for(int var = ModelBlock->Block_List[blck].Nb_Recursives; var < ModelBlock->Block_List[blck].Size; var++)
{
int varr = ModelBlock->Block_List[blck].Variable[var];
NodeID d1 = equations[eqr]->getChaineRuleDerivative(recursive_variables, varr, lag);
NodeID d1 = equations[eqr]->getChainRuleDerivative(recursive_variables, varr, lag);
if (d1 == Zero)
continue;
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = d1;
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, lag))] = d1;
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(eqr, eq), make_pair(varr, make_pair(var, lag))));
}
}
}*/
@ -2818,7 +2858,7 @@ DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_BACKWARD_COMPLETE or ModelBlock->Block_List[blck].Simulation_Type==SOLVE_FORWARD_COMPLETE)
{
//cout << "SOLVE_FORWARD_SIMPLE \n";
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->clear();
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->clear();
for(int i = 0; i < ModelBlock->Block_List[blck].Nb_Recursives; i++)
{
if (ModelBlock->Block_List[blck].Equation_Type[i] == E_EVALUATE_S)
@ -2835,8 +2875,8 @@ DynamicModel::computeChaineRuleJacobian(Model_Block *ModelBlock)
NodeID d1 = equations[eqr]->getChainRuleDerivative(getDerivID(symbol_table.getID(eEndogenous, varr), 0), recursive_variables);
if (d1 == Zero)
continue;
first_chaine_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
ModelBlock->Block_List[blck].Chaine_Rule_Derivatives->push_back(make_pair( make_pair(eq, eqr), make_pair(var, make_pair(varr, 0))));
first_chain_rule_derivatives[make_pair(eqr, make_pair(varr, 0))] = d1;
ModelBlock->Block_List[blck].Chain_Rule_Derivatives->push_back(make_pair( make_pair(0, make_pair(eq, var)), make_pair(eqr, varr)));
}
}
}
@ -2928,13 +2968,15 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
paramsDerivsFile.close();
}
void
DynamicModel::writeChaineRuleDerivative(ostream &output, int eq, int var, int lag,
DynamicModel::writeChainRuleDerivative(ostream &output, int eqr, int varr, int lag,
ExprNodeOutputType output_type,
const temporary_terms_type &temporary_terms) const
{
map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chaine_rule_derivatives.find(make_pair(eq, make_pair(var, lag)));
if (it != first_chaine_rule_derivatives.end())
map<pair<int, pair<int, int> >, NodeID>::const_iterator it = first_chain_rule_derivatives.find(make_pair(eqr, make_pair(varr, lag)));
if (it != first_chain_rule_derivatives.end())
(it->second)->writeOutput(output, output_type, temporary_terms);
else
output << 0;

14
preprocessor/DynamicModel.hh
View File

@ -77,8 +77,8 @@ private:
//! Temporary terms for the file containing parameters dervicatives
temporary_terms_type params_derivs_temporary_terms;
typedef map< pair< int, pair< int, int> >, NodeID> first_chaine_rule_derivatives_type;
first_chaine_rule_derivatives_type first_chaine_rule_derivatives;
typedef map< pair< int, pair< int, int> >, NodeID> first_chain_rule_derivatives_type;
first_chain_rule_derivatives_type first_chain_rule_derivatives;
//! Writes dynamic model file (Matlab version)
@ -110,6 +110,8 @@ private:
void computeTemporaryTermsOrdered(Model_Block *ModelBlock);
//! Write derivative code of an equation w.r. to a variable
void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &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_type &map_idx) const;
virtual int computeDerivID(int symb_id, int lag);
//! Get the type corresponding to a derivation ID
@ -120,8 +122,8 @@ private:
int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
//! Compute the column indices of the dynamic Jacobian
void computeDynJacobianCols(bool jacobianExo);
//! Computes chaine rule derivatives of the Jacobian w.r. to endogenous variables
void computeChaineRuleJacobian(Model_Block *ModelBlock);
//! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
void computeChainRuleJacobian(Model_Block *ModelBlock);
//! Computes derivatives of the Jacobian w.r. to parameters
void computeParamsDerivatives();
//! Computes temporary terms for the file containing parameters derivatives
@ -137,8 +139,8 @@ private:
/*! Writes either (i+1,j+1) or [i+j*NNZDerivatives[1]] */
void hessianHelper(ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const;
//! Write chaine rule derivative of a recursive equation w.r. to a variable
void writeChaineRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
//! Write chain rule derivative of a recursive equation w.r. to a variable
void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
public:

12
tests/ferhat/fs2000.mod
View File

@ -32,8 +32,8 @@ psi = 0.787;
del = 0.02;
toto = [2 3];
//model(sparse_dll,cutoff=1e-17);
model(sparse, cutoff=0);
model(sparse_dll);
//model(sparse);
//model;
/*0*/ exp(gam+e_a) = dA ;
/*1*/ log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
@ -91,9 +91,9 @@ var e_m; stderr 0.005;
end;
options_.solve_tolf=1e-10;
options_.maxit_=100;
steady;
//options_.solve_tolf=1e-10;
options_.maxit_=10;
steady(block_mfs);
model_info;
//check;
shocks;
@ -105,7 +105,7 @@ end;
disp(toto(1,2));
simul(periods=200, method=lu);
simul(periods=20, method=lu);
//stoch_simul(periods=200,order=1);
rplot y;

6
tests/ferhat/ireland.mod
View File

@ -23,9 +23,9 @@ scy = 0.0040;
shy = 0.0015;
shc = 0.0010;
//model(sparse_dll);
//model(sparse,cutoff=0);
model(sparse);
//model(sparse_dll, cutoff=0);
model(sparse,cutoff=0);
//model(sparse);
//model;
exp(y) = exp(a)*exp(k(-1))^theta*exp(h)^(1-theta);
a = (1-rho)*aa+rho*a(-1)+e;

6
tests/ferhat/ls2003.mod
View File

@ -17,8 +17,8 @@ rho_ys = 0.9;
rho_pies = 0.7;
//model(sparse_dll);
model(sparse);
model(sparse_dll, markowitz=0, cutoff=0);
//model(sparse);
//model;
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
@ -93,7 +93,7 @@ values 0.5;
end;
//simul(periods=200,method=bicgstab);
simul(periods=200);
simul(periods=20);
rplot vv;
rplot ww;
rplot A;

33
tests/ferhat/multimod.mod
View File

@ -852,8 +852,8 @@ W0906=0.0800069594276;
W0907=0.147854375051;
W0908=0.206834342322;
W0909=-1;
//model(SPARSE_DLL,markowitz=2.0);
model(SPARSE,markowitz=2.0);
model(SPARSE_DLL,markowitz=2/*0.5*//*2.0*/);
//model(SPARSE,markowitz=2.0);
//model;
( log(US_CPI)-(log(US_CPI(-1)))) = US_CPI1*( log(US_PIM)-(log(US_PIM(-1))))+US_CPI2*( log(US_PGNP)-(log(US_PGNP(-1))))+(1-US_CPI1-US_CPI2)*log(US_CPI(-1)/US_CPI(-2))+RES_US_CPI ;
US_UNR_A = US_UNR_FE+US_UNR_1*100*log(US_GDP/US_GDP_FE)+US_UNR_2*(US_UNR(-1)-US_UNR_FE(-1))+RES_US_UNR_A ;
@ -911,7 +911,8 @@ model(SPARSE,markowitz=2.0);
US_PIM = (S0101*US_PXM+S0201*JA_PXM*JA_ER/JA_E96+S0301*GR_PXM*GR_ER/GR_E96+S0401*FR_PXM*FR_ER/FR_E96+S0501*IT_PXM*IT_ER/IT_E96+S0601*UK_PXM*UK_ER/UK_E96+S0701*CA_PXM*CA_ER/CA_E96+S0801*SI_PXM*SI_ER/SI_E96+S0901*RW_PXM*RW_ER/RW_E96)/(US_ER/US_E96)*(1+RES_US_PIM) ;
US_PIMA = US_PIM+T01*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/US_ER/US_IM ;
US_PIT = (US_IM*US_PIMA+US_IOIL*POIL/US_ER*US_E96+US_ICOM*PCOM/US_ER*US_E96)/US_IT ;
log(US_TFP_FE) = RES_US_TFP_FE ;
//log(US_TFP_FE) = RES_US_TFP_FE ;
US_TFP_FE = exp(RES_US_TFP_FE) ;
US_GDP_FE = US_TFP_FE*US_K^US_BETA*((1-US_UNR_FE/100)*US_LF)^(1-US_BETA) ;
US_LF = US_POP*US_PART/(1+US_DEM3) ;
US_CU = 100*US_GDP/US_GDP_FE ;
@ -980,7 +981,8 @@ model(SPARSE,markowitz=2.0);
JA_PIM = (S0102*US_PXM+S0202*JA_PXM*JA_ER/JA_E96+S0302*GR_PXM*GR_ER/GR_E96+S0402*FR_PXM*FR_ER/FR_E96+S0502*IT_PXM*IT_ER/IT_E96+S0602*UK_PXM*UK_ER/UK_E96+S0702*CA_PXM*CA_ER/CA_E96+S0802*SI_PXM*SI_ER/SI_E96+S0902*RW_PXM*RW_ER/RW_E96)/(JA_ER/JA_E96)*(1+RES_JA_PIM) ;
JA_PIMA = JA_PIM+T02*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/JA_ER/JA_IM ;
JA_PIT = (JA_IM*JA_PIMA+JA_IOIL*POIL/JA_ER*JA_E96+JA_ICOM*PCOM/JA_ER*JA_E96)/JA_IT ;
log(JA_TFP_FE) = RES_JA_TFP_FE ;
//log(JA_TFP_FE) = RES_JA_TFP_FE ;
JA_TFP_FE = exp(RES_JA_TFP_FE );
JA_GDP_FE = JA_TFP_FE*JA_K^JA_BETA*((1-JA_UNR_FE/100)*JA_LF)^(1-JA_BETA) ;
JA_LF = JA_POP*JA_PART/(1+JA_DEM3) ;
JA_CU = 100*JA_GDP/JA_GDP_FE ;
@ -1048,7 +1050,8 @@ model(SPARSE,markowitz=2.0);
GR_PIM = (S0103*US_PXM+S0203*JA_PXM*JA_ER/JA_E96+S0303*GR_PXM*GR_ER/GR_E96+S0403*FR_PXM*FR_ER/FR_E96+S0503*IT_PXM*IT_ER/IT_E96+S0603*UK_PXM*UK_ER/UK_E96+S0703*CA_PXM*CA_ER/CA_E96+S0803*SI_PXM*SI_ER/SI_E96+S0903*RW_PXM*RW_ER/RW_E96)/(GR_ER/GR_E96)*(1+RES_GR_PIM) ;
GR_PIMA = GR_PIM+T03*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/GR_ER/GR_IM ;
GR_PIT = (GR_IM*GR_PIMA+GR_IOIL*POIL/GR_ER*GR_E96+GR_ICOM*PCOM/GR_ER*GR_E96)/GR_IT ;
log(GR_TFP_FE) = RES_GR_TFP_FE ;
//log(GR_TFP_FE) = RES_GR_TFP_FE ;
GR_TFP_FE = exp(RES_GR_TFP_FE) ;
GR_GDP_FE = GR_TFP_FE*GR_K^GR_BETA*((1-GR_UNR_FE/100)*GR_LF)^(1-GR_BETA) ;
GR_LF = GR_POP*GR_PART/(1+GR_DEM3) ;
GR_CU = 100*GR_GDP/GR_GDP_FE ;
@ -1116,7 +1119,8 @@ model(SPARSE,markowitz=2.0);
FR_PIM = (S0104*US_PXM+S0204*JA_PXM*JA_ER/JA_E96+S0304*GR_PXM*GR_ER/GR_E96+S0404*FR_PXM*FR_ER/FR_E96+S0504*IT_PXM*IT_ER/IT_E96+S0604*UK_PXM*UK_ER/UK_E96+S0704*CA_PXM*CA_ER/CA_E96+S0804*SI_PXM*SI_ER/SI_E96+S0904*RW_PXM*RW_ER/RW_E96)/(FR_ER/FR_E96)*(1+RES_FR_PIM) ;
FR_PIMA = FR_PIM+T04*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/FR_ER/FR_IM ;
FR_PIT = (FR_IM*FR_PIMA+FR_IOIL*POIL/FR_ER*FR_E96+FR_ICOM*PCOM/FR_ER*FR_E96)/FR_IT ;
log(FR_TFP_FE) = RES_FR_TFP_FE ;
//log(FR_TFP_FE) = RES_FR_TFP_FE ;
FR_TFP_FE = exp(RES_FR_TFP_FE) ;
FR_GDP_FE = FR_TFP_FE*FR_K^FR_BETA*((1-FR_UNR_FE/100)*FR_LF)^(1-FR_BETA) ;
FR_LF = FR_POP*FR_PART/(1+FR_DEM3) ;
FR_CU = 100*FR_GDP/FR_GDP_FE ;
@ -1184,7 +1188,8 @@ model(SPARSE,markowitz=2.0);
IT_PIM = (S0105*US_PXM+S0205*JA_PXM*JA_ER/JA_E96+S0305*GR_PXM*GR_ER/GR_E96+S0405*FR_PXM*FR_ER/FR_E96+S0505*IT_PXM*IT_ER/IT_E96+S0605*UK_PXM*UK_ER/UK_E96+S0705*CA_PXM*CA_ER/CA_E96+S0805*SI_PXM*SI_ER/SI_E96+S0905*RW_PXM*RW_ER/RW_E96)/(IT_ER/IT_E96)*(1+RES_IT_PIM) ;
IT_PIMA = IT_PIM+T05*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/IT_ER/IT_IM ;
IT_PIT = (IT_IM*IT_PIMA+IT_IOIL*POIL/IT_ER*IT_E96+IT_ICOM*PCOM/IT_ER*IT_E96)/IT_IT ;
log(IT_TFP_FE) = RES_IT_TFP_FE ;
//log(IT_TFP_FE) = RES_IT_TFP_FE ;
IT_TFP_FE = exp(RES_IT_TFP_FE);
IT_GDP_FE = IT_TFP_FE*IT_K^IT_BETA*((1-IT_UNR_FE/100)*IT_LF)^(1-IT_BETA) ;
IT_LF = IT_POP*IT_PART/(1+IT_DEM3) ;
IT_CU = 100*IT_GDP/IT_GDP_FE ;
@ -1252,7 +1257,8 @@ model(SPARSE,markowitz=2.0);
UK_PIM = (S0106*US_PXM+S0206*JA_PXM*JA_ER/JA_E96+S0306*GR_PXM*GR_ER/GR_E96+S0406*FR_PXM*FR_ER/FR_E96+S0506*IT_PXM*IT_ER/IT_E96+S0606*UK_PXM*UK_ER/UK_E96+S0706*CA_PXM*CA_ER/CA_E96+S0806*SI_PXM*SI_ER/SI_E96+S0906*RW_PXM*RW_ER/RW_E96)/(UK_ER/UK_E96)*(1+RES_UK_PIM) ;
UK_PIMA = UK_PIM+T06*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/UK_ER/UK_IM ;
UK_PIT = (UK_IM*UK_PIMA+UK_IOIL*POIL/UK_ER*UK_E96+UK_ICOM*PCOM/UK_ER*UK_E96)/UK_IT ;
log(UK_TFP_FE) = RES_UK_TFP_FE ;
//log(UK_TFP_FE) = RES_UK_TFP_FE ;
UK_TFP_FE = exp(RES_UK_TFP_FE);
UK_GDP_FE = UK_TFP_FE*UK_K^UK_BETA*((1-UK_UNR_FE/100)*UK_LF)^(1-UK_BETA) ;
UK_LF = UK_POP*UK_PART/(1+UK_DEM3) ;
UK_CU = 100*UK_GDP/UK_GDP_FE ;
@ -1320,7 +1326,8 @@ model(SPARSE,markowitz=2.0);
CA_PIM = (S0107*US_PXM+S0207*JA_PXM*JA_ER/JA_E96+S0307*GR_PXM*GR_ER/GR_E96+S0407*FR_PXM*FR_ER/FR_E96+S0507*IT_PXM*IT_ER/IT_E96+S0607*UK_PXM*UK_ER/UK_E96+S0707*CA_PXM*CA_ER/CA_E96+S0807*SI_PXM*SI_ER/SI_E96+S0907*RW_PXM*RW_ER/RW_E96)/(CA_ER/CA_E96)*(1+RES_CA_PIM) ;
CA_PIMA = CA_PIM+T07*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/CA_ER/CA_IM ;
CA_PIT = (CA_IM*CA_PIMA+CA_IOIL*POIL/CA_ER*CA_E96+CA_ICOM*PCOM/CA_ER*CA_E96)/CA_IT ;
log(CA_TFP_FE) = RES_CA_TFP_FE ;
//log(CA_TFP_FE) = RES_CA_TFP_FE ;
CA_TFP_FE = exp(RES_CA_TFP_FE) ;
CA_GDP_FE = CA_TFP_FE*CA_K^CA_BETA*((1-CA_UNR_FE/100)*CA_LF)^(1-CA_BETA) ;
CA_LF = CA_POP*CA_PART/(1+CA_DEM3) ;
CA_CU = 100*CA_GDP/CA_GDP_FE ;
@ -1388,7 +1395,8 @@ model(SPARSE,markowitz=2.0);
SI_PIM = (S0108*US_PXM+S0208*JA_PXM*JA_ER/JA_E96+S0308*GR_PXM*GR_ER/GR_E96+S0408*FR_PXM*FR_ER/FR_E96+S0508*IT_PXM*IT_ER/IT_E96+S0608*UK_PXM*UK_ER/UK_E96+S0708*CA_PXM*CA_ER/CA_E96+S0808*SI_PXM*SI_ER/SI_E96+S0908*RW_PXM*RW_ER/RW_E96)/(SI_ER/SI_E96)*(1+RES_SI_PIM) ;
SI_PIMA = SI_PIM+T08*(WTRADE-TRDE*GREAL^TME*US_INFL^TME)/SI_ER/SI_IM ;
SI_PIT = (SI_IM*SI_PIMA+SI_IOIL*POIL/SI_ER*SI_E96+SI_ICOM*PCOM/SI_ER*SI_E96)/SI_IT ;
log(SI_TFP_FE) = RES_SI_TFP_FE ;
//log(SI_TFP_FE) = RES_SI_TFP_FE ;
SI_TFP_FE = exp(RES_SI_TFP_FE);
SI_GDP_FE = SI_TFP_FE*SI_K^SI_BETA*((1-SI_UNR_FE/100)*SI_LF)^(1-SI_BETA) ;
SI_LF = SI_POP*SI_PART/(1+SI_DEM3) ;
SI_CU = 100*SI_GDP/SI_GDP_FE ;
@ -2742,7 +2750,7 @@ end;
options_.slowc = 1.0;
options_.dynatol = 1e-4;
options_.maxit_ = 100;
options_.maxit_ = 5;
simul(periods=80,datafile=mark3);
@ -2758,7 +2766,7 @@ var US_G;
periods 1;
values 4330.714737;
end;
simul(periods=80);
/*simul(periods=10);
oo_.endo_simul_sim=oo_.endo_simul;
@ -2766,3 +2774,4 @@ oo_.endo_simul=(oo_.endo_simul_sim./oo_.endo_simul_ss-1)*100;
rplot WTRADER;
rplot US_GDP;
*/

4
tests/ferhat/ramst.mod
View File

@ -10,8 +10,8 @@ bet=0.05;
aa=0.5;
model(sparse);
//model(sparse_dll);
//model(sparse);
model(sparse_dll);
//model;
//s = aa*x*k(-1)^alph - c;
c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);// + 0.00000001*s;