piecewise (constant) functions in ode's

[tom....@uhasselt.be]

Hey everyone,

With Stan's new ability to use ode's (great job!), we were interested in implementing our model in Stan. We sadly ran into a few snags.

Our ode has a piecewise constant function inside (based on model parameters), which we currently implemented using a simple loop. This works great while simulating with "generated quantities", 

but when running the sampler it hangs forever. I can see the loop and a piecewise constant function posing a problem with the sensitivities, but I don't really see how to get around it.

We could maybe change to a parameterized function if that would avoid the problem.

It's currently coded as follows:

idx <- 1;
while(idx <= M && intervals[idx] < t) {
    idx <- idx + 1;
}

beta_t <- (transmissionrate + theta[9+idx]) / (prop/gamma + (1-prop)/alpha);

dydt[] using beta_t

We would also like to lower the error tolerances to 1e-9, based on simulations we did in Matlab. It seems like this feature was discussed previously, but I can't seem to find any documentation on how to do it.

Thanks in advance.

Kind regards,

Tom

[Bob Carpenter]

What we do with discontinuities due to dosing is to is to run the
integrator between dosing events. I hope somebody else can comment
on issues your system may bring up. It might help to see more of it,
or maybe what you sent is enough for someone to help out.

We don't yet have a way to control tolerances, but we'll
be adding that. It's not much implementation, but it's a whole
bunch of testing and doc.

The current integrator is indeed prone to hang,

http://headmyshoulder.github.io/odeint-v2/doc/boost_numeric_odeint/odeint_in_detail/steppers.html

http://en.wikipedia.org/wiki/Dormand–Prince_method

As far as we can tell, what's going on is that the parameters are moving
into spaces where the current integrator requires enough steps to reach
tolerance that it effectively hangs. We haven't figured out an easy way
to limit the max number of steps and also Michael's warning us that may just
bias the results.

We'll be introducing a stiff integrator, hopefully soon, but we're
pretty backed up on issues and feature requests and I'm not sure how
that one's being prioritized.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[tom....@uhasselt.be]

Hey Bob,

I figured that doing the integration in multiple intervals would be really slow, so I tried to avoid that option.

Right now I've simplified the model down to bear minimum, replaced the piecewise function with a simple parameterized function and simulated the data. But still it hangs on me.

It's really strange how it behaves: it starts at the exact parameters from which the data is generated and immediately ventures off to regions in 1e+214.

I've looked over the model a bunch of times, but I'm clean out of ideas.

Hope someone can help.

Thanks,

Tom

[Daniel Lee]

Hi Tom,

On Tue, Dec 2, 2014 at 8:53 AM, <tom....@uhasselt.be> wrote:

I figured that doing the integration in multiple intervals would be really slow, so I tried to avoid that option.

Why do you think that? (that's not a rhetorical question.... do you have some intuition behind this?)

From what I've seen, it's not the number of intervals that will make Stan run slow. The current ode integrator is not a stiff solver and will break down if the parameter values get into stiff areas.

NUTS / HMC will definitely break down if there are discontinuities in the gradients of the log probability function, which you're doing explicitly. I know that for a fact. It'll still "run" -- just not efficiently.

 

Right now I've simplified the model down to bear minimum, replaced the piecewise function with a simple parameterized function and simulated the data. But still it hangs on me.

It's really strange how it behaves: it starts at the exact parameters from which the data is generated and immediately ventures off to regions in 1e+214.

Is there any way to double-check the math? If it's running to infinity, then I would think something is wrong with the model. Could you put in prints and make sure the solutions are correct?

I would try Bob's suggestion: run the integration in multiple intervals. I think that's the sensible solution at this point. Hopefully we can get this model fit and then we can talk about speed after that.

 

I've looked over the model a bunch of times, but I'm clean out of ideas.

Hope someone can help.

Thanks,

Tom

--

[tom....@uhasselt.be]

The integrator is already making stops according to the stepsize, but it tries to take large steps when the error tolerances allow it. Integrating in intervals will put a limit on the maximum stepsize and might force much smaller steps than it would otherwise take. Obviously, a piecewise constant function contains discontinuities and the integrator should take small steps around the discontinuity.
However when making long term projections in a multiple dosing regimen, the solution trends to a kind of steady-state and the discontinuities are relatively small. This allows the integrator to still take big steps instead of stopping for the daily dose.

This is just one specific case and does not hold in general. Besides that, there is an overhead in setting up the integrator. I'm not 100% sure how the gradients are computed but this might also give extra costs for multiple integrations. (There was some discussion about adding sensitivities to the ODE system, but from a brief look at the code it seems like reverse-mode autodiff of the integrator)

That said, in the previously attached model I replaced the piecewise constant function with a nice parametrized one. This avoids that problem entirely. Sadly, this "easier" model is also giving troubles indicating additional problems. I've added some print statements to the model (thanks for the tip, didn't know that stan had printing functionality. Much easier than using a debugger) and it seems ok, giving the correct values and log-posterior. But some parameters immediately jump to infinity. I guess the gradient must be incorrect. Is there a way to easily check that, besides the debugger again? :)
The ode integrator breaking down for some parameters will still be an issue in the end because the dynamic system just becomes unstable for some parameters. Finding the correct constraints on the parameters is however no-obvious. But my hope is that if you start somewhere reasonable, there's no need to jump to unstable parameter regions.

