diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index 056db63a5..0a72135ac 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -66,7 +66,7 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i % along with Dynare. If not, see . - global oo_ M_; + global oo_ M_ options_; Blck_size=size(y_index_eq,2); g2 = []; g3 = []; @@ -90,9 +90,9 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i g1=spalloc( Blck_size, Blck_size, nze); while ~(cvg==1 | iter>maxit_), if(is_dynamic) - [r, g1] = feval(fname, y, x, params, it_, 0, g1, g2, g3); + [r, g1, g2, g3, b] = feval(fname, y, x, params, it_, 0, g1, g2, g3); else - [r, g1] = feval(fname, y, x, params, 0); + [r, g1, g2, g3, b] = feval(fname, y, x, params, 0); end; if(~isreal(r)) max_res=(-(max(max(abs(r))))^2)^0.5; @@ -100,10 +100,16 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i max_res=max(max(abs(r))); end; if(verbose==1) - disp(['iteration : ' int2str(iter) ' => ' num2str(max_res)]); - disp([M_.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]); + disp(['iteration : ' int2str(iter) ' => ' num2str(max_res) ' time = ' int2str(it_)]); + if(is_dynamic) + disp([M_.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])]); + else + disp([M_.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])]); + end; end; - if(is_linear) + if(~isreal(max_res) | isnan(max_res)) + cvg = 0; + elseif(is_linear & iter>0) cvg = 1; else cvg=(max_res 0 + info = 0; + else + info = 1; + end + elseif(~is_dynamic & options_.solve_algo==2) + lambda=1; + stpmx = 100 ; + stpmax = stpmx*max([sqrt(y'*y);size(y_index_eq,2)]); + nn=1:size(y_index_eq,2); + g = (r'*g1)'; + f = 0.5*r'*r; + p = -g1\r ; + [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, 0); + elseif(~is_dynamic & options_.solve_algo==3) + [yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1); + y(y_index_eq) = yn; + elseif((simulation_method==0 & is_dynamic) | (~is_dynamic & options_.solve_algo==1)), + dx = g1\r; + ya = ya - lambda*dx; if(is_dynamic) y(it_,y_index_eq) = ya'; else y(y_index_eq) = ya; end; - elseif(simulation_method==2), + elseif(simulation_method==2 & is_dynamic), flag1=1; while(flag1>0) [L1, U1]=luinc(g1,luinc_tol); [za,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1); + %[za,flag1] = gmres(-g1,b',Blck_size,1e-6,Blck_size,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]); @@ -204,11 +241,12 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i end; end; end; - elseif(simulation_method==3), + elseif(simulation_method==3 & is_dynamic), flag1=1; while(flag1>0) [L1, U1]=luinc(g1,luinc_tol); [za,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1); + %[za,flag1] = bicgstab(-g1,b',1e-7,Blck_size,L1,U1); if (flag1>0 | reduced) if(flag1==1) disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]); @@ -229,6 +267,14 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i end; end; end; + else + disp('unknown option : '); + if(is_dynamic) + disp(['options_.simulation_method = ' num2str(simulation_method) ' not implemented']); + else + disp(['options_.solve_algo = ' num2str(options_.solve_algo) ' not implemented']); + end; + return; end; iter=iter+1; end @@ -254,8 +300,15 @@ function y = solve_one_boundary(fname, y, x, params, y_index_eq, nze, periods, i oo_.deterministic_simulation.status = 1; oo_.deterministic_simulation.error = max_res; oo_.deterministic_simulation.iterations = iter; - oo_.deterministic_simulation.block(Block_Num).status = 1;% Convergency failed. + oo_.deterministic_simulation.block(Block_Num).status = 1; oo_.deterministic_simulation.block(Block_Num).error = max_res; oo_.deterministic_simulation.block(Block_Num).iterations = iter; end; - return; \ No newline at end of file + return; + +function [err, G]=local_fname(yl, x, params, y, y_index_eq, fname, is_csolve) + y(y_index_eq) = yl; + [err, G] = feval(fname, y, x, params, 0); + if(is_csolve) + G = full(G); + end; diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 1c324aef5..a652caa03 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -32,7 +32,7 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ % - 3 BicGStab % % OUTPUTS -% y [matrix] All endogenous variables of the model +% y [matrix] All endogenous variables of the model % % ALGORITHM % Newton with LU or GMRES or BicGstab @@ -72,7 +72,6 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); reduced = 0; while ~(cvg==1 | iter>maxit_), - %fname [r, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, g1, g2, g3, y_kmin, Blck_size); g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size); b = b' -g1(:, 1+(y_kmin-y_kmin_l)*Blck_size:y_kmin*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)'-g1(:, (periods+y_kmin)*Blck_size+1:(periods+y_kmin+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)'; @@ -81,113 +80,117 @@ function y = solve_two_boundaries(fname, y, x, params, y_index, nze, periods, y_ else max_res=max(max(abs(r))); end; - if(iter>0) - if(~isreal(max_res) | isnan(max_res) | (max_resa1)) - if(isnan(max_res)) - detJ=det(g1aa); - if(abs(detJ)<1e-7) - max_factor=max(max(abs(g1aa))); - ze_elem=sum(diag(g1aa)1e-8) - lambda=lambda/2; - reduced = 1; - disp(['reducing the path length: lambda=' num2str(lambda,'%f')]); - y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)'; - continue; - else - if(cutoff == 0) - fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_".\n',Block_Num, iter); - else - fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_" or set "cutoff=0" in model options.\n',Block_Num, iter); - end; - oo_.deterministic_simulation.status = 0; - oo_.deterministic_simulation.error = max_res; - oo_.deterministic_simulation.iterations = iter; - oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed. - oo_.deterministic_simulation.block(Block_Num).error = max_res; - oo_.deterministic_simulation.block(Block_Num).iterations = iter; - return; - end; - else - if(lambda<1) - lambda=max(lambda*2, 1); - end; - end; - end; - ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)'; - ya_save=ya; - g1aa=g1a; - ba=b; - max_resa=max_res; - if(simulation_method==0), - dx = g1a\b- ya; - ya = ya + lambda*dx; - y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; - elseif(simulation_method==2), - flag1=1; - while(flag1>0) - [L1, U1]=luinc(g1a,luinc_tol); - [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1); - if (flag1>0 | reduced) - if(flag1==1) - disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); - elseif(flag1==2) - disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); - elseif(flag1==3) - disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); - end; - luinc_tol = luinc_tol/10; - reduced = 0; - else - dx = za - ya; - ya = ya + lambda*dx; - y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; - end; - end; - elseif(simulation_method==3), - flag1=1; - while(flag1>0) - [L1, U1]=luinc(g1a,luinc_tol); - [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1); - if (flag1>0 | reduced) - if(flag1==1) - disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); - elseif(flag1==2) - disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); - elseif(flag1==3) - disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); - end; - luinc_tol = luinc_tol/10; - reduced = 0; - else - dx = za - ya; - ya = ya + lambda*dx; - y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; - end; - end; - end; - if(is_linear) - cvg = 1; + if(~isreal(max_res) | isnan(max_res)) + cvg = 0; + elseif(is_linear & iter>0) + cvg = 1; else - cvg=(max_res0) + if(~isreal(max_res) | isnan(max_res) | (max_resa1)) + if(isnan(max_res)) + detJ=det(g1aa); + if(abs(detJ)<1e-7) + max_factor=max(max(abs(g1aa))); + ze_elem=sum(diag(g1aa)1e-8) + lambda=lambda/2; + reduced = 1; + disp(['reducing the path length: lambda=' num2str(lambda,'%f')]); + y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)'; + continue; + else + if(cutoff == 0) + fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_".\n',Block_Num, iter); + else + fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.maxit_" or set "cutoff=0" in model options.\n',Block_Num, iter); + end; + oo_.deterministic_simulation.status = 0; + oo_.deterministic_simulation.error = max_res; + oo_.deterministic_simulation.iterations = iter; + oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed. + oo_.deterministic_simulation.block(Block_Num).error = max_res; + oo_.deterministic_simulation.block(Block_Num).iterations = iter; + return; + end; + else + if(lambda<1) + lambda=max(lambda*2, 1); + end; + end; + end; + ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)'; + ya_save=ya; + g1aa=g1a; + ba=b; + max_resa=max_res; + if(simulation_method==0), + dx = g1a\b- ya; + ya = ya + lambda*dx; + y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; + elseif(simulation_method==2), + flag1=1; + while(flag1>0) + [L1, U1]=luinc(g1a,luinc_tol); + [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1); + if (flag1>0 | reduced) + if(flag1==1) + disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); + elseif(flag1==2) + disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); + elseif(flag1==3) + disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); + end; + luinc_tol = luinc_tol/10; + reduced = 0; + else + dx = za - ya; + ya = ya + lambda*dx; + y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; + end; + end; + elseif(simulation_method==3), + flag1=1; + while(flag1>0) + [L1, U1]=luinc(g1a,luinc_tol); + [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1); + if (flag1>0 | reduced) + if(flag1==1) + disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block' num2str(Block_Size,'%3d')]); + elseif(flag1==2) + disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Size,'%3d')]); + elseif(flag1==3) + disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]); + end; + luinc_tol = luinc_tol/10; + reduced = 0; + else + dx = za - ya; + ya = ya + lambda*dx; + y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)'; + end; + end; + end; + iter=iter+1; + disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]); + end + end; if (iter>maxit_) disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]); oo_.deterministic_simulation.status = 0; diff --git a/mex/2007a/simulate.dll b/mex/2007a/simulate.dll index bb3b5e74a..fd6ea7d91 100755 Binary files a/mex/2007a/simulate.dll and b/mex/2007a/simulate.dll differ diff --git a/mex/2007b/simulate.mexw32 b/mex/2007b/simulate.mexw32 index 466daec81..a1d968e8a 100755 Binary files a/mex/2007b/simulate.mexw32 and b/mex/2007b/simulate.mexw32 differ diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc index 43cb9e4ec..36f63c042 100644 --- a/mex/sources/simulate/Interpreter.cc +++ b/mex/sources/simulate/Interpreter.cc @@ -106,7 +106,7 @@ Interpreter::compute_block_time() /*throw(EvalException)*/ var=get_code_int lag=get_code_int #ifdef DEBUGC - if(var==6) + if(var==650 or var==643 or var==628) { mexPrintf(" FLD x[%d, time=%d, var=%d, lag=%d]=%f\n",it_+lag+var*nb_row_x,it_,var,lag,x[it_+lag+var*nb_row_x]); mexEvalString("drawnow;"); @@ -194,13 +194,13 @@ Interpreter::compute_block_time() /*throw(EvalException)*/ #endif y[(it_+lag)*y_size+var] = Stack.top(); #ifdef DEBUGC - if(var==153) + if(var==557 || var==558) { mexPrintf(" FSTP y[var=%d,time=%d,lag=%d,%d]=%f\n",var,it_,lag,(it_+lag)*y_size+var,y[(it_+lag)*y_size+var]); mexEvalString("drawnow;"); } - mexPrintf("%f\n",y[(it_+lag)*y_size+var]); - mexEvalString("drawnow;"); + /*mexPrintf("%f\n",y[(it_+lag)*y_size+var]); + mexEvalString("drawnow;");*/ #endif Stack.pop(); break; @@ -686,6 +686,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba set_code_pointer(begining); Per_y_=it_*y_size; compute_block_time(); + mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_SIMPLE : OK\n"); + mexEvalString("drawnow;"); y[Per_y_+Block_Contain[0].Variable] += -r[0]/g1[0]; double rr; rr=r[0]/(1+y[Per_y_+Block_Contain[0].Variable]); @@ -705,7 +707,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mxFree(g1); mxFree(r); break; - case SOLVE_TWO_BOUNDARIES_SIMPLE : + /*case SOLVE_TWO_BOUNDARIES_SIMPLE : #ifdef DEBUGC mexPrintf("SOLVE_TWO_BOUNDARIES_SIMPLE\n"); #endif @@ -777,17 +779,17 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba cvg=(max_res in SOLVE_FORWARD_COMPLETE : OK\n"); + mexEvalString("drawnow;"); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); res2=0; res1=0; @@ -833,8 +837,17 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { set_code_pointer(begining); Per_y_=it_*y_size; + iter = 0; + res1=res2=max_res=0; + /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE before compute_block_time OK\n"); + mexEvalString("drawnow;");*/ compute_block_time(); - Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + /*mexPrintf("Compute_block_time=> in SOLVE_FORWARD_COMPLETE : %d OK\n",it_); + mexEvalString("drawnow;");*/ + //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + //Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, /*true*/false, cvg, iter); + //simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter); } } mxFree(g1); @@ -866,6 +879,8 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba { set_code_pointer(begining); compute_block_time(); + mexPrintf("Compute_block_time=> in SOLVE_BACKWARD_COMPLETE : OK\n"); + mexEvalString("drawnow;"); Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); res2=0; res1=0; @@ -897,13 +912,14 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba set_code_pointer(begining); Per_y_=it_*y_size; compute_block_time(); - Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 0, false, iter); + Direct_Simulate(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, 1, false, iter); } } mxFree(g1); mxFree(r); mxFree(u); break; + case SOLVE_TWO_BOUNDARIES_SIMPLE : case SOLVE_TWO_BOUNDARIES_COMPLETE: #ifdef DEBUGC mexPrintf("SOLVE_TWO_BOUNDARIES_COMPLETE\n"); @@ -948,11 +964,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba mexPrintf("u=%x\n",u); #endif Read_SparseMatrix(bin_basename, size, periods, y_kmin, y_kmax); - //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)); + //g1=(double*)mxMalloc(size*size*sizeof(double)); r=(double*)mxMalloc(size*sizeof(double)); y_save=(double*)mxMalloc(y_size*sizeof(double)*(periods+y_kmax+y_kmin)); #ifdef DEBUGC @@ -1093,8 +1108,9 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba //mexErrMsgTxt("End of simulate"); #endif - mxFree(g1); + //mxFree(g1); mxFree(r); + mxFree(y_save); mxFree(u); mxFree(index_vara); memset(direction,0,size_of_direction); @@ -1113,8 +1129,6 @@ Interpreter::compute_blocks(string file_name, string bin_basename) { ifstream CompiledCode; int Code_Size, var; - /*mexPrintf("compute_blocks %s\n",file_name.c_str()); - mexEvalString("drawnow;");*/ //First read and store inn memory the code CompiledCode.open((file_name + ".cod").c_str(),std::ios::in | std::ios::binary| std::ios::ate); if (!CompiledCode.is_open()) @@ -1132,8 +1146,6 @@ Interpreter::compute_blocks(string file_name, string bin_basename) mexEvalString("drawnow;"); #endif CompiledCode.seekg(std::ios::beg); - /*Code_Size=CompiledCode.tellg(); - mexPrintf("Code_Size=%d\n",Code_Size);*/ Code=(char*)mxMalloc(Code_Size); CompiledCode.seekg(0); CompiledCode.read(reinterpret_cast(Code), Code_Size); diff --git a/mex/sources/simulate/Mem_Mngr.cc b/mex/sources/simulate/Mem_Mngr.cc index 41e7d6c5f..a3c2a2e3d 100644 --- a/mex/sources/simulate/Mem_Mngr.cc +++ b/mex/sources/simulate/Mem_Mngr.cc @@ -61,25 +61,32 @@ NonZeroElem* Mem_Mngr::mxMalloc_NZE() { int i; - if (!Chunk_Stack.empty()) + if (!Chunk_Stack.empty()) /*An unused block of memory available inside the heap*/ { NonZeroElem* p1 = Chunk_Stack.back(); Chunk_Stack.pop_back(); return(p1); } - else if (CHUNK_heap_pos=ti_y_kmin) + if (lag<=ti_y_kmax && lag>=ti_y_kmin) /*Build the index for sparse matrix containing the jacobian : u*/ { //mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag); var+=Size*t; @@ -592,7 +592,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map< mexPrintf("=> "); #endif } - else + else /*Build the additive terms ooutside the simulation periods related to the first lags and the las leads...*/ { #ifdef PRINT_OUT mexPrintf("nn "); @@ -604,7 +604,7 @@ void SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, std::map< //mexPrintf(" u[%d](%f)*y[%d](%f)=%f",it4->second+u_count_init*t,u[it4->second+u_count_init*t],index_vara[var+Size*(y_kmin+t)],y[index_vara[var+Size*(y_kmin+t)]],u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]]); } } - else + else /* ...and store it in the u vector*/ { #ifdef PRINT_OUT mexPrintf(""); @@ -806,11 +806,11 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s; int *diff1, *diff2; #ifdef MEM_ALLOC_CHK - mexPrintf("diff1=(int*)mxMalloc(%d*sizeof(int))\n",nop); + mexPrintf("diff1=(int*)mxMalloc(%f)\n",double(nop)*double(sizeof(int))/1024); #endif diff1=(int*)mxMalloc(nop*sizeof(int)); #ifdef MEM_ALLOC_CHK - mexPrintf("diff1=(int*)mxMalloc(%d*sizeof(int))\n",nop); + mexPrintf("diff2=(int*)mxMalloc(%f)\n",double(nop)*double(sizeof(int))/1024); #endif diff2=(int*)mxMalloc(nop*sizeof(int)); #ifdef MEM_ALLOC_CHK @@ -1010,7 +1010,9 @@ SparseMatrix::compare( int *save_op, int *save_opa, int *save_opaa, int beg_t, i mexErrMsgTxt(filename.c_str()); #endif } + //mexPrintf("mxFree(diff1)\n"); mxFree(diff1); + //mexPrintf("mxFree(diff2)\n"); mxFree(diff2); return OK; } @@ -1373,6 +1375,7 @@ SparseMatrix::bksub( int tbreak, int last_period, int Size, double slowc_l #endif yy=-(yy+y[eq]+u[b[pos]]); direction[eq]=yy; + //mexPrintf("direction[%d] = %f\n",eq,yy); y[eq] += slowc_l*yy; #ifdef PRINT_OUT_y1 mexPrintf("=%f (%f)\n",double(yy),double(y[eq])); @@ -1409,6 +1412,7 @@ SparseMatrix::Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_k double uu, yy; //char tmp_s[150]; //mexPrintf("period=%d\n",period); + //mexPrintf("Direct_Simulate\n"); #ifdef PRINT_OUT for (j = 0;j < it_ -y_kmin;j++) { @@ -1990,7 +1994,10 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax tdelete1=0; tdelete2=0; tdelete21=0; tdelete22=0; tdelete221=0; tdelete222=0; #endif if (iter>0) - mexPrintf("Sim : %f ms\n",(1000.0*(double(clock())-double(time00)))/double(CLOCKS_PER_SEC)); + { + mexPrintf("Sim : %f ms\n",(1000.0*(double(clock())-double(time00)))/double(CLOCKS_PER_SEC)); + mexEvalString("drawnow;"); + } #ifdef MEMORY_LEAKS mexEvalString("feature('memstats');"); #endif @@ -2037,12 +2044,15 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax } } res1a=res1; - mexPrintf("-----------------------------------\n"); - mexPrintf(" Simulate iteration° %d \n",iter+1); - mexPrintf(" max. error=%.10e \n",double(max_res)); - mexPrintf(" sqr. error=%.10e \n",double(res2)); - mexPrintf(" abs. error=%.10e \n",double(res1)); - mexPrintf("-----------------------------------\n"); + if(print_it) + { + mexPrintf("-----------------------------------\n"); + mexPrintf(" Simulate iteration° %d \n",iter+1); + mexPrintf(" max. error=%.10e \n",double(max_res)); + mexPrintf(" sqr. error=%.10e \n",double(res2)); + mexPrintf(" abs. error=%.10e \n",double(res1)); + mexPrintf("-----------------------------------\n"); + } if (cvg) return(0); #ifdef PRINT_OUT diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc index 336334625..d0cffcd41 100644 --- a/preprocessor/BlockTriangular.cc +++ b/preprocessor/BlockTriangular.cc @@ -103,12 +103,13 @@ BlockTriangular::Prologue_Epilogue(bool* IM, int* prologue, int* epilogue, int n void BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, BlockType type, Model_Block * ModelBlock) { - int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1; - int *tmp_size, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_exo, tmp_nb_exo, nb_lead_lag_endo; + int i, j, k, l, ls, m, i_1, Lead, Lag, first_count_equ, i1, li; + int *tmp_size, *tmp_size_other_endo, *tmp_size_exo, *tmp_var, *tmp_endo, *tmp_other_endo, *tmp_exo, tmp_nb_other_endo, tmp_nb_exo, nb_lead_lag_endo; + bool *tmp_variable_evaluated; bool *Cur_IM; bool *IM, OK; ModelBlock->Periods = periods; - int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo; + int Lag_Endo, Lead_Endo, Lag_Exo, Lead_Exo, Lag_Other_Endo, Lead_Other_Endo; if ((type == PROLOGUE) || (type == EPILOGUE)) { for(i = 0;i < size;i++) @@ -120,15 +121,24 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Temporary_terms=new temporary_terms_type (); ModelBlock->Block_List[*count_Block].Temporary_terms->clear(); 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_size_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_var = (int*)malloc(sizeof(int)); tmp_size_exo = (int*)malloc((incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1) * sizeof(int)); + memset(tmp_size_exo, 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)); - nb_lead_lag_endo = Lead = Lag = 0; - Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0; + memset(tmp_size_other_endo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + memset(tmp_other_endo, 0, symbol_table.endo_nbr*sizeof(int)); + + nb_lead_lag_endo = Lead = Lag = 0; + Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = Lag_Other_Endo = Lead_Other_Endo = 0; + + tmp_variable_evaluated = (bool*)malloc(symbol_table.endo_nbr*sizeof(bool)); + memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool)); for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++) { Cur_IM = incidencematrix.Get_IM(k, eEndogenous); @@ -139,6 +149,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { if(Cur_IM[i_1 + Index_Var_IM[*count_Equ].index]) { + tmp_variable_evaluated[Index_Var_IM[*count_Equ].index] = true; nb_lead_lag_endo++; tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if(k > Lead) @@ -149,6 +160,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { if(Cur_IM[i_1 + Index_Var_IM[*count_Equ].index]) { + tmp_variable_evaluated[Index_Var_IM[*count_Equ].index] = true; tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; nb_lead_lag_endo++; if(-k > Lag) @@ -160,6 +172,45 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc Lag_Endo = Lag; Lead_Endo = Lead; + + tmp_nb_other_endo = 0; + for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++) + { + Cur_IM = incidencematrix.Get_IM(k, eEndogenous); + if(Cur_IM) + { + i_1 = Index_Equ_IM[*count_Equ].index * symbol_table.endo_nbr; + for(j = 0;j < symbol_table.endo_nbr;j++) + { + int ij = Index_Var_IM[j].index; + if(Cur_IM[i_1 + ij]) + { + if(!tmp_variable_evaluated[ij]) + { + if(!tmp_other_endo[ij]) + { + tmp_other_endo[ij] = 1; + tmp_nb_other_endo++; + } + if(k>0 && k>Lead_Other_Endo) + Lead_Other_Endo = k; + else if(k<0 && (-k)>Lag_Other_Endo) + Lag_Other_Endo = -k; + if(k>0 && k>Lead) + Lead = k; + else if(k<0 && (-k)>Lag) + Lag = -k; + tmp_size_other_endo[k+incidencematrix.Model_Max_Lag_Endo]++; + } + } + } + } + } + ModelBlock->Block_List[*count_Block].nb_other_endo = tmp_nb_other_endo; + ModelBlock->Block_List[*count_Block].Other_Endogenous = (int*)malloc(tmp_nb_other_endo * sizeof(int)); + + + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); tmp_nb_exo = 0; @@ -214,6 +265,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Own_Derivative = (int*)malloc(sizeof(int)); ModelBlock->Block_List[*count_Block].Equation[0] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].Variable[0] = Index_Var_IM[*count_Equ].index; + + + if ((Lead > 0) && (Lag > 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_SIMPLE; else if((Lead > 0) && (Lag == 0)) @@ -267,6 +321,16 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].size_other_endo = tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].nb_other_endo = tmp_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + li] * sizeof(int)); + + + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_init = l; IM = incidencematrix.Get_IM(li - Lag, eEndogenous); if(IM) @@ -290,10 +354,24 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var[m] = 0; ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index[m] = Index_Equ_IM[*count_Equ].index; ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index[m] = Index_Var_IM[*count_Equ].index; + tmp_variable_evaluated[Index_Var_IM[*count_Equ].index] = true; l++; m++; } ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_finish = l - 1; + + m = 0; + for(k = 0;k < symbol_table.endo_nbr;k++) + if((!tmp_variable_evaluated[Index_Var_IM[k].index]) && IM[Index_Var_IM[k].index + i_1]) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].u_other_endo[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_other_endo[m] = 0; //j - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_other_endo[m] = k ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Equ_Index_other_endo[m] = Index_Equ_IM[*count_Equ].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[li].Var_Index_other_endo[m] = Index_Var_IM[k].index; + l++; + m++; + } } } else @@ -333,6 +411,9 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc free(tmp_endo); free(tmp_exo); free(tmp_var); + free(tmp_size_other_endo); + free(tmp_other_endo); + free(tmp_variable_evaluated); } } else @@ -350,14 +431,21 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc 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_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)); memset(tmp_size_exo, 0, (incidencematrix.Model_Max_Lead + incidencematrix.Model_Max_Lag + 1)*sizeof(int)); + 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)); nb_lead_lag_endo = 0; - Lag_Endo = Lead_Endo = Lag_Exo = Lead_Exo = 0; + Lag_Endo = Lead_Endo = Lag_Other_Endo = Lead_Other_Endo = Lag_Exo = Lead_Exo = 0; + //Variable by variable looking for all leads and lags its occurence in each equation of the block + tmp_variable_evaluated = (bool*)malloc(symbol_table.endo_nbr*sizeof(bool)); + memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool)); for(i = 0;i < size;i++) { ModelBlock->Block_List[*count_Block].Equation[i] = Index_Equ_IM[*count_Equ].index; @@ -375,6 +463,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc { if(Cur_IM[i_1 + Index_Equ_IM[first_count_equ + j].index*symbol_table.endo_nbr]) { + tmp_variable_evaluated[i_1] = true; tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if (!OK) { @@ -396,6 +485,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc tmp_size[incidencematrix.Model_Max_Lag_Endo + k]++; if (!OK) { + tmp_variable_evaluated[i_1] = true; tmp_endo[incidencematrix.Model_Max_Lag + k]++; nb_lead_lag_endo++; OK = true; @@ -409,8 +499,6 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc } (*count_Equ)++; } - - if ((Lag > 0) && (Lead > 0)) ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_TWO_BOUNDARIES_COMPLETE; else if(size > 1) @@ -427,12 +515,45 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc else ModelBlock->Block_List[*count_Block].Simulation_Type = SOLVE_FORWARD_SIMPLE; } - - Lag_Endo = Lag; Lead_Endo = Lead; + + tmp_nb_other_endo = 0; + for(i = 0;i < size;i++) + { + for(k = -incidencematrix.Model_Max_Lag_Endo; k<=incidencematrix.Model_Max_Lead_Endo; k++) + { + Cur_IM = incidencematrix.Get_IM(k, eEndogenous); + if(Cur_IM) + { + i_1 = Index_Equ_IM[first_count_equ+i].index * symbol_table.endo_nbr; + for(j = 0;j < symbol_table.endo_nbr;j++) + if(Cur_IM[i_1 + j]) + { + if(!tmp_variable_evaluated[j]) + { + tmp_other_endo[j] = 1; + tmp_nb_other_endo++; + } + if(k>0 && k>Lead_Other_Endo) + Lead_Other_Endo = k; + else if(k<0 && (-k)>Lag_Other_Endo) + Lag_Other_Endo = -k; + if(k>0 && k>Lead) + Lead = k; + else if(k<0 && (-k)>Lag) + Lag = -k; + tmp_size_other_endo[k+incidencematrix.Model_Max_Lag_Endo]++; + } + } + } + } + ModelBlock->Block_List[*count_Block].nb_other_endo = tmp_nb_other_endo; + ModelBlock->Block_List[*count_Block].Other_Endogenous = (int*)malloc(tmp_nb_other_endo * sizeof(int)); + + tmp_exo = (int*)malloc(symbol_table.exo_nbr * sizeof(int)); - memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); + memset(tmp_exo, 0, symbol_table.exo_nbr * sizeof(int)); tmp_nb_exo = 0; for(i = 0;i < size;i++) { @@ -481,10 +602,12 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].Max_Lead = Lead; ModelBlock->Block_List[*count_Block].Max_Lag_Endo = Lag_Endo; ModelBlock->Block_List[*count_Block].Max_Lead_Endo = Lead_Endo; + ModelBlock->Block_List[*count_Block].Max_Lag_Other_Endo = Lag_Other_Endo; + ModelBlock->Block_List[*count_Block].Max_Lead_Other_Endo = Lead_Other_Endo; ModelBlock->Block_List[*count_Block].Max_Lag_Exo = Lag_Exo; ModelBlock->Block_List[*count_Block].Max_Lead_Exo = Lead_Exo; ModelBlock->Block_List[*count_Block].IM_lead_lag = (IM_compact*)malloc((Lead + Lag + 1) * sizeof(IM_compact)); - ls = l = size; + ls = l = li = size; i1 = 0; ModelBlock->Block_List[*count_Block].Nb_Lead_Lag_Endo = nb_lead_lag_endo; for(i = 0;i < Lead + Lag + 1;i++) @@ -499,6 +622,14 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index = (int*)malloc(tmp_size[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_other_endo = tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].nb_other_endo = tmp_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i]; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index_other_endo = (int*)malloc(tmp_size_other_endo[incidencematrix.Model_Max_Lag_Endo - Lag + i] * sizeof(int)); } else ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size = 0; @@ -513,6 +644,7 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc else ModelBlock->Block_List[*count_Block].IM_lead_lag[i].size_exo = 0; ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_init = l; + memset(tmp_variable_evaluated, 0, symbol_table.endo_nbr*sizeof(bool)); IM = incidencematrix.Get_IM(i - Lag, eEndogenous); if(IM) { @@ -541,16 +673,34 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc ModelBlock->Block_List[*count_Block].IM_lead_lag[i].us[m] = ls; ls++; } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u[m] = li; ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ[m] = j - first_count_equ; ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var[m] = k - first_count_equ; ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index[m] = Index_Equ_IM[j].index; ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index[m] = Index_Var_IM[k].index; + tmp_variable_evaluated[Index_Var_IM[k].index] = true; + l++; + m++; + li++; + } + } + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = li - 1; + m = 0; + for(j = first_count_equ;j < size + first_count_equ;j++) + { + i_1 = Index_Equ_IM[j].index * symbol_table.endo_nbr; + for(k = 0;k < symbol_table.endo_nbr;k++) + if((!tmp_variable_evaluated[Index_Var_IM[k].index]) && IM[Index_Var_IM[k].index + i_1]) + { + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_other_endo[m] = l; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_other_endo[m] = j - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_other_endo[m] = k - first_count_equ; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Equ_Index_other_endo[m] = Index_Equ_IM[j].index; + ModelBlock->Block_List[*count_Block].IM_lead_lag[i].Var_Index_other_endo[m] = Index_Var_IM[k].index; l++; m++; } } - ModelBlock->Block_List[*count_Block].IM_lead_lag[i].u_finish = l - 1; } IM = incidencematrix.Get_IM(i - Lag, eExogenous); if(IM) @@ -575,10 +725,13 @@ BlockTriangular::Allocate_Block(int size, int *count_Equ, int *count_Block, Bloc } (*count_Block)++; free(tmp_size); + free(tmp_size_other_endo); free(tmp_size_exo); free(tmp_endo); + free(tmp_other_endo); free(tmp_exo); free(tmp_var); + free(tmp_variable_evaluated); } } @@ -712,7 +865,7 @@ BlockTriangular::Normalize_and_BlockDecompose(bool* IM, Model_Block* ModelBlock, SIM0[iter->first.first*n+iter->first.second]=1; if(!IM_0[iter->first.first*n+iter->first.second]) { - cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << "\n"; + cout << "Error nothing at IM_0[" << iter->first.first << ", " << iter->first.second << "]=" << IM_0[iter->first.first*n+iter->first.second] << " " << iter->second << "\n"; } } else diff --git a/preprocessor/DynareMain2.cc b/preprocessor/DynareMain2.cc index e6111e910..5ae76a02f 100644 --- a/preprocessor/DynareMain2.cc +++ b/preprocessor/DynareMain2.cc @@ -35,6 +35,9 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool no_tm // Run checking pass mod_file->checkPass(); + // Evaluate parameters initialization, initval, endval and pounds + mod_file->evalAllExpressions(); + // Do computations mod_file->computingPass(no_tmp_terms); diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc index 8b4bf6163..352eb7b73 100644 --- a/preprocessor/ExprNode.cc +++ b/preprocessor/ExprNode.cc @@ -278,7 +278,10 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type, case eModelLocalVariable: case eModFileLocalVariable: - output << datatree.symbol_table.getNameByID(type, symb_id); + if(type==oMatlabDynamicModelSparse || type==oMatlabStaticModelSparse) + datatree.local_variables_table[symb_id]->writeOutput(output, output_type,temporary_terms); + else + output << datatree.symbol_table.getNameByID(type, symb_id); break; case eEndogenous: @@ -406,9 +409,14 @@ VariableNode::eval(const eval_context_type &eval_context) const throw (EvalExcep // ModelTree::evaluateJacobian need to have the initval values applied to lead/lagged variables also /*if (lag != 0) throw EvalException();*/ + /*if(type==eModelLocalVariable) + cout << "eModelLocalVariable = " << symb_id << "\n";*/ eval_context_type::const_iterator it = eval_context.find(make_pair(symb_id, type)); if (it == eval_context.end()) - throw EvalException(); + { + //cout << "unknonw variable type = " << type << " simb_id = " << symb_id << "\n"; + throw EvalException(); + } return it->second; } @@ -464,10 +472,10 @@ VariableNode::compile(ofstream &CompileCode, bool lhs_rhs, ExprNodeOutputType ou break; case eModelLocalVariable: case eModFileLocalVariable: - cerr << "VariableNode::compile: unhandled variable type" << endl; - exit(EXIT_FAILURE); + datatree.local_variables_table[symb_id]->compile(CompileCode, lhs_rhs, output_type, temporary_terms, map_idx); + break; case eUnknownFunction: - cerr << "Impossible case" << endl; + cerr << "Impossible case: eUnknownFuncion" << endl; exit(EXIT_FAILURE); } } @@ -477,6 +485,8 @@ VariableNode::collectEndogenous(set > &result) const { if (type == eEndogenous) result.insert(make_pair(symb_id, lag)); + else if (type == eModelLocalVariable) + datatree.local_variables_table[symb_id]->collectEndogenous(result); } void @@ -484,6 +494,8 @@ VariableNode::collectExogenous(set > &result) const { if (type == eExogenous) result.insert(make_pair(symb_id, lag)); + else if (type == eModelLocalVariable) + datatree.local_variables_table[symb_id]->collectExogenous(result); } diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index dbd8eb35a..9bef8e80e 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -20,7 +20,7 @@ #include #include #include - +#include #include "ModFile.hh" ModFile::ModFile() : expressions_tree(symbol_table, num_constants), @@ -36,6 +36,116 @@ ModFile::~ModFile() delete (*it); } + +void +ModFile::evalAllExpressions() +{ + //Evaluate Parameters + cout << "Evaluating expressions ..."; + InitParamStatement *it; + int j=0; + for(vector::const_iterator it1=statements.begin();it1!=statements.end(); it1++) + { + it=dynamic_cast(*it1); + if(it) + { + try + { + const NodeID expression = it->get_expression(); + double val = expression->eval(global_eval_context); + int symb_id = symbol_table.getID(it->get_name()); + global_eval_context[make_pair(symb_id, eParameter)] = val; + j++; + } + catch(ExprNode::EvalException &e) + { + cout << "error in evaluation of param\n"; + } + } + } + if (j!=symbol_table.parameter_nbr) + { + cout << "Warning: Uninitialized parameters: \n"; + for(j=0;j first; + const NodeID expression = it->second; + SymbolType type = symbol_table.getType(name); + double val = expression->eval(global_eval_context); + int symb_id = symbol_table.getID(name); + global_eval_context[make_pair(symb_id, type)] = val; + } + catch(ExprNode::EvalException &e) + { + cout << "error in evaluation of variable\n"; + } + } + if(init_values.size()!=symbol_table.endo_nbr+symbol_table.exo_nbr+symbol_table.exo_det_nbr) + { + cout << "\nWarning: Uninitialized variable: \n"; + cout << "Endogenous\n"; + for(j=0;j ::const_iterator it = model_tree.local_variables_table.begin(); it !=model_tree.local_variables_table.end(); it++) + { + try + { + const NodeID expression = it->second; + double val = expression->eval(global_eval_context); + //cout << it->first << " " << symbol_table.getNameByID(eModelLocalVariable, it->first) << " = " << val << "\n"; + global_eval_context[make_pair(it->first, eModelLocalVariable)] = val; + } + catch(ExprNode::EvalException &e) + { + cout << "error in evaluation of pound\n"; + } + } + if(model_tree.local_variables_table.size()!=symbol_table.model_local_variable_nbr+symbol_table.modfile_local_variable_nbr) + { + cout << "Warning: Unitilialized pound: \n"; + cout << "Local variable in a model\n"; + for(j=0;j ::iterator it = statements.begin(); it != statements.end(); it++) (*it)->computingPass(); + //evalAllExpressions(); } void diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index a947dea17..8f10ca81e 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -390,7 +390,7 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & ostringstream Uf[symbol_table.endo_nbr]; map reference_count; int prev_Simulation_Type=-1, count_derivates=0; - int jacobian_max_endo_col; + int jacobian_max_endo_col, jacobian_max_exo_col; ofstream output; temporary_terms_type::const_iterator it_temp=temporary_terms.begin(); //---------------------------------------------------------------------- @@ -404,7 +404,7 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & else tmp_output << " "; - (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); + (*it)->writeOutput(tmp_output, oMatlabDynamicModelSparse, temporary_terms); } if(tmp_output.str().length()) @@ -553,8 +553,8 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & break; case SOLVE_BACKWARD_SIMPLE: case SOLVE_FORWARD_SIMPLE: - output << sps << "residual(" << i+1 << ") = ("; - goto end; + /*output << sps << "residual(" << i+1 << ") = ("; + goto end;*/ case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = residual(" << i+1 << ")"; @@ -584,8 +584,6 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & output << " " << sps << "% Jacobian " << endl << " if jacobian_eval" << endl; switch(ModelBlock->Block_List[j].Simulation_Type) { - case SOLVE_BACKWARD_SIMPLE: - case SOLVE_FORWARD_SIMPLE: case EVALUATE_BACKWARD: case EVALUATE_FORWARD: case EVALUATE_BACKWARD_R: @@ -619,13 +617,32 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_X_Index[i]; int var=ModelBlock->Block_List[j].IM_lead_lag[m].Exogenous_Index[i]; output << " g1(" << eq+1 << ", " - << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr/*ModelBlock->Block_List[j].nb_exo*/ << ") = "; + << jacobian_max_endo_col+var+1+(m+variable_table.max_exo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.exo_nbr << ") = "; writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eExogenous); output << "; % variable=" << symbol_table.getNameByID(eExogenous, var) << "(" << k << ") " << var+1 << ", equation=" << eq+1 << endl; } } + jacobian_max_exo_col=(variable_table.max_exo_lag+variable_table.max_exo_lead+1)*symbol_table.exo_nbr; + 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; + if(block_triangular.incidencematrix.Model_Max_Lag_Endo - ModelBlock->Block_List[j].Max_Lag +m >=0) + { + for(i=0;iBlock_List[j].IM_lead_lag[m].size_other_endo;i++) + { + int eq=ModelBlock->Block_List[j].IM_lead_lag[m].Equ_Index_other_endo[i]; + int var=ModelBlock->Block_List[j].IM_lead_lag[m].Var_Index_other_endo[i]; + output << " g1(" << eq+1 << ", " + << jacobian_max_endo_col+jacobian_max_exo_col+var+1+(m+variable_table.max_endo_lag-ModelBlock->Block_List[j].Max_Lag)*symbol_table.endo_nbr << ") = "; + writeDerivative(output, eq, var, k, oMatlabDynamicModelSparse, temporary_terms, eEndogenous); + output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) + << "(" << k << ") " << var+1 + << ", equation=" << eq+1 << endl; + } + } + } if (ModelBlock->Block_List[j].Simulation_Type==SOLVE_BACKWARD_SIMPLE || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) { @@ -645,6 +662,8 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & || ModelBlock->Block_List[j].Simulation_Type==SOLVE_FORWARD_SIMPLE) output << " end;" << endl; break; + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: count_derivates++; @@ -719,7 +738,6 @@ ModelTree::writeModelEquationsOrdered_M( Model_Block *ModelBlock, const string & Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_+" << k-1 << "))*y(it_+" << k << ", " << var+1 << ")"; else if (k<0) Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+y_size*(it_" << k-1 << "))*y(it_" << k << ", " << var+1 << ")"; - //output << " u(" << u+1 << "+Per_u_) = "; if(k==0) output << " g1(" << eqr+1 << "+Per_J_, " << varr+1 << "+Per_K_) = "; else if(k==1) @@ -997,10 +1015,12 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str lhs->writeOutput(output, oMatlabStaticModelSparse, temporary_terms); output << ";\n"; break; + case SOLVE_BACKWARD_SIMPLE: + case SOLVE_FORWARD_SIMPLE: case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: case SOLVE_TWO_BOUNDARIES_COMPLETE: - Uf[ModelBlock->Block_List[j].Equation[i]] << " b(" << i+1 << ") = - residual(" << i+1 << ")"; + Uf[ModelBlock->Block_List[j].Equation[i]] << "b(" << i+1 << ") = residual(" << i+1 << ")"; goto end; default: end: @@ -1034,33 +1054,8 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str break; case SOLVE_BACKWARD_SIMPLE: case SOLVE_FORWARD_SIMPLE: - output << " g1(1)="; - writeDerivative(output, ModelBlock->Block_List[j].Equation[0], ModelBlock->Block_List[j].Variable[0], 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, ModelBlock->Block_List[j].Variable[0]) - << "(" << variable_table.getLag(variable_table.getSymbolID(ModelBlock->Block_List[j].Variable[0])) - << ") " << ModelBlock->Block_List[j].Variable[0]+1 - << ", equation=" << ModelBlock->Block_List[j].Equation[0]+1 << endl; - break; case SOLVE_BACKWARD_COMPLETE: case SOLVE_FORWARD_COMPLETE: - output << " g2=0;g3=0;\n"; - m=ModelBlock->Block_List[j].Max_Lag; - for(i=0;iBlock_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]; - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y(" << var+1 << ")"; - output << " g1(" << eqr+1 << ", " << varr+1 << ") = "; - writeDerivative(output, eq, var, 0, oMatlabStaticModelSparse, temporary_terms, eEndogenous); - output << "; % variable=" << symbol_table.getNameByID(eEndogenous, var) - << "(" << variable_table.getLag(variable_table.getSymbolID(var)) << ") " << var+1 - << ", equation=" << eq+1 << endl; - } - for(i = 0;i < ModelBlock->Block_List[j].Size;i++) - output << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; - break; case SOLVE_TWO_BOUNDARIES_COMPLETE: case SOLVE_TWO_BOUNDARIES_SIMPLE: output << " g2=0;g3=0;\n"; @@ -1075,9 +1070,9 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str int varr=ModelBlock->Block_List[j].IM_lead_lag[m].Var[i]; if(!IM[eqr*ModelBlock->Block_List[j].Size+varr]) { - Uf[ModelBlock->Block_List[j].Equation[eqr]] << "+g1(" << eqr+1 + Uf[ModelBlock->Block_List[j].Equation[eqr]] << "-g1(" << eqr+1 << ", " << varr+1 << ")*y( " << var+1 << ")"; - IM[eqr*ModelBlock->Block_List[j].Size+varr]=1; + IM[eqr*ModelBlock->Block_List[j].Size+varr]=true; } output << " g1(" << eqr+1 << ", " << varr+1 << ") = g1(" << eqr+1 << ", " << varr+1 << ") + "; writeDerivative(output, eq, var, k, oMatlabStaticModelSparse, temporary_terms, eEndogenous); @@ -1090,7 +1085,6 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str #endif } } - output << " if(~jacobian_eval)\n"; for(i = 0;i < ModelBlock->Block_List[j].Size;i++) { output << " " << Uf[ModelBlock->Block_List[j].Equation[i]].str() << ";\n"; @@ -1099,7 +1093,7 @@ ModelTree::writeModelStaticEquationsOrdered_M(Model_Block *ModelBlock, const str output << " condition(" << i+1 << ")=u(" << i+1 << "+Per_u_);\n"; #endif } - output << " end\n"; + //output << " end\n"; #ifdef CONDITION for(m=0;m<=ModelBlock->Block_List[j].Max_Lead+ModelBlock->Block_List[j].Max_Lag;m++) { @@ -1909,6 +1903,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); u_count_int++; } + for(j=0;jBlock_List[num].Size;j++) { int varr=block_triangular.ModelBlock->Block_List[num].Variable[j]; @@ -1947,11 +1942,43 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b mStaticModelFile << "%/\n"; mStaticModelFile << "function [varargout] = " << static_basename << "(varargin)\n"; mStaticModelFile << " global oo_ M_ options_ ys0_ ;\n"; + bool OK=true; + ostringstream tmp_output; + for(temporary_terms_type::const_iterator it = temporary_terms.begin(); + it != temporary_terms.end(); it++) + { + if (OK) + OK=false; + else + tmp_output << " "; + (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); + } + if (tmp_output.str().length()>0) + mStaticModelFile << " global " << tmp_output.str() << " M_ ;\n"; + mStaticModelFile << " T_init=0;\n"; + tmp_output.str(""); + for(temporary_terms_type::const_iterator it = temporary_terms.begin(); + it != temporary_terms.end(); it++) + { + tmp_output << " "; + (*it)->writeOutput(tmp_output, oMatlabDynamicModel, temporary_terms); + tmp_output << "=T_init;\n"; + } + if (tmp_output.str().length()>0) + mStaticModelFile << tmp_output.str(); + mStaticModelFile << " y_kmin=M_.maximum_lag;\n"; mStaticModelFile << " y_kmax=M_.maximum_lead;\n"; mStaticModelFile << " y_size=M_.endo_nbr;\n"; + + /*tmp_output.str(""); + writeModelLocalVariables(tmp_output, oMatlabDynamicModel); + if (tmp_output.str().length()>0) + mStaticModelFile << tmp_output.str() << "\n";*/ + + mStaticModelFile << " if(length(varargin)>0)\n"; - mStaticModelFile << " %it is a simple evaluation of the dynamic model for time _it\n"; + mStaticModelFile << " %A simple evaluation of the static model\n"; //mStaticModelFile << " global it_;\n"; mStaticModelFile << " y=varargin{1}(:);\n"; mStaticModelFile << " ys=y;\n"; @@ -1994,21 +2021,6 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b } } - - /*mStaticModelFile << " y_index=["; - for(int ik=0;ikBlock_List[i].Size;ik++) - { - mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; - } - mStaticModelFile << " ];\n"; - mStaticModelFile << " y_index_eq=["; - for(int ik=0;ikBlock_List[i].Size;ik++) - { - mStaticModelFile << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; - } - mStaticModelFile << " ];\n"; - k=block_triangular.ModelBlock->Block_List[i].Simulation_Type;*/ - if (BlockTriangular::BlockSim(prev_Simulation_Type)==BlockTriangular::BlockSim(k) && (k==EVALUATE_FORWARD || k==EVALUATE_BACKWARD || k==EVALUATE_FORWARD_R || k==EVALUATE_BACKWARD_R)) skip_head=true; @@ -2030,8 +2042,6 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b } else ga_index++; - /*mStaticModelFile << " residual(y_index)=ys(y_index)-y(y_index);\n"; - mStaticModelFile << " g1(y_index_eq, y_index) = ga(" << ga_index << ", " << ga_index << ");\n";*/ break; case SOLVE_FORWARD_COMPLETE: case SOLVE_BACKWARD_COMPLETE: @@ -2055,7 +2065,7 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b mStaticModelFile << " varargout{2}=g1;\n"; mStaticModelFile << " return;\n"; mStaticModelFile << " end;\n"; - mStaticModelFile << " %it is the deterministic simulation of the block decomposed static model\n"; + mStaticModelFile << " %The deterministic simulation of the block decomposed static model\n"; mStaticModelFile << " periods=options_.periods;\n"; mStaticModelFile << " maxit_=options_.maxit_;\n"; mStaticModelFile << " solve_tolf=options_.solve_tolf;\n"; @@ -2088,32 +2098,6 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b } open_par=false; } - /*else if ((k == SOLVE_FORWARD_SIMPLE || k == SOLVE_BACKWARD_SIMPLE) && (block_triangular.ModelBlock->Block_List[i].Size)) - { - if (open_par) - { - mStaticModelFile << " end\n"; - } - open_par=false; - mStaticModelFile << " g1=0;\n"; - mStaticModelFile << " r=0;\n"; - mStaticModelFile << " cvg=0;\n"; - mStaticModelFile << " iter=0;\n"; - mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mStaticModelFile << " [r, g1] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n"; - mStaticModelFile << " y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ") = y(" << block_triangular.ModelBlock->Block_List[i].Variable[0]+1 << ")-r/g1;\n"; - mStaticModelFile << " cvg=((r*r)Block_List[i].Size)) { if (open_par) @@ -2134,46 +2118,11 @@ ModelTree::writeSparseStaticMFile(const string &static_basename, const string &b int nze, m; 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; - /*mStaticModelFile << " if(isfield(oo_.deterministic_simulation,'block'))\n"; - mStaticModelFile << " blck_num = length(oo_.deterministic_simulation.block)+1;\n"; - mStaticModelFile << " else\n"; - mStaticModelFile << " blck_num = 1;\n"; - mStaticModelFile << " end;\n";*/ mStaticModelFile << " y = solve_one_boundary('" << static_basename << "_" << i + 1 << "'" << ", y, x, params, y_index, " << nze << - ", y_kmin, " << block_triangular.ModelBlock->Block_List[i].is_linear << + ", 1, " << block_triangular.ModelBlock->Block_List[i].is_linear << ", " << Blck_Num << ", y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 0, 0);\n"; - /*mStaticModelFile << " r=0;\n"; - mStaticModelFile << " cvg=0;\n"; - mStaticModelFile << " iter=0;\n"; - mStaticModelFile << " lambda=1;\n"; - mStaticModelFile << " stpmx = 100 ;\n"; - mStaticModelFile << " stpmax = stpmx*max([sqrt(y'*y);size(y_index,2)]);\n"; - mStaticModelFile << " nn=1:size(y_index,2);\n"; - mStaticModelFile << " while ~(cvg==1 | iter>maxit_),\n"; - mStaticModelFile << " [r, g1, g2, g3, b] = " << static_basename << "_" << i + 1 << "(y, x, 0);\n"; - mStaticModelFile << " max_res=max(abs(r));\n"; - mStaticModelFile << " cvg=(max_resBlock_List[i].Size)) { @@ -2461,12 +2408,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string open_par=false; mDynamicModelFile << " g1=0;\n"; mDynamicModelFile << " r=0;\n"; - tmp_eq.str(""); + tmp.str(""); for(int ik=0;ikBlock_List[i].Size;ik++) { - tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; } - mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; + mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; int nze, m; 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; @@ -2476,7 +2423,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " blck_num = 1;\n"; mDynamicModelFile << " end;\n"; mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" << - ", y, x, params, y_index_eq, " << nze << + ", y, x, params, y_index, " << nze << ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear << ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n"; @@ -2488,12 +2435,12 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string open_par=false; mDynamicModelFile << " g1=0;\n"; mDynamicModelFile << " r=0;\n"; - tmp_eq.str(""); + tmp.str(""); for(int ik=0;ikBlock_List[i].Size;ik++) { - tmp_eq << " " << block_triangular.ModelBlock->Block_List[i].Equation[ik]+1; + tmp << " " << block_triangular.ModelBlock->Block_List[i].Variable[ik]+1; } - mDynamicModelFile << " y_index_eq = [" << tmp_eq.str() << "];\n"; + mDynamicModelFile << " y_index = [" << tmp.str() << "];\n"; int nze, m; 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; @@ -2503,7 +2450,7 @@ ModelTree::writeSparseDynamicMFile(const string &dynamic_basename, const string mDynamicModelFile << " blck_num = 1;\n"; mDynamicModelFile << " end;\n"; mDynamicModelFile << " y = solve_one_boundary('" << dynamic_basename << "_" << i + 1 << "'" << - ", y, x, params, y_index_eq, " << nze << + ", y, x, params, y_index, " << nze << ", options_.periods, " << block_triangular.ModelBlock->Block_List[i].is_linear << ", blck_num, y_kmin, options_.maxit_, options_.solve_tolf, options_.slowc, options_.cutoff, options_.simulation_method, 1, 1, 0);\n"; } @@ -3072,7 +3019,9 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_ } catch(ExprNode::EvalException &e) { - cerr << "ModelTree::evaluateJacobian: evaluation of Jacobian failed!" << endl; + cout << "evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getNameByID(eEndogenous, variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ") [" << variable_table.getSymbolID(it->first.second) << "] !" << endl; + Id->writeOutput(cout, oMatlabDynamicModelSparse, temporary_terms);cout << "\n"; + cerr << "ModelTree::evaluateJacobian: evaluation of Jacobian failed for equation " << it->first.first+1 << " and variable " << symbol_table.getNameByID(eEndogenous, variable_table.getSymbolID(it->first.second)) << "(" << variable_table.getLag(it->first.second) << ")!" << endl; } int eq=it->first.first; int var=variable_table.getSymbolID(it->first.second); diff --git a/preprocessor/NumericalInitialization.cc b/preprocessor/NumericalInitialization.cc index cf5bf6005..456d7defd 100644 --- a/preprocessor/NumericalInitialization.cc +++ b/preprocessor/NumericalInitialization.cc @@ -38,6 +38,19 @@ InitParamStatement::writeOutput(ostream &output, const string &basename) const output << param_name << " = M_.params( " << id << " );\n"; } +NodeID +InitParamStatement::get_expression() const +{ + return(param_value); +} + +string +InitParamStatement::get_name() const +{ + return(param_name); +} + + InitOrEndValStatement::InitOrEndValStatement(const init_values_type &init_values_arg, const SymbolTable &symbol_table_arg) : init_values(init_values_arg), diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc index c1109e1e7..d8652cec7 100644 --- a/preprocessor/ParsingDriver.cc +++ b/preprocessor/ParsingDriver.cc @@ -74,8 +74,8 @@ ParsingDriver::parse(istream &in, bool debug) /* Deleting filename in DynareFlex::yyterminate() is prematurate, because if there is an unexpected end of file, the call to ParsingDriver::error() needs filename */ - if (location.begin.filename) - delete location.begin.filename; + /*if (location.begin.filename) + delete location.begin.filename;*/ return mod_file; } @@ -179,13 +179,6 @@ ParsingDriver::add_model_variable(string *name, string *olag) NodeID id = model_tree->AddVariable(*name, lag); - /*if (model_tree->mode == eSparseDLLMode || model_tree->mode == eSparseMode) - { - if (type == eEndogenous) - model_tree->block_triangular.fill_IM(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); - if (type == eExogenous) - model_tree->block_triangular.fill_IM_X(model_tree->equation_number(), mod_file->symbol_table.getID(*name), lag); - }*/ delete name; delete olag; return id; @@ -250,6 +243,7 @@ ParsingDriver::dsample(string *arg1, string *arg2) delete arg2; } + void ParsingDriver::init_param(string *name, NodeID rhs) { @@ -260,7 +254,7 @@ ParsingDriver::init_param(string *name, NodeID rhs) mod_file->addStatement(new InitParamStatement(*name, rhs, mod_file->symbol_table)); // Update global eval context - try + /*try { double val = rhs->eval(mod_file->global_eval_context); int symb_id = mod_file->symbol_table.getID(*name); @@ -269,7 +263,7 @@ ParsingDriver::init_param(string *name, NodeID rhs) catch(ExprNode::EvalException &e) { } - + */ delete name; } @@ -283,11 +277,11 @@ ParsingDriver::init_val(string *name, NodeID rhs) && type != eExogenous && type != eExogenousDet) error("initval/endval: " + *name + " should be an endogenous or exogenous variable"); - - init_values.push_back(make_pair(*name, rhs)); - + //cout << "mod_file->init_values = " << mod_file->init_values << "\n"; + mod_file->init_values.push_back(make_pair(*name, rhs)); + //cout << "init_val " << *name << " mod_file->init_values.size()=" << mod_file->init_values.size() << "\n"; // Update global evaluation context - try + /*try { double val = rhs->eval(mod_file->global_eval_context); int symb_id = mod_file->symbol_table.getID(*name); @@ -296,7 +290,7 @@ ParsingDriver::init_val(string *name, NodeID rhs) catch(ExprNode::EvalException &e) { } - + */ delete name; } @@ -379,15 +373,17 @@ ParsingDriver::sparse() void ParsingDriver::end_initval() { - mod_file->addStatement(new InitValStatement(init_values, mod_file->symbol_table)); - init_values.clear(); + mod_file->addStatement(new InitValStatement(mod_file->init_values, mod_file->symbol_table)); + //mod_file->init_values.clear(); + //cout << "mod_file->init_values.clear() in end_initval()\n"; } void ParsingDriver::end_endval() { - mod_file->addStatement(new EndValStatement(init_values, mod_file->symbol_table)); - init_values.clear(); + mod_file->addStatement(new EndValStatement(mod_file->init_values, mod_file->symbol_table)); + //mod_file->init_values.clear(); + //cout << "mod_file->init_values.clear() in end_endval()\n"; } void diff --git a/preprocessor/include/DataTree.hh b/preprocessor/include/DataTree.hh index 3ee186794..e2029553a 100644 --- a/preprocessor/include/DataTree.hh +++ b/preprocessor/include/DataTree.hh @@ -56,8 +56,6 @@ protected: //! A counter for filling ExprNode's idx field int node_counter; - //! Stores local variables value - map local_variables_table; typedef map num_const_node_map_type; num_const_node_map_type num_const_node_map; @@ -81,6 +79,8 @@ public: //! The variable table VariableTable variable_table; NodeID Zero, One, MinusOne; + //! Stores local variables value + map local_variables_table; //! Raised when a local parameter is declared twice class LocalParameterException diff --git a/preprocessor/include/ExprNode.hh b/preprocessor/include/ExprNode.hh index d3d24e351..7dab87718 100644 --- a/preprocessor/include/ExprNode.hh +++ b/preprocessor/include/ExprNode.hh @@ -321,27 +321,31 @@ public: //! For one lead/lag of one block, stores mapping of information between original model and block-decomposed model struct IM_compact { - int size, u_init, u_finish, nb_endo, size_exo; + int size, u_init, u_finish, nb_endo, nb_other_endo, size_exo, size_other_endo; int *u, *us, *Var, *Equ, *Var_Index, *Equ_Index, *Exogenous, *Exogenous_Index, *Equ_X, *Equ_X_Index; + int *u_other_endo, *Var_other_endo, *Equ_other_endo, *Var_Index_other_endo, *Equ_Index_other_endo; }; //! One block of the model struct Block { - int Size, Sized, nb_exo, nb_exo_det; + int Size, Sized, nb_exo, nb_exo_det, nb_other_endo; BlockType Type; BlockSimulationType Simulation_Type; int Max_Lead, Max_Lag, Nb_Lead_Lag_Endo; int Max_Lag_Endo, Max_Lead_Endo; + int Max_Lag_Other_Endo, Max_Lead_Other_Endo; int Max_Lag_Exo, Max_Lead_Exo; bool is_linear; int *Equation, *Own_Derivative; - int *Variable, *Exogenous; + int *Variable, *Other_Endogenous, *Exogenous; temporary_terms_type *Temporary_terms; IM_compact *IM_lead_lag; int Code_Start, Code_Length; }; + + //! The set of all blocks of the model struct Model_Block { diff --git a/preprocessor/include/ModFile.hh b/preprocessor/include/ModFile.hh index 1a3f3a3b2..3dfea1c96 100644 --- a/preprocessor/include/ModFile.hh +++ b/preprocessor/include/ModFile.hh @@ -26,6 +26,7 @@ using namespace std; #include "SymbolTable.hh" #include "NumericalConstants.hh" +#include "NumericalInitialization.hh" #include "ModelTree.hh" #include "VariableTable.hh" #include "Statement.hh" @@ -49,6 +50,8 @@ public: //! Global evaluation context /*! Filled using initval blocks and parameters initializations */ eval_context_type global_eval_context; + //! Temporary storage for initval/endval blocks + InitOrEndValStatement::init_values_type init_values; private: //! List of statements @@ -59,6 +62,8 @@ private: public: //! Add a statement void addStatement(Statement *st); + //! Evaluate all the statements + void evalAllExpressions(); //! Do some checking and fills mod_file_struct /*! \todo add check for number of equations and endogenous if ramsey_policy is present */ void checkPass(); diff --git a/preprocessor/include/NumericalInitialization.hh b/preprocessor/include/NumericalInitialization.hh index ff9e4ed31..d9fd9360c 100644 --- a/preprocessor/include/NumericalInitialization.hh +++ b/preprocessor/include/NumericalInitialization.hh @@ -39,6 +39,8 @@ public: InitParamStatement(const string ¶m_name_arg, const NodeID param_value_arg, const SymbolTable &symbol_table_arg); virtual void writeOutput(ostream &output, const string &basename) const; + NodeID get_expression() const; + string get_name() const; }; class InitOrEndValStatement : public Statement diff --git a/preprocessor/include/ParsingDriver.hh b/preprocessor/include/ParsingDriver.hh index c9f60996e..67faa0700 100644 --- a/preprocessor/include/ParsingDriver.hh +++ b/preprocessor/include/ParsingDriver.hh @@ -125,8 +125,6 @@ private: SigmaeStatement::row_type sigmae_row; //! Temporary storage for Sigma_e matrix SigmaeStatement::matrix_type sigmae_matrix; - //! Temporary storage for initval/endval blocks - InitOrEndValStatement::init_values_type init_values; //! Temporary storage for histval blocks HistValStatement::hist_values_type hist_values; //! Temporary storage for homotopy_setup blocks diff --git a/tests/ferhat/fs2000.mod b/tests/ferhat/fs2000.mod index 673598ac8..879d991ee 100644 --- a/tests/ferhat/fs2000.mod +++ b/tests/ferhat/fs2000.mod @@ -31,8 +31,8 @@ rho = 0.7; psi = 0.787; del = 0.02; -//model(sparse_dll,no_compiler,cutoff=1e-17); -model(sparse); +model(sparse_dll,cutoff=1e-17); +//model(sparse); //model; dA = exp(gam+e_a); log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; @@ -93,7 +93,7 @@ options_.solve_tolf=1e-10; options_.maxit_=100; steady; model_info; -check; +//check; shocks; var e_a; periods 1; @@ -102,7 +102,7 @@ end; simul(periods=200, method=lu); -stoch_simul(periods=200,order=1); +//stoch_simul(periods=200,order=1); rplot y; rplot k; diff --git a/tests/ferhat/ireland.mod b/tests/ferhat/ireland.mod index c51cba745..ccb62f713 100644 --- a/tests/ferhat/ireland.mod +++ b/tests/ferhat/ireland.mod @@ -23,8 +23,9 @@ scy = 0.0040; shy = 0.0015; shc = 0.0010; -//model(sparse_dll,no_compiler); -model(sparse); +model(sparse_dll); +//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; @@ -70,15 +71,18 @@ eoh= 0; oy= 7.99331700544506; oc= 7.83132048725623; oh= 5.34253084908048; - +k=k+0.000001; end; + options_.dynatol=1e-12; -options_.maxit_=500; +options_.maxit_=5; options_.slowc=1; -steady(solve_algo=3); -options_.dynatol=4e-8; -check; +steady(solve_algo=2); +//steady; +options_.dynatol=4e-5; + +//check; shocks; @@ -89,7 +93,7 @@ end; options_.maxit_=20; model_info; -simul(periods=2000, method=/*LU*/GMRES/*bicgstab*/); +simul(periods=2000, method=LU/*GMRES*//*bicgstab*/); rplot y; rplot k; diff --git a/tests/ferhat/ls2003.mod b/tests/ferhat/ls2003.mod index 8a8cb8238..ef4366157 100644 --- a/tests/ferhat/ls2003.mod +++ b/tests/ferhat/ls2003.mod @@ -82,6 +82,7 @@ estimation(datafile=data_ca1,first_obs=8,nobs=79,mh_nblocks=10,prefilter=1,mh_js options_.maxit_=100; steady; + model_info; check; @@ -91,7 +92,10 @@ periods 1; values 0.5; end; -simul(periods=200,method=bicgstab); +//simul(periods=200,method=bicgstab); +simul(periods=200); +rplot vv; +rplot ww; rplot A; rplot pie; diff --git a/tests/ferhat/multimod.mod b/tests/ferhat/multimod.mod index 9ad14d376..56cc4fd64 100644 --- a/tests/ferhat/multimod.mod +++ b/tests/ferhat/multimod.mod @@ -852,8 +852,8 @@ W0906=0.0800069594276; W0907=0.147854375051; W0908=0.206834342322; W0909=-1; -//model(SPARSE_DLL,markowitz=2.0,no_compiler); -model(SPARSE); +model(SPARSE_DLL,markowitz=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 ; diff --git a/tests/ferhat/ramst.mod b/tests/ferhat/ramst.mod index 8e64dbf04..32c088065 100644 --- a/tests/ferhat/ramst.mod +++ b/tests/ferhat/ramst.mod @@ -9,9 +9,9 @@ bet=0.05; aa=0.5; -model(sparse); +//model(sparse); //model(sparse_dll); -//model; +model; c + k - aa*x*k(-1)^alph - (1-delt)*k(-1); c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam); end; @@ -24,7 +24,7 @@ end; steady; -check; +//check; model_info; shocks; var x; diff --git a/tests/ferhat/ramst_a.mod b/tests/ferhat/ramst_a.mod index 4a81935c3..4789f3549 100644 --- a/tests/ferhat/ramst_a.mod +++ b/tests/ferhat/ramst_a.mod @@ -10,7 +10,7 @@ bet=0.05; aa=0.5; -model; +model(sparse); c + k - aa*x*k(-1)^alph - (1-delt)*k(-1); c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam); end;