if we’re using your doc, can you be much more pedantic about the math? You know a lot more about this than I do and there are steps missing where I’m not following.
I can follow Michael’s doc and that’s what we implemented for boost / odeint.
Can you share the tex? And we don’t need to talk about autodiff in the math. Just zeros. Zeros are important.
wds15 Developer
October 18
Hi Daniel!Have a look at what I wrote here:
https://groups.google.com/group/stan-dev/attach/5b162997cf789/ode_system.pdf?part=0.1&authuser=0
Basically, the sensitivity states are crammed into a matrix S. Stheta is a NxM matrix, i.e. per parameter theta_m for which we take the sensitivity there is a column of N rows in that matrix (one per state).
In that notation the sensitivity ODE RHS is given by
• for the initials as
dSy0/dt = Jy * Sy0• for the parameters as
dStheta/dt = Jy * Stheta + JthetaThat means for the 3 cases with any sensitivites we endup with:
• v,d = varying initials, fixed parameters
• states (y,Sy0)
• Jy needed only
• d,v = fixed initials, varying parameters• states (y, Stheta) • Jy AND Jtheta needed• v,v = both are varying
• states (y, Sy0, Stheta) • Jy AND Jtheta neededThe code you refer to buries the math behind this - while I think it would be good to take advantage of the nice structure of matrix-vector products which will easy readability of the code a lot (at least to me).
I was saying that we have one adapter per integrator, yes. However, I do think that we should design things such that we do not need 4 different specializations underneath for each adapter. I think that can be done and it would make maintenance easier in going forward.
I’m OK with one specialization as long as it’s not a conditional/case-based
thing. Otherwise, I’d prefer multiple specializations with shared code
pulled out into subclasses or utility functions so that the objects actually
being used are simpler to reason about and test.
That is, I don’t want something like this
struct big_config_object {
int object_type;
that stores what type its for and then conditions behavior as in:
if (object_type == 1) ... do something ...
else if (object_type == 2) ... do something ...
Hi!
@ Bob: In my mind we introduce an ode_rhs_senstivity object which you can see as a utility object for the adapters. This object would exist in the 4 different specialisations we have and provide the things needed specific to each case in a general enough way so that we can feed it easily into the different data-formats of each integrator. That will save a lot of code-duplication and can include calculation of the sensitivitiy RHS and the decoupling operation. I am not sue, but maybe CRTP is the way to do it, but last time I tried this, I failed to setup CRTP in a meaningful way.
@ Daniel: p_i denotes any sensitivity parameter. So that can be y0, or/and theta. Where have I left out steps? I have written this up as markdown and can share it later, sure. To be clear - what I have written is exactly the same as Michael has put down; it’s just written differently. If you go through the for loops in the coupled_ode_system you will find that the expressions give the same as “mine”. There is only a minor difference in how the initial values of the sensitivity system are setup. The way I do it will avoid the shifting operation which Michael introduced whenever initials vary.
Sebastian
wds15 Developer
October 19
Hi!@ Bob: In my mind we introduce an ode_rhs_senstivity object which you can see as a utility object for the adapters. This object would exist in the 4 different specialisations we have and provide the things needed specific to each case in a general enough way so that we can feed it easily into the different data-formats of each integrator. That will save a lot of code-duplication and can include calculation of the sensitivitiy RHS and the decoupling operation. I am not sue, but maybe CRTP is the way to do it, but last time I tried this, I failed to setup CRTP in a meaningful way.
I think we’re saying the same thing. I wouldn’t bother with CRTP for
four cases—just use the C++ concept directly.
The way to save code duplication is always to write reusable
functions—it’s just a matter of where to package the code
that gets reused: as simple functions, as (abstract) OO base
classes, as CRTP base classes, or some combination thereof.
@ Daniel: p_i denotes any sensitivity parameter. So that can be y0, or/and theta. Where have I left out steps? I have written this up as markdown and can share it later, sure. To be clear - what I have written is exactly the same as Michael has put down; it’s just written differently. If you go through the for loops in the coupled_ode_system you will find that the expressions give the same as “mine”. There is only a minor difference in how the initial values of the sensitivity system are setup. The way I do it will avoid the shifting operation which Michael introduced whenever initials vary.
Sounds good.
Maybe Michael will have some suggestions. I don’t understand
all of this well enough to have much of an opinion other than on
API organization at a general level.
HI Sebastian,
That seems like a really big deal. Perhaps an explanation why the initial values are set up differently? That would help.
Hi Daniel,
if you do setup the initials like I suggest, then you do not need to perform the shifting operation. The sensitivities as they get calculated can be taken directly as they are. The CVODES integrator works exactly this way and the sensitivities which get from it match the one from odeint.
Sebastian