diff --git a/matlab/simul.m b/matlab/simul.m index 0f3e4f531..6c58af226 100644 --- a/matlab/simul.m +++ b/matlab/simul.m @@ -78,11 +78,15 @@ if(options_.block) eval([M_.fname '_dynamic']); end; else - if M_.maximum_endo_lag ==1 & M_.maximum_endo_lead <= 1 - sim1 ; + if(options_.bytecode) + oo_.endo_simul=bytecode('dynamic'); else - simk ; - end + if M_.maximum_endo_lag ==1 & M_.maximum_endo_lead <= 1 + sim1 ; + else + simk ; + end; + end; end; dyn2vec; diff --git a/matlab/steady_.m b/matlab/steady_.m index 9eb72303d..7e9422178 100644 --- a/matlab/steady_.m +++ b/matlab/steady_.m @@ -118,7 +118,7 @@ elseif options_.block && ~options_.bytecode [oo_.exo_steady_state; ... oo_.exo_det_steady_state], M_.params); end -elseif options_.block && options_.bytecode +elseif options_.bytecode [oo_.steady_state,check] = bytecode('static'); else [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],... diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index 9b13454a4..5505cfbac 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -21,6 +21,7 @@ #include "Interpreter.hh" #define BIG 1.0e+8; #define SMALL 1.0e-5; +//#define DEBUG Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, double *x_arg, double *steady_y_arg, double *steady_x_arg, double *direction_arg, int y_size_arg, @@ -169,6 +170,11 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) Stack.push(y[var]); #ifdef DEBUG tmp_out << " y[" << var << "](" << y[var] << ")"; + if(var<0 || var>= y_size) + { + mexPrintf("y[%d]=",var); + mexErrMsgTxt("End of bytecode"); + } #endif break; case eExogenous: @@ -238,6 +244,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) //load a temporary variable in the processor var = ((FLDT_ *) it_code->second)->get_pos(); #ifdef DEBUG + mexPrintf("T[it_=%d var=%d, y_kmin=%d, y_kmax=%d == %d]=>%f\n", it_, var, y_kmin, y_kmax, var*(periods+y_kmin+y_kmax)+it_, var); tmp_out << " T[" << it_ << ", " << var << "](" << T[var*(periods+y_kmin+y_kmax)+it_] << ")"; #endif Stack.push(T[var*(periods+y_kmin+y_kmax)+it_]); @@ -274,7 +281,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) break; case FLDZ: //load 0 in the processor - Stack.push(0); + Stack.push(0.0); #ifdef DEBUG tmp_out << " 0"; #endif @@ -397,6 +404,7 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) var += Per_u_; u[var] = Stack.top(); #ifdef DEBUG + mexPrintf("FSTPU\n"); tmp_out << "=>"; mexPrintf(" u[%d](%f)=%s\n", var, u[var], tmp_out.str().c_str()); tmp_out.str(""); @@ -406,6 +414,10 @@ Interpreter::compute_block_time(int Per_u_, bool evaluate, int block_num) case FSTPSU: //store in u variable from the processor var = ((FSTPSU_ *) it_code->second)->get_pos(); +#ifdef DEBUG + if(var>=u_count_alloc || var<0) + mexPrintf("Erreur var=%d\n",var); +#endif u[var] = Stack.top(); #ifdef DEBUG tmp_out << "=>"; @@ -1288,6 +1300,7 @@ Interpreter::simulate_a_block(const int size, const int type, string file_name, Per_u_ = (it_-y_kmin)*u_count_int; Per_y_ = it_*y_size; it_code = begining; + error_not_printed = true; compute_block_time(Per_u_, false, block_num); if (isnan(res1) || isinf(res1)) { diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc index 7917ebb86..7f67fad06 100644 --- a/mex/sources/bytecode/SparseMatrix.cc +++ b/mex/sources/bytecode/SparseMatrix.cc @@ -326,20 +326,12 @@ SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, i for (j = 0; j < Size; j++) SaveCode.read(reinterpret_cast(&index_vara[j]), sizeof(*index_vara)); if (periods+y_kmin+y_kmax > 1) - { - for (i = 1; i < periods+y_kmin+y_kmax; i++) - { - for (j = 0; j < Size; j++) - { - index_vara[j+Size*i] = index_vara[j+Size*(i-1)]+y_size; - } - } - } + for (i = 1; i < periods+y_kmin+y_kmax; i++) + for (j = 0; j < Size; j++) + index_vara[j+Size*i] = index_vara[j+Size*(i-1)]+y_size; index_equa = (int *) mxMalloc(Size*sizeof(int)); for (j = 0; j < Size; j++) - { - SaveCode.read(reinterpret_cast(&index_equa[j]), sizeof(*index_equa)); - } + SaveCode.read(reinterpret_cast(&index_equa[j]), sizeof(*index_equa)); } void @@ -367,13 +359,11 @@ SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, mapc_index]+cal_y]*u[first->u_index]; - tmp_out << "y[" << index_vara[first->c_index]+cal_y << "]" << "(" << y[index_vara[first->c_index]+cal_y] << ")*u[" << first->u_index << "](" << u[first->u_index] << ")+"; first = first->NZE_R_N; } yy = -(yy+y[eq]+u[b[pos]]); @@ -994,13 +981,9 @@ SparseMatrix::simple_bksub(int it_, int Size, double slowc_l) nb_var--; int eq = index_vara[i]; yy = 0; - ostringstream tmp_out; - tmp_out.clear(); - tmp_out << "|"; for (k = 0; k < nb_var; k++) { yy += y[index_vara[first->c_index]+it_*y_size]*u[first->u_index]; - tmp_out << "y[" << index_vara[first->c_index]+it_*y_size << "](" << y[index_vara[first->c_index]+it_*y_size] << ")*u[" << first->u_index << "](" << u[first->u_index] << ")+"; first = first->NZE_R_N; } yy = -(yy+y[eq+it_*y_size]+u[b[pos]]); @@ -1046,7 +1029,7 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, else mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); } - mexPrintf("res1=%5.25\n", res1); + mexPrintf("res1=%5.10f\n", res1); mexPrintf("The initial values of endogenous variables are too far from the solution.\n"); mexPrintf("Change them!\n"); mexEvalString("drawnow;"); @@ -1063,35 +1046,40 @@ SparseMatrix::simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, mexErrMsgTxt(filename.c_str()); } } - if (slowc_save < 1e-8) + if (fabs(slowc_save) < 1e-8) { - for (j = 0; j < y_size; j++) - { - bool select = false; - for (int i = 0; i < Size; i++) - if (j == index_vara[i]) - { - select = true; - break; - } - if (select) - mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - else - mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); - } - mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); - mexEvalString("drawnow;"); - mxFree(piv_v); - mxFree(pivj_v); - mxFree(pivk_v); - mxFree(NR); - if (steady_state) - return false; + if (slowc_save>0) + slowc_save = -slowc; else { - mexEvalString("st=fclose('all');clear all;"); - filename += " stopped"; - mexErrMsgTxt(filename.c_str()); + for (j = 0; j < y_size; j++) + { + bool select = false; + for (int i = 0; i < Size; i++) + if (j == index_vara[i]) + { + select = true; + break; + } + if (select) + mexPrintf("-> variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); + else + mexPrintf(" variable %d at time %d = %f direction = %f\n", j+1, it_, y[j+it_*y_size], direction[j+it_*y_size]); + } + mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); + mexEvalString("drawnow;"); + mxFree(piv_v); + mxFree(pivj_v); + mxFree(pivk_v); + mxFree(NR); + if (steady_state) + return false; + else + { + mexEvalString("st=fclose('all');clear all;"); + filename += " stopped"; + mexErrMsgTxt(filename.c_str()); + } } } slowc_save /= 2; @@ -1484,8 +1472,8 @@ SparseMatrix::simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax { mexPrintf("slowc_save=%g\n", slowc_save); for (j = 0; j < y_size; j++) - mexPrintf("variable %d at time %d = %f\n", j+1, it_, y[j+it_*y_size]); - mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d)\n", blck+1, it_+1, max_res_idx); + mexPrintf("variable %d at time %d = %f direction = %f variable at last step = %f b = %f\n", j+1, it_+1, y[j+it_*y_size], direction[j+it_*y_size], ya[j+it_*y_size], u[pivot[j+it_*y_size]]); + mexPrintf("Dynare cannot improve the simulation in block %d at time %d (variable %d max_res = %f, res1 = %f)\n", blck+1, it_+1, max_res_idx, max_res, res1); mexEvalString("drawnow;"); mexEvalString("st=fclose('all');clear all;"); filename += " stopped"; diff --git a/mex/sources/bytecode/SparseMatrix.hh b/mex/sources/bytecode/SparseMatrix.hh index cc93a568f..fd24b22f3 100644 --- a/mex/sources/bytecode/SparseMatrix.hh +++ b/mex/sources/bytecode/SparseMatrix.hh @@ -165,12 +165,12 @@ private: int *g_save_op; int first_count_loop; int g_nop_all; - int u_count_alloc, u_count_alloc_save; double markowitz_c_s; double res1a; long int nop_all, nop1, nop2; map, int>, int> IM_i; protected: + int u_count_alloc, u_count_alloc_save; double *u, *y, *ya; double res1, res2, max_res, max_res_idx; double slowc, slowc_save, markowitz_c; diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc index 40d2b8e41..cef457ec4 100644 --- a/mex/sources/bytecode/bytecode.cc +++ b/mex/sources/bytecode/bytecode.cc @@ -88,8 +88,9 @@ main(int argc, const char *argv[]) } } fid = fopen(tmp_out.str().c_str(), "r"); - int periods; - fscanf(fid, "%d", &periods); + int periods = 1; + if (!steady_state) + fscanf(fid, "%d", &periods); int maxit_; fscanf(fid, "%d", &maxit_); fscanf(fid, "%f", &f_tmp); diff --git a/mex/sources/bytecode/testing/simulate_debug.m b/mex/sources/bytecode/testing/simulate_debug.m index 1550585d8..f38aeca98 100644 --- a/mex/sources/bytecode/testing/simulate_debug.m +++ b/mex/sources/bytecode/testing/simulate_debug.m @@ -1,4 +1,4 @@ -function simulate_debug() +function simulate_debug(steady_state) global M_ oo_ options_; fid = fopen([M_.fname '_options.txt'],'wt'); fprintf(fid,'%d\n',options_.periods); @@ -14,10 +14,19 @@ fprintf(fid,'%d\n',M_.maximum_lag); fprintf(fid,'%d\n',M_.maximum_lead); fprintf(fid,'%d\n',M_.maximum_endo_lag); fprintf(fid,'%d\n',M_.param_nbr); -fprintf(fid,'%d\n',size(oo_.exo_simul, 1)); -fprintf(fid,'%d\n',size(oo_.exo_simul, 2)); +if steady_state==1 + fprintf(fid,'%d\n',size(oo_.exo_steady_state, 1)); + fprintf(fid,'%d\n',size(oo_.exo_steady_state, 2)); +else + fprintf(fid,'%d\n',size(oo_.exo_simul, 1)); + fprintf(fid,'%d\n',size(oo_.exo_simul, 2)); +end; fprintf(fid,'%d\n',M_.endo_nbr); -fprintf(fid,'%d\n',size(oo_.endo_simul, 2)); +if steady_state==1 + fprintf(fid,'%d\n',size(oo_.steady_state, 2)); +else + fprintf(fid,'%d\n',size(oo_.endo_simul, 2)); +end; fprintf(fid,'%d\n',M_.exo_det_nbr); fprintf(fid,'%d\n',size(oo_.steady_state,1)); @@ -29,8 +38,13 @@ fprintf(fid,'%6.20f\n',M_.params); fclose(fid); fid = fopen([M_.fname '_oo.txt'],'wt'); -fprintf(fid,'%6.20f\n',oo_.endo_simul); -fprintf(fid,'%6.20f\n',oo_.exo_simul); +if steady_state==1 + fprintf(fid,'%6.20f\n',oo_.steady_state); + fprintf(fid,'%6.20f\n',oo_.exo_steady_state); +else + fprintf(fid,'%6.20f\n',oo_.endo_simul); + fprintf(fid,'%6.20f\n',oo_.exo_simul); +end; fprintf(fid,'%6.20f\n',oo_.steady_state); fprintf(fid,'%6.20f\n',oo_.exo_steady_state); fclose(fid); diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh index 90774f6a6..e86665814 100644 --- a/preprocessor/CodeInterpreter.hh +++ b/preprocessor/CodeInterpreter.hh @@ -82,6 +82,7 @@ enum Tags FSTPR, //!< Loads a residual from the stack - 12 FSTPG, //!< Loads a derivative from the stack - 13 + FSTPG2, //!< Loads a derivative matrix from the stack - 13 FUNARY, //!< A Unary operator - 14 FBINARY, //!< A binary operator - 15 @@ -533,6 +534,28 @@ public: }; }; +class FSTPG2_ : public TagWithTwoArguments +{ +public: + inline FSTPG2_() : TagWithTwoArguments::TagWithTwoArguments(FSTPG2, 0, 0) + { + }; + inline FSTPG2_(const unsigned int pos_arg1, const unsigned int pos_arg2) : TagWithTwoArguments::TagWithTwoArguments(FSTPG2, pos_arg1, pos_arg2) + { + }; + inline unsigned int + get_row() + { + return arg1; + }; + inline unsigned int + get_col() + { + return arg2; + }; +}; + + class FUNARY_ : public TagWithOneArgument { public: @@ -728,7 +751,7 @@ public: op_code = FBEGINBLOCK; size = 0; type = UNKNOWN; /*variable = NULL; equation = NULL;*/ is_linear = false; endo_nbr = 0; Max_Lag = 0; Max_Lead = 0; u_count_int = 0; }; - inline FBEGINBLOCK_(unsigned int &size_arg, BlockSimulationType &type_arg, int unsigned first_element, int unsigned &block_size, + inline FBEGINBLOCK_(unsigned int size_arg, BlockSimulationType type_arg, int unsigned first_element, int unsigned block_size, const vector &variable_arg, const vector &equation_arg, bool is_linear_arg, int endo_nbr_arg, int Max_Lag_arg, int Max_Lead_arg, int &u_count_int_arg) { diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 20118f1eb..ff26e0bf2 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -878,7 +878,7 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct) void PlannerObjectiveStatement::computingPass() { - model_tree->computingPass(eval_context_type(), false, true, false); + model_tree->computingPass(eval_context_type(), false, true, false, false); } void diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 9e39502b8..21058288f 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -81,6 +81,18 @@ DynamicModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, } } +void +DynamicModel::initializeVariablesAndEquations() +{ + for(int j=0; jsecond->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); v_temporary_terms_inuse[block] = temporary_terms_in_use; } - // Add a mapping form node ID to temporary terms order - int j = 0; - for (temporary_terms_type::const_iterator it = temporary_terms.begin(); - it != temporary_terms.end(); it++) - map_idx[(*it)->idx] = j++; + computeTemporaryTermsMapping(); } } +void +DynamicModel::computeTemporaryTermsMapping() +{ + // Add a mapping form node ID to temporary terms order + int j = 0; + for (temporary_terms_type::const_iterator it = temporary_terms.begin(); + it != temporary_terms.end(); it++) + map_idx[(*it)->idx] = j++; +} + + void DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const { @@ -710,7 +729,118 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const } void -DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const +DynamicModel::writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_type map_idx) const +{ + ostringstream tmp_output; + ofstream code_file; + bool file_open = false; + string main_name = file_name; + + main_name += ".cod"; + code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); + if (!code_file.is_open()) + { + cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + exit(EXIT_FAILURE); + } + + + int count_u; + int u_count_int = 0; + BlockSimulationType simulation_type; + if ((max_endo_lag > 0) && (max_endo_lead > 0)) + simulation_type = SOLVE_TWO_BOUNDARIES_COMPLETE; + else if ((max_endo_lag >= 0) && (max_endo_lead == 0)) + simulation_type = SOLVE_FORWARD_COMPLETE; + else + simulation_type = SOLVE_BACKWARD_COMPLETE; + + Write_Inf_To_Bin_File(file_name, u_count_int, file_open, simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE, symbol_table.endo_nbr() ); + file_open = true; + + //Temporary variables declaration + FDIMT_ fdimt(temporary_terms.size()); + fdimt.write(code_file); + + FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(), + simulation_type, + 0, + symbol_table.endo_nbr(), + variable_reordered, + equation_reordered, + false, + symbol_table.endo_nbr(), + 0, + 0, + u_count_int + ); + fbeginblock.write(code_file); + + compileTemporaryTerms(code_file, temporary_terms, map_idx, true, false); + + compileModelEquations(code_file, temporary_terms, map_idx, true, false); + + FENDEQU_ fendequ; + fendequ.write(code_file); + vector, int > > > derivatives; + derivatives.resize(symbol_table.endo_nbr()); + count_u = symbol_table.endo_nbr(); + for (first_derivatives_type::const_iterator it = first_derivatives.begin(); + it != first_derivatives.end(); it++) + { + int deriv_id = it->first.second; + if (getTypeByDerivID(deriv_id) == eEndogenous) + { + NodeID d1 = it->second; + unsigned int eq = it->first.first; + int symb = getSymbIDByDerivID(deriv_id); + unsigned int var = symbol_table.getTypeSpecificID(symb); + int lag = getLagByDerivID(deriv_id); + if (!derivatives[eq].size()) + derivatives[eq].clear(); + derivatives[eq].push_back(make_pair(make_pair(var, lag), count_u)); + d1->compile(code_file, false, temporary_terms, map_idx, true, false); + + FSTPU_ fstpu(count_u); + fstpu.write(code_file); + count_u++; + } + } + for (int i = 0; i < symbol_table.endo_nbr(); i++) + { + FLDR_ fldr(i); + fldr.write(code_file); + for(vector, int> >::const_iterator it = derivatives[i].begin(); + it != derivatives[i].end(); it++) + { + FLDU_ fldu(it->second); + fldu.write(code_file); + FLDV_ fldv(eEndogenous, it->first.first, it->first.second); + fldv.write(code_file); + FBINARY_ fbinary(oTimes); + fbinary.write(code_file); + if (it != derivatives[i].begin()) + { + FBINARY_ fbinary(oPlus); + fbinary.write(code_file); + } + } + FBINARY_ fbinary(oMinus); + fbinary.write(code_file); + FSTPU_ fstpu(i); + fstpu.write(code_file); + } + FENDBLOCK_ fendblock; + fendblock.write(code_file); + FEND_ fend; + fend.write(code_file); + code_file.close(); +} + + + +void +DynamicModel::writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_type map_idx) const { struct Uff_l { @@ -766,7 +896,7 @@ DynamicModel::writeModelEquationsCodeOrdered(const string file_name, const strin if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE) { - Write_Inf_To_Bin_File(file_name, bin_basename, block, u_count_int, file_open, + Write_Inf_To_Bin_File_Block(file_name, bin_basename, block, u_count_int, file_open, simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE); file_open = true; } @@ -1109,7 +1239,7 @@ DynamicModel::reform(const string name1) const } void -DynamicModel::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, const int &num, +DynamicModel::Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename, const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const { int j; @@ -1940,7 +2070,7 @@ DynamicModel::collect_first_order_derivatives_endogenous() void DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives, - const eval_context_type &eval_context, bool no_tmp_terms, bool block, bool use_dll) + const eval_context_type &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode) { assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivatives)); @@ -2025,7 +2155,11 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative } else if (!no_tmp_terms) - computeTemporaryTerms(!use_dll); + { + computeTemporaryTerms(!use_dll); + if (bytecode) + computeTemporaryTermsMapping(); + } } map >, pair >, int> @@ -2274,7 +2408,9 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode { int r; if (block && bytecode) - writeModelEquationsCodeOrdered(basename + "_dynamic", basename, map_idx); + writeModelEquationsCode_Block(basename + "_dynamic", basename, map_idx); + else if (!block && bytecode) + writeModelEquationsCode(basename + "_dynamic", basename, map_idx); else if (block && !bytecode) { #ifdef _WIN32 diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh index 32da867a9..0c40f5dfb 100644 --- a/preprocessor/DynamicModel.hh +++ b/preprocessor/DynamicModel.hh @@ -96,7 +96,10 @@ private: //! Writes the Block reordred structure of the model in M output void writeModelEquationsOrdered_M(const string &dynamic_basename) const; //! Writes the code of the Block reordred structure of the model in virtual machine bytecode - void writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const; + void writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_type map_idx) const; + //! Writes the code of the model in virtual machine bytecode + void writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_type map_idx) const; + //! Computes jacobian and prepares for equation normalization /*! Using values from initval/endval blocks and parameter initializations: - computes the jacobian for the model w.r. to contemporaneous variables @@ -112,7 +115,11 @@ private: string reform(string name) const; map_idx_type map_idx; + //! sorts the temporary terms in the blocks order void computeTemporaryTermsOrdered(); + + //! creates a mapping from the index of temporary terms to a natural index + void computeTemporaryTermsMapping(); //! Write derivative code of an equation w.r. to a variable void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const; //! Write chain rule derivative code of an equation w.r. to a variable @@ -217,12 +224,12 @@ public: \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */ void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives, - const eval_context_type &eval_context, bool no_tmp_terms, bool block, bool use_dll); + const eval_context_type &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode); //! Writes model initialization and lead/lag incidence matrix to output void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll) const; //! Adds informations for simulation in a binary file - void Write_Inf_To_Bin_File(const string &dynamic_basename, const string &bin_basename, + void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename, const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const; //! Writes dynamic model file void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll) const; @@ -235,6 +242,9 @@ public: //! Writes LaTeX file with the equations of the dynamic model void writeLatexFile(const string &basename) const; + //! Initialize equation_reordered & variable_reordered + void initializeVariablesAndEquations(); + virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException); diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh index f305c6111..4403e24c5 100644 --- a/preprocessor/ExprNode.hh +++ b/preprocessor/ExprNode.hh @@ -108,6 +108,7 @@ class ExprNode friend class DataTree; friend class DynamicModel; friend class StaticModel; + friend class ModelTree; friend class ExprNodeLess; friend class NumConstNode; friend class VariableNode; diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index 0597b125a..74016ef17 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -129,11 +129,6 @@ ModFile::checkPass() exit(EXIT_FAILURE); } - if (byte_code && !block) - { - cerr << "ERROR: In 'model' block, can't use option 'bytecode' without option 'block'" << endl; - exit(EXIT_FAILURE); - } if ( (stochastic_statement_present || mod_file_struct.check_present || mod_file_struct.steady_present) && no_static) { cerr << "no_static option is incompatible with stochastic simulation, estimation, optimal policy, steady or check command" << endl; @@ -193,12 +188,18 @@ ModFile::computingPass(bool no_tmp_terms) // Compute static model and its derivatives dynamic_model.toStatic(static_model); if(!no_static) - static_model.computingPass(global_eval_context, no_tmp_terms, false, block); + { + static_model.initializeVariablesAndEquations(); + static_model.computingPass(global_eval_context, no_tmp_terms, false, block, byte_code); + } // Set things to compute for dynamic model if (dynamic_model_needed) { if (mod_file_struct.simul_present) - dynamic_model.computingPass(false, false, false, false, global_eval_context, no_tmp_terms, block, use_dll); + { + dynamic_model.initializeVariablesAndEquations(); + dynamic_model.computingPass(false, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code); + } else { if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3) @@ -209,11 +210,11 @@ ModFile::computingPass(bool no_tmp_terms) bool hessian = mod_file_struct.order_option >= 2 || mod_file_struct.identification_present; bool thirdDerivatives = mod_file_struct.order_option == 3; bool paramsDerivatives = mod_file_struct.identification_present; - dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivatives, global_eval_context, no_tmp_terms, false, use_dll); + dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivatives, global_eval_context, no_tmp_terms, false, use_dll, byte_code); } } else - dynamic_model.computingPass(true, true, false, false, global_eval_context, no_tmp_terms, false, false); + dynamic_model.computingPass(true, true, false, false, global_eval_context, no_tmp_terms, false, false, byte_code); } for (vector::iterator it = statements.begin(); diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index d7d6f1e28..4490a3385 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -1014,6 +1014,32 @@ ModelTree::writeTemporaryTerms(const temporary_terms_type &tt, ostream &output, output << ";" << endl; } +void +ModelTree::compileTemporaryTerms(ostream &code_file, const temporary_terms_type &tt, map_idx_type map_idx, bool dynamic, bool steady_dynamic) const +{ + // Local var used to keep track of temp nodes already written + temporary_terms_type tt2; + for (temporary_terms_type::const_iterator it = tt.begin(); + it != tt.end(); it++) + { + (*it)->compile(code_file, false, tt2, map_idx, dynamic, steady_dynamic); + if (dynamic) + { + FSTPT_ fstpt((int)(map_idx.find((*it)->idx)->second)); + fstpt.write(code_file); + } + else + { + FSTPST_ fstpst((int)(map_idx.find((*it)->idx)->second)); + fstpst.write(code_file); + } + // Insert current node into tt2 + tt2.insert(*it); + } +} + + + void ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const { @@ -1079,6 +1105,89 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type) } } +void +ModelTree::compileModelEquations(ostream &code_file, const temporary_terms_type &tt, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const +{ + for (int eq = 0; eq < (int) equations.size(); eq++) + { + BinaryOpNode *eq_node = equations[eq]; + NodeID lhs = eq_node->get_arg1(); + NodeID rhs = eq_node->get_arg2(); + + // Test if the right hand side of the equation is empty. + double vrhs = 1.0; + try + { + vrhs = rhs->eval(eval_context_type()); + } + catch (ExprNode::EvalException &e) + { + } + + if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs; + { + lhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic); + rhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic); + + FBINARY_ fbinary(oMinus); + fbinary.write(code_file); + + FSTPR_ fstpr(eq); + fstpr.write(code_file); + } + else // The right hand side of the equation is empty ==> residual=lhs; + { + lhs->compile(code_file, false, temporary_terms, map_idx, dynamic, steady_dynamic); + FSTPR_ fstpr(eq); + fstpr.write(code_file); + } + } +} + +void +ModelTree::Write_Inf_To_Bin_File(const string &basename, + int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const +{ + int j; + std::ofstream SaveCode; + const string bin_basename = basename + ".bin"; + if (file_open) + SaveCode.open(bin_basename.c_str(), ios::out | ios::in | ios::binary | ios::ate); + else + SaveCode.open(bin_basename.c_str(), ios::out | ios::binary); + if (!SaveCode.is_open()) + { + cout << "Error : Can't open file \"" << bin_basename << "\" for writing\n"; + exit(EXIT_FAILURE); + } + u_count_int = 0; + for (first_derivatives_type::const_iterator it = first_derivatives.begin();it != first_derivatives.end(); it++) + { + int deriv_id = it->first.second; + if (getTypeByDerivID(deriv_id) == eEndogenous) + { + int eq = it->first.first; + int symb = getSymbIDByDerivID(deriv_id); + int var = symbol_table.getTypeSpecificID(symb); + int lag = getLagByDerivID(deriv_id); + SaveCode.write(reinterpret_cast(&eq), sizeof(eq)); + int varr = var + lag * block_mfs; + SaveCode.write(reinterpret_cast(&varr), sizeof(varr)); + SaveCode.write(reinterpret_cast(&lag), sizeof(lag)); + int u = u_count_int + block_mfs; + SaveCode.write(reinterpret_cast(&u), sizeof(u)); + u_count_int++; + } + } + if (is_two_boundaries) + u_count_int += symbol_table.endo_nbr(); + for (j = 0; j < (int) symbol_table.endo_nbr(); j++) + SaveCode.write(reinterpret_cast(&j), sizeof(j)); + for (j = 0; j < (int) symbol_table.endo_nbr(); j++) + SaveCode.write(reinterpret_cast(&j), sizeof(j)); + SaveCode.close(); +} + void ModelTree::writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const { diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh index 4614a12c6..4bd6a0efa 100644 --- a/preprocessor/ModelTree.hh +++ b/preprocessor/ModelTree.hh @@ -109,11 +109,18 @@ protected: void computeTemporaryTerms(bool is_matlab); //! Writes temporary terms void writeTemporaryTerms(const temporary_terms_type &tt, ostream &output, ExprNodeOutputType output_type) const; + //! Compiles temporary terms + void compileTemporaryTerms(ostream &code_file, const temporary_terms_type &tt, map_idx_type map_idx, bool dynamic, bool steady_dynamic) const; + //! Adds informations for simulation in a binary file + void Write_Inf_To_Bin_File(const string &basename, int &u_count_int, bool &file_open, bool is_two_boundaries, int block_mfs) const; + //! Writes model local variables /*! No temporary term is used in the output, so that local parameters declarations can be safely put before temporary terms declaration in the output files */ void writeModelLocalVariables(ostream &output, ExprNodeOutputType output_type) const; //! Writes model equations void writeModelEquations(ostream &output, ExprNodeOutputType output_type) const; + //! Compiles model equations + void compileModelEquations(ostream &code_file, const temporary_terms_type &tt, map_idx_type &map_idx, bool dynamic, bool steady_dynamic) const; //! Writes LaTeX model file void writeLatexModelFile(const string &filename, ExprNodeOutputType output_type) const; diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 6bce4e2b1..581f85326 100644 --- a/preprocessor/StaticModel.cc +++ b/preprocessor/StaticModel.cc @@ -70,6 +70,16 @@ StaticModel::compileChainRuleDerivative(ofstream &code_file, int eqr, int varr, } } +void +StaticModel::initializeVariablesAndEquations() +{ + for(int j = 0; j < equation_number(); j++) + { + equation_reordered.push_back(j); + variable_reordered.push_back(j); + } +} + void StaticModel::computeTemporaryTermsOrdered() { @@ -173,11 +183,17 @@ StaticModel::computeTemporaryTermsOrdered() (*it)->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block); v_temporary_terms_inuse[block] = temporary_terms_in_use; } + computeTemporaryTermsMapping(); } +} + +void +StaticModel::computeTemporaryTermsMapping() +{ // Add a mapping form node ID to temporary terms order int j = 0; for (temporary_terms_type::const_iterator it = temporary_terms.begin(); - it != temporary_terms.end(); it++) + it != temporary_terms.end(); it++) map_idx[(*it)->idx] = j++; } @@ -388,7 +404,115 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const } void -StaticModel::writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const +StaticModel::writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_type map_idx) const +{ + + ostringstream tmp_output; + ofstream code_file; + bool file_open = false; + + string main_name = file_name; + main_name += ".cod"; + code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); + if (!code_file.is_open()) + { + cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + exit(EXIT_FAILURE); + } + int count_u; + int u_count_int = 0; + + Write_Inf_To_Bin_File(file_name, u_count_int, file_open, false, symbol_table.endo_nbr()); + file_open = true; + + //Temporary variables declaration + FDIMT_ fdimt(temporary_terms.size()); + fdimt.write(code_file); + + FBEGINBLOCK_ fbeginblock(symbol_table.endo_nbr(), + SOLVE_FORWARD_COMPLETE, + 0, + symbol_table.endo_nbr(), + variable_reordered, + equation_reordered, + false, + symbol_table.endo_nbr(), + 0, + 0, + u_count_int + ); + fbeginblock.write(code_file); + + + // Add a mapping form node ID to temporary terms order + int j = 0; + for (temporary_terms_type::const_iterator it = temporary_terms.begin(); + it != temporary_terms.end(); it++) + map_idx[(*it)->idx] = j++; + compileTemporaryTerms(code_file, temporary_terms, map_idx, false, false); + + compileModelEquations(code_file, temporary_terms, map_idx, false, false); + + FENDEQU_ fendequ; + fendequ.write(code_file); + + vector > > derivatives; + derivatives.resize(symbol_table.endo_nbr()); + count_u = symbol_table.endo_nbr(); + for (first_derivatives_type::const_iterator it = first_derivatives.begin(); + it != first_derivatives.end(); it++) + { + int deriv_id = it->first.second; + if (getTypeByDerivID(deriv_id) == eEndogenous) + { + NodeID d1 = it->second; + unsigned int eq = it->first.first; + int symb = getSymbIDByDerivID(deriv_id); + unsigned int var = symbol_table.getTypeSpecificID(symb); + if (!derivatives[eq].size()) + derivatives[eq].clear(); + derivatives[eq].push_back(make_pair(var, count_u)); + + d1->compile(code_file, false, temporary_terms, map_idx, false, false); + + FSTPSU_ fstpsu(count_u); + fstpsu.write(code_file); + count_u++; + } + } + for (int i = 0; i < symbol_table.endo_nbr(); i++) + { + FLDR_ fldr(i); + fldr.write(code_file); + for(vector >::const_iterator it = derivatives[i].begin(); + it != derivatives[i].end(); it++) + { + FLDSU_ fldsu(it->second); + fldsu.write(code_file); + FLDSV_ fldsv(eEndogenous, it->first); + fldsv.write(code_file); + FBINARY_ fbinary(oTimes); + fbinary.write(code_file); + if (it != derivatives[i].begin()) + { + FBINARY_ fbinary(oPlus); + fbinary.write(code_file); + } + } + FBINARY_ fbinary(oMinus); + fbinary.write(code_file); + FSTPSU_ fstpsu(i); + fstpsu.write(code_file); + } + FENDBLOCK_ fendblock; + fendblock.write(code_file); + FEND_ fend; + fend.write(code_file); + code_file.close(); +} + +void +StaticModel::writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_type map_idx) const { struct Uff_l { @@ -443,7 +567,7 @@ StaticModel::writeModelEquationsCodeOrdered(const string file_name, const string if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE) { - Write_Inf_To_Bin_File(file_name, bin_basename, block, u_count_int, file_open); + Write_Inf_To_Bin_File_Block(file_name, bin_basename, block, u_count_int, file_open); file_open = true; } @@ -629,7 +753,7 @@ StaticModel::writeModelEquationsCodeOrdered(const string file_name, const string } void -StaticModel::Write_Inf_To_Bin_File(const string &static_basename, const string &bin_basename, const int &num, +StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num, int &u_count_int, bool &file_open) const { int j; @@ -697,7 +821,7 @@ StaticModel::collect_first_order_derivatives_endogenous() } void -StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_terms, bool hessian, bool block) +StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode) { // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic set vars; @@ -753,11 +877,16 @@ StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_te global_temporary_terms = true; if (!no_tmp_terms) computeTemporaryTermsOrdered(); - } else - if (!no_tmp_terms) - computeTemporaryTerms(true); + { + if (!no_tmp_terms) + { + computeTemporaryTerms(true); + if (bytecode) + computeTemporaryTermsMapping(); + } + } } void @@ -894,7 +1023,9 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode) exit(EXIT_FAILURE); } if (block && bytecode) - writeModelEquationsCodeOrdered(basename + "_static", basename, map_idx); + writeModelEquationsCode_Block(basename + "_static", basename, map_idx); + else if (!block && bytecode) + writeModelEquationsCode(basename + "_static", basename, map_idx); else if (block && !bytecode) { chdir(basename.c_str()); diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh index 1dbfb07d1..8aee42427 100644 --- a/preprocessor/StaticModel.hh +++ b/preprocessor/StaticModel.hh @@ -61,7 +61,11 @@ private: void writeModelEquationsOrdered_M(const string &dynamic_basename) const; //! Writes the code of the Block reordred structure of the model in virtual machine bytecode - void writeModelEquationsCodeOrdered(const string file_name, const string bin_basename, map_idx_type map_idx) const; + void writeModelEquationsCode_Block(const string file_name, const string bin_basename, map_idx_type map_idx) const; + + //! Writes the code of the model in virtual machine bytecode + void writeModelEquationsCode(const string file_name, const string bin_basename, map_idx_type map_idx) const; + //! Computes jacobian and prepares for equation normalization /*! Using values from initval/endval blocks and parameter initializations: @@ -72,7 +76,11 @@ private: map_idx_type map_idx; + //! sorts the temporary terms in the blocks order void computeTemporaryTermsOrdered(); + //! creates a mapping from the index of temporary terms to a natural index + void computeTemporaryTermsMapping(); + //! Write derivative code of an equation w.r. to a variable void compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const; //! Write chain rule derivative code of an equation w.r. to a variable @@ -167,10 +175,10 @@ public: \param eval_context evaluation context for normalization \param no_tmp_terms if true, no temporary terms will be computed in the static files */ - void computingPass(const eval_context_type &eval_context, bool no_tmp_terms, bool hessian, bool block); + void computingPass(const eval_context_type &eval_context, bool no_tmp_terms, bool hessian, bool block, bool bytecode); - //! Adds informations for simulation in a binary file - void Write_Inf_To_Bin_File(const string &static_basename, const string &bin_basename, const int &num, + //! Adds informations for simulation in a binary file for a block decomposed model + void Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num, int &u_count_int, bool &file_open) const; //! Writes static model file @@ -182,6 +190,9 @@ public: //! Writes initializations in oo_.steady_state for the auxiliary variables void writeAuxVarInitval(ostream &output) const; + //! Initialize equation_reordered & variable_reordered + void initializeVariablesAndEquations(); + virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); //! Return the number of blocks