cvodes_integrator feedback?

[Sebastian Weber]

Hi there!

It would be great if someone would take the time to have a quick look at the design of the cvodes_integrator:

https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/mat/functor/cvodes_integrator.hpp

The cvodes_integrator should now be in the state as we discussed last week in the meeting and other important bit, the ode_system, also got discussed such that I believe we are very far with this pull, I guess.

From the cvodes pull request the objects which would still require a review and tests are

https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/mat/functor/ode_system.hpp

and

https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/arr/fun/decouple_ode_states.hpp

The entry function

https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/mat/functor/integrate_ode_bdf.hpp

already got reviewed.

As I have tomorrow a full day time to work on these things, I could get this in a state where it would be ready to be merged (at least I would hope so). For that to happen it would be great if you could give me some input on the testing strategy. My understanding of this would be that I use the objects in the intended way with my tests and check for correctness.

The pull request with all comments leading to the current state is in
https://github.com/stan-dev/math/pull/283

Thanks a lot!

Best,
Sebastian

(wishful thinking ... getting this out with 2.10 could still be an option ...)

[Michael Betancourt]

Firstly, a reminder that moving forward we should use more descriptive branch names,
at least not having mysterious numerical labels! Something like “refactor_integrate_ode”
would be appropriate here, and the branch name can be changed so don’t hesitate to
do so before a pull request.

ode_system.hpp

Line 14: No reason to have abbreviations in comments.

Line 15: “S” and “P” aren’t used in this implementation, no need to define them.

Line 21: If it’s got methods, it’s a class. Just style.

Line 57, 76, etc: Qualitifers/types on same line. So “inline void” instead of “inline\nvoid”

integrate_ode_bdf.hpp

Line 32: No abbreviations. This is worth a few sentences.

Line 87: I recommend breaking this line to


