Stiff ode solver

[Daniel Lee]

I'm looking at the stiff ode solvers again.

Should I look into this C implementation of lsoda?
http://lh3lh3.users.sourceforge.net/download/lsoda.c

Or sundial's cvode?

If there's no preference, I'll try to get the easiest one in first.



Daniel

[Bob Carpenter]

You can go back over all the old posts we got with advice,
but I think lsoda was the one that deSolve uses in R (albeit
the original Fortran implementation). I believe it computes
Jacobians with finite differences dynamically on an as-needed
basis.

My only concern with lsoda is the translation quality from Fortran.
I know other work by the translator, so have every reason to believe
it'll be reasonable.

I don't really know anything about cvode. How does it take
Jacobians? There's the much longer-about way that Michael
worked out the math for and requires higher-order autodiff (which
we conveniently now have). I'm not sure if that'll be faster in
the end, but probably not if it has to compute the Jacobian every
step, other than for problems where that's required.

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

[Tom HABER]

I did some work implementing CVODE as an ODE solver in Stan. CVODE is a really great solver, the implementation has been around for many years and based on even older fortran code.
It can deal with stiff ODEs and with failure cases when there are instabilities.

I'm collaborating with statisticians in Pharma as well as academics with lots of ODE models that don't fit in the box of nice models that RK45 can deal with. Most of the work up till now has been using custom build software, but I also would like to contribute the knowledge gained with the wonderful effort of Stan.

Unfortunately, CVODE is obviously a little bit more heavy for simpler models such as the harmonic oscillator or Lorenz. But there is definitely a gain with more complex models as it can deal with stiff models and usually take large time steps. I'm now trying to get the authorization from my collaborators to release some of the models and the data that we are working with to the Stan community. Those definitely require a more advanced solver.

I'm also working on an implementation of a more advanced technique to compute the gradients: there are cases when the sensitivity approach is too expensive. But this might take a little while longer to complete.

Would you be interested in this work or discussing it further?

As for lsoda/cvoda, both solvers are based on the same algorithms (BDF or Adams-Moulton). Both have also been around for a long long time and started out as Fortran code. I believe CVODE is a reimplementation (also geared towards parallelism) whereas lsoda is an automatic translation to C.
They are basically the de-facto standard for any ODE solving (except for the easy cases). I have a slight preference for CVODE, but I would say there is no real difference.
As for the jacobians, both packages can compute the jacobians using finite differencing but the user can also provide it as a function in dense or sparse format. I don't think that the Adams-Moulton algorithm uses Jacobians, since it's basically a Runge-Kutta + correction. For the sensitivity solution, the jacobian might involve higher-order derivation of the righ-hand side but FD might be fine. But there are different ways to go about this.

Tom

[Bob Carpenter]

That would be a big "yes". We'd love to have help
adding stiff solvers to Stan. It's our highest
priority to-do item for new language features and also
the one with which we have the least experience.

We've been doing applied work on some PK/PD and
also some soil-carbon models based on ODEs and other
users have been asking us for stiff solvers.

We're particularly interested in sorting out finite diffs
vs. autodiff gradients. As you may have seen, we have built
a coupled system for the non-stiff case that lets us compute
the first-order gradients we need for Euclidean HMC and optimization.

Michael Betancourt worked out the math for the coupled system
with Jacobians we'd need for implementing our own Jacobians
using autodiff, but we're interested also in the finite diffs
built in if they're not too cumbersome. The problem is that
doing the Jacobians would require higher-order autodiff on the
coupled system, and we're just getting the kinks worked out
of higher-order autodiff and it hasn't been very highly optimized
yet. We're also working on alternative autodiff we plan to plug
in at some point.

Other people with ODE solver experience have echoed
your sentiments about CVODE (and Sundials in general) and lsoda.
It's great they both have compatible licenses.

What did you mean that CVODE is heavier for simpler
models? More memory or time? If it's just integration headache,
we have to take that on anyway.

- Bob

[Tom HABER]

We've been working on ODE problems exclusively for over a year now. Mostly focused on parallelism and performance aspects. So sounds like a great match :)
The parallelism bit is also something I would love to contribute, but that is probably something more long term.

I have a working prototype which needs some more optimizations and linking needs to be cleaned up. The current coupled system uses std::vector a lot which doesn't fit nicely with CVODE's C interface, so at the moment it does some extra copying which isn't really required. Additionally it uses CVODE's built-in FD to compute Jacobians, which might still needs to be replaced as it would make the method much more efficient. Without resorting to higher-order autodiff, applying FD outside of CVODE could be beneficial because of the structure of the jacobian. Since CVODE has no idea about this structure it needs to do tons of extra work. I'll look into this next week.

For the simple models such as the oscillator, a stiff solver (such as BDF) is just overkill and something like Runge-Kutta or Adams-Moulton is perfect. But this can easily be an options for the user. In the next two weeks, we are hopefully able to release some of the models + data.

Tom

[Daniel Lee]

Hi Tom,

Let us know where we can be of help. Where is this work being done? Is it on a branch in a fork somewhere?

