Convert ModelSolution to Eigen

mex/sources/estimation/ModelSolution.cc

@@ -30,18 +30,19 @@
 /**
  * compute the steady state (2nd stage), and computes first order approximation
  */
-ModelSolution::ModelSolution(const std::string &dynamicDllFile,  size_t n_endo_arg, size_t n_exo_arg, const std::vector<size_t> &zeta_fwrd_arg,
-                             const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
-                             const std::vector<size_t> &zeta_static_arg, double INqz_criterium) :
+ModelSolution::ModelSolution(const std::string &dynamicDllFile, ptrdiff_t n_endo_arg,
+                             ptrdiff_t n_exo_arg, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
+                             const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
+                             const std::vector<ptrdiff_t> &zeta_static_arg, double INqz_criterium) :
   n_endo(n_endo_arg), n_exo(n_exo_arg),  // n_jcols = Num of Jacobian columns = nStat+2*nPred+3*nBoth+2*nForw+nExog
   n_jcols(n_exo+n_endo+ zeta_back_arg.size() /*nsPred*/ + zeta_fwrd_arg.size() /*nsForw*/ +2*zeta_mixed_arg.size()),
   jacobian(n_endo, n_jcols), residual(n_endo), Mx(1, n_exo),
   decisionRules(n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, INqz_criterium),
