diff --git a/mex/build/estimation.am b/mex/build/estimation.am index 30c848542..36b1837f4 100644 --- a/mex/build/estimation.am +++ b/mex/build/estimation.am @@ -1,9 +1,9 @@ noinst_PROGRAMS = logposterior logMHMCMCposterior # We use shared flags so that automake does not compile things two times -AM_CPPFLAGS += -I$(top_srcdir)/../../sources/estimation/libmat -I$(top_srcdir)/../../sources/estimation/utils $(CPPFLAGS_MATIO) $(BOOST_CPPFLAGS) -AM_LDFLAGS += $(LDFLAGS_MATIO) $(BOOST_LDFLAGS) -LDADD = $(LIBADD_DLOPEN) $(LIBADD_MATIO) +AM_CPPFLAGS += -I$(top_srcdir)/../../sources/estimation/libmat -I$(top_srcdir)/../../sources/estimation/utils $(CPPFLAGS_MATIO) $(BOOST_CPPFLAGS) $(GSL_CPPFLAGS) +AM_LDFLAGS += $(LDFLAGS_MATIO) $(BOOST_LDFLAGS) $(GSL_LDFLAGS) +LDADD = $(LIBADD_DLOPEN) $(LIBADD_MATIO) $(GSL_LIBS) TOPDIR = $(top_srcdir)/../../sources/estimation @@ -52,8 +52,12 @@ COMMON_SRCS = \ $(TOPDIR)/ModelSolution.hh \ $(TOPDIR)/Prior.cc \ $(TOPDIR)/Prior.hh \ + $(TOPDIR)/SteadyStateSolver.cc \ + $(TOPDIR)/SteadyStateSolver.hh \ $(TOPDIR)/utils/dynamic_dll.cc \ - $(TOPDIR)/utils/dynamic_dll.hh + $(TOPDIR)/utils/dynamic_dll.hh \ + $(TOPDIR)/utils/static_dll.cc \ + $(TOPDIR)/utils/static_dll.hh nodist_logposterior_SOURCES = \ $(COMMON_SRCS) \ diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am index 971ff8ff8..afb9403b1 100644 --- a/mex/build/matlab/Makefile.am +++ b/mex/build/matlab/Makefile.am @@ -3,14 +3,14 @@ ACLOCAL_AMFLAGS = -I ../../../m4 # libdynare++ must come before gensylv, k_order_perturbation, dynare_simul_ if DO_SOMETHING -SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv estimation block_kalman_filter sobol local_state_space_iterations +SUBDIRS = mjdgges kronecker bytecode libdynare++ gensylv block_kalman_filter sobol local_state_space_iterations if HAVE_MATIO SUBDIRS += k_order_perturbation dynare_simul_ endif if HAVE_GSL -SUBDIRS += ms_sbvar +SUBDIRS += ms_sbvar estimation endif if HAVE_SLICOT diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am index 1d39ab697..3b88c5910 100644 --- a/mex/build/octave/Makefile.am +++ b/mex/build/octave/Makefile.am @@ -10,7 +10,7 @@ endif if HAVE_GSL if HAVE_MATIO -SUBDIRS += ms_sbvar +SUBDIRS += ms_sbvar estimation endif endif @@ -18,7 +18,4 @@ if HAVE_SLICOT SUBDIRS += kalman_steady_state endif -if HAVE_MATIO -SUBDIRS += estimation -endif endif diff --git a/mex/build/octave/estimation/Makefile.am b/mex/build/octave/estimation/Makefile.am index bc9c58ab0..8c2e282e7 100644 --- a/mex/build/octave/estimation/Makefile.am +++ b/mex/build/octave/estimation/Makefile.am @@ -1,5 +1,3 @@ EXEEXT = .mex include ../mex.am include ../../estimation.am - -logMHMCMCposterior_LDADD = $(LIBADD_MATIO) diff --git a/mex/sources/estimation/Makefile.am b/mex/sources/estimation/Makefile.am index f1b7ece1d..fd42d7bff 100644 --- a/mex/sources/estimation/Makefile.am +++ b/mex/sources/estimation/Makefile.am @@ -38,6 +38,10 @@ EXTRA_DIST = \ Proposal.cc \ Proposal.hh \ RandomWalkMetropolisHastings.hh \ + SteadyStateSolver.cc \ + SteadyStateSolver.hh \ utils/dynamic_dll.cc \ utils/dynamic_dll.hh \ + utils/static_dll.cc \ + utils/static_dll.hh \ utils/ts_exception.h diff --git a/mex/sources/estimation/ModelSolution.cc b/mex/sources/estimation/ModelSolution.cc index 81450bc48..64543bf09 100644 --- a/mex/sources/estimation/ModelSolution.cc +++ b/mex/sources/estimation/ModelSolution.cc @@ -38,6 +38,7 @@ ModelSolution::ModelSolution(const std::string &basename, size_t n_endo_arg, si 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(basename), + steadyStateSolver(basename, n_endo), llXsteadyState(n_jcols-n_exo) { Mx.setAll(0.0); diff --git a/mex/sources/estimation/ModelSolution.hh b/mex/sources/estimation/ModelSolution.hh index 504690661..92e03a544 100644 --- a/mex/sources/estimation/ModelSolution.hh +++ b/mex/sources/estimation/ModelSolution.hh @@ -27,6 +27,7 @@ #define ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_ #include "DecisionRules.hh" +#include "SteadyStateSolver.hh" #include "dynamic_dll.hh" /** @@ -43,10 +44,10 @@ public: const std::vector &zeta_static_arg, double qz_criterium); virtual ~ModelSolution() {}; template - void compute(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &ghx, Mat2 &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException) + void compute(Vec1 &steadyState, const Vec2 &deepParams, Mat1 &ghx, Mat2 &ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException, SteadyStateSolver::SteadyStateException) { // compute Steady State - ComputeSteadyState(steadyState, deepParams); + steadyStateSolver.compute(steadyState, Mx, deepParams); // then get jacobian and @@ -64,6 +65,7 @@ private: Matrix Mx; DecisionRules decisionRules; DynamicModelDLL dynamicDLLp; + SteadyStateSolver steadyStateSolver; Vector llXsteadyState; //Matrix jacobian; template @@ -88,13 +90,6 @@ private: //compute rules decisionRules.compute(jacobian, ghx, ghu); } - template - void ComputeSteadyState(Vec1 &steadyState, const Vec2 &deepParams) - { - // does nothig for time being. - } - - }; #endif // !defined(5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_) diff --git a/mex/sources/estimation/SteadyStateSolver.cc b/mex/sources/estimation/SteadyStateSolver.cc new file mode 100644 index 000000000..8f9fbf12c --- /dev/null +++ b/mex/sources/estimation/SteadyStateSolver.cc @@ -0,0 +1,78 @@ +/* + * Copyright (C) 2013 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 . + */ + +#include "SteadyStateSolver.hh" + +SteadyStateSolver::SteadyStateSolver(const std::string &basename, size_t n_endo) + : static_dll(basename), residual(n_endo), g1(n_endo) +{ +} + +int +SteadyStateSolver::static_f(const gsl_vector *yy, void *p, gsl_vector *F) +{ + params *pp = (params *) p; + VectorConstView deepParams(pp->deepParams, pp->n_params, 1); + MatrixConstView x(pp->x, 1, pp->n_exo, 1); + + VectorView y(yy->data, yy->size, yy->stride); + VectorView residual(F->data, F->size, F->stride); + + pp->static_dll->eval(y, x, deepParams, residual, NULL, NULL); + + return GSL_SUCCESS; +} + +int +SteadyStateSolver::static_df(const gsl_vector *yy, void *p, gsl_matrix *J) +{ + params *pp = (params *) p; + VectorConstView deepParams(pp->deepParams, pp->n_params, 1); + MatrixConstView x(pp->x, 1, pp->n_exo, 1); + + VectorView y(yy->data, yy->size, yy->stride); + + pp->static_dll->eval(y, x, deepParams, *pp->residual, pp->g1, NULL); + + assert(J->size1 == J->size2 && J->size1 == J->tda); + MatrixView g1t(J->data, J->size1, J->size2, J->tda); + mat::transpose(g1t, *pp->g1); // GSL wants row-major order + + return GSL_SUCCESS; +} + +int +SteadyStateSolver::static_fdf(const gsl_vector *yy, void *p, gsl_vector *F, gsl_matrix *J) +{ + params *pp = (params *) p; + VectorConstView deepParams(pp->deepParams, pp->n_params, 1); + MatrixConstView x(pp->x, 1, pp->n_exo, 1); + + VectorView y(yy->data, yy->size, yy->stride); + VectorView residual(F->data, F->size, F->stride); + + pp->static_dll->eval(y, x, deepParams, residual, pp->g1, NULL); + + assert(J->size1 == J->size2 && J->size1 == J->tda); + MatrixView g1t(J->data, J->size1, J->size2, J->tda); + mat::transpose(g1t, *pp->g1); // GSL wants row-major order + + return GSL_SUCCESS; +} + diff --git a/mex/sources/estimation/SteadyStateSolver.hh b/mex/sources/estimation/SteadyStateSolver.hh new file mode 100644 index 000000000..c3127aa3a --- /dev/null +++ b/mex/sources/estimation/SteadyStateSolver.hh @@ -0,0 +1,102 @@ +/* + * Copyright (C) 2013 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 . + */ + +#include "Vector.hh" +#include "static_dll.hh" + +#include +#include + +class SteadyStateSolver +{ +private: + StaticModelDLL static_dll; + Vector residual; // Will be discarded, only used by df() + Matrix g1; // Temporary buffer for computing transpose + + struct params + { + StaticModelDLL *static_dll; + const double *deepParams; + size_t n_params; + const double *x; + size_t n_exo; + Vector *residual; + Matrix *g1; + }; + + static int static_f(const gsl_vector *yy, void *p, gsl_vector *F); + static int static_df(const gsl_vector *yy, void *p, gsl_matrix *J); + static int static_fdf(const gsl_vector *yy, void *p, gsl_vector *F, gsl_matrix *J); + + const static double tolerance = 1e-5; + const static size_t max_iterations = 1000; +public: + class SteadyStateException + { + }; + + SteadyStateSolver(const std::string &basename, size_t n_endo); + + template + void compute(Vec1 &steadyState, const Mat Mx, const Vec2 &deepParams) throw (SteadyStateException) + { + assert(steadyState.getStride() == 1); + assert(deepParams.getStride() == 1); + + const size_t n = steadyState.getSize(); + + gsl_vector_view ss = gsl_vector_view_array(steadyState.getData(), n); + + params p = { &static_dll, deepParams.getData(), deepParams.getSize(), Mx.getData(), Mx.getCols(), &residual, &g1 }; + + gsl_multiroot_function_fdf f = {&static_f, &static_df, &static_fdf, + n, &p}; + + const gsl_multiroot_fdfsolver_type *T = gsl_multiroot_fdfsolver_gnewton; + gsl_multiroot_fdfsolver *s = gsl_multiroot_fdfsolver_alloc(T, n); + + gsl_multiroot_fdfsolver_set(s, &f, &ss.vector); + + int status; + size_t iter = 0; + + do + { + iter++; + + status = gsl_multiroot_fdfsolver_iterate(s); + + if (status) + break; + + status = gsl_multiroot_test_residual(s->f, tolerance); + } + while(status == GSL_CONTINUE && iter < max_iterations); + + if (status != GSL_SUCCESS) + throw SteadyStateException(); + + gsl_vector_memcpy(&ss.vector, gsl_multiroot_fdfsolver_root(s)); + + gsl_multiroot_fdfsolver_free(s); + } +}; + + diff --git a/mex/sources/estimation/utils/static_dll.cc b/mex/sources/estimation/utils/static_dll.cc new file mode 100644 index 000000000..fe3366398 --- /dev/null +++ b/mex/sources/estimation/utils/static_dll.cc @@ -0,0 +1,90 @@ +/* + * Copyright (C) 2010-2013 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 . + */ + +#include "static_dll.hh" + +#include + +using namespace std; + +StaticModelDLL::StaticModelDLL(const std::string &basename) throw (TSException) +{ + std::string fName; +#if !defined(__CYGWIN32__) && !defined(_WIN32) + if (basename[0] != '/') + fName = "./"; +#endif + fName += basename + "_static" + MEXEXT; + + try + { +#if defined(__CYGWIN32__) || defined(_WIN32) + staticHinstance = LoadLibrary(fName.c_str()); + if (staticHinstance == NULL) + throw 1; + Static = (StaticFn) GetProcAddress(staticHinstance, "Static"); + if (Static == NULL) + { + FreeLibrary(staticHinstance); // Free the library + throw 2; + } +#else // Linux or Mac + staticHinstance = dlopen(fName.c_str(), RTLD_NOW); + if ((staticHinstance == NULL) || dlerror()) + { + cerr << dlerror() << endl; + throw 1; + } + Static = (StaticFn) dlsym(staticHinstance, "Static"); + if ((Static == NULL) || dlerror()) + { + dlclose(staticHinstance); // Free the library + cerr << dlerror() << endl; + throw 2; + } +#endif + + } + catch (int i) + { + std::ostringstream msg; + msg << "Error when loading " << fName << " ("; + if (i == 1) + msg << "can't dynamically load the file"; + if (i == 2) + msg << "can't locate the 'Static' symbol"; + msg << ")"; + throw TSException(__FILE__, __LINE__, msg.str()); + } + catch (...) + { + throw TSException(__FILE__, __LINE__, std::string("Can't find Static function in ") + fName); + } +} + +StaticModelDLL::~StaticModelDLL() +{ +#if defined(__CYGWIN32__) || defined(_WIN32) + bool result = FreeLibrary(staticHinstance); + if (result == 0) + throw TSException(__FILE__, __LINE__, std::string("Can't free the *_static DLL")); +#else + dlclose(staticHinstance); +#endif +} diff --git a/mex/sources/estimation/utils/static_dll.hh b/mex/sources/estimation/utils/static_dll.hh new file mode 100644 index 000000000..77f183890 --- /dev/null +++ b/mex/sources/estimation/utils/static_dll.hh @@ -0,0 +1,72 @@ +/* + * Copyright (C) 2010-2013 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 . + */ + +#if defined(_WIN32) || defined(__CYGWIN32__) +# ifndef NOMINMAX +# define NOMINMAX // Do not define "min" and "max" macros +# endif +# include +#else +# include // unix/linux DLL (.so) handling routines +#endif + +#include +#include "Matrix.hh" + +#include "ts_exception.h" + +// Pointer to the Static function in the MEX +typedef void (*StaticFn) +(const double *y, const double *x, int nb_row_x, const double *params, double *residual, double *g1, double *v2); + +/** + * creates pointer to Dynamic function inside _static.dll + * and handles calls to it. + **/ +class StaticModelDLL +{ +private: + StaticFn Static; // pointer to the Dynamic function in DLL +#if defined(_WIN32) || defined(__CYGWIN32__) + HINSTANCE staticHinstance; // DLL instance pointer in Windows +#else + void *staticHinstance; // and in Linux or Mac +#endif + +public: + // construct and load Static model DLL + StaticModelDLL(const std::string &basename) throw (TSException); + virtual ~StaticModelDLL(); + + //! evaluate Static model DLL + template + void eval(const Vec1 &y, const Mat1 &x, const Vec2 &modParams, + Vec3 &residual, Matrix *g1, Matrix *v2) throw (TSException) + { + assert(y.getStride() == 1); + assert(x.getLd() == x.getRows()); + assert(modParams.getStride() == 1); + assert(residual.getStride() == 1); + assert(g1->getLd() == g1->getRows()); + assert(v2->getLd() == v2->getRows()); + + Static(y.getData(), x.getData(), 1, modParams.getData(), residual.getData(), + g1 == NULL ? NULL : g1->getData(), v2 == NULL ? NULL : v2->getData()); + }; +};