Are you using Sundials? Were you able to link a C library into Stan?

Daniel

Tom

[Tom HABER]

Hey,

I've put my latest version on my github fork: https://github.com/tomhaber/stan/tree/cvode

I just added the code providing jacobians to CVODE as discussed previously: instead of CVODE doing finite differencing, this code does finite differencing internally for computing the hessians. It can exploit the structure of the jacobian to compute it much faster: should be in order of #parameters faster.
When/if higher-order autodiff becomes available, this internal finite differencing can then be easily replaced.

The simple unit tests seem to work, but I should run a test on my other models as well.

The linking would be something that you guys could definitely help with. I'm currently using an ugly hack to completely avoid it. The header-only thing makes it a little hard. When I briefly looked at it, I could only see a possibility of adding it into the R binding code.

For now, it is more of a prototype. I use a wrapper that I made to hide any ode solver, but it might not make sense now (Unless you consider other ode solvers at some point). The implementation can be cleaned up/simplified otherwise.

Tom

[Tom HABER]

Has anyone had a chance to look at any of this?
Can anyone help out with the linking bit? I could of course look into this myself.

Tom

[Bob Carpenter]

I hope Daniel can help you out with the linking. I'm going to
be out of regular contact for the next couple of weeks due to traveling.

- Bob

[Daniel Lee]

Hi Tom,

Sorry about the delay. We're trying to do a major refactor of the code while there's a little down-time. I'm going to take a look at your branch soon. Unfortunately, I don't know that much about linking libraries, but hopefully I'll know enough to get it working.


Daniel

[Daniel Lee]

Hi Tom,

I'm starting to dig into this. How did you get over the hurdle of installation of CVODE? I'm having a little trouble there -- were you just expecting it to be a user issue?

Daniel

[Bob Carpenter]

Did you take a look at the package Sebastian sent?

https://github.com/vladimir-ch/odeity/

They got back to us and said it's all MIT or BSD
licensed now (other than the Sundials code itself).

- Bob

[Daniel Lee]

Not in depth. I'll try that today.

[Tom Haber]

Hey Daniel,

What seems to be the problem with the installation? Should be as simple
as "make && make install".

My problem was just when to do the linking correctly according to stan
rules. I guess rstan should do it, but all other derived packages should
know about CVODE as well. I couldn't link it into the stan library
(which now also doesn't exist anymore it seems).

As a side-note, definitely not all ODEs are stiff and the real problem
might be in the model itself (but stan shouldn't silently fail here).
Runga-kutta 45 is a good solver for a lot of problems, but the
implementation in boost is just not stable. CVODE on the other hand is
tried and proven. If people worry about the cost of an implicit solver,
CVODE also includes an explicit one and it is possible for an implicit
solver to take huge steps and still outperform a simple explicit solver.
Just my 2 cents...

Regards,
Tom

> +unsub...@googlegroups.com.

> > For more options, visit
> https://groups.google.com/d/optout.
>
> --
> 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

> +unsub...@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/ZG0tP6FxkPo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to

[Daniel Lee]

Hi Tom,

I'm trying to link against it using one of the examples provided by CVODE.

Although it looks like the installation looks like it succeeded, after compiling an example, I run into this issue:
> ./example
dyld: Library not loaded: libsundials_cvode.1.dylib
  Referenced from: /Users/daniel/dev/cvode-test-source/./example
  Reason: image not found
Trace/BPT trap: 5

and I'm lost at this point. I managed to install it correctly in some other place, but I'm uncertain of what options I used at that time. Any suggestions?

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Daniel Lee]

I just figured out the install of cvode, at least for my mac. I'll try it out on a windows box next.

[Daniel Lee]

Hi, I'm resurrecting this thread.

Tom, I finally got CVODE to install on both Mac and Windows; and unfortunately, it wasn't as simple as "make && make install" on both platforms.

Is the latest branch you're working on still here? https://github.com/tomhaber/stan/tree/cvode

Michael mentioned that you're presenting at Bayes-pharma. Great! Is it at a point where it should be integrated into Stan proper?

Daniel

[Michael Betancourt]

Had a chance to talk to Tom — we’ll spend plenty more time talking over

the next few days.  Firstly, the code in Tom’s repo should be up to date

with the current implementation.  What’s great is that CVODE could be

called from the C++ with a simple extern which bodes well for integration.

He also finite-differenced the autodiff gradients to approximate the Hessian

to compute the Jacobian instead of letting CVODE compute it implicitly.

It might be worth using the block-diagonal approximation to test integration

without having to wait for the higher-order autodiff.

There seem to be some stability issues but hopefully we can figure that

out when we work through the implementation.

[Tom Haber]

Awesome news Daniel! What can I do to help?

[Daniel Lee]

I have questions, then things I think we should do before we get it into the code base.

Questions:

1) Is the branch on your fork the latest version?

2) Could we move that to the main repo (github.com/stan-dev/stan)? I'd be happy to give you permissions to the main repo (if you don't already).

3) Have you dealt with conditional builds? I haven't thought too much about it, but we definitely want Stan to build even without having installed CVODE.

