Jacobian for system of ODEs


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.
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),

  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());

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_;


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.


I doubt this was ever an officially supported feature :P, but I’d bet you can get this hacked back into a working situation by:

  1. Write out your full model and put in a placeholder ODE function (just make it compile-able – don’t worry about putting an actual ODE in it)

  2. Feed the model to stanc. That’ll produce the cpp version of your model that would have been compiled into a binary had you used the regular make model thing

  3. Find the section of code that looks like the stuff you already have and figure out what changed

Hope that helps :D. Custom Jacobians are definitely on my Stan list of things I want (think some other folks want them as well too), but it’s still at the wishful thinking stage.


I am guilty of having mentioned this “feature” which was never supported (and I never said something else).

In principle things should still work for 2.17.1… but for the forthcoming 2.18 the approach of passing in a ode_system template specialization is going to fail as we refactored that away (giving 10-20% speed increase).

It’s not too hard to adapt to 2.18 this approach, but honestly… if you can: Just use a few more CPUs and solve you speed problem in a brute force manner. Threading or MPI is in develop (MPI is almost complete).

(but, of course, I would really like to see this custom Jacobian support in Stan sooner than later, but it’s hard)


Thanks for the suggestions. I wanted to compare with variational inference which didn’t work on 2.14. My system solves 5 days with jacobian supplied and I tried without - didn’t have patience to wait.

I forgot to mention - the system is published in book chapter - one more stan publication: Stamatis, S.D., Mockus, L., Kirsch, L.E., Reklaitis, G.V., “Bayesian Hierarchical modeling of Gabapentin absorption and disposition with application to dosing regimen assessment” in Quantitative Systems Pharmacology - Models and Model-based systems with applications, Comp Aid. Chem. Eng Vol 42, D. Manca (ed), Elsevier, (in press, 2018)

The fit I am doing right now is for another book chapter - yet another stan success story in PBPK domain.

In reality making it work in 2.17 should be fairly simple - I just lack required C++ knowledge.

stan/lib/stan_math/stan/math/rev/mat/functor/cvodes_ode_data.hpp:96:9: error: no match for call to ‘(const stan::math::ode_system<pbpkautojac_model_namespace::pbpk_functor__>) (double&, const std::vector&, Eigen::Map<Eigen::Matrix<double, -1, 1>, 0, Eigen::Stride<0, 0> >&)’
ode_system_(t, y_vec, dy_dt_eig);

My code has ode_system(const F& f,
const std::vector theta,
const std::vector& x,
const std::vector& x_int,
std::ostream* msgs)

