From a4e4034028c582cf271fb97fc362b9548d579e51 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien.villemot@ens.fr>
Date: Mon, 13 Aug 2012 15:01:42 +0200
Subject: [PATCH] Convert ModelSolution to Eigen
---
mex/sources/estimation/ModelSolution.cc | 13 ++--
mex/sources/estimation/ModelSolution.hh | 48 +++++++------
mex/sources/estimation/tests/Makefile.am | 2 +-
mex/sources/estimation/tests/model1.mod | 64 +++++++++++++++++
.../estimation/tests/testModelSolution.cc | 71 +++++++++----------
mex/sources/estimation/utils/dynamic_dll.cc | 4 +-
mex/sources/estimation/utils/dynamic_dll.hh | 20 +++---
7 files changed, 144 insertions(+), 78 deletions(-)
create mode 100644 mex/sources/estimation/tests/model1.mod
diff --git a/mex/sources/estimation/ModelSolution.cc b/mex/sources/estimation/ModelSolution.cc
index 73066a795d..a8cae1c581 100644
--- a/mex/sources/estimation/ModelSolution.cc
+++ b/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(),
diff --git a/mex/sources/estimation/ModelSolution.hh b/mex/sources/estimation/ModelSolution.hh
index 3e690ffa66..2c7e63ff12 100644
--- a/mex/sources/estimation/ModelSolution.hh
+++ b/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.
}
diff --git a/mex/sources/estimation/tests/Makefile.am b/mex/sources/estimation/tests/Makefile.am
index fed293826c..d3f5faab24 100644
--- a/mex/sources/estimation/tests/Makefile.am
+++ b/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
diff --git a/mex/sources/estimation/tests/model1.mod b/mex/sources/estimation/tests/model1.mod
new file mode 100644
index 0000000000..896916fff7
--- /dev/null
+++ b/mex/sources/estimation/tests/model1.mod
@@ -0,0 +1,64 @@
+/* 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,:)
diff --git a/mex/sources/estimation/tests/testModelSolution.cc b/mex/sources/estimation/tests/testModelSolution.cc
index 45b20d8311..8b5efcfe58 100644
--- a/mex/sources/estimation/tests/testModelSolution.cc
+++ b/mex/sources/estimation/tests/testModelSolution.cc
@@ -1,5 +1,5 @@
/*
- * 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 [] = {
+ Matrix2d vCov;
+ vCov <<
+ 0.1960e-3, 0.0,
+ 0.0, 0.0250e-3;
+
+ VectorXd deepParams(npar);
+ deepParams << 0.3300,
+ 0.9900,
+ 0.0030,
+ 1.0110,
+ 0.7000,
+ 0.7870,
+ 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
- };
- double 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;
-
- double dparams[] = { 0.3300,
- 0.9900,
- 0.0030,
- 1.0110,
- 0.7000,
- 0.7870,
- 0.0200};
+ 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);
diff --git a/mex/sources/estimation/utils/dynamic_dll.cc b/mex/sources/estimation/utils/dynamic_dll.cc
index 495e0bae49..17e92c33af 100644
--- a/mex/sources/estimation/utils/dynamic_dll.cc
+++ b/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)
diff --git a/mex/sources/estimation/utils/dynamic_dll.hh b/mex/sources/estimation/utils/dynamic_dll.hh
index dfb4d8b19f..5b7a9a8e1a 100644
--- a/mex/sources/estimation/utils/dynamic_dll.hh
+++ b/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 ¶ms, 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());
};
};
--
GitLab