Putting draft pharmacometric functions on the Stan Github
Charles Margossian
The original post was in an e-mail chain, but here is a summary, and answer to some of the questions that have been asked.
My team is writing new functions for Stan with applications in Pharmacometrics. These functions can be called directly inside a stan file (so there is no dependence on whether the user uses cmdStan, PyStan, etc.)
We call this function library Torsten, after Torsten Teorell, who is regarded as the father of Pharmacokinetics.
The current code is compatible with stan/math-2.10.0.
The functions are coded in C++ and would be located in a directory under stan/math. So I expect to upload them in the stan-dev/math branch. The major advantage in doing so is that it will be easier to work on the code and get help on it, since it will be accessible to the community.
To expose some of these functions, I modified rev/mat.hpp and function_signatures.h. Other functions however are higher-order functions (to the same effect as integrate_ode_bdf and integrate_ode_rk45), and cannot simply be exposed in function_signatures.h. So I edited a few grammar files for the sole purpose of exposing these higher-order functions. I expect to share these modifications using the stan-dev/stan branch.
These higher-order functions have 14 arguments, which means we will need to increase the limit on the number of arguments (which is currently 12).
I still have to write unit tests, consistent with Stan practices. (I wrote some examples and tests using cmdStan and R, for development purposes but these are likely not consistent with Stan's practices. In due time, I will upload the examples to the stan/examples repo).
I agree with Bob's suggestion: we can start with a wiki and I am happy to attend a stan meeting to discuss design.
I hope this clarifies matter, and I'm happy to answer any other question.
Bob Carpenter
going on. The new name confused me.
At this point, I think you only need to create issues for
what you're going to add explaining what the function
signatures and names are going to be.
You also need the fwd implementations to go with the rev
ones for any functions we expose in Stan (other than the
ODEs). At least up until now, we don't have any way to
signal which functions are available at which order. It's
easy enough to write templated C++ to deal with forward-mode,
but it's harder if the arguments are polymorphic or if
there are a lot of arguments for a function.
Then it'll just be stan-dev/math and stan-dev/stan pull requests.
These shouldn't be made until a complete set of tests is
in place.
- Bob
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
Sebastian Weber
In fact, I intended to start such a wiki on how to make Stan better/easier to use for PK/PD.
From my perspective we need a discussion on the design of things. I myself have come up with large Stan programs which can exactly do what I want - this is in my case to take a Nonmem data set and fit PK/PD models.
I intentionally kept things to be programmed in the Stan language only as I think that there should be no Nonmem (or PK/PD) specifics in Stan C++ and I am also not sure if we want to have PK model specifics in Stan (unless there are good arguments for it, of course).
From my experience in programming these models I figured it would probably make sense to include an event handling mechanism in Stan. In PK models it is crucial to code the dosing history of a given patient. To me it would be natural to include into Stan a facility (functions and possibly data-structures) which can handle these events for the purpose of making it easier to code and to speedup things. Such an event mechanism will probably be very useful for other domains as well aside from PK/PD models and the event handler may not be limited to ODE models, but should be able to cope with analytic models as well.
To move forward, we should definitely discuss during a Stan meeting. I guess we also need a wiki and maybe a scratch place to keep Thorsten code around and maybe my little Stan only solution. From there we should decide how to design things, what makes sense to include into Stan itself and what should live on the Stan language side.
Sebastian
On Thursday, July 7, 2016 at 3:14:31 AM UTC+2, Bob Carpenter wrote:
> Thanks for moving to the list and clearing up what's
> going on. The new name confused me.
>
> At this point, I think you only need to create issues for
> what you're going to add explaining what the function
> signatures and names are going to be.
>
> You also need the fwd implementations to go with the rev
> ones for any functions we expose in Stan (other than the
> ODEs). At least up until now, we don't have any way to
> signal which functions are available at which order. It's
> easy enough to write templated C++ to deal with forward-mode,
> but it's harder if the arguments are polymorphic or if
> there are a lot of arguments for a function.
>
> Then it'll just be stan-dev/math and stan-dev/stan pull requests.
> These shouldn't be made until a complete set of tests is
> in place.
>
> - Bob
>
>
> > On Jul 6, 2016, at 5:38 PM, Charles Margossian
> >
> >
Charles Margossian
I'll set up a wiki so we can pool ideas, and will create a gitHub repo to post Thorsten. Once I have created the unit tests and dealt with forward auto-diff, I'll proceed to add the code to stan-dev.
Some quick thoughts:
The main benefit of using built-in stan functions is that the user does not need to write extensive code in the function block. I also imagine it would be pretty involved (if at all possible?) to write in Stan some of the code I did in C++ (I used classes and functors a lot...). Built-in functions make the coding more efficient and very easy, which I think will be key to interest users.
All the functions take in NONMEM data and return amounts in each compartment at each event. There are roughly speaking 4 functions. Here are two examples:
PKTwoCptModel() which analytically computes the solution to a two compartment model with first order absorption.
genCptModel_bdf() which computes the solution to an ODE system specified by a user using the integrate_ode_bdf function.
These functions all take advantage of exactly the same event handling mechanism, which augments the data and then sequentially predicts the amounts one event at a time. What changes from one function to the other is the function that predicts the amount at an individual event.
Right now the event mechanism is designed for NONMEM data, but I agree it would be great to generalize it. Maybe we could create a stan function that let's the user specify the one-event predictive function in the function block and then runs it sequentially through all the events.
Charles
Sebastian Weber
A huge fraction of that code is completely agnostic to the specific model. In fact, the model function (be it an ODE or not) most be able to calculate the solutions for a given state of the compartments, the parameters and the initial time.
I was able to handle a 2-cmt oral dosing situation with time-changing clearance. It's tedious work to put that together and I guess things run faster if the event handling would be done by Stan.
For Stan I think we need a data-structure to represent events and a function to "execute" these events given a model. This would be very useful for many other domains as well and as such should not be Nonmem specific. Having such a data-structure would certainly allow to quickly convert Nonmem data to this structure. This is my view on this - and yes we should certainly pool our ideas on the wiki.
BTW, this event-thing could be used to include an event like "new-patient" or better "new grouping unit". This could then potentially in the long run enable a parallelized evaluation of the likelihood (given we figure out how to extend the nested AD to such a case which is doable from my understanding).
Sebastian
Andrew Gelman
A
Devin Pastoor
Have you had an interaction with Wenping Wang and his work on PHXStan (also at Novartis but I appreciate a large organization). I was able to get ahold of the source code, (previously mentioned by Andrew http://andrewgelman.com/wp-content/uploads/2015/10/poster_PMXstan_ACoP2015_30Sep2015.pdf) however since it predated much of the ODE work they also implemented their own custom implementations, as well as some analytical solutions for the traditional 1cmpt/2 cmpt models.
Right now they want to keep the source "closed" until it gets cleaned up, however I am encouraging them to open it up to everyone, maybe this effort will be enough to throw that code into the ring as well.
I am finally getting through a swath of other projects (including my thesis work :-p) and will have time to start working this (and some cloudstan ideas) and am happy to contribute in some capacity.
Playing the devil's advocate a little, in regards to incorporating directly in stan/math, unless the functions are generic (in application) I wonder as to the potentially slippery slope of industry-specific code being "dumped" into that repo. Pharmacometrics alone has so many somewhat niche (at least in terms of design) needs in terms of things like events I wonder if some sort of independent module would be more appropriate and provide an example path forward for other industries interested in substantial domain-specific code.
Charles, have you interacted much with Kyle regarding how he is managing events and such in mrgsolve. As a DSL, mrgsolve has been very pleasant for simulation and could be a nice to have a reasonably aligned API. https://github.com/mrgsolve and docs https://github.com/metrumresearchgroup/mrgsolve/wiki/gallery
Bill Gillespie
What is missing is some convenient way to distribute and install domain-specific functionality, e.g., some sort of package distribution mechanism ala CRAN or even just being able to include external files containing Stan language functions. In the absence of that, I prefer to have the PMX functions incorporated into Stan so they are easily distributed as part of the standard Stan distribution.
Cheers,
Bill
Bob Carpenter
code. It's a real tradeoff. By including it directly in stan::math,
everything stays compatible because we agree to maintain it. And it's
easy for everyone to install and access. This is convenient, but
a heavy long-term commitment. And it means everything just takes longer
in stan::lang and stan::math, like unit testing.
This one's a bit trickier because it involves modifying the parser.
I don't see how to make that extensible. Just a bunch of function
libraries would be easy to include, but parser code is pretty entangled
with the rest of the Stan language because of the recursion
through expressions (or statements).
Hmm. I'll add it as a meeting discussion topic.
- Bob
Bob Carpenter
should be pushed directly to the develop branch without a pull request and
code review.
- Bob
Bob Carpenter
by dividing the data. We'd need to build up a TensorFlow-like interface
that allows us to specify comptuation units in terms of input arguments
that returns a value with partials that can be slotted back into the global
autodiff graph. The trick will be figuring out how to deal with memory
efficiently for subcomputations.
I'm going to the autodiff workshop in September---maybe they'll have
some good ideas about parallelization techniques.
- Bob
Bob Carpenter
deals with this autogenerated code. I'd prefer to have
functions that we can call if that's at all possible.
- Bob
Bob Carpenter
> On Jul 7, 2016, at 1:31 PM, Bill Gillespie <bi...@metrumrg.com> wrote:
>
> While we (Charles and I) have built some prototype pharmacometrics (PMX) functions in C++, I'm pretty agnostic about the Stan language vs C++ issue. I agree that certain core functionality should be part of Stan math like numerical solvers for systems of ODEs and algebraic equations and some sort of event handling like that suggested by Sebastian. The more PMX-specific functions could be either C++ or Stan language. I look forward to further discussion about which approach may be optimal or most consistent with Stan design philosophy.
add more users and more developers and more interfaces and add-ons.
Dependencies are the bane of all computational work.
> What is missing is some convenient way to distribute and install domain-specific functionality, e.g., some sort of package distribution mechanism ala CRAN or even just being able to include external files containing Stan language functions. In the absence of that, I prefer to have the PMX functions incorporated into Stan so they are easily distributed as part of the standard Stan distribution.
way to extend Stan.
Unless you code analytic gradients or want to use fancy C++
features like template metaprogramming, it's just as efficient
to code in Stan as in C++. But we need a general way to extend
Stan programs, too. There's a placeholder in R, but we want
to make a general include mechanism.
And double exactly on making it easy on users to install. The
more hurdles users have to jump through, the less likely they
will be to install or upgrade.
- Bob
Sebastian Weber
Do you have in mind a simple include facility for Stan programs and C++? Or are you even opting for a well-defined API on the C++ side which is intended to stay stable over time? This would be very ambitious, I think.
If I get your words right, then my Stan only solution should performance wise be very similar to Thorsten. Would be interesting to benchmark this.
My main thought of going a Stan only solution was simple - you made a commitment to the Stan language and as such I can be sure that my programs will keep working for the foreseeable time. Anything under the hood appears to be under constant refactoring which is how a good software evolves over time.
Sebastian
On Thursday, July 7, 2016 at 8:31:33 PM UTC+2, Bob Carpenter wrote:
> > On Jul 7, 2016, at 1:31 PM, Bill Gillespie
> >
Bob Carpenter
> On Jul 8, 2016, at 2:48 AM, Sebastian Weber <sdw....@gmail.com> wrote:
>
> Hi Bob!
>
> Do you have in mind a simple include facility for Stan programs and C++? Or are you even opting for a well-defined API on the C++ side which is intended to stay stable over time? This would be very ambitious, I think.
keep the include functionality simple. I don't want to
write an entire preprocessor.
> If I get your words right, then my Stan only solution should performance wise be very similar to Thorsten. Would be interesting to benchmark this.
comes down to memory locality and reusing shared computations
efficiently.
> My main thought of going a Stan only solution was simple - you made a commitment to the Stan language and as such I can be sure that my programs will keep working for the foreseeable time. Anything under the hood appears to be under constant refactoring which is how a good software evolves over time.
I think the math library's structure will be stable longer term
and so will stan-dev/stan after Daniel finishes the algorithm/command
refactor and we clean up the model class interface some more.
- Bob
Charles Margossian
Put the code on a branch on stan-dev/stan and stan-dev/math. Nothing
should be pushed directly to the develop branch without a pull request and
code review.
I intended to start such a wiki on how to make Stan better/easier to use for PK/PD.
Charles, have you interacted much with Kyle regarding how he is managing events and such in mrgsolve.
- Bob
Ben Goodrich
For Stan I think we need a data-structure to represent events and a function to "execute" these events given a model. This would be very useful for many other domains as well and as such should not be Nonmem specific. Having such a data-structure would certainly allow to quickly convert Nonmem data to this structure. This is my view on this - and yes we should certainly pool our ideas on the wiki.
I don't remember the NONMEM data structure but I think a lot of problems like these could be solved if the Stan language had a tuple object. Then we could probably "typedef" a NONMEM structure as
vector[N] tuple{int, int, real, ...}
or whatever types the NONMEM fields are in reality.
Ben
Sebastian Weber
+1
YES, a tuple object would have made my life in coding this stuff a lot easier!
Sebastian
Sebastian Weber
Thanks for sharing Torsten. After a first quick look I have the impression that what you did is conceptually very close to what I came up myself. My language to describe things is a bit different (more oriented at the math behind it probably), but the design is the same. I am wondering if we should make the wiki around "complex ODE" models instead of being domain specific on Pharmacometrics. While Pharmacometrics has special needs, these should probably be phrased in more general terms such that other disciplines benefit as well. Given that stan is a general purpose software this would make it a better fit. Another domain where events are used is systems biology and I am sure there are more disciplines.
Is there an example program I can try with Torsten?
One more comment: Out of my own experience in contributing CVODES to the Stan code base, I would advise you to think at the moment about tests, but hold off for the moment to implement all of them. I am saying this as I have myself written lots of tests before there was a definite agreement on design and before I knew where to put things within the Stan source system. As a result I had to rewrite my tests a few times until I learned where to put stuff and until the design was settled. There were many design aspects which I was not aware of when starting this. This process was at times a bit frustrating for me (no blame on anyone here, I think we are just learning this process), so keep this in mind.
Best,
Sebastian
Devin Pastoor
For example, I wanted to open a question about the currently implementation of one/two compartment models is very clearly coded around oral dosing, with CL/V/K. For example, https://github.com/charlesm93/math/blob/feature/issue-314-torsten-pharmacometrics/stan/math/torsten/PKModel/PKModel_class.hpp#L54 shows onecmpt having two compartments, which for oral given the depot compartment makes complete sense, however when managing IV you don't need that depot. This goes in the other direction, like if you have multiple depot compartments for complex absorption models
Do you have any spec docs or otherwise to expand on the thought process behind some of these implementation details, as well as get a better understanding of the overall implementation without tracing through the code itself. For example, I was trying to find some more detail about things like handling parameterizations (physiological CL/Ka/V vs rate constants).
Really excited about this implementation - great work so far!
Charles Margossian
I am wondering if we should make the wiki around "complex ODE" models instead of being domain specific on Pharmacometrics.
Is there an example program I can try with Torsten?
I would advise you to think at the moment about tests, but hold off for the moment to implement all of them.
Charles Margossian
Is the listserve the best place to discuss implementation details (I don't think so but don't know where?)
This goes in the other direction, like if you have multiple depot compartments for complex absorption models