[WIP] Adding iterators

[pwl]

Hi, I have been playing around with iterator version of oderkf and came up with my own implementation. My idea was what follows: make oderkf(F,x0,tspan::Vector) behave as previously and add oderkf(F,x0,tstart::Number) which creates an iterator. I reimplemented oderkf(F,x0,tspan::Vector) using the iterator version of oderkf, so that both implementations are not redundant. I also made changes to the keyword arguments accepted by oderkf

• New keyword arguments: minstep, maxstep, initstep, I decided that hmin and hmax from the previous implementation were already arbitrary, so why not allow the user to define his own problem specific values? In the previous implementation hmin and hmax were computed from maximal and minimal time, for iterator this is impossible because, in principle, we might not now the maximal time. Therefore, I set the default value of minstep to eps(T)^(1/3), where T is the type of tstart; the default value for initstep is minstep while maxstep=1/minstep. Alternatively we could set maxstep=Inf, but all these values are a matter for a discussion.
• I added a tstop=Inf optional argument to the iterator version of oderkf, naturally oderkf stops the integration exactly when t=tstop is reached.
• I changed the reaction for abs(h)<abs(hmin), now it produces a warining and does not output x (only t and h). Printing x is rather annoying when x is a large vector (the value of x is the last step before abs(h)<abs(hmin) so it is returned by the solver anyway).
• With tstop=Inf it is possible to integrate the system ad infinitum, so I had to add a test for finiteness of the resuts.

The implementation of the non-iterator version of oderkf is simply

function oderkf(F,x0,tspan::Vector,coefs...;args...)
    t0    = tspan[1]
    tstop = tspan[end]
    tout = Array(typeof(t0),0)
    xout = Array(typeof(x0),0)
    for (t,x) in oderkf(F,x0,t0,coefs...;args...,tstop=tstop)
        push!(tout,t)
        push!(xout,x)
    end
    return (tout,xout)
end

You can also use collect(ode45((t,x)->x, 1., 0., tstop = 10)).

You can now do something similar to what @ivarne had in mind in issue #27 :

ode = ode45((t,x)->x, 1., 0.)
for (t,y) in ode
    if y > 10
        @show (t,y)
        break
    end
end

With iterators we would be one step closer to implementing root finding for solutions to ODE's. I am curious if iterators could get along with dense output #40 #44 .

PS I know that API is not fixed yet, but I needed iterators for my own purposes and decided to share this draft with you.

Cheers, Paweł

[coveralls]

Coverage decreased (-1.88%) when pulling 7d6265b8ab324d7fbed3334652792519587883dd on pwl:master into 568c21332d32a301901e9ba4719deb3d417d9dff on JuliaLang:master.

[pwl]

This seems to be a hot season for PRs, any thoughts on iterator versions of solvers? I am willing to rebase this PR to the current master if anyone is interested in merging it.

[tomasaschan]

I love this API idea, so I say do it!

In the interest of keeping the API consistent, we probably want to make this style possible in other solvers, but if that seems like a daunting task there's no problem in splitting that into more PR's. Also, note that @acroy implemented the keywords for minstep etc (including a nice estimator for the initial step) in #59 so you might want to base on that branch to avoid conflicts.

[tomasaschan]

This would probably also be interesting for creating a better solution to timholy/ProgressMeter.jl#15, so there are real use cases.

[pwl]

Ok, I think I will wait until #59 is merged and rework this PR after that.

[tomasaschan]

Now that both #59 and #61 are merged, this is probably a good time for this too, if you're up for it :)

[acroy]

OTOH it might make sense to wait for @mauro3's RK implementation. Otherwise the iterator version will also count as "derived" from oderkf.

[tomasaschan]

Good point.

[pwl]

I agree with @acroy, I will wait and build on top of @mauro3's work.

[mauro3]

The pressure is on...

[berceanu]

OK, so I want to help build the iterator version, because I really need it at this point. Are the any outstanding issues with @mauro3's PR or can we merge that and move on? :)

[mauro3]

From my side, it's ready to be merged.

