Jacobian for system of ODEs

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)

1 Like

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

1 Like

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:

39,40c39,42
<                              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.

1 Like

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.

Will you be able to help to incorporate jacobian of system of ODEs for cmdstan 2.18? The posterior estimation with the system of ODEs I am working with takes quite a bit so in cmdstan 2.14 and 2.17 I have provided jacobians which sped up estimation quite a bit. My guess providing jacobian for system of ODEs will speed up things in 2.18 as compared with numerical jacobian estimation.

The problem I am working with is individualized dosing. Each subject is different so for some 100 mg IV once a day is within therapeutic window while some need 83mg once a day. Even some physicians are interested. For compartmental models this is not so complicated but for physiologically based PK models (described by a large system of ODEs) this is quite challenging. On previous versions of cmd stan we did some study with gabapentin and it worked quite well although was quite computationally intensive (book chapter is published). However, variational inference didn’t work well in those versions. Idea is to use gabapentin data to show that variational inference produces similar results to MCMC and we hope that variational inference will make calculations realistic. I think this may be a hook to decision support system.

I may have some code ready somewhere which does what you want… but apart from that it sounds that you are solving the large ODE system for multiple patients, right? If that is the case, then you should by all means use the new map_rect function and parallelize the evaluation of the log-likelihood. It is a bit of a pain to do the packing and unpacking of the data, yes, but in your case you will get almost for sure a linear speedup of your computation time with the number of cores per chain.

Yes, hierarchical model for multiple patients. If there is a way to supply jacobian for ODEs in cmdstan 2.18 I would be grateful. I would like to use map_rect - maybe you could share some example?

In general, to speed up things I was thinking about approximated inference (variational inference) and transforming ODEs to algebraic equations (via orthogonal collocation or some other techniques, perhaps you could suggest). Since stan already has AE solver it would be great to supply jacobian for system of AEs.

In reality I think it would be nice to expose C++ interface for folks who love supplying jacobians for ODEs or AEs. It is OK that it is hard to use - but sometimes one has to bite a bullet.

The manual has an example and there is a StanCon contribution forthcoming from my end which will demonstrate things. For your models you really should use map_rect. These type of problems are ideal for threading or MPI and I would expect that you get linear performance increases with the number of cores up to the point where you have 1 core per patient.

The need to pass in the Jacobian will be much reduced (and BTW, the stiff ODE solver got ~15% faster in 2.18 in examples I tried).

Looking forward for your StanCon contribution while Bob’s logistic regression example is quite self explanatory (http://discourse.mc-stan.org/t/map-rect-in-the-stan-language-mpi-is-coming/3547/2). I still suspect that even with map_rect supplying jacobian for ODEs would reduce time even more (simply because calculating jacobian numerically requires multiple solutions of ODEs). 15% from 5 days is significant but I would like to push the boundaries.

Of course, supplying the Jacobian will help on top of the parallelization… but I can tell you that with brute-force 80 cores I was able to get a 2.7 days runtime down to 1h for an ODE problem. At that point I don’t bother with Jacobians any more… but this was running on a high-performance computing cluster not all of us have available.

Cool, I don’t have such luxury. BTW, if interested I can share model/data to test on the cluster.