In stanCmd 2.14 there was a possibility to supply jacobian for system of ODEs. Unfortunately in 2.17.1 the same code didn’t work.
Thanks for the help. My guess is that naming conventions/parameters have changed when going from 2.14 to 2.17. Is it possible to get some working example?
The code is below - I wasn’t able to attach hpp files.
Essentially I generate hpp file and at the end insert a line
#include “pbpkauto_jacobianCmdStan.hpp”
The content of this file is:
#pragma once
#include <stan/math/rev/mat/functor/ode_system.hpp>
// define an analytic jacobian for the hornberg for Stan’s CVODES
// solver by partial template specialisation
namespace stan {
namespace math {
// todo: idealy, we derive the specialised ode_model from the
// general one which gives possibly syntax problem which is why we
// chose to use a typedef, maybe we need another helper... not
// sure.
template<>
struct ode_system<pbpkautojac_model_namespace::pbpk_functor__> {
typedef pbpkautojac_model_namespace::pbpk_functor__ F;
const F& f_;
const std::vector<double> theta_;
const std::vector<double>& x_;
const std::vector<int>& x_int_;
std::ostream* msgs_;
ode_system(const F& f,
const std::vector<double> theta,
const std::vector<double>& x,
const std::vector<int>& x_int,
std::ostream* msgs)
: f_(f),
theta_(theta),
x_(x),
x_int_(x_int),
msgs_(msgs)
{}
inline void operator()(const double t,
const std::vector<double>& y,
std::vector<double>& dy_dt) const {
dy_dt = f_(t, y, theta_, x_, x_int_, msgs_);
}
template <typename Derived1, typename Derived2>
inline void
jacobian(const double t,
const std::vector<double>& y,
Eigen::MatrixBase<Derived1>& fy,
Eigen::MatrixBase<Derived2>& Jy
) const {
using Eigen::VectorXd;
using std::vector;
using std::pow;
const vector<double> f = f_(t, y, theta_, x_, x_int_, msgs_);
fy = VectorXd::Map(&f[0], f.size());
Jy.setZero();
Eigen::Map<Eigen::VectorXd> J = Eigen::Map<Eigen::VectorXd>(&Jy(0,0), Jy.cols()*Jy.rows());
// dy/dt=f(y,theta)
// df/dy
#include “states.hpp”
}
template <typename Derived1, typename Derived2>
inline void
jacobian(const double t,
const std::vector<double>& y,
Eigen::MatrixBase<Derived1>& fy,
Eigen::MatrixBase<Derived2>& Jy,
Eigen::MatrixBase<Derived2>& Jtheta
) const {
using Eigen::VectorXd;
using std::vector;
using std::pow;
jacobian(t, y, fy, Jy);
//const vector<double>& parms = theta_;
Jtheta.setZero();
Eigen::Map<Eigen::VectorXd> J = Eigen::Map<Eigen::VectorXd>(&Jtheta(0,0), Jtheta.cols()*Jtheta.rows());
// df/dtheta
#include “params.hpp”
}
};
}
}
where params.hpp and states.hpp are actual jacobians.