An alternative to your problem in #71 could also be to introduce a poststep function which gets executed after each step. This could be useful for diagnostics and such as well. Less general than iterators but probably a lot easier to program.

[berceanu]

It's not just for the normalization problem that I need it. It is also useful for visualizing the convergence, so that I can see for example when the steady state is reached.

[pwl]

If @mauro3 PR is merged in time I will be able to take a look at the iterators this weekend and rework them to fit the new solvers.

[berceanu]

That would be great! Please, let me know how I can be of assistance.

[acroy]

I guess poststep is a simple example of an event function and would be very useful. The question is how to realize the callback?

[berceanu]

@pwl any progress on implementing the iterator version on top of @mauro3 's new RK code? :)

[pwl]

@berceanu Sorry for not keeping you up to date. I keep this issue at the back of my head constantly but the last week was pretty busy and currently I am at a workshop (which lasts until Friday). I gave the code a look some time ago and implementing iterators will take considerably more time than I anticipated. The reason is I think it is necessary to unify fixed step and adaptive methods to avoid wasting resources into supporting two separate iterators, and this would be a major rework of the existing source code. From what portions of code I understand it seems possible, but as I said, it might take me some time to dig into the algorithms.

If you are in a hurry and you need the iterators right now feel free to use my fork, I have been using this implementation for quite some time and it works flawlessly (although some new functionality like points=:specified or custom types are not implemented there). The API is described at the beginning of this PR. Also, I don't want to hold you back, so you could give this issue a shot yourself.

[mauro3]

One thing I noted yesterday using DASSL is that the error message for an error in the objective function did not report in which function the error happened and was thus pretty much useless. I suspect because of the task-based iterator. Is there a way around this? Or is there at least an issue in Julia tracking this?

[pwl]

Certainly, that's caused by the coroutines implementation of iterators in DASSL. I don't know if there is any way to fix this, I started a "discussion" on a Google groups [1] but there wasn't much interest in it. I chose to implement iterators via coroutines for DASSL because it was the least invasive method (essentially I just added one produce() statement to the main function) with practically no downsides, apart from the error reporting issues. So coroutines provide a sort of poor man's iterators.

I think that for ODE we should have a proper iterator implemented via specialized start, next and stop routines, that's why it might take a little bit longer to get it right.

[1] https://groups.google.com/forum/#!topic/julia-users/R9yILEz_uno

[mauro3]

Thanks for the clarification (I think you should open an issue with the test case you posted to julia-users). Doing it the other way will need some careful designing of the state variable, but that is probably worth it in the long run.

[mauro3]

Could use or take inspiration from: https://github.com/spencerlyon2/IterationManagers.jl

[pwl]

I re-implemented the iterators basing on the current master, that is I took into account the hue reimplementation of solvers by @mauro3 in PR #68. I tried my best to leave the algorithms intact but this is an initial version and there will be some bugs. I also cut out some of the stuff, which I plan to add gradually back but first I wanted to hear out some opinions. In particular I removed the negative direction of integration (this can be easily added as a wrapper to the current iterator protocol) and I removed most of the stuff regarding types, so at this point iterators probably don't support custom types. On the positive side, the dense output is still supported (as a wrapper to the actual solver)! To construct an iterator (based on e.g. bt_feuler) you can call the ODE.solver function

sol = ODE.solver(f,y0,t0;
   method = ODE.bt_feuler,
   reltol=1e-6,
   points=:specified,
   tspan=(t0.+[1,2,3,4]))

The ODE.solver integrates both the adaptive and fixed step methods and switches between them depending on what you pass to the method keyword. In the end we can use the new iterators as usual:

 for (t,y) in sol
     @show (t,y)
     if t > 3
         break
     end
 end

Once we fix the interface I could work on the removed features (arbitrary scalar-like types and reverse time integration).

As of yet, there is no standardized way of collecting the results but there has been a recent discussion of what type of result we expect from a solver so I decided to wait until that issue is resolved.

And of course tests don't pass at this point.

