dense output support

[berceanu]

I would like to know if dense output, in particular for ode45, is in the works or not.
If not, I could give it a try, however

1. i do not know if the current codebase will not change completely due to the GPL licensing (what are the plans?)
2. i do not know if dense output can be integrated nicely into the existing routines, or a rewrite is anyhow in order

[ivarne]

What do you mean by dense output?

If you want a matrix of the solution vectors you can use hcat(sol...)

[acroy]

@berceanu : AFAIK no one is working on dense output, but it would be a nice feature. It might be that we will replace the current oderkf by something evolving from #33, but I would think that a dense output routine should be very similar for both cases. Eventually, we probably want to have a separate DOPRI5 routine (which might find its way into Base and for which dense output is very efficient) and a routine like the current oderkf, which can handle different RK-triples.

@ivarne : Dense output refers to interpolation methods which give the solution between two integration steps. For example, within DOPRI5 this can be done without additional evaluations of the RHS.

[berceanu]

So are you suggesting to implement a stand-along dense output routine, which could then be called by the current oderkf and the future DOPRI5? Or will it be part of the DOPRI5 (like in the Hairer code)?

[acroy]

For DOPRI5 it might make sense to integrate the dense output into the routine, but since we don't have a dedicated DOPRI5 solver now, I would start with something that can be used with oderkf. (If DOPRI5 makes it into Base we would anyways need separate implementations.)

[tomasaschan]

