perfect_foresight_problem MEX: optimization for linear models

When the model is linear, there is no need to reevaluate the Jacobian for each
time period, since it is invariant.

Closes: #1662

mex/sources/perfect_foresight_problem/DynamicModelCaller.cc

@@ -75,8 +75,8 @@ DynamicModelDllCaller::unload_dll()
 #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},
+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 linear_arg, bool compute_jacobian_arg) :
+  DynamicModelCaller{linear_arg, 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);
@@ -94,11 +94,14 @@ DynamicModelDllCaller::eval(int it, double *resid)
     {
       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());
+
+      if (linear)
+        compute_jacobian = false; // If model is linear, no need to recompute Jacobian later
     }
 }
 
-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},
+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 linear_arg, bool compute_jacobian_arg) :
+  DynamicModelCaller{linear_arg, compute_jacobian_arg},
   basename{std::move(basename_arg)},
   T_mx{mxCreateDoubleMatrix(ntt, 1, mxREAL)},
   y_mx{mxCreateDoubleMatrix(ndynvars, 1, mxREAL)},
@@ -178,6 +181,9 @@ DynamicModelMatlabCaller::eval(int it, double *resid)
         jacobian_mx = cmplxToReal(plhs[0]);
       else
         jacobian_mx = plhs[0];
+
+      if (linear)
+        compute_jacobian = false; // If model is linear, no need to recompute Jacobian later
     }
 }
 

mex/sources/perfect_foresight_problem/DynamicModelCaller.hh

@@ -33,12 +33,13 @@
 class DynamicModelCaller
 {
 public:
-  const bool compute_jacobian;
+  const bool linear;
+  bool compute_jacobian; // Not constant, because will be changed from true to false for linear models after first evaluation
 
   // 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} {};
+  DynamicModelCaller(bool linear_arg, bool compute_jacobian_arg) : linear{linear_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;
@@ -62,7 +63,7 @@ private:
   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);
+  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 linear_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]; };
@@ -82,7 +83,7 @@ private:
      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(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 linear_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(); };

mex/sources/perfect_foresight_problem/perfect_foresight_problem.cc

@@ -86,6 +86,11 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     mexErrMsgTxt("options_.use_dll should be a logical scalar");
   bool use_dll = static_cast<bool>(mxGetScalar(use_dll_mx));
 
+  const mxArray *linear_mx = mxGetField(options_mx, 0, "linear");
+  if (!(linear_mx && mxIsLogicalScalar(linear_mx)))
+    mexErrMsgTxt("options_.linear should be a logical scalar");
+  bool linear = static_cast<bool>(mxGetScalar(linear_mx));
+
   const mxArray *threads_mx = mxGetField(options_mx, 0, "threads");
   if (!threads_mx)
     mexErrMsgTxt("Can't find field options_.threads");
@@ -211,9 +216,9 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     // 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);
+      m = std::make_unique<DynamicModelDllCaller>(ntt, nx, ny, ndynvars, exo_path, nb_row_x, params, steady_state, linear, compute_jacobian);
     else
-      m = std::make_unique<DynamicModelMatlabCaller>(basename, ntt, ndynvars, exo_path_mx, params_mx, steady_state_mx, compute_jacobian);
+      m = std::make_unique<DynamicModelMatlabCaller>(basename, ntt, ndynvars, exo_path_mx, params_mx, steady_state_mx, linear, compute_jacobian);
 
     // Main computing loop
 #pragma omp for

tests/deterministic_simulations/multiple_lead_lags/AR2_forward.mod

@@ -6,7 +6,7 @@ rho_1=0.2;
 rho_2=0.1;
 
 // Equilibrium conditions
-model;
+model(linear);
 y_backward=rho_1*y_backward(-1)+rho_2*y_backward(-2);
 dummy_var=0.9*dummy_var(+1);
 end;