EDIT: I removed the call method for the tableaus and the need to use DenseProblem explicitly, now all you need is ODE.solver with the method parameter and all the dense output parameters like points=:all/:specified and tspan.

[pwl]

Another update, I added the reverse time integration and some simple root finding option (via already existing Hermite interpolation in dense output). The root finding uses a very simple bisection algorithm. Here is an example summarizing the new features

import ODE

f(t,x)=[x[2],-x[1]]
y0  = [0.0, 1.0]
t0  = 0.0
tspan = -collect(1:5)
stopevent(t,y) = y[1] > 0.0

prob = ODE.solver(f,y0,t0;
                  reltol = 1e-7,
                  abstol = 1e-7,
                  tstop = -Inf,
                  method = ODE.bt_feuler,
                  tspan = tspan,
                  points = :specified,
                  stopevent = stopevent,
                  roottol = 1e-10)

for (t,y) in prob
    @show (t,y)
end

and the result is

(t,y) = (0.0,[0.0,1.0])
(t,y) = (-1.0,[-0.8414842895974096,0.5403108491723219])
(t,y) = (-2.0,[-0.909326182080129,-0.4161599958706816])
(t,y) = (-3.0,[-0.14112670311715536,-0.9900394570025749])
(t,y) = (-3.141592654651042,[3.9867758613623664e-11,-1.0000496742146774])

The tspan together with :specified ensure that the results will be returned at times t=-1,-2,-3,-4,-5 (plus the initial time t=0, we should decide whether this is the default behavior) but before we reach time t=-4 the root finder triggers and finds the root (the function stopevent returns true). The tolerance for the root finder is given by roottol, we find that the root is t=-pi+/-10e-9 and this is of course true because the solution is y[1]=sin(t). There are still some rough edges, for example the tstop and tspan are somewhat redundant. Apart from that the support for scalar-like types is not yet there.

And there is a bonus as well: without any modifications we can use the algorithms from DASSL! (although we need a wrapper to conform to the API of DASSL)

using DASSL
using Iterators

dassl_wrapper(f,y0,t0;kargs...)=imap(x->(x[1],x[2]),
                                     dasslIterator((t,y,dy)->f(t,y)-dy, y0, t0; kargs...))

