diff --git a/matlab/perfect-foresight-models/det_cond_forecast.m b/matlab/perfect-foresight-models/det_cond_forecast.m
index 2bef6c5a6b480a9bf812a98edcd681f884c38f3b..38d72d34594a9b29350c0cd42c305e7645c15d60 100644
--- a/matlab/perfect-foresight-models/det_cond_forecast.m
+++ b/matlab/perfect-foresight-models/det_cond_forecast.m
@@ -105,15 +105,16 @@ else
sym_dset = dset(dates(-range(1)):dates(range(range.ndat)));
periods = options_.periods + M_.maximum_lag + M_.maximum_lead;
+ total_periods = periods + range.ndat;
if isfield(oo_, 'exo_simul')
- if size(oo_.exo_simul, 1) ~= max(range.ndat + 1, periods)
- oo_.exo_simul = repmat(oo_.exo_steady_state',max(range.ndat + 1, periods),1);
+ if size(oo_.exo_simul, 1) ~= total_periods
+ oo_.exo_simul = repmat(oo_.exo_steady_state',total_periods,1);
end
else
- oo_.exo_simul = repmat(oo_.exo_steady_state',max(range.ndat + 1, periods),1);
+ oo_.exo_simul = repmat(oo_.exo_steady_state',total_periods,1);
end
- oo_.endo_simul = repmat(oo_.steady_state, 1, max(range.ndat + 1, periods));
+ oo_.endo_simul = repmat(oo_.steady_state, 1, total_periods);
for i = 1:sym_dset.vobs
iy = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.endo_names)));
@@ -168,7 +169,7 @@ else
oo_.exo_simul = exo;
end
endo = endo';
- endo_l = size(endo(1+M_.maximum_lag:end,:),1);
+ endo_l = size(endo(1+M_.maximum_lag:end,:),1);
jrng = dates(plan.date(1)):dates(plan.date(1)+endo_l);
data_set = dseries(nan(endo_l, dset.vobs), plan.date(1), dset.name);
for i = 1:length(dset.name)
diff --git a/mex/sources/bytecode/ErrorHandling.hh b/mex/sources/bytecode/ErrorHandling.hh
index 39fc01ec3acdebc2b2a286b0d0d765132d141dac..eb2c906f840efc30d776a775776d62ea218aa111 100644
--- a/mex/sources/bytecode/ErrorHandling.hh
+++ b/mex/sources/bytecode/ErrorHandling.hh
@@ -316,7 +316,7 @@ public:
double *y, *ya;
int y_size;
double *T;
- int nb_row_xd, nb_row_x, col_x;
+ int nb_row_xd, nb_row_x, col_x, col_y;
int y_kmin, y_kmax, periods;
double *x, *params;
double *u;
diff --git a/mex/sources/bytecode/Interpreter.cc b/mex/sources/bytecode/Interpreter.cc
index 4b2dccb9be93a5f7db4712c7cf52b0101567e782..2fbec06b9a441b3012176fd8ef2861d16bae3c77 100644
--- a/mex/sources/bytecode/Interpreter.cc
+++ b/mex/sources/bytecode/Interpreter.cc
@@ -35,7 +35,7 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
int maxit_arg_, double solve_tolf_arg, size_t size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
- bool steady_state_arg, bool print_it_arg, int col_x_arg
+ bool steady_state_arg, bool print_it_arg, int col_x_arg, int col_y_arg
#ifdef CUDA
, const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
#endif
@@ -73,7 +73,8 @@ Interpreter::Interpreter(double *params_arg, double *y_arg, double *ya_arg, doub
solve_algo = solve_algo_arg;
global_temporary_terms = global_temporary_terms_arg;
print = print_arg;
- col_x = col_x_arg;
+ col_x = col_x_arg;
+ col_y = col_y_arg;
GlobalTemporaryTerms = GlobalTemporaryTerms_arg;
print_error = print_error_arg;
//steady_state = steady_state_arg;
@@ -854,44 +855,33 @@ Interpreter::elastic(string str, unsigned int len, bool left)
bool
Interpreter::extended_path(string file_name, string bin_basename, bool evaluate, int block, int &nb_blocks, int nb_periods, vector<s_plan> sextended_path, vector<s_plan> sconstrained_extended_path, vector<string> dates, table_conditional_global_type table_conditional_global)
{
- CodeLoad code;
+ CodeLoad code;
+
ReadCodeFile(file_name, code);
it_code = code_liste.begin();
it_code_type Init_Code = code_liste.begin();
/*size_t size_of_direction = y_size*(periods + y_kmax + y_kmin)*sizeof(double);
double *y_save = (double *) mxMalloc(size_of_direction);
double *x_save = (double *) mxMalloc((periods + y_kmax + y_kmin) * col_x *sizeof(double));*/
- int max_periods = max(periods, nb_periods);
- size_t size_of_direction = y_size*(/*max_periods*/ nb_periods + y_kmax + y_kmin)*sizeof(double);
- double *y_save = (double *) mxMalloc(size_of_direction);
- test_mxMalloc(y_save, __LINE__, __FILE__, __func__, size_of_direction);
-
-
- double *x_save = (double *) mxMalloc((/*max_periods*/nb_periods + y_kmax + y_kmin) * col_x *sizeof(double));
- test_mxMalloc(x_save, __LINE__, __FILE__, __func__, (/*max_periods*/nb_periods + y_kmax + y_kmin) * col_x *sizeof(double));
+ size_t size_of_direction = y_size*col_y*sizeof(double);
+ double *y_save = (double *) mxMalloc(size_of_direction);
+ test_mxMalloc(y_save, __LINE__, __FILE__, __func__, size_of_direction);
+ double *x_save = (double *) mxMalloc(nb_row_x * col_x *sizeof(double));
+ test_mxMalloc(x_save, __LINE__, __FILE__, __func__, nb_row_x * col_x *sizeof(double));
vector_table_conditional_local_type vector_table_conditional_local;
vector_table_conditional_local.clear();
- int endo_name_length_l = endo_name_length;
+ int endo_name_length_l = endo_name_length;
for (int j = 0; j < col_x* nb_row_x; j++)
{
x_save[j] = x[j];
x[j] = 0;
- }
- for (int j = 0; j < col_x; j++)
- for (int i = nb_row_x; i < nb_periods; i++)
- x_save[j * (nb_periods + y_kmin + y_kmax) + i] = 0;//x[(j * nb_row_x) + 1];
+ }
for (int j = 0; j < col_x; j++)
- x[y_kmin + j * nb_row_x] = x_save[1 + y_kmin + j * nb_row_x];
- //for (int i = 0; i < (y_size*(periods + y_kmax + y_kmin)); i++)
- for (int j = 0; j < y_size; j++)
- for (int i = 0; i < nb_periods; i++)
- if (i < periods)
- y_save[j + (i + y_kmin) * y_size] = y[ j + (i + y_kmin) * y_size];
- else
- y_save[j + (i + y_kmin) * y_size] = y[ j + (periods + y_kmin) * y_size];
-
+ x[y_kmin + j * nb_row_x] = x_save[y_kmin + j * nb_row_x];
+ for (int i = 0; i < y_size * col_y; i++)
+ y_save[i] = y[i];
if (endo_name_length_l < 8)
endo_name_length_l = 8;
bool old_print_it = print_it;
@@ -922,9 +912,9 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
mexPrintf("|%s|",elastic(dates[t], date_length+2, false).c_str());
mexEvalString("drawnow;");
}
- for (vector<s_plan>::const_iterator it = sextended_path.begin(); it != sextended_path.end(); it++)
- {
- x[y_kmin + (it->exo_num - 1) * /*(nb_periods + y_kmax + y_kmin)*/ nb_row_x] = it->value[t];
+ for (vector<s_plan>::const_iterator it = sextended_path.begin(); it != sextended_path.end(); it++)
+ {
+ x[y_kmin + (it->exo_num - 1) * /*(nb_periods + y_kmax + y_kmin)*/ nb_row_x] = it->value[t];
}
it_code = Init_Code;
@@ -965,10 +955,10 @@ Interpreter::extended_path(string file_name, string bin_basename, bool evaluate,
{
for(int k = nb_periods; k < periods; k++)
y_save[j + (k + y_kmin) * y_size] = y[ j + ( k - (nb_periods-1) + y_kmin) * y_size];
- }*/
- for (int i = 0; i < (y_size*(nb_periods + y_kmax + y_kmin)); i++)
+ }*/
+ for (int i = 0; i < y_size * col_y; i++)
y[i] = y_save[i];
- for (int j = 0; j < col_x* nb_row_x; j++)
+ for (int j = 0; j < col_x * nb_row_x; j++)
x[j] = x_save[j];
if (Init_Code->second)
mxFree(Init_Code->second);
diff --git a/mex/sources/bytecode/Interpreter.hh b/mex/sources/bytecode/Interpreter.hh
index 8e8a00e0a8c4301a86104628b343658bddf72e61..1822c1c95b3a65af7ee30cf63d7110cd530da844 100644
--- a/mex/sources/bytecode/Interpreter.hh
+++ b/mex/sources/bytecode/Interpreter.hh
@@ -59,7 +59,7 @@ public:
int maxit_arg_, double solve_tolf_arg, size_t size_of_direction_arg, double slowc_arg, int y_decal_arg, double markowitz_c_arg,
string &filename_arg, int minimal_solving_periods_arg, int stack_solve_algo_arg, int solve_algo_arg,
bool global_temporary_terms_arg, bool print_arg, bool print_error_arg, mxArray *GlobalTemporaryTerms_arg,
- bool steady_state_arg, bool print_it_arg, int col_x_arg
+ bool steady_state_arg, bool print_it_arg, int col_x_arg, int col_y_arg
#ifdef CUDA
, const int CUDA_device, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
#endif
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index e999dd025a8d73ff83dbf2114d43bfb038a903f7..a0752b7fb7cbbe4eb775213086b64e937373f759 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -3409,10 +3409,8 @@ dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, do
int eq = index_vara[i+Size*(y_kmin)];
int flip_exo = vector_table_conditional_local[i].var_exo;
double yy = -(res[i] + x[y_kmin + flip_exo*nb_row_x]);
- direction[eq] = 0;
- //mexPrintf("solving for init x[%d]=%f\n", flip_exo*nb_row_x + y_kmin, x[flip_exo*nb_row_x + y_kmin]);
- x[flip_exo*nb_row_x + y_kmin] += slowc_l * yy;
- //mexPrintf(" ==> x[flip_exo=%d (%d)]=%f yy=%f, slowc_l=%f, res[%d]=%f, y[5]=%f\n", flip_exo, flip_exo*nb_row_x + y_kmin, x[flip_exo*nb_row_x + y_kmin], yy, slowc_l, i, res[i], y[5]);
+ direction[eq] = 0;
+ x[flip_exo*nb_row_x + y_kmin] += slowc_l * yy;
}
else
{
diff --git a/mex/sources/bytecode/bytecode.cc b/mex/sources/bytecode/bytecode.cc
index c67bc09207dd697160c0ad3d221f05340d53cf59..27a2d850e2bd79eaa2905fb38113608f2e8fc156 100644
--- a/mex/sources/bytecode/bytecode.cc
+++ b/mex/sources/bytecode/bytecode.cc
@@ -660,7 +660,7 @@ main(int nrhs, const char *prhs[])
conditional_local.constrained_value = specific_constrained_paths_[j];
table_conditional_global[constrained_int_date[j]][sconditional_extended_path[i].exo_num - 1] = conditional_local;
sconditional_extended_path[i].per_value[j] = make_pair(constrained_int_date[j], specific_constrained_paths_[j]);
- sconditional_extended_path[i].value[constrained_int_date[j]] = specific_constrained_paths_[j];
+ sconditional_extended_path[i].value[constrained_int_date[j]] = specific_constrained_paths_[j];
if (max_periods < constrained_int_date[j] + 1)
max_periods = constrained_int_date[j] + 1;
}
@@ -1060,7 +1060,6 @@ main(int nrhs, const char *prhs[])
if (stack_solve_algo == 7 && !steady_state)
DYN_MEX_FUNC_ERR_MSG_TXT("bytecode has not been compiled with CUDA option. Bytecode Can't use options_.stack_solve_algo=7\n");
#endif
-
size_t size_of_direction = col_y*row_y*sizeof(double);
double *y = (double *) mxMalloc(size_of_direction);
error_msg.test_mxMalloc(y, __LINE__, __FILE__, __func__, size_of_direction);
@@ -1069,7 +1068,7 @@ main(int nrhs, const char *prhs[])
direction = (double *) mxMalloc(size_of_direction);
error_msg.test_mxMalloc(direction, __LINE__, __FILE__, __func__, size_of_direction);
memset(direction, 0, size_of_direction);
- /*mexPrintf("col_x : %d, row_x : %d\n",col_x, row_x);*/
+ /*mexPrintf("col_x : %d, row_x : %d\n",col_x, row_x);*/
double *x = (double *) mxMalloc(col_x*row_x*sizeof(double));
error_msg.test_mxMalloc(x, __LINE__, __FILE__, __func__, col_x*row_x*sizeof(double));
for (i = 0; i < row_x*col_x; i++)
@@ -1084,32 +1083,10 @@ main(int nrhs, const char *prhs[])
size_t y_size = row_y;
size_t nb_row_x = row_x;
- if (extended_path)
- {
- if (max_periods > periods)
- {
- double* y_save = (double*) mxMalloc(y_size * (max_periods + y_kmin + y_kmax) * sizeof(double));
- for (int j = 0; j < y_size; j++)
- for (int i = 0; i < periods; i++)
- y_save[j + (i + y_kmin) * y_size] = y[ j + (i + y_kmin) * y_size];
-
- y = (double*) mxRealloc(y, y_size * (max_periods + y_kmin + y_kmax) * sizeof(double));
- x = (double*) mxRealloc(x, nb_row_x * (max_periods + y_kmin + y_kmax) * sizeof(double));
- for (int j = 0; j < nb_row_x * (max_periods + y_kmin + y_kmax); j++)
- x[j] = 0;
- for (int j = 0; j < y_size; j++)
- for (int i = 0; i < max_periods; i++)
- if (i < periods)
- y[j + (i + y_kmin) * y_size] = y_save[ j + (i + y_kmin) * y_size];
- else
- y[j + (i + y_kmin) * y_size] = y_save[ j + (periods + y_kmin) * y_size];
- mxFree(y_save);
- }
- }
clock_t t0 = clock();
Interpreter interprete(params, y, ya, x, steady_yd, steady_xd, direction, y_size, nb_row_x, nb_row_xd, periods, y_kmin, y_kmax, maxit_, solve_tolf, size_of_direction, slowc, y_decal,
markowitz_c, file_name, minimal_solving_periods, stack_solve_algo, solve_algo, global_temporary_terms, print, print_error, GlobalTemporaryTerms, steady_state,
- print_it, col_x
+ print_it, col_x, col_y
#ifdef CUDA
, CUDA_device, cublas_handle, cusparse_handle, descr
#endif
@@ -1174,29 +1151,26 @@ main(int nrhs, const char *prhs[])
pind[i] = residual[i];
}
else
- {
- int out_periods;
- if (extended_path)
- out_periods = max_periods + y_kmin + y_kmax;
- else
+ {
+ int out_periods;
+ if (extended_path)
+ out_periods = max_periods + y_kmin;
+ else
out_periods = row_y;
- plhs[1] = mxCreateDoubleMatrix(/*int(row_y)*/out_periods, int(col_y), mxREAL);
+ plhs[1] = mxCreateDoubleMatrix(out_periods, int(col_y), mxREAL);
pind = mxGetPr(plhs[1]);
- for (i = 0; i < /*row_y*/out_periods*col_y; i++)
- {
- pind[i] = y[i];
- }
-
+ for (i = 0; i < out_periods*col_y; i++)
+ pind[i] = y[i];
}
}
else
{
- int out_periods;
- if (extended_path)
- out_periods = max_periods;
- else
+ int out_periods;
+ if (extended_path)
+ out_periods = max_periods + y_kmin;
+ else
out_periods = col_y;
- plhs[1] = mxCreateDoubleMatrix(int(row_y), /*int(col_y)*/out_periods, mxREAL);
+ plhs[1] = mxCreateDoubleMatrix(int(row_y), out_periods, mxREAL);
pind = mxGetPr(plhs[1]);
if (evaluate)
{
@@ -1205,10 +1179,8 @@ main(int nrhs, const char *prhs[])
pind[i] = residual[i];
}
else
- for (i = 0; i < row_y*out_periods /** col_y*/; i++)
- {
+ for (i = 0; i < row_y*out_periods; i++)
pind[i] = y[i];
- }
}
if (nlhs > 2)
{