-  dynamicDLLp(dynamicDllFile, n_exo),
+  dynamicDLLp(dynamicDllFile),
   llXsteadyState(n_jcols-n_exo)
 {
-  Mx.setAll(0.0);
-  jacobian.setAll(0.0);
+  Mx.setZero();
+  jacobian.setZero();
 
   set_union(zeta_fwrd_arg.begin(), zeta_fwrd_arg.end(),
             zeta_mixed_arg.begin(), zeta_mixed_arg.end(),

mex/sources/estimation/ModelSolution.hh

@@ -26,6 +26,10 @@
 #if !defined(ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_)
 #define ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_
 
+#include <Eigen/Core>
+
+using namespace Eigen;
+
 #include "DecisionRules.hh"
 #include "dynamic_dll.hh"
 
@@ -38,12 +42,14 @@ class ModelSolution
 {
 
 public:
-  ModelSolution(const std::string &dynamicDllFile,  size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg,
-                const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg,
-                const std::vector<size_t> &zeta_static_arg, double qz_criterium);
+  ModelSolution(const std::string &dynamicDllFile, ptrdiff_t n_endo, ptrdiff_t n_exo, const std::vector<ptrdiff_t> &zeta_fwrd_arg,
+                const std::vector<ptrdiff_t> &zeta_back_arg, const std::vector<ptrdiff_t> &zeta_mixed_arg,
+                const std::vector<ptrdiff_t> &zeta_static_arg, double qz_criterium);
   virtual ~ModelSolution() {};
-  template <class VEC>
-  void compute(VEC &steadyState, const VectorView &deepParams,      Matrix &ghx, Matrix &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+  template <class Derived1, class Derived2>
+  void compute(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams,
+               MatrixXd &ghx, MatrixXd &ghu) throw (DecisionRules::BlanchardKahnException,
+                                                    GeneralizedSchurDecomposition::GSDException)
   {
     // compute Steady State
     ComputeSteadyState(steadyState, deepParams);
@@ -55,31 +61,31 @@ public:
   }
 
 private:
-  const size_t n_endo;
-  const size_t n_exo;
-  const size_t n_jcols; // Num of Jacobian columns
-  std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed;
-  Matrix jacobian;
-  Vector residual;
-  Matrix Mx;
+  const ptrdiff_t n_endo;
+  const ptrdiff_t n_exo;
+  const ptrdiff_t n_jcols; // Num of Jacobian columns
+  std::vector<ptrdiff_t> zeta_fwrd_mixed, zeta_back_mixed;
+  MatrixXd jacobian;
+  VectorXd residual;
+  MatrixXd Mx;
   DecisionRules decisionRules;
   DynamicModelDLL dynamicDLLp;
-  Vector llXsteadyState;
+  VectorXd llXsteadyState;
   //Matrix jacobian;
-  template <class VEC>
-  void ComputeModelSolution(VEC &steadyState, const VectorView &deepParams,         
-			    Matrix &ghx, Matrix &ghu) 
+  template <class Derived1, class Derived2>
+  void ComputeModelSolution(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams,         
+			    MatrixXd &ghx, MatrixXd &ghu) 
     throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
   {
     // set extended Steady State
 
-    for (size_t i = 0; i < zeta_back_mixed.size(); i++)
+    for (ptrdiff_t i = 0; i < (ptrdiff_t) zeta_back_mixed.size(); i++)
       llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
 
-    for (size_t i = 0; i < n_endo; i++)
+    for (ptrdiff_t i = 0; i < n_endo; i++)
       llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
 
-    for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
+    for (ptrdiff_t i = 0; i < (ptrdiff_t) zeta_fwrd_mixed.size(); i++)
       llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
 
     //get jacobian
@@ -88,8 +94,8 @@ private:
     //compute rules
     decisionRules.compute(jacobian, ghx, ghu);
   }
-  template <class VEC>
-  void ComputeSteadyState(VEC &steadyState, const VectorView &deepParams)
+  template <class Derived1, class Derived2>
+  void ComputeSteadyState(PlainObjectBase<Derived1> &steadyState, const PlainObjectBase<Derived2> &deepParams)
   {
     // does nothig for time being.
   }

mex/sources/estimation/tests/Makefile.am

@@ -4,7 +4,7 @@ test_dr_SOURCES = ../libmat/GeneralizedSchurDecomposition.cc ../libmat/QRDecompo
 test_dr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
 test_dr_CPPFLAGS = -I.. -I../libmat -I../../
 
-testModelSolution_SOURCES = ../libmat/Matrix.cc ../libmat/Vector.cc ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../libmat/LUSolver.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc testModelSolution.cc
+testModelSolution_SOURCES = ../libmat/QRDecomposition.cc ../libmat/GeneralizedSchurDecomposition.cc ../utils/dynamic_dll.cc ../DecisionRules.cc ../ModelSolution.cc testModelSolution.cc
 testModelSolution_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_DLOPEN)
 testModelSolution_CPPFLAGS = -I.. -I../libmat -I../../ -I../utils
 

mex/sources/estimation/tests/model1.mod

+/* To be used with testModelSolution; derived from fs2000.mod */
+
+var m P c e W R k d n l gy_obs gp_obs y dA ;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model (use_dll);
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+k = 6;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam); 
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.196e-3;
+var e_m; stderr 0.0250e-3;
+end;
+
+steady;
+
+stoch_simul(nograph,noprint,order=1);
+
+// Output ghx and ghu using the ordering used by the DLL
+
+[s,i]=sort(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.npred+oo_.dr.nstatic));
+oo_.dr.ghx(oo_.dr.inv_order_var,i)
+
+oo_.dr.ghu(oo_.dr.inv_order_var,:)

mex/sources/estimation/tests/testModelSolution.cc

 /*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -17,8 +17,9 @@
  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-// Test  for ModelSolution
-// Uses fs2000k2.mod and its ..._dynamic.mexw32
+// Test for ModelSolution; to be used with model1.mod
+
+#include <iostream>
 
 #include "ModelSolution.hh"
 
@@ -27,56 +28,50 @@ main(int argc, char **argv)
 {
   if (argc < 2)
     {
-      std::cerr << argv[0] << ": please provide as argument the name of the dynamic DLL generated from fs2000k2.mod (typically fs2000k2_dynamic.mex*)" << std::endl;
+      std::cerr << argv[0] << ": please provide as argument the name of the dynamic DLL generated from model1.mod (typically model1_dynamic.mex*)" << std::endl;
       exit(EXIT_FAILURE);
     }
 
   std::string modName = argv[1];
   const int npar = 7;
-  const size_t n_endo = 15, n_exo = 2;
-  std::vector<size_t> zeta_fwrd_arg;
-  std::vector<size_t> zeta_back_arg;
-  std::vector<size_t> zeta_mixed_arg;
-  std::vector<size_t> zeta_static_arg;
+  const ptrdiff_t n_endo = 15, n_exo = 2;
+  std::vector<ptrdiff_t> zeta_fwrd_arg;
+  std::vector<ptrdiff_t> zeta_back_arg;
+  std::vector<ptrdiff_t> zeta_mixed_arg;
+  std::vector<ptrdiff_t> zeta_static_arg;
   double qz_criterium = 1.0+1.0e-9;
-  Vector deepParams(npar);
 
-  double dYSparams [] = {
-    1.0110,  2.2582,  0.4477,  1.0000,
-    4.5959,  1.0212,  5.8012,  0.8494,
-    0.1872,  0.8604,  1.0030,  1.0080,
-    0.5808,  1.0030,  2.2093
-  };
-  double vcov[] = {
+  Matrix2d vCov;
+  vCov <<
     0.1960e-3, 0.0,
-    0.0, 0.0250e-3
-  };
-  int nVCVpar = 2;
-  MatrixView vCovVW(vcov, nVCVpar, nVCVpar, nVCVpar);
-  Matrix vCov(nVCVpar, nVCVpar);
-  vCov = vCovVW;
+    0.0, 0.0250e-3;
   
-  double dparams[] = { 0.3300,
+  VectorXd deepParams(npar);
+  deepParams << 0.3300,
     0.9900,
     0.0030,
     1.0110,
     0.7000,
     0.7870,
-                       0.0200};
+    0.0200;
+
+  VectorXd steadyState(n_endo);
+  steadyState <<
+    1.0110,  2.2582,  0.4477,  1.0000,
+    4.5959,  1.0212,  5.8012,  0.8494,
+    0.1872,  0.8604,  1.0030,  1.0080,
+    0.5808,  1.0030,  2.2093;
 
-  VectorView modParamsVW(dparams, npar, 1);
-  deepParams = modParamsVW;
-  VectorView steadyState(dYSparams, n_endo, 1);
   std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl;
   std::cout << "Matrix vCov: " << std::endl << vCov << std::endl;
   std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl;
 
   // Set zeta vectors [0:(n-1)] from Matlab indices [1:n]
   //order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
-  size_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
-  size_t back[] = {1, 7, 13};
-  size_t both[] = {2};
-  size_t fwd[] = { 3, 15};
+  ptrdiff_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
+  ptrdiff_t back[] = {1, 7, 13};
+  ptrdiff_t both[] = {2};
+  ptrdiff_t fwd[] = { 3, 15};
   for (int i = 0; i < 9; ++i)
     zeta_static_arg.push_back(statc[i]-1);
   for (int i = 0; i < 3; ++i)
@@ -86,8 +81,8 @@ main(int argc, char **argv)
   for (int i = 0; i < 2; ++i)
     zeta_fwrd_arg.push_back(fwd[i]-1);
 
-  Matrix ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size());
-  Matrix ghu(n_endo, n_exo);
+  MatrixXd ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size());
+  MatrixXd ghu(n_endo, n_exo);
 
   ModelSolution modelSolution(modName, n_endo, n_exo,
                               zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium);

mex/sources/estimation/utils/dynamic_dll.cc

@@ -19,12 +19,12 @@
 
 #include "dynamic_dll.hh"
 
+#include <iostream>
 #include <sstream>
 
 using namespace std;
 
-DynamicModelDLL::DynamicModelDLL(const std::string &dynamicDllFile, size_t n_exog_arg) throw (TSException) :
-  n_exog(n_exog_arg)
+DynamicModelDLL::DynamicModelDLL(const std::string &dynamicDllFile) throw (TSException)
 {
   std::string fName;
 #if !defined(__CYGWIN32__) && !defined(_WIN32)

mex/sources/estimation/utils/dynamic_dll.hh

@@ -27,7 +27,9 @@
 #endif
 
 #include <string>
-#include "Matrix.hh"
+#include <Eigen/Core>
+
+using namespace Eigen;
 
 #include "ts_exception.h"
 
@@ -44,7 +46,6 @@ class DynamicModelDLL
 {
 private:
   DynamicFn Dynamic; // pointer to the Dynamic function in DLL
-  const size_t n_exog; // no of exogenous
 #if defined(_WIN32) || defined(__CYGWIN32__)
   HINSTANCE dynamicHinstance;  // DLL instance pointer in Windows
 #else
@@ -53,17 +54,16 @@ private:
 
 public:
   // construct and load Dynamic model DLL
-  DynamicModelDLL(const std::string &dynamicDllFile, size_t n_exog_arg) throw (TSException);
+  DynamicModelDLL(const std::string &dynamicDllFile) throw (TSException);
   virtual ~DynamicModelDLL();
 
   //! evaluate Dynamic model DLL
-  void eval(const Vector &y, const Matrix &x, const Vector &params, VectorView &ySteady,
-            Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
-  template<class VEC>
-  void eval(const Vector &y, const Matrix &x, const VectorView &modParams, VEC &ySteady,
-                      Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
+  template<class Derived1, class Derived2>
+  void eval(const VectorXd &y, const MatrixXd &x, const PlainObjectBase<Derived1> &modParams,
+            PlainObjectBase<Derived2> &ySteady, VectorXd &residual, MatrixXd *g1, MatrixXd *g2,
+            MatrixXd *g3) throw (TSException)
   {
-    Dynamic(y.getData(), x.getData(), 1, modParams.getData(), ySteady.getData(), 0, residual.getData(),
-	    g1 == NULL ? NULL : g1->getData(), g2 == NULL ? NULL : g2->getData(), g3 == NULL ? NULL : g3->getData());
+    Dynamic(y.data(), x.data(), 1, modParams.data(), ySteady.data(), 0, residual.data(),
+	    g1 == NULL ? NULL : g1->data(), g2 == NULL ? NULL : g2->data(), g3 == NULL ? NULL : g3->data());
   };
 };