and inline void operator()(const double t,
const std::vector& y,
std::vector& dy_dt

what resembles pretty well the error.

Thanks for any help, unfortunately with undocumented feature of stan, hopefully it is not too complicated.


Could you attach the model, the c++ file, and the command you usually use to compile it?

It doesn’t have to be the full model if you don’t want to. Just enough that would call the external stuff. I can have a quick look tomorrow.

This sounds like a non-identified parameter or something weird. How many parameters?


Terrific. I am re-installing cmdstan 2.14 - but your help is much
appreciated. In 2.14 I was unable to run variational inference on my system
of ODEs - there are 33 of them, around 10 parameters and run time 193938
seconds (Total), 200 warm-ups and 1000 sampling iterations.

The stan model is attached. I compile it with …/bin/stanc pbpkautojac.stan
–o=pbpkautojac.hpp using stanc from cmdstan 2.14. Then I add line #include
“pbpkauto_jacobianCmdStan.hpp” at the and of pbpkautojac.hpp.The final file
is attached. All required files are also attached.

Evidently something changed from 2.14 so the problem is in
pbpkauto_jacobianCmdStan.hpp. I believe it will be quick fix for an expert.
Many thanks for the help.

pbpkautojac.stan (13.1 KB)


Whoops, my bad, this was trickier than I thought :D. Can you attach the file “pbpkauto_jacobianCmdStan.hpp”? I think I’ll need to actually compile this and see the error and the formatting is all messed up in the first post (and I don’t wanna try to fix it).


Sorry to cause so much trouble. I really appreciate your help. Working with 2.14 is bit awkward.


I think you may have forgotten to attach the .hpp :D. Either that or Discourse is hiding it


I think the extension makes google suspect virus. I renamed it with .x. Thanks a lot for your help.


I think the extension makes google suspect virus. I renamed it with .txt. Thanks a lot for your help.
pbpkautojac.stan (13.1 KB)
states.txt (24.2 KB)
params.txt (9.3 KB)
pbpkauto_jacobianCmdStan.txt (2.5 KB)


Just FYI, I enabled hpp and cpp as accepted extensions for uploads.


Thank you.


Take 1! I didn’t investigate the reasons why, but if you make these changes (replacing lines preceded by “<” with the stuff preceded by “>”) then hopefully the model will compile:

<                              std::vector<double>& dy_dt) const {
< 	dy_dt = f_(t, y, theta_, x_, x_int_, msgs_);
> 			     Eigen::Map<Eigen::Matrix<double, -1, 1>>& dy_dt) const {
> 	std::vector<double> dy_dt_ = f_(t, y, theta_, x_, x_int_, msgs_);
> 	for(size_t i = 0; i < dy_dt_.size(); i++)
> 	  dy_dt(i) = dy_dt_[i];

Lemme know if it goes anywhere!

edit: Haha, man, I hope dt_dt is allocated somewhere. I just assume it is since it’s an Eigen::Map


@syclik I have a copy of @linas’s model:

jac.hpp (2.6 KB)
pbpkautojac.stan (12.9 KB)
params.txt (9.3 KB)
states.txt (24.2 KB)

I should be able to compile this in cmdstan with:

make pbpkautojac right?

I’m getting the error:

--- Translating Stan model to C++ code ---
bin/stanc  examples/pbpkautojac.stan --o=examples/pbpkautojac.hpp
Model name=pbpkautojac_model
Input file=examples/pbpkautojac.stan
Output file=examples/pbpkautojac.hpp

could not find include file
make/models:24: recipe for target 'examples/pbpkautojac.hpp' failed
make: *** [examples/pbpkautojac.hpp] Error 255

Is there something else I should be doing to use #include ... w/ cmdstan?

(To do the testing for my previous post I just copy-pasted the g++ command line and added -include lines manually)


Thanks a lot. It seems it worked. One of the chains stopped immediately with an error:
Unrecoverable error evaluating the log probability at the initial value.
ran beyond end of program in trace()
ran beyond end of program in trace()

but other 3 chains are still running. Hopefully the error is not related to our hacking the jacobian.

I don’t know if it is related to the missing include file but to make it easier before calling make I do the following:
cd c:/Tools/cmdstan-2.17.1/tmp
…/bin/stanc pbpkautojac.stan --o=pbpkautojac.hpp
add line #include "pbpkauto_jacobianCmdStan.txt\ to the end of pbpkautojac.hpp
cd …
make tmp/pbpkautojac.exe


The diagnose mode of cmdstan should tell you quickly if your Jacobian is ok. This compares finite differencing vs autodiff. However, have a look at the output from cmdstan from 2.14 (where you know things are ok) in order to judge the output from 2.17 (or better: compile the model without the Jacobian thing with 2.17 as well). I have seen noticeable differences for ODEs for cases when it is correct, but when your Jacobian thing did not work correctly it should become apparent.


Thanks. At least in 2.14 I had trouble comparing autodiff with numerical. Every time I compared autodiff-numerical for some parameter was different and not in 0.0000001 but in 0.1. Eventually I arrived to autodiff=numerical. My suspicion that it tries to compare at different initial values. I guess I have just to set inits?

There are two jacobians - one for parameters, another one for ODE state variables. Are both jacobians being used by cmdstan?


That depends on your problem, but very likely the answer is yes. Only if you do not vary the parameters of the model, then only the derivatives wrt to the states are needed. If parameters vary (as usual), then both (states+parameter) Jacobians are used.