From 1d9260251d66931287b1d0f66adefa1dfccaa3f0 Mon Sep 17 00:00:00 2001 From: ferhat Date: Tue, 21 Jul 2009 15:50:12 +0000 Subject: [PATCH] - sparse_dll option works fine with feedback variables git-svn-id: https://www.dynare.org/svn/dynare/trunk@2851 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/solve_two_boundaries.m | 37 +- mex/sources/simulate/Interpreter.cc | 209 +++++-- mex/sources/simulate/Interpreter.hh | 6 +- mex/sources/simulate/Mem_Mngr.cc | 78 +-- mex/sources/simulate/Mem_Mngr.hh | 8 +- mex/sources/simulate/SparseMatrix.cc | 509 +++++++++++----- mex/sources/simulate/SparseMatrix.hh | 35 +- mex/sources/simulate/simulate.cc | 3 - mex/sources/simulate/simulate.hh | 81 +-- mex/sources/simulate/testing/mex_interface.cc | 15 +- mex/sources/simulate/testing/mex_interface.hh | 2 +- preprocessor/BlockTriangular.cc | 36 +- preprocessor/BlockTriangular.hh | 6 +- preprocessor/DynamicModel.cc | 552 ++++++++++-------- preprocessor/DynamicModel.hh | 14 +- tests/ferhat/fs2000.mod | 12 +- tests/ferhat/ireland.mod | 6 +- tests/ferhat/ls2003.mod | 6 +- tests/ferhat/multimod.mod | 33 +- tests/ferhat/ramst.mod | 4 +- 20 files changed, 971 insertions(+), 681 deletions(-) diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index 68e084590..b34efe742 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -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 . - 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) diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc index 16a068070..5cf8e255d 100644 --- a/mex/sources/simulate/Interpreter.cc +++ b/mex/sources/simulate/Interpreter.cc @@ -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_=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= 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_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;isolve_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=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(nop_all), sizeof(*nop_all)); SaveCode_swp.write(reinterpret_cast(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 #include #include -#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 NZE_Mem_Allocated; int swp_f_b; fstream SaveCode_swp; string filename; diff --git a/mex/sources/simulate/SparseMatrix.cc b/mex/sources/simulate/SparseMatrix.cc index 654edb709..2044db163 100644 --- a/mex/sources/simulate/SparseMatrix.cc +++ b/mex/sources/simulate/SparseMatrix.cc @@ -18,7 +18,7 @@ */ #include - +#include #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;jNZE_R_N; + } + nb_var=NbNZCol[c]; + first=FNZE_C[c]; + for(int j=0;jNZE_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_indexc_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_indexr_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;jNZE_R_N; + } + nb_var=NbNZCol[c]; + first=FNZE_C[c]; + for(int j=0;jNZE_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(&eq), sizeof(eq)); SaveCode.read(reinterpret_cast(&var), sizeof(var)); SaveCode.read(reinterpret_cast(&lag), sizeof(lag)); SaveCode.read(reinterpret_cast(&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;ifirst.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(lagsecond+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=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;tfirst+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;tlag<((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;jr_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;ic_index; + else + column = -9999999; + for(int j=0;ju_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; iepsilon) + 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 ,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; ic_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-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; i0 && 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_absNZE_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;jr_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 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 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_subr_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_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); } diff --git a/mex/sources/simulate/SparseMatrix.hh b/mex/sources/simulate/SparseMatrix.hh index 8f587dfaa..e5b57a76a 100644 --- a/mex/sources/simulate/SparseMatrix.hh +++ b/mex/sources/simulate/SparseMatrix.hh @@ -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 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 ,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; }; diff --git a/mex/sources/simulate/simulate.cc b/mex/sources/simulate/simulate.cc index 0ea7a0087..737afb8b5 100644 --- a/mex/sources/simulate/simulate.cc +++ b/mex/sources/simulate/simulate.cc @@ -178,10 +178,7 @@ main( int argc, const char* argv[] ) mxFree(ya); if(direction) mxFree(direction); - - free(params); - } #else diff --git a/mex/sources/simulate/simulate.hh b/mex/sources/simulate/simulate.hh index 928f87cfb..8e1756562 100644 --- a/mex/sources/simulate/simulate.hh +++ b/mex/sources/simulate/simulate.hh @@ -20,87 +20,20 @@ #ifndef SIMULATE_HH_INCLUDED #define SIMULATE_HH_INCLUDED -/*#include -#include -#include -#include -#include -#include */ -//#include -//#include -/*#include -#include -#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, 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; diff --git a/mex/sources/simulate/testing/mex_interface.cc b/mex/sources/simulate/testing/mex_interface.cc index 09c64f1a7..8f9ff5d20 100644 --- a/mex/sources/simulate/testing/mex_interface.cc +++ b/mex/sources/simulate/testing/mex_interface.cc @@ -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 diff --git a/mex/sources/simulate/testing/mex_interface.hh b/mex/sources/simulate/testing/mex_interface.hh index 4880f365c..3999121db 100644 --- a/mex/sources/simulate/testing/mex_interface.hh +++ b/mex/sources/simulate/testing/mex_interface.hh @@ -5,7 +5,7 @@ #include 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); diff --git a/preprocessor/BlockTriangular.cc b/preprocessor/BlockTriangular.cc index 65813d5bd..3260e5469 100644 --- a/preprocessor/BlockTriangular.cc +++ b/preprocessor/BlockTriangular.cc @@ -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 &equations, } else { + //vector > 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, int> >, int> + +map >, pair >, int> BlockTriangular::get_Derivatives(Model_Block *ModelBlock, int blck) { - map, pair, int> >, int> Derivatives; + map >, pair >, 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(eqBlock_List[blck].Nb_Recursives and varBlock_List[blck].Nb_Recursives) - {*/ bool OK = true; - map, pair, int> >, int>::const_iterator its = Derivatives.find(make_pair(make_pair(eqr, eq), make_pair(make_pair(varr, var), lag))); + map >, pair >, 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 eqBlock_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(eqBlock_List[blck].Nb_Recursives and varBlock_List[blck].Nb_Recursives)*/ if(varBlock_List[blck].Nb_Recursives) { int eqs = ModelBlock->Block_List[blck].Equation[var]; for(int vars=ModelBlock->Block_List[blck].Nb_Recursives; varsBlock_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; } } } diff --git a/preprocessor/BlockTriangular.hh b/preprocessor/BlockTriangular.hh index bd30c3c77..3d6e6dabf 100644 --- a/preprocessor/BlockTriangular.hh +++ b/preprocessor/BlockTriangular.hh @@ -44,7 +44,7 @@ typedef vector > t_vtype; typedef set temporary_terms_inuse_type; -typedef vector, pair > > > chaine_rule_derivatives_type; +typedef vector >, pair > > 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, int> >, int> get_Derivatives(Model_Block *ModelBlock, int Blck); + map >, pair >, int> get_Derivatives(Model_Block *ModelBlock, int Blck); void Normalize_and_BlockDecompose_Static_0_Model(jacob_map &j_m, vector &equations, t_etype &V_Equation_Type, map >, NodeID> &first_order_endo_derivatives); diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 0ce1534b2..27471170f 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -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 >, 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; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + for(i=0; iBlock_List[j].Chain_Rule_Derivatives->size();i++) { - pair< pair, pair > > 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 >, pair > 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; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + for(i=0; iBlock_List[j].Chain_Rule_Derivatives->size();i++) { - pair< pair, pair > > 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 >, pair > 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 (iBlock_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 (iBlock_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;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]; - 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; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + for(i=0; iBlock_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, pair > > 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 >, pair > 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;iBlock_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; iBlock_List[j].Chaine_Rule_Derivatives->size();i++) + for(i=0; iBlock_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, pair > > 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 >, pair > 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 reference_count; - //map ModelBlock_Aggregated_Size, ModelBlock_Aggregated_Number; vector 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(&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(&v),sizeof(v)); v=ModelBlock->Block_List[j].Simulation_Type; code_file.write(reinterpret_cast(&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(&ModelBlock->Block_List[j].Variable[i]),sizeof(ModelBlock->Block_List[j].Variable[i])); code_file.write(reinterpret_cast(&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(&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(&v),sizeof(v)); v=symbol_table.endo_nbr(); code_file.write(reinterpret_cast(&v),sizeof(v)); @@ -945,16 +959,12 @@ ModelBlock_Aggregated_Count = ModelBlock->Size; code_file.write(reinterpret_cast(&v),sizeof(v)); v=block_triangular.ModelBlock->Block_List[j].Max_Lead; code_file.write(reinterpret_cast(&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(&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(&v),sizeof(v)); code_file.write(&FSTPR, sizeof(FSTPR)); - code_file.write(reinterpret_cast(&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(&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; iBlock_List[j].Chain_Rule_Derivatives->size();i++) { - int u = feedback_variables.size(); - for(i=0; iBlock_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 >, pair > 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, pair > > 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(&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;iBlock_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(&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(&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(&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(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); - code_file.write(&FLDV, sizeof(FLDV)); - char vc=eEndogenous; - code_file.write(reinterpret_cast(&vc), sizeof(vc)); - code_file.write(reinterpret_cast(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var)); - int v1=0; - code_file.write(reinterpret_cast(&v1), sizeof(v1)); - code_file.write(&FBINARY, sizeof(FBINARY)); - v1=oTimes; - code_file.write(reinterpret_cast(&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(&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(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); + code_file.write(&FLDV, sizeof(FLDV)); + char vc=eEndogenous; + code_file.write(reinterpret_cast(&vc), sizeof(vc)); + code_file.write(reinterpret_cast(&Uf[v].Ufl->var), sizeof(Uf[v].Ufl->var)); + int v1=0; + code_file.write(reinterpret_cast(&v1), sizeof(v1)); + code_file.write(&FBINARY, sizeof(FBINARY)); + v1=oTimes; + code_file.write(reinterpret_cast(&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(&v), sizeof(v)); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast(&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(&v), sizeof(v)); + code_file.write(&FSTPU, sizeof(FSTPU)); + v = i - ModelBlock->Block_List[j].Nb_Recursives; + code_file.write(reinterpret_cast(&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; iBlock_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 >, pair > 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(&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;iBlock_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(&u), sizeof(u)); -#ifdef CONDITION - output << " if (fabs(condition[" << eqr << "])Block_List[j].Size;i++) { - code_file.write(&FLDR, sizeof(FLDR)); - code_file.write(reinterpret_cast(&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(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); - code_file.write(&FLDV, sizeof(FLDV)); - char vc=eEndogenous; - code_file.write(reinterpret_cast(&vc), sizeof(vc)); - int v1=Uf[v].Ufl->var; - code_file.write(reinterpret_cast(&v1), sizeof(v1)); - v1=Uf[v].Ufl->lag; - code_file.write(reinterpret_cast(&v1), sizeof(v1)); - code_file.write(&FBINARY, sizeof(FBINARY)); - v1=oTimes; - code_file.write(reinterpret_cast(&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(&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(&Uf[v].Ufl->u), sizeof(Uf[v].Ufl->u)); + code_file.write(&FLDV, sizeof(FLDV)); + char vc=eEndogenous; + code_file.write(reinterpret_cast(&vc), sizeof(vc)); + int v1=Uf[v].Ufl->var; + code_file.write(reinterpret_cast(&v1), sizeof(v1)); + v1=Uf[v].Ufl->lag; + code_file.write(reinterpret_cast(&v1), sizeof(v1)); + code_file.write(&FBINARY, sizeof(FBINARY)); + v1=oTimes; + code_file.write(reinterpret_cast(&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(&v), sizeof(v)); - code_file.write(&FSTPU, sizeof(FSTPU)); - code_file.write(reinterpret_cast(&i), sizeof(i)); -#ifdef CONDITION - output << " if (fabs(condition[" << i << "])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(&v), sizeof(v)); + code_file.write(&FSTPU, sizeof(FSTPU)); + v = i - ModelBlock->Block_List[j].Nb_Recursives; + code_file.write(reinterpret_cast(&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; iBlock_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 >, pair > 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(&v), sizeof(v)); + int var=it.first.second.second-block_triangular.ModelBlock->Block_List[num].Nb_Recursives + k * Size; + SaveCode.write(reinterpret_cast(&var), sizeof(var)); + SaveCode.write(reinterpret_cast(&k), sizeof(k)); + int u = u_count_int + Size; + SaveCode.write(reinterpret_cast(&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;jBlock_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(&u), sizeof(u)); u_count_int++; } - } + }*/ if (is_two_boundaries) { - for (j=0;jBlock_List[num].Size;j++) + for (j=0;jBlock_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(&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;jBlock_List[num].Size;j++) + for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;jBlock_List[num].Size;j++) { int varr=block_triangular.ModelBlock->Block_List[num].Variable[j]; SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); } - for (j=0;jBlock_List[num].Size;j++) + for (j=block_triangular.ModelBlock->Block_List[num].Nb_Recursives;jBlock_List[num].Size;j++) { int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j]; SaveCode.write(reinterpret_cast(&eqr1), sizeof(eqr1)); @@ -1548,7 +1590,7 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri mDynamicModelFile << tmp1.str(); tmp1.str(""); } - for (int ik=0;ikBlock_List[i].Size;ik++) + for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_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;ikBlock_List[i].Size;ik++) + for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_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;ikBlock_List[i].Size;ik++) + for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_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;ikBlock_List[i].Size;ik++) + for (int ik=block_triangular.ModelBlock->Block_List[i].Nb_Recursives ;ikBlock_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 recursive_variables; - first_chaine_rule_derivatives.clear(); + first_chain_rule_derivatives.clear(); for(int blck = 0; blckSize; 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,int> > , int > Derivatives = block_triangular.get_Derivatives(ModelBlock, blck); - for(map, pair,int> > , int >::const_iterator it = Derivatives.begin(); it != Derivatives.end(); it++) + map >, pair >, int> Derivatives = block_triangular.get_Derivatives(ModelBlock, blck); + for(map >, pair >, 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 eqBlock_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 >, NodeID>::const_iterator it = first_chaine_rule_derivatives.find(make_pair(eq, make_pair(var, lag))); - if (it != first_chaine_rule_derivatives.end()) + map >, 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; diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh index a93244141..54263fbec 100644 --- a/preprocessor/DynamicModel.hh +++ b/preprocessor/DynamicModel.hh @@ -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: diff --git a/tests/ferhat/fs2000.mod b/tests/ferhat/fs2000.mod index a3dc38fb9..73ab64bbb 100644 --- a/tests/ferhat/fs2000.mod +++ b/tests/ferhat/fs2000.mod @@ -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; diff --git a/tests/ferhat/ireland.mod b/tests/ferhat/ireland.mod index e567a81db..79875d08c 100644 --- a/tests/ferhat/ireland.mod +++ b/tests/ferhat/ireland.mod @@ -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; diff --git a/tests/ferhat/ls2003.mod b/tests/ferhat/ls2003.mod index eac804916..be9abd699 100644 --- a/tests/ferhat/ls2003.mod +++ b/tests/ferhat/ls2003.mod @@ -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; diff --git a/tests/ferhat/multimod.mod b/tests/ferhat/multimod.mod index eadae2655..b4dbc0fe5 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); -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; +*/ \ No newline at end of file diff --git a/tests/ferhat/ramst.mod b/tests/ferhat/ramst.mod index 9fe47f5b8..1b6caf92f 100644 --- a/tests/ferhat/ramst.mod +++ b/tests/ferhat/ramst.mod @@ -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;