prob = ODE.solver(f,y0,t0;
                  (...)
                  method = dassl_wrapper,
                  (...)

[pwl]

@tlycken @mauro3 would you mind taking a look at my implementation of the iterators? Holiday season is coming so I might have some time to add the tests and optimize it a bit. It could also be an opportunity to modify the output according to #82 or #80.

[mauro3]

Sorry to not have looked at this yet! I'll try over the next few days, but I suspect I may well not be able to give much input until mid-January.

[tomasaschan]

I took a quick look, and although I won't have time to prod deeper, at least on the surface it looks quite good to me.

Do we have usage docs for this library somewhere, that should be updated?

[mauro3]

I had a brief look over it: looks good. A few comments in-line, probably redundant as this is still WIP.

To test performance, maybe use https://github.com/mauro3/IVPTestSuite.jl: Threebody and HIRES should work for the explicit solvers. (I should get that package cleaned up to be more user friendly and more directly integrated...)

Random thoughts/questions:

• would multi step solvers fit into this scheme?
• would this integrate with the Sundials.jl solvers as well? It would be cool to eventually be able to call those through ODE.jl.
• could the start and done be generic and a particular class of solvers would only need to implement next?
• This design assumes that dense output can be attached to the solver afterwards. I think, for more sophisticated (higher order) dense output, the solvers and the dense-output are often more tightly coupled. Not something to implement now but something to ponder.
• there seems to be some overlap between the types DenseProblem and Problem, maybe could be refactored?
• eventually, I would like to see in-place functions to be used in the solvers (see e.g. #72). Would this work with this?

I'll look at it in more detail next year... thanks for the work!

[mauro3]

Actually, above points about "sophisticated dense output" and multistep methods may need the same solution: a state which holds more than just the last output but either several or intermediate steps. Pawel, as your DASSL code is multistep, maybe implementing that within this framework would lead to valuable insights on a general design (as well as more work!)?

[pwl]

Hi, I made some major updates to the PR and refactored the code loosely basing on the ideas by @mauro3. The changes include

• Support for the legacy(?) interfaces like ODE.odeXX, so most of the tests pass (but not all, there might be something wrong with the adaptive steppers).
• Complete rewrite of the lower level structures, we now have ExplicitODE to hold the ode data, Options to keep all the options in one place and anAbstractStepper with the algorithms and the associated iterator structure. Everything is held together by the Solution{AbstractStepper} type. Again, the names can be changed if their meaning is not clear enough.
• Given ode=ExplicitODE(...) and stepper=SomeStepper(...) one can construct an iterator by calling sol=solve(ode, stepper; kargs...)
• The modified Rosenbrock method ode23s is now supported
• The dense output wrapper is now a DenseStepper and can be used as any other stepper, to get the dense output use dense(sol) with sol as above. The dense stepper takes its options from the Options component of the Solution
• All implemented algorithms support the in-place formulation (i.e. I switched from dy=F(t,y) to F!(t,y,dy)). In fact this is the default way of defining an ODE now and the higher level interface converts F to F!, and does the same for the Jacobian.
• Removed the ugly implementation of call for the Tableau.

Future plans

• Implement the multistep methods (I think this is the only remaining unimplemented big chunk of the existing code).
• Add tests for the lower level interface and the iterator protocol.
• Incorporate DASSL?
• Change the default output type of the higher level interface to an array of tuples [(t_n,y_n)] instead of the current tuple of arrays ([t_n],[y_n]), this way we could get rid of the useless conversions of the results coming from the iterator.
• So far all the dense output uses a hard-coded 3rd order interpolation, as @mauro3 suggested some steppers could provide a more accurate and efficient implementations.
• Change the top-level interface so that it can accept the in place F! and jac!? Maybe have ODE.odeXX!(...) and ODE.odeXX(...).

[pwl]

I'm pretty sure there is a bug in the RK stepper (I don't get any convergence with bt_dopri5), if someone spots it please let me know.

EDIT: Found the culprit, in my tests I didn't account for overshooting. Fixed now.

[mauro3]

Cool! I'll look at this Friday in a week (currently I'm in a pre-conference frenzy).

[mauro3]

Thanks @pwl, this is coming together nicely! Big PR, big review ;-) I didn't run much code, apart from the tests. So this is just a reading-the-code review.

In fact, can you add an example to example/ on how to use all levels of solver API? That would also make reviewing a lot easier as I could just run those.

My internet was patchy when reviewing this (on the train). So instead I made the in-line comments in the code and did a PR against @pwl's fork: https://github.com/pwl/ODE.jl/pull/1/files . I hope this is ok.

We need to solve the problem with scalar types! Otherwise this PR will be to big a regression. I for one would be happy to make scalar types not a performance priority. But I would hate to see them go. What about just wrapping it in an array, why does that not work anymore? This is what used to be in runge-kutta.jl:

function oderk_fixed(fn, y0, tspan, btab::TableauRKExplicit)
    # Non-arrays y0 treat as scalar
    fn_(t, y) = [fn(t, y[1])]
    t,y = oderk_fixed(fn_, [y0], tspan, btab)
    return t, vcat_nosplat(y)
end

Does that not work anymore?

What are the proposed semantics of the returned values in the iteration? From a performance point of view, it would make sense that y is what is contained in the state itself and then would need copying before storage. This wold allow preallocated storage to be used effectively:

yout = similar(y0, length(y0), length(tspan))
tout = similar(y0, length(tspan))
for (i,(t,y)) in enumerate(oderkf(F,y0,t0,coefs...;args...,tstop=tstop))
    tout[i] = t
    xout[:,i] = x
end

I think this is what is currently implemented, but we should probably make a conscious decision on that.

This needs benchmarking. I or a GSoC student will resurrect the IVPtestsuite.jl and try to run this on there. But I'm mostly away from the computer for about three weeks, so it might be after that.

make_consistent_types is not used anymore, is that on the todo or to be scrapped?