Btw: a stiff ode solver might still break down if the system is unstable. I ran into this problem many times and it seems like the sampler should be capable of dealing with these cases if the parameters cannot be constrained somehow.

[Bob Carpenter]

I'll have to try to get my head around how one of these
stepwise integrators can handle discontinuities. It might
be more expensive to compute stepwise, though my guess would
be the opposite --- by explicitly telling it where the
discontinuities are, it could then use better stepping between
them. But I don't have much experience with these algorithms
and none with discontinuous systems.

The integrator has sensitivities added via a nested
reverse mode autodiff. Our original prototype just
autodiffed the integrator.

You can control the initial jumping into unstable regions
by reducing the initial step size. You can keep the
jumps small (and thus more curvature respecting) by increasing
the target acceptance rate. I've found 0.001 and 0.95 to be
pretty stable across ODE cases I've tried so far.

- Bob

[tom....@uhasselt.be]

Alright! Thanks for the information.

I've made some progress by added (fake) constraints to the parameters to stop them from jumping into unstable regions and used the suggested step_size and adapt_delta.

Now it's doing something, but sadly my machine runs out of memory before finishing 1000 sampling iterations. The machine has 96GB of memory.

Do you already have some idea of where all this memory is going or should I investigate this further?

Optimizing gives some good results on the model, but also has the memory problem.

Thanks in advance,

Tom

[Bob Carpenter]

For the previous thread, I had a discussion with Michael Betancourt
yesterday and he assures me that the integrator should be able to
handle discontinuous derivatives as long as they're not pure delta
functions.

Which version and interface of Stan are you using?

You may not need 1000 sampling iterations for inference, but this
sounds like a problem.

I don't know if the integrator can wind up taking an unbounded
amount of memory --- there aren't any other memory leaks we're aware
of.

- Bob

[tom....@uhasselt.be]

I agree that the integrator should and have no doubts that it does deal with discontinuities. I was mostly wondering about the jumps to infinity which could indicate a problem with the model or the gradients.

I'm using rstan 2.5.0, but can also try another one.
My guess would be a memory leak somewhere in the integrator. I'll have a closer look when I get the time.

Thanks for the help. I'd love to get this working and then maybe start contributing. We're mostly looking at improving the performance of these models through optimisation and parallelism.

Kind regards,
Tom

[tom....@uhasselt.be]

I looked into the memory problem and it seems that there is a leak in the reverse-mode autodiff: I believe the precomputed_gradients_vari's allocated in precomputed_gradients are never released.

My knowledge of the code is at the moment very rudimentary so this might need to be double checked by someone else.

I can continue to look at the problem/code to see if I can come up with a patch if this hasn't been fixed already.

Tom

[Bob Carpenter]

Thanks for tracking this down. Do you have an example (preferably
simplified) that would let us reproduce the problem?

I created an issue and prioritized it for the next release:

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

- Bob