[R-sig-dyn-mod] Skinny square waves (John Harrold)

[Tim Keitt]

I'll pipe-in, because this was the problem that led me to use ODEINT [1]
and eventually create odeinter [2]. I was running very large ode systems
with external driving and it was dramatically slow using solutions I could
come up with in R (mostly interpolation of the external signal). I ended up
using a two-scale solution where I used a piecewise-linear driver updated
at regular intervals and continuous integration between these time points
[3]. To do this efficiently, I needed all steps in compiled code. The
ODEINT library facilitates flexible implementations because it gives you
direct access to integration steps. It also can overshoot and interpolate
back to a given time point. I found if you need to do something
non-standard, it is a very flexible library.

THK


[1] odeint.com
[2] https://github.com/thk686/odeintr
[3]
https://sites.cns.utexas.edu/keittlab/publications/resilience-vs-historical-contingency-microbial-responses-environmental-change


http://www.keittlab.org/

On Thu, Mar 9, 2017 at 10:11 AM, Thomas Petzoldt <
thomas.petzoldt at tu-dresden.de> wrote:

> Hello,
>
> thanks for this interesting discussion. A model that skips external
> forcings was indeed my first problem report regarding "odesolve" to Woody
> Setzer long time ago and a beginning of a great adventure. We know now
> different ways how to handle this, all having their pros (+) and cons (-).
> Here an incomplete list of what we do in our lab:
>
> 1) set times to match external forcings
> - slows down simulation,
> - produce lots of output data
>
> 2) set hmax to limit the internal time step
> - slows down simulation
> + much better than (1)
> + output data still under control (better than 1)
>
> 3) use another solver, e.g. Runge-Kutta
> - less efficient than any of the odepack solvers
> + steps are known before (incl. intermediate steps),
>    so that data can be supplied in correct resolution
>
> 4) make the external signal more smooth, e.g. by using (steep) sigmoidal
> instead of rectangular signals.
>
> In all these cases, forcing data needs to be interpolated, either in the
> model function (1-4) or beforehand (3). It is also possible, to pass a
> "hook function" to the derivative function for interpolation. This can
> enhance efficiency for simulations with big amount of external data,
> especially if the interpolator has "knowledge" about the data structure.
>
> It is also good to know that most solvers are able to interpolate the
> states. In case of doubt, we set defaults towards robustness instead of
> maximum efficiency. There are still some options for fine-tuning, and in
> case of doubt, set hmax.
>
> 5) always check precision (atol, rtol) especially if your states occur at
> very small or large scales.
>
> 6) use events!
> + elegant to implement (thanks to Karline)
> - slows down simulation (but not as much as (1) or (2)), because the
> solver does "internal stop-and-go" between integration steps.
> + can be used in combination with 1-3
>
> Events are good if external manipulations are very short, to be considered
> as "instantaneous change of state" while processes are quite slow during
> the rest of the time, e.g. if you want to follow a patient over months that
> gets an 1min injection every morning.
>
> 7) use "external stop-and-go" mode by chaining calls to ode()
> + external forcing data can be kept constant like parameters during each
> step
> + optimum performance of the solver (no intermediate steps)
> - little bit more programming (e.g. loops)
> - some computational overhead, due to repeated invocation of ode()
>
> ... see also the vignette of package "rodeo" about this:
> https://cran.r-project.org/web/packages/rodeo/vignettes/rode
> oVignette.pdf#14
>
> To sum up: R is a full programming language and many roads lead to Rome.
>
> Thomas
>
>
> Am 09.03.2017 um 14:45 schrieb Setzer, Woodrow:
>
>> Can you arrange to code this using deSolve's "events"? You would have to
>> code the model so that it is state variables that change abruptly.
>>
>> Woody
>>
>> -----Original Message-----
>> From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-b
>> ounces at r-project.org] On Behalf Of Karline Soetaert
>> Sent: Thursday, March 09, 2017 2:22 AM
>> To: Special Interest Group for Dynamic Simulation Models in R <
>> r-sig-dynamic-models at r-project.org>
>> Subject: Re: [R-sig-dyn-mod] Skinny square waves (John Harrold)
>>
>> Hi John,
>>
>> Usually in such cases, the problem is that the default solver is too
>> efficient and thus takes too large timesteps. You could either try a less
>> efficient solver that by default takes smaller time steps or restrict the
>> maximal time step (argument hmax).
>>
>> Karline
>>
>> -----Original Message-----
>> From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-b
>> ounces at r-project.org] On Behalf Of John Harrold
>> Sent: donderdag 9 maart 2017 2:56
>> To: Special Interest Group for Dynamic Simulation Models in R <
>> r-sig-dynamic-models at r-project.org>
>> Subject: Re: [R-sig-dyn-mod] Skinny square waves (John Harrold)
>>
>> Howdy Steve,
>>
>> I've run into this problem in the past in Matlab. I tested it in R at the
>> time and I didn't seem to have the same problem. But it seems to be
>> somewhat situational :).
>>
>> I had a similar solution in Matlab to what you did. The R-equivalent of
>> adding events at each of the switching times. So as a test I added events
>> at each of the switching times. Basically I told ode to add 0 to the first
>> state at each of those times, and that seemed to do it. These events should
>> force the ode solver to stop and restart at each time point. Now I just
>> need to automate this in my simulation routines.
>>
>> I think it would be good if we could just provide the ode function a list
>> of times and have it do this automatically. This way it would make it
>> easier for me :). And it would also allow this to be optimized better than
>> my poor coding skills will probably allow. Is this the appropriate list to
>> make these suggestions?
>>
>> John
>>
>> On Wed, Mar 8, 2017 at 8:33 PM, Stephen Paul Ellner <spe2 at cornell.edu>
>> wrote:
>>
>> John,
>>>
>>> the efficient deSolve solvers can't deal with such rapid on-off
>>> changes (even if you make the changes continuous). And even with start
>>> and stop times in the output times, the solvers may well leap over the
>>> 30 minute periods when the forcing is 'on' and give you interpolated
>>> output values for output times during the 'on' periods.
>>>
>>> The way I've always dealt with such situations is to solve the 'on'
>>> and 'off' periods separately.
>>>
>>> for(j in 1:whatever) {
>>>         solve from stop[j] to start[j+1] with forcing=0
>>>         solve from start[j+1] to stop[j+1] with forcing = value }
>>> storing the results along the way in one big data frame or matrix.
>>>
>>> Steve
>>>
>>> Stephen P. Ellner
>>> Department of Ecology & Evolutionary Biology, E339 Corson Hall Cornell
>>> University, Ithaca NY 14853-2701
>>>
>>
>
>> --
>> -------------------------------------
>> John M. Harrold
>> john.m.harrold _at_gmail
>> -------------------------------------
>>
>
> --
> Dr. Thomas Petzoldt
> Technische Universitaet Dresden
> Faculty of Environmental Sciences
> Institute of Hydrobiology
> 01062 Dresden, Germany
>
> E-Mail: thomas.petzoldt at tu-dresden.de
> http://tu-dresden.de/Members/thomas.petzoldt
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>

	[[alternative HTML version deleted]]