Things we need to do:

1) tests: unit tests, tests to make sure all generated ODEs work

2) doc: at the class level, in the documentation

3) evaluations: how do the different solvers (within and outside of CVODE) compare. Does it actually do a reasonable job with a stiff solver on a non-stiff problem? Should we be looking at ARKode or LSODA instead?

Daniel

[Daniel Lee]

Thanks, Tom. All great points. More inline.

On Thu, May 21, 2015 at 10:27 AM, Tom Haber <tom....@uhasselt.be> wrote:

On Thu, 2015-05-21 at 09:58 -0400, Daniel Lee wrote:
> I have questions, then things I think we should do before we get it
> into the code base.
>
>
> Questions:
>
> 1) Is the branch on your fork the latest version?
>

It is my latest/final version, but it is based on an older version of
stan/rstan. Since the big refactoring of the stan code, things seem to
have moved around quite a bit.

Ok; good to know that's the final version.

Yes, things have moved, but it's still tractable.

> 2) Could we move that to the main repo (github.com/stan-dev/stan)? I'd
> be happy to give you permissions to the main repo (if you don't
> already).
>

I could do that. I don't know how to deal with the dependence on an old
version of rstan however.

Ok. I'll let you know when it's on a branch in Stan. I don't see the connection to RStan?

 

> 3) Have you dealt with conditional builds? I haven't thought too much
> about it, but we definitely want Stan to build even without having
> installed CVODE.
>

I have done this many times and it's easy to do. Not sure how to do it
in the stan build system however. If this is not available in the
current build system, I'm afraid you'll need to add this feature. I
usually work with autotools or cmake.

We don't have it in yet and we don't use autotools or cmake. Do you know the mechanism that does it? C macros? Something else? The problem is that CVODE doesn't include enough headers to compile unless you run their install instructions. That means a user needs to install to make all of Stan compile. That isn't what we want.

 

> Things we need to do:
>
> 1) tests: unit tests, tests to make sure all generated ODEs work
>

Since my implementation is a drop-in replacement for the existing
solver, I'm just using the existing unit tests. I did need to change
some of the thresholds as I'm asking for such an accurate solution from
CVODE that the tests fail because the FD gradients are not good as a
reference. This was all checked in matlab with the true analytic
solutions. Maybe it's a good idea to add this as a unit test as well,
using the true analytic solution as a reference instead of the FD
estimate.

Yes, that would be great to add those as unit tests as well.

 

> 2) doc: at the class level, in the documentation
>

The CVODE implementation is a drop-in replacement, so I kept the
available documentation. I can surely add more documentation wherever
needed.

Thanks. Let's start by moving things and we'll take inventory later.

 

> 3) evaluations: how do the different solvers (within and outside of
> CVODE) compare. Does it actually do a reasonable job with a stiff
> solver on a non-stiff problem? Should we be looking at ARKode or LSODA
> instead?
>
>

The underlying algorithms in CVODE are the same as LSODA ((BDF and
Adams-Moulton). I haven't looked at the LSODA implementation, but I
suspect that it might have some small numerical differences. An
interesting aspect of LSODA is that it has heuristics builtin to switch
automatically between the explicit and implicit methods. Comparisons are
great, but I fear the results might be too problem-specific to
generalize. I would maybe even advocate the ability for the user to
decide on the method used (as well as absolute and relative tolerances!
). Using the simulation mode in stan, the user can investigate which
method gives the best trade-off.

Thanks for the info. I didn't realize it was using a heuristic to go between the two.

[Tom Haber]

On Thu, 2015-05-21 at 09:58 -0400, Daniel Lee wrote:

> I have questions, then things I think we should do before we get it
> into the code base.
>
>
> Questions:
>
> 1) Is the branch on your fork the latest version?
>

It is my latest/final version, but it is based on an older version of

stan/rstan. Since the big refactoring of the stan code, things seem to
have moved around quite a bit.

> 2) Could we move that to the main repo (github.com/stan-dev/stan)? I'd
> be happy to give you permissions to the main repo (if you don't
> already).
>

I could do that. I don't know how to deal with the dependence on an old
version of rstan however.

> 3) Have you dealt with conditional builds? I haven't thought too much
> about it, but we definitely want Stan to build even without having
> installed CVODE.
>

I have done this many times and it's easy to do. Not sure how to do it
in the stan build system however. If this is not available in the
current build system, I'm afraid you'll need to add this feature. I
usually work with autotools or cmake.
>
>

> Things we need to do:
>
> 1) tests: unit tests, tests to make sure all generated ODEs work
>

Since my implementation is a drop-in replacement for the existing
solver, I'm just using the existing unit tests. I did need to change
some of the thresholds as I'm asking for such an accurate solution from
CVODE that the tests fail because the FD gradients are not good as a
reference. This was all checked in matlab with the true analytic
solutions. Maybe it's a good idea to add this as a unit test as well,
using the true analytic solution as a reference instead of the FD
estimate.

> 2) doc: at the class level, in the documentation
>

