From 24b2993f50b000056d234d89b45ad616c10b5d19 Mon Sep 17 00:00:00 2001 From: ferhat Date: Wed, 17 Dec 2014 09:37:43 +0100 Subject: [PATCH] Solves issues in deterministic simulation: - replaces maximum_endo_(lead|lag) by maximum_(lead|lag) to determine the maximum number of lead and lag in deterministic simulation - allows to use bytecode in solve_perfect_foresight_model.m - Adds model information in bytecode --- matlab/make_y_.m | 6 +-- matlab/perfect_foresight_solver.m | 4 +- matlab/solve_perfect_foresight_model.m | 5 +++ mex/sources/bytecode/Evaluate.cc | 11 +++-- mex/sources/bytecode/Interpreter.cc | 8 ++-- preprocessor/CodeInterpreter.hh | 28 +++++++++--- preprocessor/DynamicModel.cc | 60 +++++++++++++++++++------- 7 files changed, 89 insertions(+), 33 deletions(-) diff --git a/matlab/make_y_.m b/matlab/make_y_.m index ae8928693..788d6eb93 100644 --- a/matlab/make_y_.m +++ b/matlab/make_y_.m @@ -43,14 +43,14 @@ end if isempty(M_.endo_histval) if isempty(ys0_) - oo_.endo_simul = [oo_.steady_state*ones(1,M_.maximum_endo_lag+options_.periods+M_.maximum_endo_lead)]; + oo_.endo_simul = [oo_.steady_state*ones(1,M_.maximum_lag+options_.periods+M_.maximum_lead)]; else - oo_.endo_simul = [ys0_*ones(1,M_.maximum_endo_lag) oo_.steady_state*ones(1,options_.periods+M_.maximum_endo_lead)]; + oo_.endo_simul = [ys0_*ones(1,M_.maximum_lag) oo_.steady_state*ones(1,options_.periods+M_.maximum_lead)]; end else if ~isempty(ys0_) error('histval and endval cannot be used simultaneously') end oo_.endo_simul = [M_.endo_histval ... - oo_.steady_state*ones(1,options_.periods+M_.maximum_endo_lead)]; + oo_.steady_state*ones(1,options_.periods+M_.maximum_lead)]; end diff --git a/matlab/perfect_foresight_solver.m b/matlab/perfect_foresight_solver.m index bfeb16e4f..db0f265d6 100644 --- a/matlab/perfect_foresight_solver.m +++ b/matlab/perfect_foresight_solver.m @@ -53,7 +53,7 @@ if isoctave && options_.stack_solve_algo == 2 end -if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maximum_endo_lag+options_.periods+M_.maximum_endo_lead ]) +if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maximum_lag+options_.periods+M_.maximum_lead ]) error('PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?') end @@ -91,7 +91,7 @@ else exosim = oo_.exo_simul; exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1); endosim = oo_.endo_simul; - endoinit = repmat(oo_.steady_state, 1,M_.maximum_endo_lag+options_.periods+M_.maximum_endo_lead); + endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead); current_weight = 0; % Current weight of target point in convex combination step = 1; diff --git a/matlab/solve_perfect_foresight_model.m b/matlab/solve_perfect_foresight_model.m index 4a6464089..5de112ad8 100644 --- a/matlab/solve_perfect_foresight_model.m +++ b/matlab/solve_perfect_foresight_model.m @@ -32,6 +32,11 @@ function [flag,endo_simul,err] = solve_perfect_foresight_model(endo_simul,exo_si fprintf('\n') ; end + if pfm.use_bytecode + [flag, endo_simul]=bytecode(Y, exo_simul, pfm.params); + return; + end + z = Y(find(pfm.lead_lag_incidence')); [d1,jacobian] = model_dynamic(z,exo_simul,pfm.params,pfm.steady_state,2); diff --git a/mex/sources/bytecode/Evaluate.cc b/mex/sources/bytecode/Evaluate.cc index 76aec3b1b..60863cfef 100644 --- a/mex/sources/bytecode/Evaluate.cc +++ b/mex/sources/bytecode/Evaluate.cc @@ -322,7 +322,10 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int var = ((FLDV_ *) it_code->second)->get_pos(); lag = ((FLDV_ *) it_code->second)->get_lead_lag(); #ifdef DEBUG - mexPrintf("FLDV y[var=%d, lag=%d, it_=%d], y_size=%d evaluate=%d\n", var, lag, it_, y_size, evaluate); + if (evaluate) + mexPrintf("FLDV y[var=%d, lag=%d, it_=%d], y_size=%d evaluate=%d, ya[%d]=%f\n", var, lag, it_, y_size, evaluate, (it_+lag)*y_size+var, ya[(it_+lag)*y_size+var]); + else + mexPrintf("FLDV y[var=%d, lag=%d, it_=%d], y_size=%d evaluate=%d, y[%d]=%f\n", var, lag, it_, y_size, evaluate, (it_+lag)*y_size+var, y[(it_+lag)*y_size+var]); #endif if (evaluate) Stack.push(ya[(it_+lag)*y_size+var]); @@ -336,8 +339,8 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int var = ((FLDV_ *) it_code->second)->get_pos(); lag = ((FLDV_ *) it_code->second)->get_lead_lag(); #ifdef DEBUG - mexPrintf("FLDV x[var=%d, lag=%d, it_=%d], nb_row_x=%d evaluate=%d\n", var, lag, it_, nb_row_x, evaluate); - tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")"; + mexPrintf("FLDV x[var=%d, lag=%d, it_=%d], nb_row_x=%d evaluate=%d x[%d]=%f\n", var, lag, it_, nb_row_x, evaluate, it_+lag+var*nb_row_x, x[it_+lag+var*nb_row_x]); + //tmp_out << " x[" << it_+lag << ", " << var << "](" << x[it_+lag+var*nb_row_x] << ")"; #endif Stack.push(x[it_+lag+var*nb_row_x]); break; @@ -837,7 +840,7 @@ Evaluate::compute_block_time(const int Per_u_, const bool evaluate, /*const int case oLess: Stack.push(double (v1 < v2)); #ifdef DEBUG - tmp_out << " |" << v1 << "<" << v2 << "|"; + mexPrintf("v1=%f v2=%f v1 < v2 = %f\n",v1,v2,double(v1 < v2)); #endif break; case oGreater: diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc index fb47033c4..4b7af71f7 100644 --- a/mex/sources/bytecode/Interpreter.cc +++ b/mex/sources/bytecode/Interpreter.cc @@ -582,8 +582,8 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool evaluate #ifdef DEBUG mexPrintf("allocates jacobian_exo_block( %d, %d, mxREAL)\n", fb->get_size(), fb->get_exo_size()); #endif - jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_exo_size(), mxREAL)); - jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_det_exo_size(), mxREAL)); + jacobian_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_exo_jacob(), mxREAL)); + jacobian_det_exo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_det_exo_jacob(), mxREAL)); jacobian_other_endo_block.push_back(mxCreateDoubleMatrix(fb->get_size(), fb->get_nb_col_other_endo_jacob(), mxREAL)); } if (block >= 0) @@ -598,8 +598,8 @@ Interpreter::compute_blocks(string file_name, string bin_basename, bool evaluate else { #ifdef DEBUG - mexPrintf("endo in block=%d, type=%d, steady_state=%d, print_it=%d, Block_Count=%d, fb->get_is_linear()=%d, fb->get_endo_nbr()=%d, fb->get_Max_Lag()=%d, fb->get_Max_Lead()=%d, fb->get_u_count_int()=%d\n", - fb->get_size(), fb->get_type(), steady_state, print_it, Block_Count, fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); + mexPrintf("endo in Block_Count=%d, block=%d, type=%d, steady_state=%d, print_it=%d, Block_Count=%d, fb->get_is_linear()=%d, fb->get_endo_nbr()=%d, fb->get_Max_Lag()=%d, fb->get_Max_Lead()=%d, fb->get_u_count_int()=%d\n", + Block_Count, fb->get_size(), fb->get_type(), steady_state, print_it, Block_Count, fb->get_is_linear(), fb->get_endo_nbr(), fb->get_Max_Lag(), fb->get_Max_Lead(), fb->get_u_count_int()); #endif result = simulate_a_block(); if (result == ERROR_ON_EXIT) diff --git a/preprocessor/CodeInterpreter.hh b/preprocessor/CodeInterpreter.hh index b95908f52..b6068394a 100644 --- a/preprocessor/CodeInterpreter.hh +++ b/preprocessor/CodeInterpreter.hh @@ -20,6 +20,8 @@ #ifndef _CODEINTERPRETER_HH #define _CODEINTERPRETER_HH //#define DEBUGL +#include + #include #include #include @@ -1414,7 +1416,7 @@ private: int u_count_int; int nb_col_jacob; unsigned int det_exo_size, exo_size, other_endo_size; - unsigned int nb_col_other_endo_jacob; + unsigned int nb_col_det_exo_jacob, nb_col_exo_jacob, nb_col_other_endo_jacob; public: inline FBEGINBLOCK_() @@ -1426,7 +1428,7 @@ public: 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, int nb_col_jacob_arg, - unsigned int det_exo_size_arg, unsigned int exo_size_arg, unsigned int other_endo_size_arg, unsigned int nb_col_other_endo_jacob_arg, + unsigned int det_exo_size_arg, unsigned int nb_col_det_exo_jacob_arg, unsigned int exo_size_arg, unsigned int nb_col_exo_jacob_arg, unsigned int other_endo_size_arg, unsigned int nb_col_other_endo_jacob_arg, const vector &det_exogenous_arg, const vector &exogenous_arg, const vector &other_endogenous_arg) { op_code = FBEGINBLOCK; size = size_arg; type = type_arg; @@ -1436,8 +1438,10 @@ public: exogenous = vector(exogenous_arg); other_endogenous = vector(other_endogenous_arg); is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg; - nb_col_jacob = nb_col_jacob_arg; det_exo_size = det_exo_size_arg; exo_size = exo_size_arg; other_endo_size = other_endo_size_arg; - nb_col_other_endo_jacob = nb_col_other_endo_jacob_arg; + nb_col_jacob = nb_col_jacob_arg; + det_exo_size = det_exo_size_arg; nb_col_det_exo_jacob = nb_col_det_exo_jacob_arg; + exo_size = exo_size_arg; nb_col_exo_jacob = nb_col_exo_jacob_arg; + other_endo_size = other_endo_size_arg; nb_col_other_endo_jacob = nb_col_other_endo_jacob_arg; }; inline FBEGINBLOCK_(unsigned int size_arg, BlockSimulationType type_arg, int unsigned first_element, int unsigned block_size, @@ -1450,7 +1454,7 @@ public: is_linear = is_linear_arg; endo_nbr = endo_nbr_arg; Max_Lag = Max_Lag_arg; Max_Lead = Max_Lead_arg; u_count_int = u_count_int_arg; nb_col_jacob = nb_col_jacob_arg; det_exo_size = 0; exo_size = 0; other_endo_size = 0; - nb_col_other_endo_jacob = 0; + nb_col_det_exo_jacob = 0;nb_col_exo_jacob = 0;nb_col_other_endo_jacob = 0; } inline unsigned int get_size() @@ -1503,11 +1507,21 @@ public: return exo_size; }; inline unsigned int + get_nb_col_exo_jacob() + { + return nb_col_exo_jacob; + }; + inline unsigned int get_det_exo_size() { return det_exo_size; }; inline unsigned int + get_nb_col_det_exo_jacob() + { + return nb_col_det_exo_jacob; + }; + inline unsigned int get_other_endo_size() { return other_endo_size; @@ -1539,7 +1553,9 @@ public: } CompileCode.write(reinterpret_cast(&nb_col_jacob), sizeof(nb_col_jacob)); CompileCode.write(reinterpret_cast(&det_exo_size), sizeof(det_exo_size)); + CompileCode.write(reinterpret_cast(&nb_col_det_exo_jacob), sizeof(nb_col_det_exo_jacob)); CompileCode.write(reinterpret_cast(&exo_size), sizeof(exo_size)); + CompileCode.write(reinterpret_cast(&nb_col_exo_jacob), sizeof(nb_col_exo_jacob)); CompileCode.write(reinterpret_cast(&other_endo_size), sizeof(other_endo_size)); CompileCode.write(reinterpret_cast(&nb_col_other_endo_jacob), sizeof(nb_col_other_endo_jacob)); @@ -1577,7 +1593,9 @@ public: } memcpy(&nb_col_jacob, code, sizeof(nb_col_jacob)); code += sizeof(nb_col_jacob); memcpy(&det_exo_size, code, sizeof(det_exo_size)); code += sizeof(det_exo_size); + memcpy(&nb_col_det_exo_jacob, code, sizeof(nb_col_det_exo_jacob)); code += sizeof(nb_col_det_exo_jacob); memcpy(&exo_size, code, sizeof(exo_size)); code += sizeof(exo_size); + memcpy(&nb_col_exo_jacob, code, sizeof(nb_col_exo_jacob)); code += sizeof(nb_col_exo_jacob); memcpy(&other_endo_size, code, sizeof(other_endo_size)); code += sizeof(other_endo_size); memcpy(&nb_col_other_endo_jacob, code, sizeof(nb_col_other_endo_jacob)); code += sizeof(nb_col_other_endo_jacob); diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 0043b76d7..09390ca96 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -868,6 +868,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen prev_lag = -999999999; int prev_type = -1; int count_col_exo = 0; + int count_col_det_exo = 0; for (map, pair >, expr_t>::const_iterator it = first_derivatives_reordered_exo.begin(); it != first_derivatives_reordered_exo.end(); it++) @@ -880,7 +881,10 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen prev_var = var; prev_lag = lag; prev_type = type; - count_col_exo++; + if (type == eExogenous) + count_col_exo++; + else if (type == eExogenousDet) + count_col_det_exo++; } } @@ -892,13 +896,15 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen equation_reordered, false, symbol_table.endo_nbr(), - 0, - 0, + max_endo_lag, + max_endo_lead, u_count_int, count_col_endo, symbol_table.exo_det_nbr(), + count_col_det_exo, + symbol_table.exo_nbr(), count_col_exo, - other_endo_size, + 0, 0, exo_det, exo, @@ -1101,6 +1107,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin unsigned int block_exo_det_size = exo_det_block[block].size(); unsigned int block_other_endo_size = other_endo_block[block].size(); int block_max_lag = max_leadlag_block[block].first; + int block_max_lead = max_leadlag_block[block].second; if (simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE || simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE) @@ -1135,19 +1142,40 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin count_col_endo++; } } + unsigned int count_col_det_exo = 0; vector exo_det; for (lag_var_t::const_iterator it = exo_det_block[block].begin(); it != exo_det_block[block].end(); it++) - exo_det.push_back(it->first); + for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++) + { + count_col_det_exo++; + if (find (exo_det.begin(), exo_det.end(), *it1) == exo_det.end()) + exo_det.push_back(*it1); + } + + unsigned int count_col_exo = 0; vector exo; for (lag_var_t::const_iterator it = exo_block[block].begin(); it != exo_block[block].end(); it++) - exo.push_back(it->first); + for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++) + { + count_col_exo++; + if (find (exo.begin(), exo.end(), *it1) == exo.end()) + exo.push_back(*it1); + } + vector other_endo; unsigned int count_col_other_endo = 0; for (lag_var_t::const_iterator it = other_endo_block[block].begin(); it != other_endo_block[block].end(); it++) - { - other_endo.push_back(it->first); - count_col_other_endo += it->second.size(); - } + for (var_t::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++) + { + count_col_other_endo++; + if (find (exo_det.begin(), exo_det.end(), *it1) == exo_det.end()) + { + //other_endo.push_back(it->first); + //count_col_other_endo += it->second.size(); + other_endo.push_back(*it1); + } + } + FBEGINBLOCK_ fbeginblock(block_mfs, simulation_type, getBlockFirstEquation(block), @@ -1157,19 +1185,21 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin blocks_linear[block], symbol_table.endo_nbr(), block_max_lag, - block_max_lag, + block_max_lead, u_count_int, count_col_endo, - block_exo_det_size, + exo_det.size(), + count_col_det_exo, + exo.size(), getBlockExoColSize(block), - block_other_endo_size, + other_endo.size(), count_col_other_endo, exo_det, exo, other_endo ); fbeginblock.write(code_file, instruction_number); - + // The equations for (i = 0; i < (int) block_size; i++) { @@ -1412,7 +1442,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin } prev_var = -1; prev_lag = -999999999; - int count_col_exo = 0; + count_col_exo = 0; for (map >, expr_t>::const_iterator it = tmp_exo_derivative.begin(); it != tmp_exo_derivative.end(); it++) { int lag = it->first.first;