Complex ODE Based Models wiki: suggestions

[Daniel Lee]

Charles, Bill, Sebastian, Andrew, and Bob: thanks for getting on a video call last week to help facilitate my understanding of what you're wanting to do with ODEs. It really helped me understand what's going on.

I didn't want to edit the page without some agreement, so I'll start with edits here. Please feel free to address all, some, or none of it.

Could we rename the wiki page? When I google "complex ode," it leads me to a wiki page on odes where the states are complex variables. I don't think that's what this page is about. Perhaps something like "Advanced techniques for models involving ODEs"?

One additional remark about renaming the wiki page: it's not focused on the complexity of ODEs at all. It abstracts that away. It's focused on techniques once you have a "solution operator." In the cases you're envisioning, the Stan program may not even specify an ODE. The dyanmic system may indeed have an ODE, but that may have been analytically solved so the representation is no longer an ODE from within the Stan language. (Honestly, this was my confusion going into the discussion and when I left the discussion, I felt much more comfortable about your plans.)

A shift in what we've worked on in the past from what you're suggesting is before we treated the ODE function as the thing that matters and is what the user is interested in specifying. Here, the focus has shifted to the "solution operator." I think we should clearly define this, make it more prominent, and that will make it a lot easier to follow. For example, we should start by defining a "solution operator" (which isn't defined anywhere on the wiki as far as I can tell):

A *solution operator* is a function that accepts a start time, a start state, and an end time, and provides the end state. The solution is either a numeric integration of an ODE system or an analytic solution (where the ODE is not specified as a Stan function).

Should we change the arguments of integrate_ode_*() to accept a start time, a start state, an end time (instead of vector of times), and the parameters and data? This would play much nicer into what you're trying to do. It would also make life a lot simpler in terms of returning just a single state. 

Is "solution operator" the right name to give this thing? Does the community give this a name? Google isn't showing anything useful -- it shows hits related to linear differential operators.

If I understand a solution operator, it should have a function signature like this:

vector solution_operator(real t0, vector state0, real t) {

  ...  // return a vector the same size as state0

}

And what you want to do with this is feed this into an event handler. So what function signature does the event handler have? It would help if this were more clearly defined. (If we had global scope or some way of currying arguments, we could call integrate_ode_*() inside there; for now, we can augment that with roughly the same signature as integrate_ode_*() with the parameters and data and we'd be able to move forward).

For the event handler, is there a real desire to use different solution operators between different times? (It was mentioned in our meeting.) Is there any real use case for that? How would you specify it? If we don't have that requirement, this would be a lot easier.

It helps to have concrete language examples / prototypes. It helps me understand what you're trying to do. If you can provide that for:

- solution operators

- event handling

- DAEs

I would have some hope of catching up and providing additional information.

Is there any concern with misspecification of the solution operator?

What is an "event handler"? Does it just handle the event times? Does it do more? Why does it need a solution operator? I can't tell if an event operator requires a solution operator in order or if the event handler and a solution operator could be used in tandem. And if it's a design decision, why is it described the first way?

----

Summary / short list of concrete suggestions:

- rename the wiki page to deemphasize ODE and emphasize solutions once you have a solution operator.

- Give definitions of key terms / concepts up front. Things to define: solution operator, event handler. Please be as clear as possible.

- Refocus the page so the unnecessary pieces are away from the key parts of the discussion. The discussion of stiff / non-stiff / not-implemented solvers is not at all relevant on this page; unless I'm missing something (in which case, it should be clear why we should care).

Thanks for putting up the wiki and walking me through your thoughts.

Daniel

[Michael Betancourt]

“Solution operator” really isn’t the terminology, it’s more of

an “evolution” operator.  In any case, yes the abstraction is

to different ways of evolving a state between time points.

Events correspond to analytic transformations of the state

at discrete times and the modifications of the system for 

later times, such as zeroing out the states or elevating

them to a specific value or turning on an infusion.

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

[Bob Carpenter]

> On Aug 16, 2016, at 5:41 AM, Daniel Lee <bea...@alum.mit.edu> wrote:
>
> Charles, Bill, Sebastian, Andrew, and Bob: thanks for getting on a video call last week to help facilitate my understanding of what you're wanting to do with ODEs. It really helped me understand what's going on.
>
> I didn't want to edit the page without some agreement, so I'll start with edits here. Please feel free to address all, some, or none of it.
>
> Could we rename the wiki page? When I google "complex ode," it leads me to a wiki page on odes where the states are complex variables. I don't think that's what this page is about. Perhaps something like "Advanced techniques for models involving ODEs"?

Title doesn't matter, but -1 for both "advanced" and "complex". This
is all relatively basic and simple!

> One additional remark about renaming the wiki page: it's not focused on the complexity of ODEs at all. It abstracts that away. It's focused on techniques once you have a "solution operator." In the cases you're envisioning, the Stan program may not even specify an ODE. The dyanmic system may indeed have an ODE, but that may have been analytically solved so the representation is no longer an ODE from within the Stan language. (Honestly, this was my confusion going into the discussion and when I left the discussion, I felt much more comfortable about your plans.)
>
> A shift in what we've worked on in the past from what you're suggesting is before we treated the ODE function as the thing that matters and is what the user is interested in specifying. Here, the focus has shifted to the "solution operator." I think we should clearly define this, make it more prominent, and that will make it a lot easier to follow. For example, we should start by defining a "solution operator" (which isn't defined anywhere on the wiki as far as I can tell):
> A *solution operator* is a function that accepts a start time, a start state, and an end time, and provides the end state. The solution is either a numeric integration of an ODE system or an analytic solution (where the ODE is not specified as a Stan function).
>
> Should we change the arguments of integrate_ode_*() to accept a start time, a start state, an end time (instead of vector of times), and the parameters and data? This would play much nicer into what you're trying to do. It would also make life a lot simpler in terms of returning just a single state.

That would be less efficient than solving for multiple states at once
because of the way the integrators interpolate.

> Is "solution operator" the right name to give this thing? Does the community give this a name? Google isn't showing anything useful -- it shows hits related to linear differential operators.

It's not *a* solution operator, because in some cases, I believe
you still do have an ODE in there.

> If I understand a solution operator, it should have a function signature like this:
> vector solution_operator(real t0, vector state0, real t) {
> ... // return a vector the same size as state0
> }
> And what you want to do with this is feed this into an event handler. So what function signature does the event handler have? It would help if this were more clearly defined. (If we had global scope or some way of currying arguments, we could call integrate_ode_*() inside there; for now, we can augment that with roughly the same signature as integrate_ode_*() with the parameters and data and we'd be able to move forward).
>
> For the event handler, is there a real desire to use different solution operators between different times? (It was mentioned in our meeting.) Is there any real use case for that? How would you specify it? If we don't have that requirement, this would be a lot easier.
>
> It helps to have concrete language examples / prototypes. It helps me understand what you're trying to do. If you can provide that for:
> - solution operators
> - event handling
> - DAEs
> I would have some hope of catching up and providing additional information.
>
> Is there any concern with misspecification of the solution operator?

What I'm missing here is how the system is specified in such a way that
it'll work with analytical, numerical-analytical, and numerical-stepwise
solutions.

> What is an "event handler"? Does it just handle the event times? Does it do more? Why does it need a solution operator? I can't tell if an event operator requires a solution operator in order or if the event handler and a solution operator could be used in tandem. And if it's a design decision, why is it described the first way?

I think they mean this to be like NONMEM's notion of "event", so
I believe that encompasses:

* drug inputs
- state delta ("bolus" for PK/PD), with rate and compartment
- state continuous input (output?) ("infusion" for PK/PD), with rate and compartment
- 1st order abs (from PMXStan, no idea what this is)

The problem I have is that I don't really understand NONMEM's
notion of event or how the solutions work when there's not an
explicit numerical integration of an ODE system.

* desired solutions
- states and times (not necessarily every state in system)
- steady state (is this related to matrix exponentials?)

Sebastian mentioned we don't need derivatives with respect to solutions
if the solutions aren't used (that is, only one component of the state
might be used in the likelihood).

Another thing I didn't understand was how the system was going to be specified
if not as an ODE system. Do we have a fixed number of named systems that
we have explicit solvers for? What about these matrix exponential things
that somehow provide alternative solutions---is that just steady state?

Maybe the PMXStan package from Novartis corporate would shed some light?

http://andrewgelman.com/wp-content/uploads/2015/10/poster_PMXstan_ACoP2015_30Sep2015.pdf

- Bob

[Charles Margossian]

I'll start by answering the things I have a better grasp on, and extend later/ let Bill and Sebastian jump in. Will probably break things up in more than one message. 

“Solution operator” really isn’t the terminology, it’s more of
an “evolution” operator.

ok. 

I encountered the term "solution operator" in Sebastian's ReadMe on his prototype PKPD library. In pharmacometrics, this term makes sense, because the equations describing the evolution of a system are in ODE form, and said ODE needs to be solved. From the perspective of the user, the operator solves the ODE system. From the perspective of the developer, the operator may or may not involve an ODE solver. In particular it doesn't if we code in an analytical solution. 

What's more, not every time-evolving system is described by an ODE system, and requires a solver, even from the user's perspective. So I'm happy to roll with "evolution operator" and will adopt it for now, but remain open to other name suggestions. 

Proposed definition: 

evolution operator - Takes in the state of a system at time t0 and returns the state of that system at time t0 + t, provided that:

(1) knowing the state at time t0 fully defines the state at finite times

(2) between t0 and t0 + t, the system is isolated

Events correspond to analytic transformations of the state
at discrete times and the modifications of the system for 
later times, such as zeroing out the states or elevating
them to a specific value or turning on an infusion.

We can also have observation events. These events do not alter the system, but indicate we want some data at a certain time. So I like to distinguish between:

(1) observation

(2) state changers

The tricky thing is that, if we follow NONMEM conventions, a certain event can correspond to multiple events, and the event handler needs to interpret it as such. For example, 5 doses administered at regular time intervals will in the data set correspond to one event, when in fact these really are 5 events. 

--
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+unsubscribe@googlegroups.com.

[Charles Margossian]

It helps to have concrete language examples / prototypes.

I'll post example code from Torsten on the wiki -- as promised during the meeting!

Summary / short list of concrete suggestions:

- rename the wiki page to deemphasize ODE and emphasize solutions once you have a solution operator.

I'm ok with not exclusively focusing on ODEs, but I think it's fair to mention a lot of these models are ODE based.  

 

- Give definitions of key terms / concepts up front. Things to define: solution operator, event handler. Please be as clear as possible.

Will do. 

 

- Refocus the page so the unnecessary pieces are away from the key parts of the discussion. The discussion of stiff / non-stiff / not-implemented solvers is not at all relevant on this page; unless I'm missing something (in which case, it should be clear why we should care).

Solvers play a key role in the construction of evolution operators for ODE based models. I could link to the wiki page on ODE integrators  but the latter page would have to be updated. 

Title doesn't matter, but -1 for both "advanced" and "complex".  This
is all relatively basic and simple!

Ok with dropping "advanced" and "complex".

It's not *a* solution operator, because in some cases, I believe
you still do have an ODE in there.

True. But a function that returns integrate_ode_*(ODE, t0, t, ....) would be an evolution operator, according to the definition. 

The problem I have is that I don't really understand NONMEM's
notion of event or how the solutions work when there's not an
explicit numerical integration of an ODE system.

Will define event. 

We do not use a numerical integration if we can code an analytical solution. In pharmacometrics, there are certain recurrent models for which an analytical solution has been worked out (coding these is a pain, many many lines of code). In these cases where we want the user to be able to call built-in evolution operators. I'll point to examples in Torsten. 

1st order abs (from PMXStan, no idea what this is)

This has to do with model. It states that the drug is not directly dosed into the central compartment (often times blood) but that it first gets to the gut, and from there gets absorbed into the central compartment. An orally administered drug will go through the gut, whereas an intravenously dosed drug goes directly into the central compartment. 

What about these matrix exponential things
that somehow provide alternative solutions---is that just steady state?

 

Matrix exponentials give us the solutions to a system of linear ODEs, and as far as I can tell are not employed for Steady State approximations. I need to brush up on Steady States, but I'll produce an explanation on how to solve this problem. 

Ok - will edit the wiki, and will try to answer more of the raised questions. 

[Charles Margossian]

Updates on wiki:

-  Terminology section at the top. I admit I didn't feel I was really adding anything that wasn't already on the wiki page -- so please let me know if these definitions are unclear. 

- Implementation Example: Torsten  Looks at an example and shows how Torsten uses an event handler and an evolution operator

- Improving from Torsten's design These are notes I took during Friday's meeting on what a general event handler in Stan might look like. 

[Sebastian Weber]

Hi!

Much has already been said, so I will be brief:

- looking into my books on ODEs, I found "time evolution operator", so let's stick with that

- I would suggest to say that the time evolution operator function is responsible to calculate the states as long as no event is applied. As such we should leave the integrate_ode_* functions just as they are such that they take as arguments a scalar t0, initial state y0 and a vector of output times t. Otherwise things would get quite a lot more inefficient due to repeated adaptation performed by the ODE integrators.

- How we are going to define events will be a challenge! A definition must be very flexible and could be "the system is not event free evolving". Following this thought it could make sense to define another function which is the event handler function. This function is called upon the event times (data) and receives as arguments the current state y at the time of the event, some data and some parameters. It would return the modified state.

- Having defined an event handler as outlined would allow to include the "steady-state" thing into the event handler function as the function would simply set all states to the steady-state which may be pre-computed or being computed in that function.

- As events can be very different from case to case and even within a given model, it would be beneficial if the data structures being passed around are flexible - maybe tuples would be good to have.

- I am not sure if we want to deal with the circumstance that some events have an extended time-frame or if we decide to let the user handle things like "infusion is on at t1 and off at t2".

- I would like to point out the opportunity to define a "reset" or "next patient event". Having such an event would allow for openMP based parallelization of evaluating things which would yield a major speedup. I understand that this complicates things, but keeping it in mind would be potentially very useful. In case the design can accommodate this, then we would end up having a facility which can be used to parallelize the evaluation of the likelihood in hierarchical models.

Best,
Sebastian

[Devin Pastoor]

+1 to all Sebastian said. Extending the request for a rest or next patient event, I would hope that tied to such events would be an easy mechanism to output internally calculated "outcomes" - for example, a every common example is to calculate time after (last) dose, where the pseudocode logic is

last_dose_time = 0
if (dosing_event) {
last_dose_time = TIME
}

TAD = TIME - last_dose_time

export TAD

A reset event is also very handy for situations like titration designs where you are adaptively determining the next dose based on some calculation, like the final trough calculation.

Very exciting stuff!