There are many types with fields which have some non-concrete type. This makes me a bit nervous as this can introduce type instability, however, it may not matter. Maybe comment where you use non-concrete that their use there are justified?

Including DASSL would be great, but leave that to another PR, this is huge already.

I'm not a fan of the spaces around :: (in particular not in function signatures) and is not usual julia style. Maybe change that?

I think (now) that copy!(x1,x2) is clearer than x1[:] = x2. Maybe it's even faster.

I don't see why the meta-programming in interfaces.jl is needed. I think this can be done otherwise.

Do you think it would be possible to refactor the adaptive time stepper (i.e. the next-method) such that the same can be used for ode23s and RK? No need to do this now, just wondering.

Documentation is lacking, but I'm sure you're aware of that. Would make reviewing easier though!

In general, this looks good, but probably needs a few more iterations until it's ready. Which brings me to the planning: I'll be around from mid May onwards and should be able to review code within a few days of posting it. If you got time then too, we could try and make some progress on this.

Another thing on the planning side: there is going to be a GSoC student (@obiajulu) working on ODE.jl this summer (mid May onward). It would be great if this gets to a state that the student could base his work on this PR from sometime June onwards. It is probably unlikely that this PR will be merged in this time-frame but if it starts to stabilize the student could base his work off your branch.

Tnx for all the work!

[pwl]

The comments look as thorough as the PR itself:-). I will take a look at them today and get back to you.

[mauro3]

Cool! Well, I'm off on holidays now, so maybe no response from me for three weeks.

[pwl]

The top level interface (interfaces.jl) is still a mess, but it looks slightly better now. I really dislike the way the initstep is handled, just having it there requires a lot of additional code in other places. Any ideas how to improve it?

@mauro3 The scalars were never left out, it's just that the lower-level interface doesn't support them (steppers). Perhaps there is a way to implement it without changing the steppers but it would probably be another wrapper. And I don't think someone willing to use the lower-level API would be disappointed to see that the scalars are not supported. But if there is a nice solution I'm all for it.

As for the make_consistent_types I don't think it's necessary any more, I would just convert the Butcher tableau to T=typeof(t0).

[mauro3]

Tnx for all the work. I don't quite see the ugliness of the initstep but maybe I'm missing something. One thing to add is that dtinit should probably dispatch on the stepper type too.

Good to hear that the scalars are still working, sorry about the confusion.

Concerning the make_consistent_types: you're right, I probably over-engineered that. There are two issues:

(1) question is whether we need to provide for a situation where typeof(y)!=typeof(y./t). (This is suggested in the API docs.) But I can't think of any time where that would be the case. (And this is what's in make_consistent_types.)

(2) One reason to make the make_consistent_types function was to allow calls like ode((t,y)->2t, 0., [0,1]) where eltype(tspan)==Int. This should be handled graciously by the at least the high-level interface (probably should add a test for this). I think this PR may fail on this, but I haven't checked. Edit: see PR #96

Anyway, maybe make_consistent_types is a bit over the top and could be trimmed.

[pwl]

Now for some reason I can't make it work with BigFloat:

julia> ODE.ode4((t,x)->x,BigFloat[2,3],[1,2])
ERROR: MethodError: `__ode#7__` has no method matching __ode#7__(::Function, ::Float64, ::Float64, ::BigFloat, ::Array{Any,1}, ::Function, ::Array{BigFloat,1}, ::Array{Float64,1}, ::ODE.RKStepper{:fixed,:bt_rk4,Float64})
Closest candidates are:
  __ode#7__{T}(::Function, ::T, ::T, ::T, ::Any, ::Any, ::Any, ::AbstractArray{T,1}, ::ODE.AbstractStepper)
 in ode_conv at /home/pawel/.julia/v0.4/ODE/src/interfaces.jl:82
 in ode4 at /home/pawel/.julia/v0.4/ODE/src/interfaces.jl:70

