Commit 82930ee2 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Implementation of steady state computation in estim DLL, using the GSL

Does not seem very robust (fails on fs2000), need to investigate why.
parent 871cbbe4
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) \
......
......@@ -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
......
......@@ -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
EXEEXT = .mex
include ../mex.am
include ../../estimation.am
logMHMCMCposterior_LDADD = $(LIBADD_MATIO)
......@@ -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
......@@ -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);
......
......@@ -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<size_t> &zeta_static_arg, double qz_criterium);
virtual ~ModelSolution() {};
template <class Vec1, class Vec2, class Mat1, class Mat2>
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 <class Vec1, class Vec2, class Mat1, class Mat2>
......@@ -88,13 +90,6 @@ private:
//compute rules
decisionRules.compute(jacobian, ghx, ghu);
}
template <class Vec1, class Vec2>
void ComputeSteadyState(Vec1 &steadyState, const Vec2 &deepParams)
{
// does nothig for time being.
}
};
#endif // !defined(5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_)
/*
* 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 <http://www.gnu.org/licenses/>.
*/
#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;
}
/*
* 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 <http://www.gnu.org/licenses/>.
*/
#include "Vector.hh"
#include "static_dll.hh"
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
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 <class Vec1, class Mat, class Vec2>
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);
}
};
/*
* 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 <http://www.gnu.org/licenses/>.
*/
#include "static_dll.hh"
#include <sstream>
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
}
/*
* 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 <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 <string>
#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 <model>_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<class Vec1, class Vec2, class Vec3, class Mat1>
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());
};
};
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment