Feature/issue 310 coupled ode system refactor

[wds15]

Submisison Checklist

• Run unit tests: ./runTests.py test/unit
• Run cpplint: make cpplint
• Declare copyright holder and open-source license: see below

Summary:

Refactor coupled_ode_system object to take advantage of ode_system and decoupled_ode_states which have been introduced lately.

Intended Effect:

1. Easier to maintain code base (about 500 lines less of code)
2. Make the non-stiff ODE RK45 integrator use the ode_system object which will enable to provide analytic Jacobians to the system (gives almost 2x speedup on simple pharmacokinetic model; the speedup will depend on the ODE specifics, but likely improve with larger system size).

How to Verify:

Run the ODE related tests which have almost all been moved now to test/unit/math/rev/mat/functor/.

Side Effects:

Minor internal only change to intergate_ode_rk45 and minor generalization of the ode_system.jacobian member function. The generalization of the ode_system.jacobian function may require to write tests for the generalized functionality which I will provide once design is agreed on to go into Stan as is.

Documentation:

Doxygen comments are updated and should be complete.

Reviewer Suggestions:

@bob-carpenter @syclik @betanalpha

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company): Sebastian Weber

By submitting this pull request, the copyright holder is requiring to license the submitted work under the following licenses:

• Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[syclik]

Hi @wds15, sorry it's taken so long to review. Since it's changing things around, I had to get a pen and paper to track through the changes.

[syclik]

The biggest design issue is that integrate_ode_rk45 is now a rev functor. It used to be in prim as a base template with a template specialization in rev.

I think if you changed the code to work that way and moved integrate_ode_rk45 back to prim, I'd be ok with it. I counted the lines of code and this pull request has ~250 less lines of code, but that's because the base template class was removed. I think if it were to be added back in, we'd see close to the same number of lines of code.

[syclik]

These rhs_sens() functions are what are forcing this class to be in rev instead of in prim. I would suggest moving the bulk of the class back to prim and templated appropriately. Then have a template instantiation in rev like it was done before.

[syclik]