(Isn't there a long-time goal to get all the functionality in this package into Base, once the API is stable and the code is well-tested?)

[ivarne]

It should at least have a fair chance of being included in the standard distribution as a default package.

[berceanu]

@acroy I guess that by "something that can be used with oderkf" you mean changing oderkf itself, since the dense output would need partial results from the runge kutta integration at each time step.

[acroy]

Yes, I guess you will have to modify oderkf. Maybe one can nevertheless put
most of the code into a separate function? And it should ideally work for
different RK-triples.

On Friday, June 20, 2014, berceanu [email protected] wrote:

@acroy https://github.com/acroy I guess that by "something that can be
used with oderkf" you mean changing oderkf itself, since the dense output
would need partial results from the runge kutta integration at each time
step.

—
Reply to this email directly or view it on GitHub
#40 (comment).

[berceanu]

I wrote some notes on implementing dense output here: http://nbviewer.ipython.org/github/berceanu/notebooks/blob/master/julia/dense_output.ipynb
Unfortunately as there is no support for MathJax in Github (at least as far as I know), I could not do a PR on the wiki. Anyway, looking forward for your feedback.

[berceanu]

For the API details on dense output, see #44

[tomasaschan]

Having read your notebook again, I'm beginning to get the feeling that a solver with dense output is something a user would want when the workflow is to first solve the ODE problem, and then find the solution at some specific points, which are not necessarily known at the time of solving. In other words, the output from a solver with dense output support, whatever type of object(s) we return, somehow needs to support a posteriori evaluation anywhere on the timeline.

Since this is extra work which we don't want to do if the user doesn't explicitly ask for it, I think we should probably do it entirely separate from the current oderkf implementation. (We want the dense output to be calculated "along the way", and the solution returned as a "dense type", which will break type stability).

Since dopri5 does intrinsically support this, maybe we could start by adding an implementation of that, with some smart output format? Then we could generalize from there, e.g. to other RK methods and/or cases that need extra rhs evaluations.

[berceanu]

dopri5 is nothing but our beloved oderkf using the Dormand-Prince set of coefficients, and with dense output integrated. i see no convincing reason to port it at this point.

about the usage cases, dense output can for instance be useful in case of an event function. imagine one want to integrate until one has g(t, x(t)) = 0 for some user-defined g. an example could be, when solving an periodic astrophysics problem, stopping after one complete period.
for me personally, its useful because i would like to do an FFT of the returned x(t) values, so they have to be on an evenly spaced grid in time. also, i do not want to store all the xout values, for each step, since they will be complex matrices of dimension 2x1024x1024.

about breaking type stability, im not sure i get why that would be. for theta = 1 (using the notation in my notebook) one recoveres the original time-steps taken by the solver. for theta between 0 and 1 one gets more output, in between those original points, but the return type would always be the same as the type of x0. am i missing something?

[tomasaschan]

I might have misunderstood how dense output usually works, but my impression of how it is as useful as possible, is if it is possible to do something like this:

tout, yout = odedense(F, x0, tspan)

and then evaluate the solution in any point (t,y) simply by evaluating yout[t], for any t between extrema(tspan), and not just for any t in tout (or, as it is today, an index). For example, this means that I can first solve the problem, then examine some part of the solution, and then decide to evaluate the solution at some other place, which was impossible to determine before the call to the solver.

If this functionality is to go in oderkf, one of two things must happen: either we always calculate the interpolation for dense output, and return a type that supports it (meaning extra work and lost performance for most use cases) or we calculate dense output only when the user has explicitly requested it, but then we have to sacrifice type stability, since regular arrays can never support dense output as I've outlined it above.

If, on the other hand, the API is to be more or less the same as for some step-limiting dense output-like feature, requiring the user to specify all the points for dense output up front when calling the solver, then this won't be a problem. But I imagine dense output is much less useful that way...

[berceanu]

Ok, so what I had in mind was an API like matlab's, see in particular the refine keywork: http://stackoverflow.com/questions/20898040/using-refine-option-in-matlabs-ode45
That would mean always calculating the interpolation, with an extra cost in performance. How much this cost is remains to be seen in practice, but if we limit ourselves to the case when no extra steps (evaluations of the rhs) are needed, then it should not be significant.

I agree that having the possibility of evaluating the solution at any point in the given interval would be very useful, but to my undestanding neither Hairer's dopri5 or matlab's ode45 provide that.
How would one implement it? Suppose i do tout, xout = odedense(F, x0, [0:1]). Then xout can no longer be an array-like type, but a function of t, right? However, the interpolating polynomials in the dense output interpolate the solution between two consecutive time steps, and not on the whole integration interval.

[tomasaschan]

I would implement it by letting xout (we should really decide if we're solving x' = F(t,x) or y' = F(t,y), by the way - it gets so confusing when one throws in the y' = F(x,y) option...) be of a special type, e.g. ODEDenseSolution <: ODESolution. This type would save all the interpolation coefficients for the piecewise interpolations between solution points, as well as the time and function values at the points (i.e. tout and xout in current oderkf). By implementing a method for getindex(s::ODEDenseSolution, i::FloatingPoint) we could allow the same syntax as today to access individual function values - x_at_t_equals_3_point_6 = xout[3.6] (given that tspan[1] == 1 - we'd have to discuss the exact mapping used). Returning the value from the arbitrary location would be a simple matter of 1) finding the two solution points which the requested points lie between, and then 2) use the interpolation coefficients for the corresponding interval to evaluate the solution.

This is admittedly a lot more work than just adding dense output to the existing function, but I do think it adds quite a lot of value. It's very possible that we could make ODEProblem <: AbstractArray, to make the solution object work in virtually every context where the current code works.

[berceanu]

I gave your suggestion a bit more thought. Let's say we want the solution for any time tau in the interval tspan.

We define

type MyDense
    t::FloatingPoint
    h::FloatingPoint
    x::typeof(x0)
    k::AbstractVector{typeof(x0)}
end

And the function odedense would be doing this

function odedense(F, x0, tspan)
    ..same as oderkf but store also the steps taken and the coefficients k..
    ..return an AbstractVector{MyDense} of length tot_no_of_steps..
end

We call it to get the array of times, solutions, steps and coefficients xout = odedense(F, x0, tspan). Finally to get the solution at the specified time, we need something like

function getsol(tau::FloatingPoint, xout::AbstractVector{MyDense})
    ..calculate step number n and required theta to get to tau, knowing that tau = t_n + theta*h_n
    ..call interpoly(xout[n], theta)..
end

where

function interpoly(sol::MyDense, theta::float in interval(0,1))
    ..return x at time sol.t + theta*sol.h using the notebook formula for x^d_{n + theta}
end

[berceanu]

Sorry for the overlap between our last 2 comments, I posted mine before I saw yours. I guess they converge somehow?

[tomasaschan]

@berceanu Yeah, mostly =) The only difference being the exact API - I was thinking to mimic how Grid.jl does it -but that's not my prime concern at the moment. First we need to decide a few other things:

• Do we want the possibility to do dense output even from "ordinary" oderkf, requiring that the user specifies the dense t-grid in advance and filling it in along the way?
• If so, is there a clean way to save the computation time by skipping dense output when the user does not request it? Would it make a difference?
• Do we, regardless of what we decide for dense output of oderkf, want a function e.g. odedense which returns something that behaves like an interpolation object?

[berceanu]

I can try to answer for the class of problems I'm typically insterested in. That is relatively long integration times (thousands of steps) of multidimensional problems.

• As soon as the dimensionality of x0 (I call the unknown function x because this is how its called in oderkf) goes to N=2 or N=3, storing all the coefficients k for each integration step becomes a problem memory-wise. So I do want the possibility to specify a grid to oderkf and have it return solutions at those exact points, in order to do an fft afterwards, for instance.
• This depends on how its implemented. If we implement it in oderkf in the matlab (or dopri5) way, then it can be done without breaking type stability. For the difference in computation time I guess its method and problem dependent, so benchmarking would be required
• Would probably be useful, even though impractical for large problems

[tomasaschan]

OK, I see. I think I'm beginning to lean more and more heavily toward the following approach, which I hinted at in #44:

• The API for points=:specified and points=:dense are exactly the same, except for the keyword argument; in both cases, solution points are given exactly and only for the time values in tspan.
• In case of points=:specified, a check such as this one is done for every entry in tspan, to ensure that the solver actually steps exactly to the specified points.
• In case of points=:dense, an if block is added such that if the last ODE step stepped past a one or more points in tspan, then (and only then) the interpolation coefficients for the past interval is calculated, and solution values for all the specified points in the interval are added to the output.

Thus, the main difference is whether the specified points are obtained by controlling the steps or by interpolation, but the results of the two methods will be almost the same - especially, length(tout) and length(yout) should not change if I switch :dense for :specified or the other way around. (Of course, doing it by interpolation will be much more effective if the grid spacing is much smaller than the optimal step size...)

If we do it this way, both can be easily integrated in oderkf, without the need for any additional odedense.

Does this sound like a reasonable approach?

[ivarne]

That sounds like a reasonable solution.
The only question left is why we even want :specified as an option. :dense (with reduced tolerances, if you need better accuracy) is better for all situations I can imagine.

I think the difference between :specified and :dense might be considered an implementation detail for each solver.

[acroy]

I agree, @tlycken's proposal sounds good. Dense output is mainly (only?) relevant for explicit Runge Kutta methods, so I would say that the general API should only support : specified. (This is at least my impression, for example Sundials doesn't have a dense option). Our RK methods would offer :dense as an additional option which might be more efficient in some cases as Tomas wrote above.

[tomasaschan]

@ivarne I've sometimes used :specified (when implemented as above) as a way to ensure I get the accuracy I need in certain parts of the interval, by simply providing a tspan vector that isn't linearly spaced. So I think there's merit in letting :specified and :dense coexist (for solvers that can support :dense of course...). If it doesn't make sense in your applications, don't use it ;)

@acroy makes a valid point too: any solver can be made to support :specified, but :dense is only relevant for solvers where there are good interpolation schemes. (For solvers where there aren't, the user can just wrap the solution in an interpolation object from Grid.jl afterwards...)