new ODE functions
Bob Carpenter
We need to add the stiff ODE solver to the language, and so
need a name.
A. Current proposal is for four functions with two names (the
control params being optional):
integrate_ode_cvodes(...args..., ...control params...)
integrate_ode_cvodes_stiff(...args..., ...control params...);
B. An alternative, which I prefer, would be to drop "cvodes" from
the name and just go with
integrate_ode(...)
integrate_ode_stiff(...)
We'd then replace the Boost solver currently behind integrate_ode
with the CVODES non-stiff solver. Note that we don't name other
functions with the packages providing algorithms used to compute them!
But I seem to recall some issue with efficiency. Is the non-stiff
CVODES solver just overall better rthan the Runge-Kutta implementation
in Boost?
C. Add in the algorithm name, rather than implementation pacakge:
integrate_ode_adams_moulton(...)
and so on. I don't even know what these are, but they might be more
familiar to users than "cvodes".
Any preferences?
- Bob
> On Apr 12, 2016, at 2:59 AM, wds15 <notifi...@github.com> wrote:
>
> Hi Daniel,
>
> great that you look into this!
>
> On Tue, Apr 12, 2016 at 5:20 AM, Daniel Lee <notifi...@github.com>
> wrote:
>
> > Hi Sebastian, I'm finally able to look at it. I reread through the
> > existing comments. Here are mine.
> >
> > 1. stan/math/rev/arr/functor/coupled_ode_system_cvode.hpp was changed
> > to stan/math/rev/arr/functor/cvodes_integrator.hpp. What was the
> > reason for the name change?
> >
> > The logic of how things work is changed now. So before Stan was wrapping
> around the base ODE system and expanded it to become the coupled ODE
> system. Then the coupled ODE system was given to the integrator. However,
> CVODES knows itself of how to calculate the sensitivities, all what we have
> to do from the Stan side is to enable this integrator functionality and
> feed into CVODES the required input. This is why I changed to name to
> cvodes_integrator as there is not anymore a coupled ODE system on the Stan
> side, instead we feed into a sensitivity enabled integrator the ODE RHS and
> the RHS of the sensitivity system.
>
> Since coupled_ode_system_cvode.hpp included already all the code needed to
> interface with the C code, it made perfect sense to use it as a starting
> base.
>
>
> >
> > 1. This change looks like it propagated to
> > stan/math/rev/arr/functor/integrate_ode_cvode.hpp. Can you summarize
> > the change? I'm having getting my head around the change.
> >
> > Not sure what you refer to. In any case, we should settle on how the
> functionality is exposed to the language. I.e. the current proposal is to
> remove integrate_ode_cvode and replace it with 4 flavors of
> integrate_ode_cvodes (stiff/non-stiff simple/extended interface). However,
> it was left open so far to what extend the language definitions can be used
> to generate those function signatures. Maybe we quickly discuss later today.
>
> >
> > 1. We can remove the cvode code if we're going forward with this. I
> > believe we decided not to keep both CVODE and CVODES.
> >
> > Yes, CVODE is to be removed and my branch did replace CVODE with CVODES
> right away. There is really no point in keeping the old one.
>
> >
> >
> > I'm a bit tired now, so I'll try to look at it again tomorrow. I'm running
> > this command to tell what files have changed:
> >
> > > git diff --name-only develop -- stan
> > stan/math/rev/arr.hpp
> > stan/math/rev/arr/functor/coupled_ode_system_cvode.hpp
> > stan/math/rev/arr/functor/cvodes_integrator.hpp
> > stan/math/rev/arr/functor/decouple_states.hpp
> > stan/math/rev/arr/functor/integrate_ode_cvode.hpp
> > stan/math/rev/arr/functor/integrate_ode_cvodes.hpp
> > stan/math/rev/arr/functor/ode_model.hpp
>
>
Krzysztof Sakrejda
Sebastian Weber
I agree that the specific implementation (CVODES) is not so of interest to the user. I am suggesting a blend of the above:
- we have two simple go-to default functions for users:
integrate_ode => non-stiff
integrate_ode_stiff => stiff
These functions set the default abs+rel tol and maximum number of steps.
- Additionally, we have functions with more control which are named for their algorithm:
integrate_ode_rk45
integrate_ode_adams_moulton
integrate_ode_bdf
These functions would take extra arguments for abs+rel tol and maximum number of steps.
This way we would have covered expert users and normal mortals. This would also be in line with what I see in other ODE packages, i.e. they make the options named after the algorithm.
I would prefer to keep the rk45 from boost in addition to the adams_moulton from CVODES. The material is far to complex for me to exhaustively cover this with benchmark tests to decide which one is really better and I do think this will actually always depend on the exact model, etc.. RK45 is around for ages as is adams_moulton. Moreover, adams_moulton is advertised as requiring fewer calls to the ODE RHS and the sensitivity ODE RHS. From that angle adams_moulton should be better than RK45, since it avoids a few calls to AD.
If you think getting rid of boost odeint is a big advantage, then so be it.
Best,
Sebastian
Bob Carpenter
> integrate_ode => non-stiff
> integrate_ode_stiff => stiff
>
> These functions set the default abs+rel tol and maximum number of steps.
>
> - Additionally, we have functions with more control which are named for their algorithm:
>
> integrate_ode_rk45
> integrate_ode_adams_moulton
> integrate_ode_bdf
users and developers.
For users, there's much less to read and many fewer decisions to
make, so much less chance of being overwhelmed by detail or being
tempted into trying all the options. Of course, this only makes sense if expert
users don't need to make those decisions.
On the dev side, fewer functions and code dependencies means
proportionally less doc, testing, and ongoing maintenance.
- Bob
Sebastian Weber
Hmm... your arguments make sense. Also, thinking about how much code could go away if we remove the odeint RK45 integrator...that's quite a lot of code+testing code. From that perspective going with CVODES only is a big win, true. Ok, even though this won't cover everything I will fire off tomorrow a benchmark for a non-stiff problem which takes 2-3 days. That is a good stress test for the Adams-Moulton and the RK45 integrator.
And 4 variants of integrate_ode are easier to read than using algorithms name for the advanced functions, also agreed. I only think that this naming scheme is not really expandable (and only experts use it, so it's fine to have the names), but maybe that is just a minor point. As I have been exposed to these things so much, I like my details; but maybe they don't matter so much in the end. Bottom line: I can live with your option A.
Sebastian
Bob Carpenter
everyone else is OK getting rid of Boost's odeint and
just going with CVODES.
No matter what we do with function names, the details
of what's going on under the hood will be available.
I think it's critical to make the high-level code to
lower-level code path as transparent as possible for
power users --- it's the only way to write good code
in a language, be it C++ or Stan or R.
- 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
I have done now a bit of benchmarking and from this I can say that RK45 is 3-5x FASTER in comparison to the Adams-Moulton from CVODES! So I would definitely recommend to keep RK45.
In order to double check that my approaches with ode_model (used in the Adams-Moulton) etc. is not the root cause for the slow down I also went ahead an reimplemented the odeint integrator with these building blocks. The speed of this newly implemented odeint RK45 integrator is the same as the old one. The nice thing though is that I can cut down on lines of code from ~800 to ~300 and out of the ~300 about 200 come from ode_model and decouple_ode_states. Hence, I would recommend to refactor the old code with these pieces, but surely keep RK45.
I have the code written, but as it depends on pieces in the CVODES branch, can I just branch of that one?
Sebastian
(I ran the new RK45 code with the existing tests for odeint and all were good)
On Thursday, April 14, 2016 at 10:09:27 PM UTC+2, Bob Carpenter wrote:
> Let's discuss this at the meeting Tuesday and see if
> everyone else is OK getting rid of Boost's odeint and
> just going with CVODES.
>
> No matter what we do with function names, the details
> of what's going on under the hood will be available.
> I think it's critical to make the high-level code to
> lower-level code path as transparent as possible for
> power users --- it's the only way to write good code
> in a language, be it C++ or Stan or R.
>
> - Bob
>
> > On Apr 13, 2016, at 3:02 PM, Sebastian Weber
> >
Daniel Lee
Krzysztof Sakrejda
> Hi!
>
> I have done now a bit of benchmarking and from this I can say that RK45 is 3-5x FASTER in comparison to the Adams-Moulton from CVODES! So I would definitely recommend to keep RK45.
Sebastian, I'm sure this is stuff that's "known" somewhere by some experts but it's fantastic to have careful benchmarks done in-house for Stan to help us figure out how to incorporate these algorithms. Thanks! Krzysztof
Sebastian Weber
The benchmarks I did apply to my test case (2cmt oral PK model). I would recommend (in order of importance):
1. Keep integrate_ode with the RK45 for non-stiff
2. replace integrate_ode_cvode with a CVODES based stiff solver
3. Add a further integrate_ode function which allows to set abs+rel tol and maximum number of steps
4. use code from step 2 to refactor integrate_ode in order to have less code to maintain to go forward
5. consider exposing the adams_moulton of CVODES - I am sure some models will benefit from this
We still need to decide of how to name those functions; to be discussed Tuesday?
Sebastian
Sebastian Weber
Thanks - and I can tell you, this is quite some work. By now I have formalized a few tests which makes it a bit less work. However, from my experience you have to benchmark things in order to know. Moreover, there are "experts" out there, yes; but apparently there is not a golden rule or a standard. Otherwise, we would not see so many options out there. Different ODE systems my run with Adams-Moulton faster than with RK45, I am almost sure about that...
Sebastian
Krzysztof Sakrejda
> Hi Krzystof!
>
> Thanks - and I can tell you, this is quite some work. By now I have formalized a few tests which makes it a bit less work. However, from my experience you have to benchmark things in order to know. Moreover, there are "experts" out there, yes; but apparently there is not a golden rule or a standard. Otherwise, we would not see so many options out there. Different ODE systems my run with Adams-Moulton faster than with RK45, I am almost sure about that...
Agreed, and that's a great reason to formalize and keep these longer benchmarks around. Really it might be worth a repo of it's own to have performance-oriented benchmarks for different models/configurations... do you have these benchmarks on github?
Daniel Lee
Hi Daniel!
The benchmarks I did apply to my test case (2cmt oral PK model). I would recommend (in order of importance):
1. Keep integrate_ode with the RK45 for non-stiff
2. replace integrate_ode_cvode with a CVODES based stiff solver
3. Add a further integrate_ode function which allows to set abs+rel tol and maximum number of steps
4. use code from step 2 to refactor integrate_ode in order to have less code to maintain to go forward
5. consider exposing the adams_moulton of CVODES - I am sure some models will benefit from this
We still need to decide of how to name those functions; to be discussed Tuesday?
Krzysztof Sakrejda
> Hi Sebastian,
>
>
>
>
> Hi Daniel!
>
>
>
> The benchmarks I did apply to my test case (2cmt oral PK model). I would recommend (in order of importance):
>
>
>
> 1. Keep integrate_ode with the RK45 for non-stiff
>
>
>
> No action necessary.
>
>
>
> 2. replace integrate_ode_cvode with a CVODES based stiff solver
>
>
>
> This was merged yesterday.
>
>
>
> 3. Add a further integrate_ode function which allows to set abs+rel tol and maximum number of steps
>
> 4. use code from step 2 to refactor integrate_ode in order to have less code to maintain to go forward
>
> 5. consider exposing the adams_moulton of CVODES - I am sure some models will benefit from this
>
>
> We still need to decide of how to name those functions; to be discussed Tuesday?
>
>
>
> Ok.
>
>
> Sebastian, do we need to have your new `integrate_ode_cvodes` in the codebase before 2.10? Or is the action taken from #2 enough and we can work on everything else after we release.
>
>
> We really don't need to rush, especially when it's clear that the benefits of Boost vs CVODE vs CVODES are uncertain.
Sebastian Weber
To clarify what point 2: Yesterday we merge a replacement of the CVODE library with CVODES library, not more. So for stiff problems the integrate_ode_cvode function is still using the old approach of exposing the coupled_ode system to the cvode ode solver, but ignoring the extra capeabilities of CVODES. What needs to be merge is the pull request 283 on math for which I still owe you some documentation and a final review.
So currently the stiff solver in Stan is not fun to work with. In fact the non-stiff system I used for benchmarks of odeint freezes the integrate_ode_cvode function as it is right now. Only when you switch to the intgrate_ode_cvodes function it works.
So for stiff there is still a huge gain by switching to integrate_ode_cvodes.
@ Krzysztof: I don't have this benchmark in a repo yet. I guess I should clean the code and put it onto github, you are right.
Sebastian
Krzysztof Sakrejda
> Hi Daniel!
>
> To clarify what point 2: Yesterday we merge a replacement of the CVODE library with CVODES library, not more. So for stiff problems the integrate_ode_cvode function is still using the old approach of exposing the coupled_ode system to the cvode ode solver, but ignoring the extra capeabilities of CVODES. What needs to be merge is the pull request 283 on math for which I still owe you some documentation and a final review.
>
> So currently the stiff solver in Stan is not fun to work with. In fact the non-stiff system I used for benchmarks of odeint freezes the integrate_ode_cvode function as it is right now. Only when you switch to the intgrate_ode_cvodes function it works.
>
> So for stiff there is still a huge gain by switching to integrate_ode_cvodes.
>
> @ Krzysztof: I don't have this benchmark in a repo yet. I guess I should clean the code and put it onto github, you are right.
Bob Carpenter
> On Apr 16, 2016, at 5:51 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> On Saturday, April 16, 2016 at 12:12:56 PM UTC-4, Sebastian Weber wrote:
>
> I don't believe in cleaning code before you put it on a repo. Yes it would be nice to have prettified benchmarks but just usable benchmarks are a good step. :)
files, large blocks of commented out code, etc. Also, code
that's checked in should run for others and not depend on features
of a local environment. All that requires some degree of cleaning.
I always start with a repo or branch and try to keep everything
"clean" as I go. I find it's easier in the long run that way,
sort of like cleaning as you go when cooking a large meal.
We're less picky for test/benchmark code than for the Stan code
itself.
- Bob
Krzysztof Sakrejda
> > On Apr 16, 2016, at 5:51 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
> >
> > On Saturday, April 16, 2016 at 12:12:56 PM UTC-4, Sebastian Weber wrote:
> >> ...
> >> @ Krzysztof: I don't have this benchmark in a repo yet. I guess I should clean the code and put it onto github, you are right.
> >
> > I don't believe in cleaning code before you put it on a repo. Yes it would be nice to have prettified benchmarks but just usable benchmarks are a good step. :)
>
> Clean is a relative concept. We really don't want unused
> files, large blocks of commented out code, etc. Also, code
> that's checked in should run for others and not depend on features
> of a local environment. All that requires some degree of cleaning.
I don't think it has to go into one of the official repos but I refuse to temper
my enthusiasm on the subject!
>
> I always start with a repo or branch and try to keep everything
> "clean" as I go. I find it's easier in the long run that way,
> sort of like cleaning as you go when cooking a large meal.
You should see my sink.
K
Bob Carpenter
> On Apr 17, 2016, at 6:46 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> On Sunday, April 17, 2016 at 3:50:29 PM UTC-4, Bob Carpenter wrote:
>>> On Apr 16, 2016, at 5:51 PM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>>>
>>> On Saturday, April 16, 2016 at 12:12:56 PM UTC-4, Sebastian Weber wrote:
>>>> ...
>>>> @ Krzysztof: I don't have this benchmark in a repo yet. I guess I should clean the code and put it onto github, you are right.
>>>
>>> I don't believe in cleaning code before you put it on a repo. Yes it would be nice to have prettified benchmarks but just usable benchmarks are a good step. :)
>>
>> Clean is a relative concept. We really don't want unused
>> files, large blocks of commented out code, etc. Also, code
>> that's checked in should run for others and not depend on features
>> of a local environment. All that requires some degree of cleaning.
>
> I don't think it has to go into one of the official repos but I refuse to temper
> my enthusiasm on the subject!
We're more lenient on branches that aren't merged with
develop; but still, I'd like to at least keep the unused
files down.
>> I always start with a repo or branch and try to keep everything
>> "clean" as I go. I find it's easier in the long run that way,
>> sort of like cleaning as you go when cooking a large meal.
>
> You should see my sink.
- Bob
Sebastian Weber
My benchmarks tests are somewhat specialized to our computing environment. So I use a system like
1. An R script which sets up simulated data, inits and everything else what is needed to kick off a Stan run. This is the "model" file as I call it.
2. Another R script which takes the model file, parses it, creates a sub-directory, serializes the data into that directory, copies the stan source and executable into it, sets up init files and finally creates "fit.sh" and optionally "post.sh".
3. then I submit things into the SGE cluster queue
However, the script fit.sh can be called from the command line like
./fit.sh
(trivial). This will start chain id 1; if you want chain X then you do
SGE_TASK_ID=X ./fit.sh
Is that too "local" or would that suffices to linger in one of the stan repos?
Sebastian
Putting things into the stan github raises the bar a bit for me - I don't always go forward with clean code (although I put quite some tests in my code) and before I share code some polishing should be done...