The CVODE implementation is a drop-in replacement, so I kept the
available documentation. I can surely add more documentation wherever
needed.

> 3) evaluations: how do the different solvers (within and outside of
> CVODE) compare. Does it actually do a reasonable job with a stiff
> solver on a non-stiff problem? Should we be looking at ARKode or LSODA
> instead?
>
>

The underlying algorithms in CVODE are the same as LSODA ((BDF and
Adams-Moulton). I haven't looked at the LSODA implementation, but I
suspect that it might have some small numerical differences. An
interesting aspect of LSODA is that it has heuristics builtin to switch
automatically between the explicit and implicit methods. Comparisons are
great, but I fear the results might be too problem-specific to
generalize. I would maybe even advocate the ability for the user to
decide on the method used (as well as absolute and relative tolerances!
). Using the simulation mode in stan, the user can investigate which
method gives the best trade-off.

Tom

> 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

[Bob Carpenter]

Thanks! This looks great.

All we have to do to add different solvers and different parameters
is figure out what we want to do. I see we want at least this:

solver, absolute-tolerance, relative-tolerance

It'd also be nice to allow one or both of the data vectors to be empty.

It's a trivial update (no thinking, just coding I mean), but it touches the
parser and code generator and user-facing interface and thus doc. And the
more permutations there are, the more work it is now and to maintain going
forward.

- Bob

[Daniel Lee]

Tom, do you have install instructions for your branch? I've gotten it merged to develop, but I can't run the test cause I'm not sure how to configure cvode so it'll work with the changes to the makefile:

LDLIBS = `sundials-config -l c -m cvode -t s -s libs`

I'm not sure what that means; before I apply the fixes I made to the installation, I'd like to run your cvode branch first.

Daniel

[Daniel Lee]

Tom, I was able to merge your code and my install code and got the test running on a branch on Stan: feature/issue-1127-cvode

Installation instructions:

> make install-cvode

Then you should be able to run Stan as you normally would. For example, you can run the test:

> ./runTests.py src/test/unit/math/rev/arr/functor/integrate_ode_test.cpp

We need to do a lot more evaluation before this gets into Stan. On the unit tests (non-stiff, simple), I'm seeing a timings of:

develop (boost's solver):
harmonic_oscillator_finite_diff: 414 ms
lorenz_finite_diff: 380 ms

feature/issue-1127-cvode:

harmonic_oscillator_finite_diff: 2113 ms
lorenz_finite_diff: 2940 ms

(This is on my Mac; I haven't done any real evaluation on a problem that's stiff)

Are you seeing the same sort of timing differences on this simple test?

Daniel

[Tom Haber]

On Tue, 2015-05-26 at 16:53 -0400, Daniel Lee wrote:
> Tom, I was able to merge your code and my install code and got the
> test running on a branch on Stan: feature/issue-1127-cvode
>

Great work! Regarding the timings, I remember seeing something which was
in the same ballpark. Did you change the tolerances? I believe the
original version was using 1e-6 (abs. and rel.) whereas I'm using 1e-12
and 1e-9. This might explain the difference.
But remember it is a heavier solver and in some cases RK45 will be
faster.

Regards,
Tom

> +unsub...@googlegroups.com.

> > For more options, visit
> https://groups.google.com/d/optout.
>
> --
> 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

> +unsub...@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/ZG0tP6FxkPo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to

[Sebastian Weber]

Hi!

I got the code up and running on a fairly large example:

- Installation went semi-smoothly due to the yet experimental state of things, but install-cvode worked (after a cmake install)

- My model compiled with the develop cmdstan

- It starts to run and things look promising

- ...but after 2h or so the progress appears not be faster than with the RK integrator (not 100% sure though on this speed thing) *and* memory consumption continuously increases into the GB range (while the model starts off with a few 100 MB).

It looks to me that there is a memory leak and possibly further issues. Not sure.

Looking forward to the progress here, I really hope this gets ODEs going faster.

Cheers,
Sebastian

[Tom Haber]

There is indeed a small leak in the CVODE code that I fixed a while
back. In the destructor of CVodeIntegrator there should be an
additional call to "N_VDestroy_Serial(nvCurrentState_);"

As for the speed comments, keep in mind the tolerances used (see my
previous post)

Regards,
Tom

[Daniel Lee]

Hi Tom,

Do you have push permissions to the github.com/stan-dev/stan repo? If not, I'll give you access and we can continue to work on the branch there so we can eventually get it into Stan.

I think there's still a lot more to do before it gets in -- I'm still trying to get it up and running in Windows. I thought I had the installation working, but I found out that it was compiling CVODE into an incompatible object.

Daniel

[Tom Haber]

I don't think I have access :)
A colleague of mine has CVODE working in windows. I think he just used
the supplied CMake files and had it working pretty quickly. If you want
I can ask the details. He's using visual studio for this btw.

Tom

> 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/ZG0tP6FxkPo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to

[Sebastian Weber]

Concerning speed: What confused me is that you, Tom, suggest smaller tolerances which I would expect to make things slower, no?

In any case - I agree that having these fairly crucial parameters configurable by the user is indeed important.

Out of curiosity: To my knowledge cvode needs Jacobians as input - is the Stan autodiff used to get this?

I am happy to throw my example again at Stan-cvode...

On Saturday, May 30, 2015 at 5:47:57 PM UTC+2, Tom HABER wrote:
> I don't think I have access :)
> A colleague of mine has CVODE working in windows. I think he just used
> the supplied CMake files and had it working pretty quickly. If you want
> I can ask the details. He's using visual studio for this btw.
>
> Tom
>
> On Sat, 2015-05-30 at 11:37 -0400, Daniel Lee wrote:
> > Hi Tom,
> >
> >
> > Do you have push permissions to the github.com/stan-dev/stan repo? If
> > not, I'll give you access and we can continue to work on the branch
> > there so we can eventually get it into Stan.
> >
> >
> > I think there's still a lot more to do before it gets in -- I'm still
> > trying to get it up and running in Windows. I thought I had the
> > installation working, but I found out that it was compiling CVODE into
> > an incompatible object.
> >
> >
> >
> > Daniel
> >
> >
> >
> > On Sat, May 30, 2015 at 11:34 AM,

[Bob Carpenter]

> On May 30, 2015, at 12:56 PM, Sebastian Weber <sdw....@gmail.com> wrote:
>
> Concerning speed: What confused me is that you, Tom, suggest smaller tolerances which I would expect to make things slower, no?

Tom used smaller tolerances. When we compare speeds,
the tolerances should be the same.

I have no idea how much these tolerances matter. We
need to experiment more.

> In any case - I agree that having these fairly crucial parameters configurable by the user is indeed important.

We can add control for tolerances as one or two more arguments
to integrate_ode. Both absolute and relative?

Stan functions don't have optional parameters, so I'm always
reluctant to add params until we're sure we'll need them.

We could also do this through the compiler if we think we
just need to find the right global setting. That is, they
wouldn't be part of a model, but part of how it's compiled:

program <- stan_model(file="foo.stan", ode_tol_absolute=1e-6, ode_tol_relative=1e-6, ...);

> Out of curiosity: To my knowledge cvode needs Jacobians as input - is the Stan autodiff used to get this?

Nope --- not Stan autodiff yet --- I believe it just uses finite diffs
internally.

Michael's worked out the math so we can supply the Jacobians
with the sensitivities --- it's going to require higher-order
autodiff, though, with multiple passes per evaluation to get the
Jacobian, so it's not necessarily going to be any faster.

- Bob

[Tom Haber]

I have a model which can be a little sensitive to errors in the ODE
solution, so I made the tolerances much lower and thus slower.

The choice of the tolerances depends entirely on the model used and what
is expected from it. Like most tolerances and thresholds, it's not an
easy matter.

CVODE does require jacobians, but it has a default option to compute
them using Finite Differencing. However, in my implementation, stan's
autodiff is used to supply the jacobians to CVODE. In the sensitivities
case, you need the hessian of the ODE which I currently computed using
FD on top of autodiff. The stability needs some further investigation,
but speed-wise it is better than FD-ing a huge, very sparse matrix.

Tom

[Sebastian Weber]

@Bob: I think making these tolerances (absolute+relative) configurable just like the sampler options is the right thing to do. I think this is quite the same in the sense as it affects how the background numerics of Stan is configured and there is no need to ever adjust these tolerances during a stan run (i.e. for different parts of the data different tolerances is not needed). In the meantime - where can I set these thresholds manually, i.e. what do I need to "grep" for in the hpp files?

@Tom: Sounds like what you did made sense - is that how the current Stan cvode branch is setup to work? Also: Great that you already pinpointed the cvode memory issue - did you sent te patch to the cvode folks?

Best,
Sebastian

On Saturday, May 30, 2015 at 8:00:09 PM UTC+2, Tom HABER wrote:
> I have a model which can be a little sensitive to errors in the ODE
> solution, so I made the tolerances much lower and thus slower.
>
> The choice of the tolerances depends entirely on the model used and what
> is expected from it. Like most tolerances and thresholds, it's not an
> easy matter.
>
> CVODE does require jacobians, but it has a default option to compute
> them using Finite Differencing. However, in my implementation, stan's
> autodiff is used to supply the jacobians to CVODE. In the sensitivities
> case, you need the hessian of the ODE which I currently computed using
> FD on top of autodiff. The stability needs some further investigation,
> but speed-wise it is better than FD-ing a huge, very sparse matrix.
>
> Tom
>
> On Sat, 2015-05-30 at 13:34 -0400, Bob Carpenter wrote:

> > > On May 30, 2015, at 12:56 PM, Sebastian Weberw rote:

[Daniel Lee]

Unfortunately, that isn't going to work with Stan. For Windows, we need to get CVODE to install under RTools g++. The alternative is to get Stan to compile using MSVC. Last time we tried, the templating was beyond the capability of MSVC.

If you could ask your colleague for help, that would be great. 

[Ben Goodrich]

On Sunday, May 31, 2015 at 9:20:23 AM UTC-4, Daniel Lee wrote:

Unfortunately, that isn't going to work with Stan. For Windows, we need to get CVODE to install under RTools g++. The alternative is to get Stan to compile using MSVC. Last time we tried, the templating was beyond the capability of MSVC.

Someone on stan-users got Stan to install several months ago with Boost 1.56 and a recent version of MSVC.

Ben

[Ben Goodrich]

On Sunday, May 31, 2015 at 9:20:23 AM UTC-4, Daniel Lee wrote:

Unfortunately, that isn't going to work with Stan. For Windows, we need to get CVODE to install under RTools g++.

There is an R package that embeds sundials

http://cran.r-project.org/web/packages/Rsundials/

which seems to compile on Windows. Their Makevars file just has

PKG_CFLAGS=-w

Ben

[Bob Carpenter]

~/cmdstan(develop)$ grep -R "1e-6" stan/src/stan
...
stan/src/stan/math/prim/arr/functor/integrate_ode.hpp: const double absolute_tolerance = 1e-6;
stan/src/stan/math/prim/arr/functor/integrate_ode.hpp: const double relative_tolerance = 1e-6;


- Bob

[Sebastian Weber]

Hi!

Thanks Bob - that was obvious...

@ Tom: I was trying to find the respective code which needs to be fixed in cvode, but was not able to pinpoint it. Could you be more specific which file in cvode needs this fix? I would really like to try this out once again without a memory leak.

... also, I found no pointers to cvode on your github.

Many thanks!

Best,
Sebastian

On Sunday, May 31, 2015 at 11:22:56 PM UTC+2, Bob Carpenter wrote:
> ~/cmdstan(develop)$ grep -R "1e-6" stan/src/stan
> ...
> stan/src/stan/math/prim/arr/functor/integrate_ode.hpp: const double absolute_tolerance = 1e-6;
> stan/src/stan/math/prim/arr/functor/integrate_ode.hpp: const double relative_tolerance = 1e-6;
>
>
> - Bob
>

[Bill Gillespie]

Hi Bob,

I would prefer to see both relative and absolute tolerances as arguments, so that different calls to the solver within the same model could set different values. That would be a relatively rare requirement, but I prefer to have that flexibility.

What about the idea of "standard" and "expert" forms of the integrator function? The standard version would have fewer arguments and set the integrator parameters to defaults while the expert version would have a loner argument list that allows users to set those parameters.

Thanks,
Bill

[Bob Carpenter]

I think that makes most sense. Ones with and without a stream
of final config parameters.

Are there other bits of config that need to go in other than
absolute and relative tolerances?

- Bob

[Sebastian Weber]

Ok, fine with me. One more point to consider when goong forward is that we need this general dosing/event thing for which michael has a feature request. How can this be included or aligned with this approach?

Moreover you could consider adding a solver option which switches between possibly future choices of ode solvers

Sebastian

[Bill Gillespie]

Hi Bob,

There's quite a large collection of possible arguments for CVODE, so I'm struggling to decide what to make available and what to leave as defaults. Also I'm new to CVODE having been a regular LSODA user. Based on an initial read of the CVODE docs I recommend the following as arguments in the "expert" function. I nay suggest more on further reading.

* Arguments of CVodeCreate(lmm, iter)
** lmm = linear multistep method: CV ADAMS or CV BDF
** iter = type of nonlinear solver iteration: CV NEWTON
or CV FUNCTIONAL.
** You could restrict this to the pairs (CV ADAMS, CV FUNCTIONAL) for nonstiff
problems and (CV BDF, CV NEWTON) for stiff problems.

* Linear solver used for Newton iterations: cvdense, cvband, cvklu, cvsuperlumt, cvdiag, cvspgmr, cvspbcg, or cvsptfqmr.
** I could live with the default (CVDENSE), but models involving large systems with specific matrix structures could benefit from the choice.

Cheers,
Bill


On Monday, June 1, 2015 at 2:55:42 PM UTC-4, Bob Carpenter wrote:

[Bob Carpenter]

We haven't even decided on CVODE vs. other solvers in Sundials or
some reimplementation/translation of LSODA in C(++).

The expert function would have to set all of the options provided,
or we'd have to invent a new syntax just for optional args to
ODEs --- there are no optional arguments in Stan itself as is.

- Bob

[Tom Haber]

Hey,

In attachment I have a patch for that memory leak.
I'm having problems finding a version of rstan that works with this
version of stan. I would love some help with this.

Tom

On Sat, 2015-05-30 at 11:37 -0400, Daniel Lee wrote:

> 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/ZG0tP6FxkPo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to

[Tom Haber]

I would love to be able to switch between ADAMS and BDF in CVODE in
addition to the maximum number of iterations allowed and maybe the
initial step size (if you have expert knowledge about the ODE system,
you can avoid some computations).

Tom

[Tom Haber]

There are two checks in integrate_ode that aren't necessary when using
CVODE:

stan::math::check_ordered("integrate_ode", "times", ts);
stan::math::check_less("integrate_ode", "initial time", t0, ts[0]);

Since CVODE can go backwards in time (but this tends to be a harder
problem to solve) and some of the problems I come across have
measurements associated with the initial condition.

How do you feel about removing those?

Kind regards,
Tom

On Sat, 2015-05-30 at 11:37 -0400, Daniel Lee wrote:

> 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/ZG0tP6FxkPo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to

[Sebastian Weber]

Hi Tom!

Aweseome. Thanks for sharing the patch to cvode. Unfortunatley, I am not able to find cmdstan/stan combinations which run together. I am getting

icc -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -I src -I /home/weberse2/libs/cmdstan/stan/src -isystem /home/weberse2/libs/cmdstan/stan/lib/stan_math_2.6.3 -isystem /home/weberse2/libs/cmdstan/stan/lib/eigen_3.2.4 -isystem /home/weberse2/libs/cmdstan/stan/lib/boost_1.58.0 -Wall -pipe -DEIGEN_NO_DEBUG -DEIGEN_USE_MKL_ALL -DNDEBUG -isystem bin/cvode-inst/include -c -O0 -o bin/stan/lang/grammars/functions_grammar_inst.o stan/src/stan/lang/grammars/functions_grammar_inst.cpp
/home/weberse2/libs/cmdstan/stan/lib/boost_1.58.0/boost/proto/detail/preprocessed/make_expr_.hpp(73): error: class
"boost::proto::transform<PrimitiveTransform, X>::result<boost::phoenix::phoenix_generator (boost::proto::exprns_::basic_expr<boost::proto::tagns_::tag::assign, boost::proto::argsns_::list2<boost::proto::exprns_::basic_expr<boost::proto::tagns_::tag::terminal, boost::proto::argsns_::term<boost::phoenix::argument<3>>, 0L>, boost::spirit::_pass_type>, 2L>)> [with PrimitiveTransform=boost::proto::switch_<boost::phoenix::phoenix_generator, boost::proto::tag_of<boost::proto::_> ()>, X=void]"
has no member "type"
typedef typename proto_generator::template result<proto_generator(expr_type)>::type result_type;
^

but I think we should wait for the 2.7.0 release to be over before continuing this branch here (I even tried to merge the stan/develop into the cvode branch).

Best,
Sebastian

On Tuesday, June 16, 2015 at 11:32:44 AM UTC+2, Tom HABER wrote:
> There are two checks in integrate_ode that aren't necessary when using
> CVODE:
>
> stan::math::check_ordered("integrate_ode", "times", ts);
> stan::math::check_less("integrate_ode", "initial time", t0, ts[0]);
>
> Since CVODE can go backwards in time (but this tends to be a harder
> problem to solve) and some of the problems I come across have
> measurements associated with the initial condition.
>
> How do you feel about removing those?
>
> Kind regards,
> Tom
>
> On Sat, 2015-05-30 at 11:37 -0400, Daniel Lee wrote:
> > Hi Tom,
> >
> >
> > Do you have push permissions to the github.com/stan-dev/stan repo? If
> > not, I'll give you access and we can continue to work on the branch
> > there so we can eventually get it into Stan.
> >
> >
> > I think there's still a lot more to do before it gets in -- I'm still
> > trying to get it up and running in Windows. I thought I had the
> > installation working, but I found out that it was compiling CVODE into
> > an incompatible object.
> >
> >
> >
> > Daniel
> >
> >
> >
> > On Sat, May 30, 2015 at 11:34 AM,

[Bob Carpenter]

Thanks for the input. We'll probably have two calls: config-everythnig
and use-the-defaults.

We have to decide whether to make the functions in Stan solver-dependent.
Right now, it's just "solve_ode()" to call our non-stiff solver. But are
we going to turn that into:

solve_ode_sundials_adams(...., ...optional config params for sundials adams...)

???

- Bob

[Sebastian Weber]

Hi!

I finally got the Stan cvode branch up and running with the patch from Tom. Thanks again. I can confirm that the memory leak is gone. From the initial message of how long a single gradient evaluation takes it looks as if cvode takes 2-3 times as long as the boost solver. I have increased the tolerances to 1e-6 for relative and absolute error.

I am right now benchmarking the code on a large problem, let's see what happens, but it looks to me that we still need some tuning here. I will report once over - with boost my test case takes 3days per chain...

Best,
Sebastian

On Tuesday, June 16, 2015 at 9:50:00 PM UTC+2, Bob Carpenter wrote:
> Thanks for the input. We'll probably have two calls: config-everythnig
> and use-the-defaults.
>
> We have to decide whether to make the functions in Stan solver-dependent.
> Right now, it's just "solve_ode()" to call our non-stiff solver. But are
> we going to turn that into:
>
> solve_ode_sundials_adams(...., ...optional config params for sundials adams...)
>
> ???
>
> - Bob
>
> > On Jun 16, 2015, at 5:26 AM, Tom Haber
> >

> > I would love to be able to switch between ADAMS and BDF in CVODE in
> > addition to the maximum number of iterations allowed and maybe the
> > initial step size (if you have expert knowledge about the ODE system,
> > you can avoid some computations).
> >
> > Tom
> >
> > On Tue, 2015-06-02 at 12:58 -0400, Bob Carpenter wrote:
> >> We haven't even decided on CVODE vs. other solvers in Sundials or
> >> some reimplementation/translation of LSODA in C(++).
> >>
> >> The expert function would have to set all of the options provided,
> >> or we'd have to invent a new syntax just for optional args to
> >> ODEs --- there are no optional arguments in Stan itself as is.
> >>
> >> - Bob
> >>
> >>
> >>> On Jun 2, 2015, at 9:45 AM, Bill Gillespie
> >>>

[Bob Carpenter]

Thanks for looking into this. It'd be nice to develop a
suite of benchmarks. I have a couple other compartment models
for soil carbon respiration we could use.

- Bob

[Sebastian Weber]

Yes, a test-suite would make a lot of sense. To report back: the cvode runs are far slower, but they seem to run more consistent, i.e. all chains progress at a very similar pace. The runs are not yet completed, I will report back once the runs are over.

This puzzles me a bit, I hope cvode can be run in a faster way (I already went to 1e-6 absolute and relative tolerances). @Tom: Any suggestions? What parameters are important to take care of (Stan wise and cvode wise)?

Out of curiosity: Is it thinkable to get

http://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/tutorial/parallel_computation_with_openmp_and_mpi.html

up and running in Stan with the good old boost odeint integrator? The instructions sound very straightforward or am I overlooking something which prevents openMP for ODE solving in Stan fundamentally? Usually I get to hear "autodiff" is the killer, but isn't this now separate enough from autodiff?

Best,
Sebastian

On Friday, June 19, 2015 at 8:26:04 PM UTC+2, Bob Carpenter wrote:
> Thanks for looking into this. It'd be nice to develop a
> suite of benchmarks. I have a couple other compartment models
> for soil carbon respiration we could use.
>
> - Bob
>
> > On Jun 19, 2015, at 9:53 AM, Sebastian Weber
> >

[Bob Carpenter]

> On Jun 22, 2015, at 5:34 AM, Sebastian Weber <sdw....@gmail.com> wrote:
>
> Yes, a test-suite would make a lot of sense. To report back: the cvode runs are far slower, but they seem to run more consistent, i.e. all chains progress at a very similar pace. The runs are not yet completed, I will report back once the runs are over.
>
> This puzzles me a bit, I hope cvode can be run in a faster way (I already went to 1e-6 absolute and relative tolerances).

I'd have only expected it to be faster if there was substantial
stiffness in the equations. Otherwise, the Boost thing's doing far
less work per iteration.

> @Tom: Any suggestions? What parameters are important to take care of (Stan wise and cvode wise)?
>
> Out of curiosity: Is it thinkable to get
>
> http://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/tutorial/parallel_computation_with_openmp_and_mpi.html
>
> up and running in Stan with the good old boost odeint integrator?

> The instructions sound very straightforward or am I overlooking something which prevents openMP for ODE solving in Stan fundamentally? Usually I get to hear "autodiff" is the killer, but isn't this now separate enough from autodiff?

This should be possible. Our interface now just calls the ODE solver
with double values, so there's no race conditions with autodiff variables.
It may be possible to do the same with CVODE or other approaches, too.

I created a feature request and assigned it to Daniel, since this
is mostly going to be ODE and make-related (sorry, Daniel!):

https://github.com/stan-dev/stan/issues/1511

We haven't been too worried about parallelism up until now because most
of our users fit modestly sized models on desktops/notebooks and have
at most four cores, which are fully occupied running four chains in
parallel. But now that we're up against bigger models, putting each chain out
to a multiprocessor makes sense if we can manage.

- Bob

[Sebastian Weber]

Uhh... I did not intend to put this on daniels plate, but this parallelism will be a very important step - i run stan on machines with 20 cpus and if you give me a mpi enabled stan than i scale this up into the hundreds on infiniband wired clusters. That would be awesome!

Sebastian

[Joshua N Pritikin]

Has anybody looked at taking advantage of OpenCL?

--
Joshua N. Pritikin
Department of Psychology
University of Virginia
485 McCormick Rd, Gilmer Hall Room 102
Charlottesville, VA 22904
http://people.virginia.edu/~jnp3bc

[Sebastian Weber]

Hi!

I have no clue how much this will bring us in terms of speed (stuff needs to be copied around from RAM to the GPU and back), but according to the docs:

http://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/tutorial/using_opencl_via_vexcl.html

This also looks very easy to do. The only drawback is the C++11 from the VexCL library.

Best,
Sebastian

[Bob Carpenter]

I don't think Stan's ready for C++11 to be turned on
by default or even not by default because it'll break
CRAN compatibility. (Ben should know more). And it'll
break everything forced to use older compilers, including
R and any other project that forces older compiler versions.

But there's no reason we can't play around with it in experimental
branches. Or provide some way to modify.

The problem is that we don't want to just keep pasting on
dependencies and ifdefs, and whatnot --- the dependencies are
already killing us.

- Bob