@wds15, I don't think we've been clear about how our math library works. We want there to be versions that work with double only and do not require var or fvar. These functions live in prim. If we include stan/math/prim/*.hpp, we should only pull in doubles and not require our reverse mode autodiff. (I think there might be some dependencies in there somewhere.)

This pull request is removing functionality from prim and moving only to rev, which is also not going to work as we are trying to make fvars work also.

The code looks solid, but I believe the way we need to move forward with this is to:

• add back a stan/math/prim/mat/functor/coupled_ode_system.hpp (I'm ok with the move from arr to mat)
• make the class in stan/math/prim/mat/functor/coupled_ode_system.hpp be the general templated case; this needs to be generic enough to work for double and fvar
• make the class in stan/math/rev/mat/functor/coupled_ode_system.hpp be the template specialization for var
• move the tests from rev back to prim: they served a purpose

I think that's the way forward with this pull request, if that makes sense. Definitely open to having a design discussion if you think there's a different way to do this.

[wds15]

@syclik Thanks for taking the time and sorry to confuse things. I wasn't aware of the need to have something working in "prim". BTW, will fvar ever work with ODEs?

As I understand the main point is to restore the functionality in "prim" for the double only case (along with moving tests back). From what I understand I need to

1. Add a prim version for the double double case
2. Move tests for the prim case

I hope I got that right, we can discuss tomorrow as needed.

BTW, I hope that the refactoring makes the ODE code easier to follow, but the main advantage is that execution flows through the ode_system which makes the ODE subsystem amenable to providing analytic Jacobians.

Thanks for your time.

[wds15]

@syclik I hope I understood things right as to how the design is intended. As such I went ahead and enabled a prim functionality for coupled_ode_system, decoupled_ode_states and integrate_ode_rk45. This should get integrate_ode_rk45 and all its utils in line with what you described to me.

I am not sure where the merge conflicts come from, I need to follow those up (possible one of the new merges).

[wds15]

Ok, now things should be good. So what I did essentially is to place under prim a forward declaration of the objects needed along with the template specialization for double, double. Then in rev I give a definition of the "general" template which uses internally is_var meta programs (and this rhs_sens logic thing) in order to provide the correct implementation for any of the cases whenever at least one var is involved.

I also moved around tests, etc.. So I hope this is fine now as it is and brings back the old design idea of integrate_ode_rk45 living in prim and having special templates in rev for var types.

[syclik]

• The class defined in prim should be the general template class. The general template class shouldn't really be in rev. At least, that's how the rest of Stan math works. We could change that, but I think being consistent is better here. So, instead of forward declaring the templated class, the general templated class should live here.
• The place where the coupled_ode_system class is poorly designed is here because there are three rhs_sens() functions: 1, 2, 3. Only one can be instantiated by the due to the template types. To make this clearer, I think splitting these out into template specializations is the right thing. It'll make it easier to extend this as we move towards forward-mode. I'm thinking a little further along with the design than just this.

So, go back to what I suggested a few comments ago. Those are still the correct comments. In prim, there needs to be the general template class. In rev, you probably want 3 specializations of the class. This is going to make the code clearer and more extensible. Right now, you've got one class that does the role of 3 completely different things.

Does that make sense?

[wds15]

This is somewhat against my intuition of abstracting code - but I certainly don't have the bigger picture in mind. So what you suggest applies to

• coupled_ode_system
• decouple_ode_states

Right?

[bob-carpenter]

On Jul 21, 2016, at 12:28 PM, wds15 notifications@github.com wrote:

This is somewhat against my intuition of abstracting code - but I certainly don't have the bigger picture in mind. So what you suggest applies to

• coupled_ode_system
• decouple_ode_states
Right?

I want to make sure we're all on the same page here and I'm not
sure what your intuition is. What we're aiming for is modularity
where dependencies run in only one direction:

prim <- fwd <- mixed
<- rev <- mixed

Wherever possible, we want prim-only versions of functions.
This reduces dependencies, which is the particular abstraction
we're aiming for here. The goal isn't to minimize the size
of the code, but the overall complexity. Branches are considered
high complexity because of the overhead they add to testing
and because they make code hard to read.

In general, we prefer C++ code that doesn't have things like branching
on types, arguments, etc. This goes against the R style of having
single entry point functions with a lot of functionality.

• Bob

[wds15]

the approach to have a prim version with few dependencies makes sense, sure.

I figured doing the same logic as for CVODES stuff would be fine in the first place, but it looks like we now return to the very original design for good reasons.

Anyway, I hope that things are now as you expect them to be.

[syclik]

@wds15, there's a branch conflict. Mind looking into that?

[syclik]

@wds15, I'm taking another look at this now.

Additional questions / comments:

• I think the decoupled_ode_states() function should move back to being a member method of the coupled_ode_system class.... (edit -- I got through and found out where else it was being used)... if we're going down this path, I think we need to look at where these things are used. This is all really confusing and even more difficult to follow.
• I'm having a lot of trouble following the changes for this pull request. It looks like the only real substantial difference is using stan/math/rev/mat/functor/ode_system.hpp's jacobian(...) function instead of having it be in stan/math/rev/mat/functor/coupled_ode_system.hpp's operator() call. Is that right? If so, what is causing so many of these changes to the files?
• I see why you've split out decoupled_ode_states(), but this all seems to be factored into a way that's more complicated than it needs to be. Can we go back to thinking about what we need for integrate_ode_bdf() and integrate_ode_rk45() and come up with a coherent design? Right now, I see these as being in question:
  • class coupled_ode_system
  • function decouple_ode_system()
  • class ode_system -- in particular the jacobian() member methods.
    Maybe if you could help explain your design or explain what's required from the integrate_ode_* methods would help? As is, I see the changes made in this pull request and it's making things much more complicated by shifting things into objects that don't have clear function. Most of this problem isn't due to you -- it was clearer before with only one of the integrate_ode_* functions, but now we have a new use case.

I just need help with outlining what the objects do and why you've split the functionality across 3 different places instead of keeping them within one object where everything is coherent.

[wds15]

@syclik Thanks for taking a look.

The original and main intent of this pull is to route all ODE integrators through the jacobin function of the ode_system object, yes. This abstracts out the Jacobian such that we can provide partial template specialisations to the ODE framework and enable analytic Jacobians (this gives 2-3x speedups, depending on the ODE system vs using AD).

As I learned from you that integrate_ode_rk45 was for good reason under prim, I went ahead and moved it back there. The integrate_ode_bdf function must reside in rev, since the stiff solver always needs a jacobian (even for the double, double aka prim case).

The decouple_ode_states function had been factored out from the CVODES codes as I was anticipating that this can be reused for the coupled_ode_system object; just like I did.

The functionality which is needed is about:

• A jacobian function (wrt to parameters and wrt to states) just like in ode_system. The intent is to enable partial template specialization's such that we can have fully analytic Jacobians - supplying a single ode_system template specialization should enable analytic Jacobians for rk45 and bdf.
• The Jacobian is needed to create the respective ODE sensitivity RHS (what operator() from the coupled_ode_system did). The difference between CVODES stuff and coupled_ode_system is that the ODE RHS is separated from the sensitivity system. If you are willing to pay some efficiency, then we could make CVODES use the coupled_ode_system (but you would still have to give the jacobin wrt to states to CVODES as it is needed for the stiff integrator).
• As each integrator essentially produces the same output object, I introduced the decouple_ode_states function. All it does is take the raw output from the integrators and create a var instance as needed within Stan
• All this functionality is needed in a flexible manner, i.e. that is why I used the Eigen::MatrixBase stuff which allows flexible assignment to the different objects.

I have to say, that I find the design nice and straightforward... maybe I need to explain some more bits and pieces to you such that it makes more sense?

[syclik]

Here's where I'm confused. We now have two entry points into the ode systems (we used to only have one in the past):

1. stan/math/prim/mat/functor/integrate_ode_rk45.hpp
2. stan/math/rev/mat/functor/integrate_ode_bdf.hpp

---

In the first, integrate_ode_rk45, the abstraction makes sense, I think. Given:

• f, system of equations
• y0, initial conditions
• t0, initial time
• ts, times
• theta, variables
• x, data
• x_int, int data

We construct a coupled_ode_system with those things. This object provides these methods:

• initial_state()
• operator() for boost's integrator to use
• decouple_states() to go from a coupled ode system back to the original ode system

We have 4 ways of instantiating coupled_ode_system: T_initial can either be {double, var} and T_param can either be {double, var}. In the existing code we have 4 implementations where some of the implementation is shared.

---

The second use case, integrate_ode_bdf, is where I get lost.

Instead of the abstraction above, we now have a separate cvodes_ode_data structure, which is really close, but not the same as coupled_ode_system. This structure provides static functions, ode_rhs and ode_rhs_jacobian, required by CVODES.

I'm having trouble following the rest of the implementation here. Why isn't essentially all of cvodes_ode_data just living here? Why can't we used coupled_ode_system? Why do we need cvodes_ode_data? I don't think exposing decouple_ode_states() as a general function makes much sense.

---

Regarding what the integrator needs -- it's not clear to me why there are 3 different classes floating around when perhaps it should be 1? Or 0? Or 2? I'm having trouble keeping it straight. Perhaps this would be done quicker over a video chat. (Those are really questions; I don't know what's driving the decisions to split things into multiple classes or merge them.)

Regarding Eigen::MatrixBase -- that's fine as long as it works.

What happens to the integrate_ode_bdf code if it's instantiated with doubles?

[syclik]

So an example of where I'm really confused is inside: stan/math/rev/mat/functor/coupled_ode_system.hpp -- line 97, inside the operator() method.

Why doesn't ode_system_.jacobian(t, y_base, dz_dt_eig, Jy, dZ_dt_sens); take care of everything? Why doesn't the action of the jacobian() call just move into this method? I don't understand why the operator() method has to both set up the call to Jacobian and then also have to do things outside the Jacobian call. Why doesn't .jacobian() take care of things like

        dZ_dt_sens.leftCols(N_).setZero();
        dZ_dt_sens += Jy * Z_sens;

inside? Each instance of operator() calls a different instance of jacobian(), so no one function is called multiple times.

[wds15]

Maybe a hangout is really easier, but let me try in order:

• having the ode_system object with the jacobian functions is a design decision for two reasons:
  1. This isolates Jacobian calculations wrt to parameters and states. As such we can easily insert analytic expressions into the system. This gives large speedups for ODE systems. So I consider this a feature. I have already tested this and it is just an awesome speedup for ODE matters with Stan.
  2. So far we needed the Jacobians only to construct the coupled ODE system. With the introduction of the stiff solver we now also need the Jacobian wrt to states for the stiff ODE integrator which relies on this information to deal with stiffness. Hence, the need for Jacobians of the ODE RHS pops up at multiple places now (not only in the coupled_ode_system).
• The integrate_ode_bdf needs the cvodes_ode_data because this class has the required static call-backs needed by CVODES.
• cvodes_ode_data is also written with max. efficiency in mind when it comes to gluing things together.
• cvodes_ode_data could be written using the coupled_ode_system if we want to pay the price of not gluing things together with maximal efficiency. Doing so would very likely eliminate the need for decouple_ode_states as well.
• if you give integrate_ode_bdf only double+double then it will use ode_system,jacobian to get the Jacobian wrt to states in order to do stiff integration. No coupled system will be formed, just the base system is being integrated.

Here is a proposal to move forward, let's split things into two pull requests:

1. I rewrite coupled_ode_system to use the ode_system object. Doing so will enable great speed gains for the rk45 and the bdf integrator (whenever the user specifies the analytic Jacobians, I have working examples of that). Everything else will stay untouched.
2. refactor cvodes_ode_data to use the coupled_ode_system. Do so in a way such that we remove the decouple_ode_states function which will be part of coupled_ode_system again.

The first pull is the critical one to get ODEs faster in Stan, the second may cost us a little bit of performance (I don't know how much, probably not too much), but it will shrink the code base a bit.

[syclik]

Maybe the two pull separate pull requests will help me see what you're trying to do. When is the user able to specify analytic Jacobians? That isn't possible now, right? Maybe that's where I'm missing a key use of these classes.

(I still don't understand what's changing and for what reason. If you can create a more narrow pull request that just addresses speed, that'd be great. And a separate one to change the code, I hope it's a lot easier to understand.)

[wds15]

The user has to provide C++ code which is a template specialisation of ode_system. This specialisation then has to provide the analytic jacobians. I have examples which demonstrate this for the case of the bdf integrator. Doing my proposed pull 1 will enable the same speedup for the RK45 integrator.

So the speedup is not for free - the extra C++ code which is needed has to be provided. In the long run it will make sense to have scripts generating this.

Would such a working example help? And if yes, where should I put these files?

[syclik]

Yes, a working example would help. Where depends on what it is. I'm not sure where. What does it touch? The math library? The stan library? CmdStan? A generated model?

Are you creating a handful of specialized systems that have analytic Jacobians so you can call that from Stan math with some sort of name? If so, then the bulk would be a math change and putting that on a branch in math is the right thing.

[wds15]

The main point of the design is that providing the analytic Jacobians requires very little C++ code, i.e. just a single specialization of the ode_system object.

I create a branch which focuses on this aspect and I add a test which will run integrate_ode_rk45 without AD. This is more of an example rather than a real test, but you will see immediately how this is intended to work.

[syclik]

Thanks. Can you send a message on the stan-dev mailing list and I'll take a
look? That seems like a good use case. I just need to see what you're
intending to get my head around it.

On Thu, Aug 4, 2016 at 1:33 PM, wds15 notifications@github.com wrote:

The main point of the design is that providing the analytic Jacobians
requires very little C++ code, i.e. just a single specialization of the
ode_system object.

I create a branch which focuses on this aspect and I add a test which will
run integrate_ode_rk45 without AD. This is more of an example rather than a
real test, but you will see immediately how this is intended to work.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#312 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AAZ_F3OppotuhHSLbeatT3ENWWHK2M32ks5qciJvgaJpZM4JGyYc
.