I have never seen an error like that before (look at the mangled name of the ode function). Something similar was already reported https://github.com/JuliaLang/julia/issues/13919 and it seems it has something to do with varargs. If you use BigFloat with the third argument it hangs for a while and crashes Julia (0.4.5).

[mauro3]

I can reproduce the bug. I suspect it's a Julia bug. Removing the initstep from the signature makes it work again.

function ode{T}(F, y0, tspan::AbstractVector{T}, stepper::AbstractStepper;
                jacobian::Function  = (t,y)->fdjacobian(F, t, y),
                # we need these options explicitly for the dtinit
                reltol::T   = eps(T)^T(1//3)/10,
                abstol::T   = eps(T)^T(1//2)/10,
                kargs...)

    initstep::T = dtinit(F, y0, tspan, reltol, abstol; order=order(stepper))
...

[mauro3]

I'm not sure this is related, are you aware of the BigFloat deepcopy issue: https://github.com/JuliaLang/julia/issues/16667? In fact, if I remove the::T from the initstep:

function ode{T}(F, y0, tspan::AbstractVector{T}, stepper::AbstractStepper;
                jacobian::Function  = (t,y)->fdjacobian(F, t, y),
                # we need these options explicitly for the dtinit
                reltol::T   = eps(T)^T(1//3)/10,
                abstol::T   = eps(T)^T(1//2)/10,
                initstep = dtinit(F, y0, tspan, reltol, abstol; order=order(stepper)),
                kargs...)

then I get a signal (6): Aborted fault and Julia-crash as in above issue.

[mauro3]

Ok, so it is two separate issues:

• if the ::T is left off then it works but hits the BigFloat-deepcopy issue JuliaLang/julia#16667. (tested by replacing all deepcopy with copy)
• if the ::T is used, then there is this odd methods finding error. (which persists with the deepcopy->copy replacement)

[mauro3]

Ok, I found it, it's not a Julia-bug after all. The call dtinit(F, y0, tspan, reltol, abstol; order=order(stepper)) returns a BigFloat but T==Float64, thus it fails. The strange mangled function name comes from internals of the keyword function machinery. Yes, there is this issue: https://github.com/JuliaLang/julia/issues/13599

[pwl]

Ok, I added a conversion to dtinit and it solved the issue with y0, but having a BigFloat tspan still crashes Julia 0.4.5.

julia> ODE.ode4((t,x)->x,[2,3],BigFloat[1,2])
(BigFloat[1.000000000000000000000000000000000000000000000000000000000000000000000000000000,get_str.c:153: MPFR assertion failed: size_s1 >= m

I guess we should leave it as is and wait for the backport of the fix. I don't have julia 0.5, but I guess it should work with https://github.com/JuliaLang/julia/pull/16999.

[pwl]

I also get a crash from

julia> ODE.ode4((t,y)->2y, Number[1,1.1,BigInt(1)], [0,1])

[pwl]

To me this "feature" seems to be causing more trouble than good and interfaces.jl now looks like a patchwork and is difficult to follow. How about we simply print an error if somebody tries to use Array{Integer} or Integer? Seems much simpler to me than trying to figure out what the user had in mind. And it's trivial to fix on the user side. On the other hand if somebody tries to use Array{Any} or something similar we just print a warning and let them do whatever they like.

EDIT: we could add a function isexact to test against Integer and others (perhaps Rational).

[mauro3]

But ODE.ode4((t,y)->2y, Any[1,1.1,BigInt(1)], [0,1]) crashes too. How about just making the warning an error?

[pwl]

So we end up supporting only leaf types and arrays of thereof? I'm all for it. But then we can perhaps do it with method dispatch. Is there a type test that would filter these out on dispatch?

EDIT: how about you set up a gitter channel for ODE?

[mauro3]

No, I don't think so. But the check is simple and can give an informative error.

[mauro3]

BTW, are you planning to work on this during tomorrow's hackathon? I could be around (remote) from about 2pm-5pm.

[pwl]

That's the plan, although I believe that this is the last serious issue left out to fix and it's almost fixed now (apart from putting back the old methods). I was also planning to work with @obiajulu to port his algorithm to the new API.