cvodes_integrator<F, T_initial, T_param>
integrator(f, y0, t0,…

to free up lots of vertical space.

cvodes_integrator.hpp

Line 44: Why is ts_ a member variable? It’s used only in integrate and would be better
left as an argument to that method.

Line 45 - 47: N, M, and S need to be doc’d.

Line 177: Split this up into two if/then and move the common code in between, as in


if (solver == 0)
cvode_mem_ = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
else
cvode_mem_ = CVodeCreate(CV_BDF, CV_NEWTON);

if (cvode_mem_ == NULL)
throw std::runtime_error("CVodeCreate failed to allocate memory”);

set_cvode_options(rel_tol, abs_tol, max_num_steps);

if (solver != 0) {
// for the stiff solvers we need to reserve additional
// memory and provide a Jacobian function call
check_flag(CVDense(cvode_mem_, N_), "CVDense");
check_flag(CVDlsSetDenseJacFn(cvode_mem_, &solver::dense_jacobian),
"CVDlsSetDenseJacFn");
}

if (solver == 2)
check_flag(CVodeSetStabLimDet(cvode_mem_, 1), "CVodeSetStabLimDet”);

Note sure if "set_cvode_options(rel_tol, abs_tol, max_num_steps);” needs to be
factored into its own function as it would be called only once, but not a big deal
either way.

Line 214: I don’t think you need the braces, but not critical to remove them.

Line 252: Why cast to a solver const pointer? Did you mean const solver?

Line 312+: Can’t all of these methods be moved into ode_system, and then
pass a pointer to the ode_system implementation in CVodeSetUserData?
I think that’s the more natural factoring here.

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

Hi!

Thanks for your comments. Let me answer inline…

> Am 08.05.2016 um 17:07 schrieb Michael Betancourt <betan...@gmail.com>:
>
> Firstly, a reminder that moving forward we should use more descriptive branch names,
> at least not having mysterious numerical labels! Something like “refactor_integrate_ode”
> would be appropriate here, and the branch name can be changed so don’t hesitate to
> do so before a pull request.

Well, the CVODES stuff is getting a (long) history. I can rename, but then we loose the current pull request

https://github.com/stan-dev/math/pull/283

And the comments there were good, I think. Let’s quickly touch base tomorrow if we want to start from scratch with a new name and a new pull.

>
> ode_system.hpp
>
> Line 14: No reason to have abbreviations in comments.
>
> Line 15: “S” and “P” aren’t used in this implementation, no need to define them.
>
> Line 21: If it’s got methods, it’s a class. Just style.
>
> Line 57, 76, etc: Qualitifers/types on same line. So “inline void” instead of “inline\nvoid”

Thanks. Some comments still described old design (the abbreviations); I adapted this and tried to follow those style rules as I assume these are standard in Stan.

>
> integrate_ode_bdf.hpp
>
> Line 32: No abbreviations. This is worth a few sentences.
>
> Line 87: I recommend breaking this line to
>
>
> cvodes_integrator<F, T_initial, T_param>
> integrator(f, y0, t0,…
>
> to free up lots of vertical space.

Ok, done.

>
> cvodes_integrator.hpp
>
> Line 44: Why is ts_ a member variable? It’s used only in integrate and would be better
> left as an argument to that method.

I would prefer to leave it where it is. It belongs together with y0 and t0, I think. Moreover, we may want to revisit at some point the initial step size of cvode - it sets the initial step size based an a scary guess of y’’ and it is nowhere documented how this is done. So maybe we would like to set the initial step-size when configuring cvode in the constructor based on a scan through ts at some point.

>
> Line 45 - 47: N, M, and S need to be doc’d.
>

Sure. I add some notes for other bookkeeping ints as well.

> Line 177: Split this up into two if/then and move the common code in between, as in
>
>
> if (solver == 0)
> cvode_mem_ = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
> else
> cvode_mem_ = CVodeCreate(CV_BDF, CV_NEWTON);
>
> if (cvode_mem_ == NULL)
> throw std::runtime_error("CVodeCreate failed to allocate memory”);
>
> set_cvode_options(rel_tol, abs_tol, max_num_steps);
>
> if (solver != 0) {
> // for the stiff solvers we need to reserve additional
> // memory and provide a Jacobian function call
> check_flag(CVDense(cvode_mem_, N_), "CVDense");
> check_flag(CVDlsSetDenseJacFn(cvode_mem_, &solver::dense_jacobian),
> "CVDlsSetDenseJacFn");
> }
>
> if (solver == 2)
> check_flag(CVodeSetStabLimDet(cvode_mem_, 1), "CVodeSetStabLimDet”);
>
> Note sure if "set_cvode_options(rel_tol, abs_tol, max_num_steps);” needs to be
> factored into its own function as it would be called only once, but not a big deal
> either way.
>

Your proposal is very close to what I had before. Daniel and Bob found it to be preferable to have a single if/then/else - as such a set_cvode_option function makes a lot of sense. I am leaving it for the moment as is, since this way of writing was seen as a step forward.

> Line 214: I don’t think you need the braces, but not critical to remove them.

Yes, they are not needed, but actually I prefer them to write them one time too often rather than forgetting.

>
> Line 252: Why cast to a solver const pointer? Did you mean const solver?

const applies to what is on the left as I learned recently (and only to the right if nothing is left of it); I changed to const solver if that is more natural to read.

>
> Line 312+: Can’t all of these methods be moved into ode_system, and then
> pass a pointer to the ode_system implementation in CVodeSetUserData?
> I think that’s the more natural factoring here.

The intent of the design is to keep all C related code inside cvodes_integrator such that we have a single point where this interface to the C library happens. Moreover, having ode_system like it is now will allow us to (i) refactor the coupled_ode_system object using ode_system and (ii) have a straightforward way to supply analytic Jacobians to our ODE solvers using partial template specialization. I hope this makes sense to you, I believe Bob and Daniel came to agree to this approach.

In the meantime I added tests for the decouple_ode_states and the ode_system object. I hope those tests are inline with the overall testing strategy.

Thanks again for looking.

Best,
Sebastian

>
> On May 5, 2016, at 3:03 PM, Sebastian Weber <sdw....@gmail.com> wrote:
>
>> Hi there!
>>
>> It would be great if someone would take the time to have a quick look at the design of the cvodes_integrator:
>>
>> https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/mat/functor/cvodes_integrator.hpp
>>
>> The cvodes_integrator should now be in the state as we discussed last week in the meeting and other important bit, the ode_system, also got discussed such that I believe we are very far with this pull, I guess.
>>
>> From the cvodes pull request the objects which would still require a review and tests are
>>
>> https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/mat/functor/ode_system.hpp
>>
>> and
>>
>> https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/arr/fun/decouple_ode_states.hpp
>>
>> The entry function
>>
>> https://github.com/stan-dev/math/blob/feature/issue-262-CVODES-2/stan/math/rev/mat/functor/integrate_ode_bdf.hpp
>>
>> already got reviewed.
>>
>> As I have tomorrow a full day time to work on these things, I could get this in a state where it would be ready to be merged (at least I would hope so). For that to happen it would be great if you could give me some input on the testing strategy. My understanding of this would be that I use the objects in the intended way with my tests and check for correctness.
>>
>> The pull request with all comments leading to the current state is in
>> https://github.com/stan-dev/math/pull/283
>>
>> Thanks a lot!
>>
>> Best,
>> Sebastian
>>
>>
>> (wishful thinking ... getting this out with 2.10 could still be an option ...)
>>
>> --
>> 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.
>
> --

> You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/--YRRtMqFhE/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

[Bob Carpenter]

> On May 9, 2016, at 4:19 PM, Sebastian Weber <sdw....@gmail.com> wrote:
>
> Hi!
>
> Thanks for your comments. Let me answer inline…
>
>> Am 08.05.2016 um 17:07 schrieb Michael Betancourt <betan...@gmail.com>:
>>
>> Firstly, a reminder that moving forward we should use more descriptive branch names,
>> at least not having mysterious numerical labels! Something like “refactor_integrate_ode”
>> would be appropriate here, and the branch name can be changed so don’t hesitate to
>> do so before a pull request.
>
> Well, the CVODES stuff is getting a (long) history. I can rename, but then we loose the current pull request

You can just rename the pull request. It's no big deal either way.

>
> https://github.com/stan-dev/math/pull/283
>
> And the comments there were good, I think. Let’s quickly touch base tomorrow if we want to start from scratch with a new name and a new pull.
>
>>
>> ode_system.hpp
>>
>> Line 14: No reason to have abbreviations in comments.
>>
>> Line 15: “S” and “P” aren’t used in this implementation, no need to define them.
>>
>> Line 21: If it’s got methods, it’s a class. Just style.
>>
>> Line 57, 76, etc: Qualitifers/types on same line. So “inline void” instead of “inline\nvoid”
>
> Thanks. Some comments still described old design (the abbreviations); I adapted this and tried to follow those style rules as I assume these are standard in Stan.
>
>>
>> integrate_ode_bdf.hpp
>>
>> Line 32: No abbreviations. This is worth a few sentences.
>>
>> Line 87: I recommend breaking this line to
>>
>>
>> cvodes_integrator<F, T_initial, T_param>
>> integrator(f, y0, t0,…
>>
>> to free up lots of vertical space.
>
> Ok, done.
>
>>
>> cvodes_integrator.hpp
>>
>> Line 44: Why is ts_ a member variable? It’s used only in integrate and would be better
>> left as an argument to that method.
>
> I would prefer to leave it where it is. It belongs together with y0 and t0, I think.

I would think all of those would be arguments. I prefer to have
classes with constant members (aka immutable). They're much easier
to reason about. For instance, if you ever change any of y0, t0
or ts, you have to enforce consistency (t0 < ts[n]).

> Moreover, we may want to revisit at some point the initial step size of cvode - it sets the initial step size based an a scary guess of y’’ and it is nowhere documented how this is done. So maybe we would like to set the initial step-size when configuring cvode in the constructor based on a scan through ts at some point.

Should we leave the initial step size in as a parameter?

>> Line 312+: Can’t all of these methods be moved into ode_system, and then
>> pass a pointer to the ode_system implementation in CVodeSetUserData?
>> I think that’s the more natural factoring here.
>
> The intent of the design is to keep all C related code inside cvodes_integrator such that we have a single point where this interface to the C library happens. Moreover, having ode_system like it is now will allow us to (i) refactor the coupled_ode_system object using ode_system and (ii) have a straightforward way to supply analytic Jacobians to our ODE solvers using partial template specialization. I hope this makes sense to you, I believe Bob and Daniel came to agree to this approach.

Yes, fine by me. I like the encapsulation.

I have a very strong preference to use references rather than
pointers wherever possible. But I understand the external interface
requires pointers in some cases.

Thanks for looking into all this. I'm sure if we had a couple
more programmers, you could get two more different opinions :-)

- Bob

[Sebastian Weber]

Hi!

Sounds like we are converging as opinions (other than mine) start to even revert code.

Concerning the step_size: I don't think we should expose it as the automatic defaults appear to work very good in all instances I have seen so far. The rationale Michael noted to not expose stuff which gets adapted is a good one from my perspective.

Wrt to the ts vector: You argumentation supports the current design, right? I mean the t0, y0 and ts are stored inside cvodes_integrator in const references.

To your last point about pointers: They cannot be totally avoided in this context, but I did try to avoid them as much as possible.

As now opinions start to diverge - is there a process in place to decide when things are final/ready to be merged?

Sebastian

On Monday, May 9, 2016 at 11:33:07 PM UTC+2, Bob Carpenter wrote:
> > On May 9, 2016, at 4:19 PM, Sebastian Weber
> >

> > Hi!
> >
> > Thanks for your comments. Let me answer inline…
> >
> >> Am 08.05.2016 um 17:07 schrieb Michael Betancourt
> >>