From 5b591fac422448ed30797b4c828f7f90aa1d560e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Mon, 24 Jun 2019 17:52:09 +0200
Subject: [PATCH] New perfect_foresight_problem MEX file
It constructs the stacked residuals and jacobian of the perfect foresight
problem.
It is an almost perfect replacement for the perfect_foresight_problem.m
routine, while being much more efficient.
Note however that the DLL never return complex numbers (it instead puts NaNs at
the place where there would have been complex). This may create problems for
some MOD files; the algorithms will need to be adapted to use a more
line-search method.
---
.../perfect_foresight_problem.m | 102 -------
.../perfect_foresight_solver.m | 33 +--
.../perfect_foresight_solver_core.m | 34 +--
.../private/initialize_stacked_problem.m | 14 +-
matlab/perfect-foresight-models/sim1.m | 220 +++++---------
.../solve_stacked_problem.m | 9 +-
mex/build/matlab/Makefile.am | 2 +-
mex/build/matlab/configure.ac | 3 +-
.../perfect_foresight_problem/Makefile.am | 2 +
mex/build/octave/Makefile.am | 2 +-
mex/build/octave/configure.ac | 5 +-
.../perfect_foresight_problem/Makefile.am | 3 +
mex/build/perfect_foresight_problem.am | 15 +
mex/sources/Makefile.am | 3 +-
.../DynamicModelCaller.cc | 215 ++++++++++++++
.../DynamicModelCaller.hh | 83 ++++++
.../perfect_foresight_problem.cc | 277 ++++++++++++++++++
17 files changed, 725 insertions(+), 297 deletions(-)
delete mode 100644 matlab/perfect-foresight-models/perfect_foresight_problem.m
create mode 100644 mex/build/matlab/perfect_foresight_problem/Makefile.am
create mode 100644 mex/build/octave/perfect_foresight_problem/Makefile.am
create mode 100644 mex/build/perfect_foresight_problem.am
create mode 100644 mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
create mode 100644 mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
create mode 100644 mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
diff --git a/matlab/perfect-foresight-models/perfect_foresight_problem.m b/matlab/perfect-foresight-models/perfect_foresight_problem.m
deleted file mode 100644
index 07220cdcef..0000000000
--- a/matlab/perfect-foresight-models/perfect_foresight_problem.m
+++ /dev/null
@@ -1,102 +0,0 @@
-function [residuals,JJacobian] = perfect_foresight_problem(y, dynamic_function, Y0, YT, ...
- exo_simul, params, steady_state, ...
- maximum_lag, T, ny, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, ...
- i_cols_j, i_cols_0, i_cols_J0)
-
-% Computes the residuals and the Jacobian matrix for a perfect foresight problem over T periods.
-%
-% INPUTS
-% - y [double] N*1 array, terminal conditions for the endogenous variables
-% - dynamic_function [handle] function handle to _dynamic-file
-% - Y0 [double] N*1 array, initial conditions for the endogenous variables
-% - YT [double] N*1 array, terminal conditions for the endogenous variables
-% - exo_simul [double] nperiods*M_.exo_nbr matrix of exogenous variables (in declaration order)
-% for all simulation periods
-% - params [double] nparams*1 array, parameter values
-% - steady_state [double] endo_nbr*1 vector of steady state values
-% - maximum_lag [scalar] maximum lag present in the model
-% - T [scalar] number of simulation periods
-% - ny [scalar] number of endogenous variables
-% - i_cols [double] indices of variables appearing in M.lead_lag_incidence
-% and that need to be passed to _dynamic-file
-% - i_cols_J1 [double] indices of contemporaneous and forward looking variables
-% appearing in M.lead_lag_incidence
-% - i_cols_1 [double] indices of contemporaneous and forward looking variables in
-% M.lead_lag_incidence in dynamic Jacobian (relevant in first period)
-% - i_cols_T [double] columns of dynamic Jacobian related to contemporaneous and backward-looking
-% variables (relevant in last period)
-% - i_cols_j [double] indices of contemporaneous variables in M.lead_lag_incidence
-% in dynamic Jacobian (relevant in intermediate periods)
-% - i_cols_0 [double] indices of contemporaneous variables in M.lead_lag_incidence in dynamic
-% Jacobian (relevant in problems with periods=1)
-% - i_cols_J0 [double] indices of contemporaneous variables appearing in M.lead_lag_incidence (relevant in problems with periods=1)
-%
-% OUTPUTS
-% - residuals [double] (N*T)*1 array, residuals of the stacked problem
-% - JJacobian [double] (N*T)*(N*T) array, Jacobian of the stacked problem
-%
-% ALGORITHM
-% None
-%
-% SPECIAL REQUIREMENTS
-% None.
-
-% Copyright (C) 1996-2019 Dynare Team
-%
-% This file is part of Dynare.
-%
-% Dynare is free software: you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation, either version 3 of the License, or
-% (at your option) any later version.
-%
-% Dynare is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License
-% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
-
-
-YY = [Y0; y; YT];
-
-residuals = zeros(T*ny,1);
-if nargout == 2
- iJacobian = cell(T,1);
-end
-
-i_rows = 1:ny;
-i_cols_J = i_cols;
-offset = 0;
-
-for it = maximum_lag+(1:T)
- if nargout == 1
- residuals(i_rows) = dynamic_function(YY(i_cols), exo_simul, params, steady_state, it);
- elseif nargout == 2
- [residuals(i_rows),jacobian] = dynamic_function(YY(i_cols), exo_simul, params, steady_state, it);
- if T==1 && it==maximum_lag+1
- [rows, cols, vals] = find(jacobian(:,i_cols_0));
- iJacobian{1} = [rows, i_cols_J0(cols), vals];
- elseif it == maximum_lag+1
- [rows,cols,vals] = find(jacobian(:,i_cols_1));
- iJacobian{1} = [offset+rows, i_cols_J1(cols), vals];
- elseif it == maximum_lag+T
- [rows,cols,vals] = find(jacobian(:,i_cols_T));
- iJacobian{T} = [offset+rows, i_cols_J(i_cols_T(cols)), vals];
- else
- [rows,cols,vals] = find(jacobian(:,i_cols_j));
- iJacobian{it-maximum_lag} = [offset+rows, i_cols_J(cols), vals];
- i_cols_J = i_cols_J + ny;
- end
- offset = offset + ny;
- end
- i_rows = i_rows + ny;
- i_cols = i_cols + ny;
-end
-
-if nargout == 2
- iJacobian = cat(1,iJacobian{:});
- JJacobian = sparse(iJacobian(:,1), iJacobian(:,2), iJacobian(:,3), T*ny, T*ny);
-end
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m
index 9816d65967..a67ad34a2b 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m
@@ -176,26 +176,23 @@ if ~oo_.deterministic_simulation.status && ~options_.no_homotopy
end
-if ~isreal(oo_.endo_simul(:)) % can only happen without bytecode
- y0 = real(oo_.endo_simul(:,1));
- yT = real(oo_.endo_simul(:,periods+2));
- yy = real(oo_.endo_simul(:,2:periods+1));
- illi = M_.lead_lag_incidence';
- [i_cols,~,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,~,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
- if periods==1
- i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
- i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
+if ~isreal(oo_.endo_simul(:)) % cannot happen with bytecode or the perfect_foresight_problem DLL
+ if M_.maximum_lag > 0
+ y0 = real(oo_.endo_simul(:, M_.maximum_lag));
else
- i_cols_0 = [];
- i_cols_J0 = [];
+ y0 = NaN(ny, 1);
end
- residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
- oo_.exo_simul,M_.params,oo_.steady_state, ...
- M_.maximum_lag, periods, M_.endo_nbr, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
+ if M_.maximum_lead > 0
+ yT = real(oo_.endo_simul(:, M_.maximum_lag+periods+1));
+ else
+ yT = NaN(ny, 1);
+ end
+ yy = real(oo_.endo_simul(:,M_.maximum_lag+(1:periods)));
+ model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']);
+ [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
+
+ residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll);
+
if max(abs(residuals))< options_.dynatol.f
oo_.deterministic_simulation.status = 1;
oo_.endo_simul=real(oo_.endo_simul);
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index 24464e8f19..76a6e2d2f5 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -120,31 +120,27 @@ else
end
if nargout>1
- y0 = oo_.endo_simul(:,1);
- yT = oo_.endo_simul(:,periods+2);
- yy = oo_.endo_simul(:,2:periods+1);
- illi = M_.lead_lag_incidence';
- [i_cols,~,i_cols_j] = find(illi(:));
- illi = illi(:,2:3);
- [i_cols_J1,~,i_cols_1] = find(illi(:));
- i_cols_T = nonzeros(M_.lead_lag_incidence(1:2,:)');
- if periods==1
- i_cols_0 = nonzeros(M_.lead_lag_incidence(2,:)');
- i_cols_J0 = find(M_.lead_lag_incidence(2,:)');
- else
- i_cols_0 = [];
- i_cols_J0 = [];
- end
if options_.block && ~options_.bytecode
maxerror = oo_.deterministic_simulation.error;
else
if options_.bytecode
[chck, residuals, ~]= bytecode('dynamic','evaluate', oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, 1);
else
- residuals = perfect_foresight_problem(yy(:),str2func([M_.fname '.dynamic']), y0, yT, ...
- oo_.exo_simul,M_.params,oo_.steady_state, ...
- M_.maximum_lag, periods,M_.endo_nbr,i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
+ if M_.maximum_lag > 0
+ y0 = oo_.endo_simul(:, M_.maximum_lag);
+ else
+ y0 = NaN(ny, 1);
+ end
+ if M_.maximum_lead > 0
+ yT = oo_.endo_simul(:, M_.maximum_lag+periods+1);
+ else
+ yT = NaN(ny, 1);
+ end
+ yy = oo_.endo_simul(:,M_.maximum_lag+(1:periods));
+ model_dynamic_g1_nz = str2func([M_.fname,'.dynamic_g1_nz']);
+ [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
+
+ residuals = perfect_foresight_problem(yy(:), M_.fname, sum(M_.dynamic_tmp_nbr(1:2)), y0, yT, oo_.exo_simul, M_.params, oo_.steady_state, periods, M_.endo_nbr, M_.maximum_lag, M_.maximum_endo_lag, M_.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M_.has_external_function, options_.use_dll);
end
maxerror = max(max(abs(residuals)));
end
diff --git a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
index 55290fc9f9..5677add41d 100644
--- a/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
+++ b/matlab/perfect-foresight-models/private/initialize_stacked_problem.m
@@ -70,8 +70,16 @@ elseif (options.solve_algo == 11)
end
end
-y0 = endogenousvariables(:,M.maximum_lag);
-yT = endogenousvariables(:,M.maximum_lag+periods+1);
+if M.maximum_lag > 0
+ y0 = endogenousvariables(:, M.maximum_lag);
+else
+ y0 = NaN(ny, 1);
+end
+if M.maximum_lead > 0
+ yT = endogenousvariables(:, M.maximum_lag+periods+1);
+else
+ yT = NaN(ny, 1);
+end
z = endogenousvariables(:,M.maximum_lag+(1:periods));
illi = M.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
@@ -85,4 +93,4 @@ else
i_cols_0 = [];
i_cols_J0 = [];
end
-dynamicmodel = str2func([M.fname,'.dynamic']);
\ No newline at end of file
+dynamicmodel = str2func([M.fname,'.dynamic']);
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 1453875a78..1379ade92c 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -1,21 +1,17 @@
function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariables, steadystate, M, options)
-
-% Performs deterministic simulations with lead or lag on one period. Uses sparse matrices.
+% Performs deterministic simulations with lead or lag of one period, using
+% a basic Newton solver on sparse matrices.
+% Uses perfect_foresight_problem DLL to construct the stacked problem.
%
% INPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
+% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
% - exogenousvariables [double] T*M array, paths for the exogenous variables.
% - steadystate [double] N*1 array, steady state for the endogenous variables.
% - M [struct] contains a description of the model.
% - options [struct] contains various options.
% OUTPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
+% - endogenousvariables [double] N*(T+M.maximum_lag+M.maximum_lead) array, paths for the endogenous variables (solution of the perfect foresight model).
% - info [struct] contains informations about the results.
-% ALGORITHM
-% ...
-%
-% SPECIAL REQUIREMENTS
-% None.
% Copyright (C) 1996-2019 Dynare Team
%
@@ -36,38 +32,25 @@ function [endogenousvariables, info] = sim1(endogenousvariables, exogenousvariab
verbose = options.verbosity;
-endogenous_terminal_period = options.endogenous_terminal_period;
-vperiods = options.periods*ones(1,options.simul.maxit);
-azero = options.dynatol.f/1e7;
-
-lead_lag_incidence = M.lead_lag_incidence;
-
ny = M.endo_nbr;
+periods = options.periods;
+vperiods = periods*ones(1,options.simul.maxit);
-maximum_lag = M.maximum_lag;
-max_lag = M.maximum_endo_lag;
-
-nyp = nnz(lead_lag_incidence(1,:)) ;
-ny0 = nnz(lead_lag_incidence(2,:)) ;
-nyf = nnz(lead_lag_incidence(3,:)) ;
-
-nd = nyp+ny0+nyf;
-stop = 0 ;
+if M.maximum_lag > 0
+ y0 = endogenousvariables(:, M.maximum_lag);
+else
+ y0 = NaN(ny, 1);
+end
-periods = options.periods;
-params = M.params;
+if M.maximum_lead > 0
+ yT = endogenousvariables(:, M.maximum_lag+periods+1);
+else
+ yT = NaN(ny, 1);
+end
-i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
-i_cols_A1 = find(lead_lag_incidence(2:3,:)');
-i_cols_A1 = i_cols_A1(:);
-i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
-i_cols_0 = nonzeros(lead_lag_incidence(2,:)');
-i_cols_A0 = find(lead_lag_incidence(2,:)');
-i_cols_A0 = i_cols_A0(:);
-i_cols_j = (1:nd)';
-i_upd = maximum_lag*ny+(1:periods*ny);
+y = reshape(endogenousvariables(:, M.maximum_lag+(1:periods)), ny*periods, 1);
-Y = endogenousvariables(:);
+stop = false;
if verbose
skipline()
@@ -76,121 +59,78 @@ if verbose
skipline()
end
-model_dynamic = str2func([M.fname,'.dynamic']);
-
-res = zeros(periods*ny,1);
-
-o_periods = periods;
-
-if endogenous_terminal_period
- ZERO = zeros(length(i_upd),1);
-end
+model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']);
+[nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
-h1 = clock ;
-iA = zeros(periods*M.NNZDerivatives(1),3);
+h1 = clock;
for iter = 1:options.simul.maxit
- h2 = clock ;
- i_rows = (1:ny)';
- i_cols_A = find(lead_lag_incidence');
- i_cols_A = i_cols_A(:);
- i_cols = i_cols_A+(maximum_lag-1)*ny;
- m = 0;
- for it = (maximum_lag+1):(maximum_lag+periods)
- [d1,jacobian] = model_dynamic(Y(i_cols), exogenousvariables, params, steadystate,it);
- if periods==1 && it==maximum_lag+1
- [r,c,v] = find(jacobian(:,i_cols_0));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A0(c(:)),v(:)];
- elseif it == maximum_lag+periods
- [r,c,v] = find(jacobian(:,i_cols_T));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A(i_cols_T(c(:))),v(:)];
- elseif it == maximum_lag+1
- [r,c,v] = find(jacobian(:,i_cols_1));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A1(c(:)),v(:)];
- else
- [r,c,v] = find(jacobian(:,i_cols_j));
- iA((1:length(v))+m,:) = [i_rows(r(:)),i_cols_A(c(:)),v(:)];
- end
- m = m + length(v);
- res(i_rows) = d1;
- if endogenous_terminal_period && iter>1
- dr = max(abs(d1));
- if dr<azero
- vperiods(iter) = it;
- periods = it-maximum_lag+1;
+ h2 = clock;
+
+ [res, A] = perfect_foresight_problem(y, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll);
+
+ if options.endogenous_terminal_period && iter > 1
+ for it = 1:periods
+ if max(abs(res((it-1)*ny+(1:ny)))) < options.dynatol.f/1e7
+ if it < periods
+ res = res(1:(it*ny));
+ A = A(1:(it*ny), 1:(it*ny));
+ yT = y(it*ny+(1:ny));
+ endogenousvariables(:, M.maximum_lag+((it+1):periods)) = reshape(y(it*ny+(1:(ny*(periods-it)))), ny, periods-it);
+ y = y(1:(it*ny));
+ periods = it;
+ end
break
end
end
- i_rows = i_rows + ny;
- i_cols = i_cols + ny;
- if it > maximum_lag+1
- i_cols_A = i_cols_A + ny;
- end
+ vperiods(iter) = periods;
end
+
err = max(abs(res));
if options.debug
fprintf('\nLargest absolute residual at iteration %d: %10.3f\n',iter,err);
- if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y))
+ if any(isnan(res)) || any(isinf(res)) || any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n');
end
- if ~isreal(res) || ~isreal(Y)
- fprintf('\nWARNING: Imaginary parts detected in the residuals or endogenous variables.\n');
- end
skipline()
end
if verbose
- str = sprintf('Iter: %s,\t err. = %s, \t time = %s',num2str(iter),num2str(err), num2str(etime(clock,h2)));
- disp(str);
+ fprintf('Iter: %d,\t err. = %g,\t time = %g\n', iter, err, etime(clock,h2));
end
if err < options.dynatol.f
- stop = 1 ;
+ stop = true;
break
end
- iA = iA(1:m,:);
- A = sparse(iA(:,1),iA(:,2),iA(:,3),periods*ny,periods*ny);
- if endogenous_terminal_period && iter>1
- dy = ZERO;
- if options.simul.robust_lin_solve
- dy(1:i_rows(end)) = -lin_solve_robust( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)),verbose );
- else
- dy(1:i_rows(end)) = -lin_solve( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)), verbose );
- end
+ if options.simul.robust_lin_solve
+ dy = -lin_solve_robust(A, res, verbose);
else
- if options.simul.robust_lin_solve
- dy = -lin_solve_robust( A, res, verbose );
- else
- dy = -lin_solve( A, res, verbose );
- end
+ dy = -lin_solve(A, res, verbose);
end
- if any(~isreal(dy)) || any(isnan(dy)) || any(isinf(dy))
+ if any(isnan(dy)) || any(isinf(dy))
if verbose
display_critical_variables(reshape(dy,[ny periods])', M);
end
end
- Y(i_upd) = Y(i_upd) + dy;
+ y = y + dy;
end
-if endogenous_terminal_period
- err = evaluate_max_dynamic_residual(model_dynamic, Y, exogenousvariables, params, steadystate, o_periods, ny, max_lag, lead_lag_incidence);
- periods = o_periods;
-end
+endogenousvariables(:, M.maximum_lag+(1:periods)) = reshape(y, ny, periods);
+if options.endogenous_terminal_period
+ periods = options.periods;
+ err = evaluate_max_dynamic_residual(str2func([M.fname,'.dynamic']), endogenousvariables, exogenousvariables, M.params, steadystate, periods, ny, M.maximum_endo_lag, M.lead_lag_incidence);
+end
if stop
- if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) || ~isreal(res) || ~isreal(Y)
+ if any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
info.status = false;% NaN or Inf occurred
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
- endogenousvariables = reshape(Y,ny,periods+maximum_lag+M.maximum_lead);
if verbose
skipline()
- disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
- if ~isreal(res) || ~isreal(Y)
- disp('Simulation terminated with imaginary parts in the residuals or endogenous variables.')
- else
- disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
- end
+ fprintf('Total time of simulation: %g.\n', etime(clock,h1))
+ disp('Simulation terminated with NaN or Inf in the residuals or endogenous variables.')
display_critical_variables(reshape(dy,[ny periods])', M);
disp('There is most likely something wrong with your model. Try model_diagnostics or another simulation method.')
printline(105)
@@ -198,19 +138,18 @@ if stop
else
if verbose
skipline();
- disp(sprintf('Total time of simulation: %s', num2str(etime(clock,h1))))
+ fprintf('Total time of simulation: %g.\n', etime(clock,h1))
printline(56)
end
info.status = true;% Convergency obtained.
info.error = err;
info.iterations = iter;
info.periods = vperiods(1:iter);
- endogenousvariables = reshape(Y,ny,periods+maximum_lag+M.maximum_lead);
end
elseif ~stop
if verbose
skipline();
- disp(sprintf('Total time of simulation: %s.', num2str(etime(clock,h1))))
+ fprintf('Total time of simulation: %g.\n', etime(clock,h1))
disp('Maximum number of iterations is reached (modify option maxit).')
printline(62)
end
@@ -224,21 +163,21 @@ if verbose
skipline();
end
-function x = lin_solve( A, b,verbose)
-if norm( b ) < sqrt( eps ) % then x = 0 is a solution
+function x = lin_solve(A, b, verbose)
+if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
return
end
x = A\b;
-x( ~isfinite( x ) ) = 0;
-relres = norm( b - A * x ) / norm( b );
+x(~isfinite(x)) = 0;
+relres = norm(b - A*x) / norm(b);
if relres > 1e-6 && verbose
- fprintf( 'WARNING : Failed to find a solution to the linear system.\n' );
+ fprintf('WARNING : Failed to find a solution to the linear system.\n');
end
-function [ x, flag, relres ] = lin_solve_robust( A, b , verbose)
-if norm( b ) < sqrt( eps ) % then x = 0 is a solution
+function [ x, flag, relres ] = lin_solve_robust(A, b ,verbose)
+if norm(b) < sqrt(eps) % then x = 0 is a solution
x = 0;
flag = 0;
relres = 0;
@@ -246,8 +185,8 @@ if norm( b ) < sqrt( eps ) % then x = 0 is a solution
end
x = A\b;
-x( ~isfinite( x ) ) = 0;
-[ x, flag, relres ] = bicgstab( A, b, [], [], [], [], x ); % returns immediately if x is a solution
+x(~isfinite(x)) = 0;
+[ x, flag, relres ] = bicgstab(A, b, [], [], [], [], x); % returns immediately if x is a solution
if flag == 0
return
end
@@ -255,38 +194,38 @@ end
disp( relres );
if verbose
- fprintf( 'Initial bicgstab failed, trying alternative start point.\n' );
+ fprintf('Initial bicgstab failed, trying alternative start point.\n');
end
old_x = x;
old_relres = relres;
-[ x, flag, relres ] = bicgstab( A, b );
+[ x, flag, relres ] = bicgstab(A, b);
if flag == 0
return
end
if verbose
- fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' );
+ fprintf('Alternative start point also failed with bicgstab, trying gmres.\n');
end
if old_relres < relres
x = old_x;
end
-[ x, flag, relres ] = gmres( A, b, [], [], [], [], [], x );
+[ x, flag, relres ] = gmres(A, b, [], [], [], [], [], x);
if flag == 0
return
end
if verbose
- fprintf( 'Initial gmres failed, trying alternative start point.\n' );
+ fprintf('Initial gmres failed, trying alternative start point.\n');
end
old_x = x;
old_relres = relres;
-[ x, flag, relres ] = gmres( A, b );
+[ x, flag, relres ] = gmres(A, b);
if flag == 0
return
end
if verbose
- fprintf( 'Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' );
+ fprintf('Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n');
end
if old_relres < relres
x = old_x;
@@ -294,15 +233,15 @@ if old_relres < relres
end
old_x = x;
old_relres = relres;
-x = pinv( full( A ) ) * b;
-relres = norm( b - A * x ) / norm( b );
+x = pinv(full(A)) * b;
+relres = norm(b - A*x) / norm(b);
if old_relres < relres
x = old_x;
relres = old_relres;
end
flag = relres > 1e-6;
if flag ~= 0 && verbose
- fprintf( 'WARNING : Failed to find a solution to the linear system\n' );
+ fprintf('WARNING : Failed to find a solution to the linear system\n');
end
function display_critical_variables(dyy, M)
@@ -321,10 +260,3 @@ if any(isinf(dyy))
fprintf('%s, ', endo_names{:}),
fprintf('\n'),
end
-if any(~isreal(dyy))
- indx = find(any(~isreal(dyy)));
- endo_names = M.endo_names(indx);
- disp('Last iteration provided complex number for the following variables:')
- fprintf('%s, ', endo_names{:}),
- fprintf('\n'),
-end
diff --git a/matlab/perfect-foresight-models/solve_stacked_problem.m b/matlab/perfect-foresight-models/solve_stacked_problem.m
index 6290d41765..797b23f9da 100644
--- a/matlab/perfect-foresight-models/solve_stacked_problem.m
+++ b/matlab/perfect-foresight-models/solve_stacked_problem.m
@@ -53,11 +53,10 @@ if (options.solve_algo == 10 || options.solve_algo == 11)% mixed complementarity
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0, ...
eq_index);
else
- [y, check] = dynare_solve(@perfect_foresight_problem,z(:),options, ...
- dynamicmodel, y0, yT, ...
- exogenousvariables, M.params, steadystate, ...
- M.maximum_lag, options.periods, M.endo_nbr, i_cols, ...
- i_cols_J1, i_cols_1, i_cols_T, i_cols_j, i_cols_0, i_cols_J0);
+ model_dynamic_g1_nz = str2func([M.fname,'.dynamic_g1_nz']);
+ [nzij_pred, nzij_current, nzij_fwrd] = model_dynamic_g1_nz();
+
+ [y, check] = dynare_solve(@perfect_foresight_problem,z(:), options, M.fname, sum(M.dynamic_tmp_nbr(1:2)), y0, yT, exogenousvariables, M.params, steadystate, options.periods, M.endo_nbr, M.maximum_lag, M.maximum_endo_lag, M.lead_lag_incidence, nzij_pred, nzij_current, nzij_fwrd, M.has_external_function, options.use_dll);
end
if all(imag(y)<.1*options.dynatol.x)
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index 0f72250826..8d555f1adc 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -1,7 +1,7 @@
ACLOCAL_AMFLAGS = -I ../../../m4
if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations
+SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if HAVE_MATIO
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index eb8f8b7b03..906aeee2e7 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -172,6 +172,7 @@ AC_CONFIG_FILES([Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
- local_state_space_iterations/Makefile])
+ local_state_space_iterations/Makefile
+ perfect_foresight_problem/Makefile])
AC_OUTPUT
diff --git a/mex/build/matlab/perfect_foresight_problem/Makefile.am b/mex/build/matlab/perfect_foresight_problem/Makefile.am
new file mode 100644
index 0000000000..d2ad7ff72c
--- /dev/null
+++ b/mex/build/matlab/perfect_foresight_problem/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../perfect_foresight_problem.am
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 704a86bf08..4703f18c61 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -1,7 +1,7 @@
ACLOCAL_AMFLAGS = -I ../../../m4
if DO_SOMETHING
-SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations
+SUBDIRS = mjdgges kronecker bytecode block_kalman_filter sobol local_state_space_iterations perfect_foresight_problem
# libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_
if HAVE_MATIO
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 478b9b22f7..7bb01f7047 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -1,6 +1,6 @@
dnl Process this file with autoconf to produce a configure script.
-dnl Copyright © 2009-2017 Dynare Team
+dnl Copyright © 2009-2019 Dynare Team
dnl
dnl This file is part of Dynare.
dnl
@@ -135,6 +135,7 @@ AC_CONFIG_FILES([Makefile
ms_sbvar/Makefile
block_kalman_filter/Makefile
sobol/Makefile
- local_state_space_iterations/Makefile])
+ local_state_space_iterations/Makefile
+ perfect_foresight_problem/Makefile])
AC_OUTPUT
diff --git a/mex/build/octave/perfect_foresight_problem/Makefile.am b/mex/build/octave/perfect_foresight_problem/Makefile.am
new file mode 100644
index 0000000000..65ceb979fa
--- /dev/null
+++ b/mex/build/octave/perfect_foresight_problem/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../perfect_foresight_problem.am
diff --git a/mex/build/perfect_foresight_problem.am b/mex/build/perfect_foresight_problem.am
new file mode 100644
index 0000000000..e295c7e252
--- /dev/null
+++ b/mex/build/perfect_foresight_problem.am
@@ -0,0 +1,15 @@
+mex_PROGRAMS = perfect_foresight_problem
+
+TOPDIR = $(top_srcdir)/../../sources/perfect_foresight_problem
+
+nodist_perfect_foresight_problem_SOURCES = perfect_foresight_problem.cc DynamicModelCaller.cc
+
+perfect_foresight_problem_CPPFLAGS = $(AM_CPPFLAGS) -I$(TOPDIR)
+perfect_foresight_problem_CXXFLAGS = $(AM_CXXFLAGS) -fopenmp
+perfect_foresight_problem_LDADD = $(LIBADD_DLOPEN)
+
+BUILT_SOURCES = $(nodist_perfect_foresight_problem_SOURCES)
+CLEANFILES = $(nodist_perfect_foresight_problem_SOURCES)
+
+%.cc: $(TOPDIR)/%.cc
+ $(LN_S) -f $< $@
diff --git a/mex/sources/Makefile.am b/mex/sources/Makefile.am
index 654100d3ab..6b469e1cfd 100644
--- a/mex/sources/Makefile.am
+++ b/mex/sources/Makefile.am
@@ -13,7 +13,8 @@ EXTRA_DIST = \
sobol \
local_state_space_iterations \
gensylv \
- dynare_simul_
+ dynare_simul_ \
+ perfect_foresight_problem
clean-local:
rm -rf `find mex/sources -name *.o`
diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
new file mode 100644
index 0000000000..2a59e0937e
--- /dev/null
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.cc
@@ -0,0 +1,215 @@
+/*
+ * Copyright © 2019 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#if defined(_WIN32) || defined(__CYGWIN32__)
+# ifndef NOMINMAX
+# define NOMINMAX // Do not define "min" and "max" macros
+# endif
+# include <windows.h>
+#else
+# include <dlfcn.h> // unix/linux DLL (.so) handling routines
+#endif
+
+#include <dynmex.h>
+
+#include "DynamicModelCaller.hh"
+
+#include <algorithm>
+
+std::string DynamicModelCaller::error_msg;
+
+void *DynamicModelDllCaller::dynamic_mex{nullptr};
+DynamicModelDllCaller::dynamic_tt_fct DynamicModelDllCaller::residual_tt_fct{nullptr}, DynamicModelDllCaller::g1_tt_fct{nullptr};
+DynamicModelDllCaller::dynamic_deriv_fct DynamicModelDllCaller::residual_fct{nullptr}, DynamicModelDllCaller::g1_fct{nullptr};
+
+void
+DynamicModelDllCaller::load_dll(const std::string &basename)
+{
+ // Load symbols from dynamic MEX
+ std::string mex_name;
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ mex_name = "./";
+#endif
+ mex_name += "+" + basename + "/dynamic" + MEXEXT;
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ dynamic_mex = dlopen(mex_name.c_str(), RTLD_NOW);
+#else
+ dynamic_mex = LoadLibrary(mex_name.c_str());
+#endif
+ if (!dynamic_mex)
+ mexErrMsgTxt("Can't load dynamic MEX file");
+
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamic_mex, "dynamic_resid_tt"));
+ residual_fct = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamic_mex, "dynamic_resid"));
+ g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamic_mex, "dynamic_g1_tt"));
+ g1_fct = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamic_mex, "dynamic_g1"));
+#else
+ residual_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamic_mex, "dynamic_resid_tt"));
+ residual_fct = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamic_mex, "dynamic_resid"));
+ g1_tt_fct = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamic_mex, "dynamic_g1_tt"));
+ g1_fct = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamic_mex, "dynamic_g1"));
+#endif
+ if (!residual_tt_fct || !residual_fct || !g1_tt_fct || !g1_fct)
+ mexErrMsgTxt("Can't load functions in dynamic MEX file");
+}
+
+void
+DynamicModelDllCaller::unload_dll()
+{
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+ dlclose(dynamic_mex);
+#else
+ FreeLibrary(dynamic_mex);
+#endif
+}
+
+DynamicModelDllCaller::DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool compute_jacobian_arg) :
+ DynamicModelCaller{compute_jacobian_arg},
+ nb_row_x{nb_row_x_arg}, x{x_arg}, params{params_arg}, steady_state{steady_state_arg}
+{
+ tt = std::make_unique<double[]>(ntt);
+ y_p = std::make_unique<double[]>(ndynvars);
+ if (compute_jacobian)
+ jacobian_p = std::make_unique<double[]>((ndynvars+nx)*ny);
+}
+
+void
+DynamicModelDllCaller::eval(int it, double *resid)
+{
+ residual_tt_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get());
+ residual_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get(), resid);
+ if (compute_jacobian)
+ {
+ g1_tt_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get());
+ g1_fct(y_p.get(), x, nb_row_x, params, steady_state, it, tt.get(), jacobian_p.get());
+ }
+}
+
+DynamicModelMatlabCaller::DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool compute_jacobian_arg) :
+ DynamicModelCaller{compute_jacobian_arg},
+ basename{std::move(basename_arg)},
+ T_mx{mxCreateDoubleMatrix(ntt, 1, mxREAL)},
+ y_mx{mxCreateDoubleMatrix(ndynvars, 1, mxREAL)},
+ it_mx{mxCreateDoubleScalar(0)},
+ T_flag_mx{mxCreateLogicalScalar(false)},
+ jacobian_mx{nullptr},
+ x_mx{mxDuplicateArray(x_mx_arg)},
+ params_mx{mxDuplicateArray(params_mx_arg)},
+ steady_state_mx{mxDuplicateArray(steady_state_mx_arg)}
+{
+}
+
+DynamicModelMatlabCaller::~DynamicModelMatlabCaller()
+{
+ mxDestroyArray(T_mx);
+ mxDestroyArray(y_mx);
+ mxDestroyArray(it_mx);
+ mxDestroyArray(T_flag_mx);
+ if (jacobian_mx)
+ mxDestroyArray(jacobian_mx);
+ mxDestroyArray(x_mx);
+ mxDestroyArray(params_mx);
+ mxDestroyArray(steady_state_mx);
+}
+
+void
+DynamicModelMatlabCaller::eval(int it, double *resid)
+{
+ *mxGetPr(it_mx) = it + 1;
+
+ {
+ // Compute temporary terms
+ std::string funcname = basename + (compute_jacobian ? ".dynamic_g1_tt" : ".dynamic_resid_tt");
+ mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx };
+
+ mxArray *exception = mexCallMATLABWithTrap(1, plhs, 6, prhs, funcname.c_str());
+ if (exception)
+ error_msg = std::string{"An error occurred when calling "} + funcname;
+
+ mxDestroyArray(T_mx);
+ T_mx = plhs[0];
+ }
+
+ {
+ // Compute residuals
+ std::string funcname = basename + ".dynamic_resid";
+ mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx, T_flag_mx };
+
+ mxArray *exception = mexCallMATLABWithTrap(1, plhs, 7, prhs, funcname.c_str());
+ if (exception)
+ error_msg = std::string{"An error occurred when calling "} + funcname;
+
+ if (mxIsComplex(plhs[0]))
+ plhs[0] = cmplxToReal(plhs[0]);
+
+ std::copy_n(mxGetPr(plhs[0]), mxGetNumberOfElements(plhs[0]), resid);
+ mxDestroyArray(plhs[0]);
+ }
+
+ if (compute_jacobian)
+ {
+ // Compute Jacobian
+ std::string funcname = basename + ".dynamic_g1";
+ mxArray *plhs[1], *prhs[] = { T_mx, y_mx, x_mx, params_mx, steady_state_mx, it_mx, T_flag_mx };
+
+ mxArray *exception = mexCallMATLABWithTrap(1, plhs, 7, prhs, funcname.c_str());
+ if (exception)
+ error_msg = std::string{"An error occurred when calling "} + funcname;
+
+ if (jacobian_mx)
+ {
+ mxDestroyArray(jacobian_mx);
+ jacobian_mx = nullptr;
+ }
+
+ if (mxIsComplex(plhs[0]))
+ jacobian_mx = cmplxToReal(plhs[0]);
+ else
+ jacobian_mx = plhs[0];
+ }
+}
+
+mxArray *
+DynamicModelMatlabCaller::cmplxToReal(mxArray *cmplx_mx)
+{
+ mxArray *real_mx = mxCreateDoubleMatrix(mxGetM(cmplx_mx), mxGetN(cmplx_mx), mxREAL);
+
+#if MX_HAS_INTERLEAVED_COMPLEX
+ mxComplexDouble *cmplx = mxGetComplexDoubles(cmplx_mx);
+#else
+ double *cmplx_real = mxGetPr(cmplx_mx);
+ double *cmplx_imag = mxGetPi(cmplx_mx);
+#endif
+ double *real = mxGetPr(real_mx);
+
+ for (size_t i = 0; i < mxGetNumberOfElements(cmplx_mx); i++)
+#if MX_HAS_INTERLEAVED_COMPLEX
+ if (cmplx[i].imag == 0.0)
+ real[i] = cmplx[i].real;
+#else
+ if (cmplx_imag[i] == 0.0)
+ real[i] = cmplx_real[i];
+#endif
+ else
+ real[i] = std::numeric_limits<double>::quiet_NaN();
+
+ mxDestroyArray(cmplx_mx);
+ return real_mx;
+}
diff --git a/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
new file mode 100644
index 0000000000..8c64938fdb
--- /dev/null
+++ b/mex/sources/perfect_foresight_problem/DynamicModelCaller.hh
@@ -0,0 +1,83 @@
+/*
+ * Copyright © 2019 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <string>
+#include <memory>
+#include <limits>
+
+class DynamicModelCaller
+{
+public:
+ const bool compute_jacobian;
+
+ // Used to store error messages (as exceptions cannot cross the OpenMP boundary)
+ static std::string error_msg;
+
+ DynamicModelCaller(bool compute_jacobian_arg) : compute_jacobian{compute_jacobian_arg} {};
+ virtual ~DynamicModelCaller() = default;
+ virtual double &y(size_t i) const = 0;
+ virtual double jacobian(size_t i) const = 0;
+ virtual void eval(int it, double *resid) = 0;
+};
+
+class DynamicModelDllCaller : public DynamicModelCaller
+{
+private:
+ static void *dynamic_mex;
+ using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);
+ using dynamic_deriv_fct = void (*) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv);
+ static dynamic_tt_fct residual_tt_fct, g1_tt_fct;
+ static dynamic_deriv_fct residual_fct, g1_fct;
+ size_t nb_row_x;
+ const double *x, *params, *steady_state;
+ std::unique_ptr<double[]> tt, y_p, jacobian_p;
+
+public:
+ DynamicModelDllCaller(size_t ntt, mwIndex nx, mwIndex ny, size_t ndynvars, const double *x_arg, size_t nb_row_x_arg, const double *params_arg, const double *steady_state_arg, bool compute_jacobian_arg);
+ virtual ~DynamicModelDllCaller() = default;
+ double &y(size_t i) const override { return y_p[i]; };
+ double jacobian(size_t i) const override { return jacobian_p[i]; };
+ void eval(int it, double *resid) override;
+ static void load_dll(const std::string &basename);
+ static void unload_dll();
+};
+
+class DynamicModelMatlabCaller : public DynamicModelCaller
+{
+private:
+ std::string basename;
+ mxArray *T_mx, *y_mx, *it_mx, *T_flag_mx, *jacobian_mx, *x_mx, *params_mx, *steady_state_mx;
+ /* Given a complex matrix, returns a real matrix of same size.
+ Real elements of the original matrix are copied as-is to the new one.
+ Complex elements are replaced by NaNs.
+ Destroys the original matrix. */
+ static mxArray *cmplxToReal(mxArray *m);
+public:
+ DynamicModelMatlabCaller(std::string basename_arg, size_t ntt, size_t ndynvars, const mxArray *x_mx_arg, const mxArray *params_mx_arg, const mxArray *steady_state_mx_arg, bool compute_jacobian_arg);
+ ~DynamicModelMatlabCaller() override;
+ double &y(size_t i) const override { return mxGetPr(y_mx)[i]; };
+ double jacobian(size_t i) const override { return jacobian_mx ? mxGetPr(jacobian_mx)[i] : std::numeric_limits<double>::quiet_NaN(); };
+ void eval(int it, double *resid) override;
+ class Exception {
+ public:
+ const std::string msg;
+ Exception(std::string msg_arg) : msg{std::move(msg_arg)} {};
+ };
+};
+
diff --git a/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
new file mode 100644
index 0000000000..f743c92b02
--- /dev/null
+++ b/mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc
@@ -0,0 +1,277 @@
+/*
+ * Copyright © 2019 Dynare Team
+ *
+ * This file is part of Dynare.
+ *
+ * Dynare is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Dynare is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <string>
+#include <memory>
+#include <algorithm>
+
+#include <dynmex.h>
+
+#include "DynamicModelCaller.hh"
+
+void
+mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+ if (nlhs < 1 || nlhs > 2 || nrhs != 18)
+ mexErrMsgTxt("Must have 18 input arguments and 1 or 2 output arguments");
+ bool compute_jacobian = nlhs == 2;
+
+ // Give explicit names to input arguments
+ const mxArray *y_mx = prhs[0];
+ const mxArray *basename_mx = prhs[1];
+ const mxArray *ntt_mx = prhs[2];
+ const mxArray *y0_mx = prhs[3];
+ const mxArray *yT_mx = prhs[4];
+ const mxArray *exo_path_mx = prhs[5];
+ const mxArray *params_mx = prhs[6];
+ const mxArray *steady_state_mx = prhs[7];
+ const mxArray *periods_mx = prhs[8];
+ const mxArray *ny_mx = prhs[9];
+ const mxArray *maximum_lag_mx = prhs[10];
+ const mxArray *maximum_endo_lag_mx = prhs[11];
+ const mxArray *lead_lag_incidence_mx = prhs[12];
+ const mxArray *nzij_pred_mx = prhs[13];
+ const mxArray *nzij_current_mx = prhs[14];
+ const mxArray *nzij_fwrd_mx = prhs[15];
+ const mxArray *has_external_function_mx = prhs[16];
+ const mxArray *use_dll_mx = prhs[17];
+
+ // Check input and map it to local variables
+ if (!(mxIsChar(basename_mx) && mxGetM(basename_mx) == 1))
+ mexErrMsgTxt("basename should be a character string");
+ char *basename = mxArrayToString(basename_mx);
+
+ if (!(mxIsScalar(ntt_mx) && mxIsNumeric(ntt_mx)))
+ mexErrMsgTxt("ntt should be a numeric scalar");
+ size_t ntt = mxGetScalar(ntt_mx);
+
+ if (!(mxIsScalar(periods_mx) && mxIsNumeric(periods_mx)))
+ mexErrMsgTxt("periods should be a numeric scalar");
+ mwIndex periods = static_cast<mwIndex>(mxGetScalar(periods_mx));
+
+ if (!(mxIsScalar(ny_mx) && mxIsNumeric(ny_mx)))
+ mexErrMsgTxt("ny should be a numeric scalar");
+ mwIndex ny = static_cast<mwIndex>(mxGetScalar(ny_mx));
+
+ if (!(mxIsScalar(maximum_lag_mx) && mxIsNumeric(maximum_lag_mx)))
+ mexErrMsgTxt("maximum_lag should be a numeric scalar");
+ mwIndex maximum_lag = static_cast<mwIndex>(mxGetScalar(maximum_lag_mx));
+
+ if (!(mxIsScalar(maximum_endo_lag_mx) && mxIsNumeric(maximum_endo_lag_mx)))
+ mexErrMsgTxt("maximum_endo_lag should be a numeric scalar");
+ mwIndex maximum_endo_lag = static_cast<mwIndex>(mxGetScalar(maximum_endo_lag_mx));
+
+ if (!(mxIsDouble(y_mx) && mxGetM(y_mx) == static_cast<size_t>(ny*periods) && mxGetN(y_mx) == 1))
+ mexErrMsgTxt("y should be a double precision column-vector of ny*periods elements");
+ const double *y = mxGetPr(y_mx);
+
+ if (!(mxIsDouble(y0_mx) && mxGetM(y0_mx) == static_cast<size_t>(ny) && mxGetN(y0_mx) == 1))
+ mexErrMsgTxt("y0 should be a double precision column-vector of ny elements");
+ const double *y0 = mxGetPr(y0_mx);
+
+ if (!(mxIsDouble(yT_mx) && mxGetM(yT_mx) == static_cast<size_t>(ny) && mxGetN(yT_mx) == 1))
+ mexErrMsgTxt("yT should be a double precision column-vector of ny elements");
+ const double *yT = mxGetPr(yT_mx);
+
+ if (!(mxIsDouble(exo_path_mx) && mxGetM(exo_path_mx) >= static_cast<size_t>(periods+maximum_lag)))
+ mexErrMsgTxt("exo_path should be a double precision matrix with at least periods+maximum_lag rows");
+ mwIndex nx = static_cast<mwIndex>(mxGetN(exo_path_mx));
+ size_t nb_row_x = mxGetM(exo_path_mx);
+ const double *exo_path = mxGetPr(exo_path_mx);
+
+ if (!(mxIsDouble(params_mx) && mxGetN(params_mx) == 1))
+ mexErrMsgTxt("params should be a double precision column-vector");
+ const double *params = mxGetPr(params_mx);
+
+ if (!(mxIsDouble(steady_state_mx) && mxGetN(steady_state_mx) == 1))
+ mexErrMsgTxt("steady_state should be a double precision column-vector");
+ const double *steady_state = mxGetPr(steady_state_mx);
+
+ if (!(mxIsDouble(lead_lag_incidence_mx) && mxGetM(lead_lag_incidence_mx) == static_cast<size_t>(2+maximum_endo_lag)
+ && mxGetN(lead_lag_incidence_mx) == static_cast<size_t>(ny)))
+ mexErrMsgTxt("lead_lag_incidence should be a double precision matrix with 2+maximum_endo_lag rows and endo_nbr columns");
+ const double *lead_lag_incidence = mxGetPr(lead_lag_incidence_mx);
+
+ if (!(mxIsInt32(nzij_pred_mx) && mxGetN(nzij_pred_mx) == 2))
+ mexErrMsgTxt("nzij_pred should be an int32 matrix with 2 columns");
+ size_t nnz_pred = mxGetM(nzij_pred_mx);
+ const int32_T *nzij_pred = static_cast<const int32_T *>(mxGetData(nzij_pred_mx));
+
+ if (!(mxIsInt32(nzij_current_mx) && mxGetN(nzij_current_mx) == 2))
+ mexErrMsgTxt("nzij_current should be an int32 matrix with 2 columns");
+ size_t nnz_current = mxGetM(nzij_current_mx);
+ const int32_T *nzij_current = static_cast<const int32_T *>(mxGetData(nzij_current_mx));
+
+ if (!(mxIsInt32(nzij_fwrd_mx) && mxGetN(nzij_fwrd_mx) == 2))
+ mexErrMsgTxt("nzij_fwrd should be an int32 matrix with 2 columns");
+ size_t nnz_fwrd = mxGetM(nzij_fwrd_mx);
+ const int32_T *nzij_fwrd = static_cast<const int32_T *>(mxGetData(nzij_fwrd_mx));
+
+ if (!(mxIsLogicalScalar(has_external_function_mx)))
+ mexErrMsgTxt("has_external_function should be a logical scalar");
+ bool has_external_function = static_cast<bool>(mxGetScalar(has_external_function_mx));
+
+ if (!(mxIsLogicalScalar(use_dll_mx)))
+ mexErrMsgTxt("use_dll should be a logical scalar");
+ bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
+
+ // Allocate output matrices
+ plhs[0] = mxCreateDoubleMatrix(periods*ny, 1, mxREAL);
+ double *stacked_residual = mxGetPr(plhs[0]);
+
+ mwIndex nzmax = periods*nnz_current+(periods-1)*(nnz_pred+nnz_fwrd);
+
+ double *stacked_jacobian = nullptr;
+ mwIndex *ir = nullptr, *jc = nullptr;
+ if (compute_jacobian)
+ {
+ plhs[1] = mxCreateSparse(periods*ny, periods*ny, nzmax, mxREAL);
+ stacked_jacobian = mxGetPr(plhs[1]);
+ ir = mxGetIr(plhs[1]);
+ jc = mxGetJc(plhs[1]);
+
+ /* Create the index vectors (IR, JC) of the sparse stacked jacobian. This
+ makes parallelization across periods possible when evaluating the model,
+ since all indices are known ex ante. */
+ mwIndex k = 0;
+ jc[0] = 0;
+ for (mwIndex T = 0; T < periods; T++)
+ {
+ size_t row_pred = 0, row_current = 0, row_fwrd = 0;
+ for (int32_T j = 0; j < static_cast<int32_T>(ny); j++)
+ {
+ if (T != 0)
+ while (row_fwrd < nnz_fwrd && nzij_fwrd[row_fwrd+nnz_fwrd]-1 == j)
+ ir[k++] = (T-1)*ny + nzij_fwrd[row_fwrd++]-1;
+ while (row_current < nnz_current && nzij_current[row_current+nnz_current]-1 == j)
+ ir[k++] = T*ny + nzij_current[row_current++]-1;
+ if (T != periods-1)
+ while (row_pred < nnz_pred && nzij_pred[row_pred+nnz_pred]-1 == j)
+ ir[k++] = (T+1)*ny + nzij_pred[row_pred++]-1;
+ jc[T*ny+j+1] = k;
+ }
+ }
+ }
+
+ size_t ndynvars = static_cast<size_t>(*std::max_element(lead_lag_incidence, lead_lag_incidence+(maximum_endo_lag+2)*ny));
+
+ if (use_dll)
+ DynamicModelDllCaller::load_dll(basename);
+
+ DynamicModelCaller::error_msg.clear();
+
+ /* Parallelize the main loop, if use_dll and no external function (to avoid
+ parallel calls to MATLAB) */
+#pragma omp parallel if (use_dll && !has_external_function)
+ {
+ // Allocate (thread-private) model evaluator (which allocates space for temporaries)
+ std::unique_ptr<DynamicModelCaller> m;
+ if (use_dll)
+ m = std::make_unique<DynamicModelDllCaller>(ntt, nx, ny, ndynvars, exo_path, nb_row_x, params, steady_state, compute_jacobian);
+ else
+ m = std::make_unique<DynamicModelMatlabCaller>(basename, ntt, ndynvars, exo_path_mx, params_mx, steady_state_mx, compute_jacobian);
+
+ // Main computing loop
+#pragma omp for
+ for (mwIndex T = 0; T < periods; T++)
+ {
+ // Fill vector of dynamic variables
+ for (mwIndex j = 0; j < maximum_endo_lag+2; j++)
+ for (mwIndex i = 0; i < ny; i++)
+ {
+ int idx = static_cast<int>(lead_lag_incidence[j+i*(2+maximum_endo_lag)])-1;
+ if (idx != -1)
+ {
+ if (T+j == maximum_endo_lag-1)
+ m->y(idx) = y0[i];
+ else if (T+j == maximum_endo_lag+periods)
+ m->y(idx) = yT[i];
+ else
+ m->y(idx) = y[i+(T+j-maximum_endo_lag)*ny];
+ }
+ }
+
+ // Compute the residual and Jacobian, and fill the stacked residual
+ m->eval(T+maximum_lag, stacked_residual+T*ny);
+
+ if (compute_jacobian)
+ {
+ // Fill the stacked jacobian
+ for (mwIndex col = T > maximum_endo_lag ? (T-maximum_endo_lag)*ny : 0; // We can't use std::max() here, because mwIndex is unsigned under MATLAB
+ col < std::min(periods*ny, (T+2)*ny); col++)
+ {
+ mwIndex k = jc[col];
+ while (k < jc[col+1])
+ {
+ if (ir[k] < T*ny)
+ {
+ k++;
+ continue;
+ }
+ if (ir[k] >= (T+1)*ny)
+ break;
+
+ mwIndex eq = ir[k]-T*ny;
+ mwIndex lli_row = col/ny-(T-maximum_endo_lag); // 0, 1 or 2
+ mwIndex lli_col = col%ny;
+ mwIndex dynvar = static_cast<mwIndex>(lead_lag_incidence[lli_row+lli_col*(2+maximum_endo_lag)])-1;
+ stacked_jacobian[k] = m->jacobian(eq+dynvar*ny);
+ k++;
+ }
+ }
+ }
+ }
+ }
+
+ /* Mimic a try/catch using a global string, since exceptions are not allowed
+ to cross OpenMP boundary */
+ if (!DynamicModelCaller::error_msg.empty())
+ mexErrMsgTxt(DynamicModelCaller::error_msg.c_str());
+
+ if (compute_jacobian)
+ {
+ /* The stacked jacobian so far constructed can still reference some zero
+ elements, because some expressions that are symbolically different from
+ zero may still evaluate to zero for some endogenous/parameter values. The
+ following code further compresses the Jacobian by removing the references
+ to those extra zeros. This makes a significant speed difference when
+ inversing the Jacobian for some large models. */
+ mwIndex k_orig = 0, k_new = 0;
+ for (mwIndex col = 0; col < periods*ny; col++)
+ {
+ while (k_orig < jc[col+1])
+ {
+ if (stacked_jacobian[k_orig] != 0.0)
+ {
+ if (k_new != k_orig)
+ {
+ stacked_jacobian[k_new] = stacked_jacobian[k_orig];
+ ir[k_new] = ir[k_orig];
+ }
+ k_new++;
+ }
+ k_orig++;
+ }
+ jc[col+1] = k_new;
+ }
+ }
+
+ if (use_dll)
+ DynamicModelDllCaller::unload_dll();
+}
--
GitLab