[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 7d6265b on pwl:master into 568c213 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.

[ChrisRackauckas]

Hey, talked with @pwl last night in depth about this change. I think that, when it comes time, we should

1. Tag the master before merging this.
2. Fork the master to another package (LegacyODE.jl?)
3. Pull this in.
4. Give it a some time to be thoroughly tested.

The reason for 2. is because the current source code has good educational purposes and is very approachable. It's pretty much what people would think the standard implementation for ODE solvers would be. On the other hand, while there is a legacy API to make it act the same (though the idea is to deprecate this: JuliaAttic/ODE.jl#25), the iterator format definitely has a bit of a barrier to entry. And DifferentialEquations.jl has a lot of extra stuff in its algorithms which are related to getting good performance which also makes it not so good in terms of readability (I definitely put performance > elegance whenever there's a tradeoff). So I think going forward that it would be good to keep this (and with the next few updates) master around as a different package which touts the legacy odexx interface and is easily hackable.

Because of this, at first I suggested that maybe this spawn off to another package, but @pwl convinced me that it's a good idea for ODE to take up this change since otherwise ODE will probably lose most of its development activity (i.e. @pwl, @mauro3, and everyone else interested in the iterator API) and its probably best to choose the path where ODE gets the most development.

Thoughts?

[ChrisRackauckas]

On a different note, I think I should mention that this branch does have a bunch of extra tests. Since I choose to wrap the iterator version of ODE from here into DifferentialEquations.jl, I have been following this branch pretty closely and have a good amount of my own tests for it. I also have some benchmarks, but since I never wrapped the current ODE, no regression tests (though if we are going to go forward with the LegacyODE idea, I can wrap the master as well and then we'll have a bunch of regression benchmarks. Otherwise it's not worth doing if that interface is going to disappear).

Remind me to share these when it comes closer to merge time.

[tkelman]

this won't work until final is released

[ChrisRackauckas]

Hey,
I was wondering if we should talk about this becoming a separate package. There's a few reasons for this:

1. This a huge breaking change, essentially a new package. It seems that almost all (if not all) of the code has been modified in some way, and the API is different. No previous code or examples (which there are quite a few examples if you Google around) will work with this library after the change. Although for this kind of library there aren't that many dependents, I think that judging the usage by dependents understates its usage as seen by Googling for examples. Besides, it already has its own specific issues, PRs, etc.
2. This is a large enough change that it basically invalidates the other PRs for this repository. Other people have some interesting contributions which could be added.
3. @tkelman has mentioned a few times that he would like to move away from the name ODE.jl. Going out of the way to ensure this other code essentially becomes ODE.jl is then clearly the wrong direction.
4. pwl already mentioned that he will not be maintaining this in the future. This is a really large change and the API is much more involved, and I'm not as comfortable with replacing all of the code in a public package with one for which the main author is no longer present.
5. This does add dependencies. Not many, but a few. This is something to think about if we want a low-dependency library, of which ODE.jl is the only one.
6. There's no problem with making a new package. We can transfer JuliaODE/ODE.jl into JuliaDiffEq with a new name (ODEIterators.jl?) and register it tonight with no problems, and it could be updated as needed. The only issue with why it needs to be so perfect for the PR is that it's re-writing the entire codebase of an already used public package, so it must cover all loose ends. Making it a new package is also less misleading to users.

The end result of making it a separate package is that, not only would it be easier on @pwl, it would also be less confusing to current users and offer a a completely different design path with different advantages and disadvantages (the upvotes for spinning the current off to LegacyODE.jl shows that there people do see upsides to the current ODE.jl architecture and want to keep it intact).

I'd like a second opinion. @ViralBShah?

[pwl]

The main reasons for replacing ODE.jl is that the new version would be easier to maintain, with less code duplication and more documentation. But I'm not opposed to the idea of making it a separate package as you suggest (see point 6. below).

1. We designed the new version to be backward compatible, there might be some failing edge cases but replacing the current ODE with our ODE shouldn't cause any code to break (at least that's what we were aiming for).
2. I agree on this front. Just a comment: some of these PRs are already implemented in our version, so there is no point to hold on to them too much.
3. I had no idea about that. Does it mean that ODE.jl will be renamed? Where did this discussion take place?
4. Even if I won't have too much time to work on it it doesn't matter that other people wouldn't be able to contribute, especially @mauro3 .
5. The only dependency that we add is ForwardDiff, which any package implementing implicit methods should use in place of finite difference approximations. And it's the first time I hear about the need for a low-dependency library for solving ODEs.
6. That is a completely valid point and an alternative solution to a somewhat stagnant situation we have now. Here is my take on it: instead of replacing ODE.jl we leave it as is but deprecate it, which would encourage people to use ODEIterators.jl. This way we could keep the old package with all the issues, feature requests and PRs and let the users gradually move to a different package. Because the API is backward compatible (to a degree) it should be as easy as changing using ODE with using ODEIterators. What do you think about that? I'm not too sure about deprecation but it might have a positive catalyzing effect and it would speed up the adoption of ODEIterators, it's certainly worth discussing.

[mauro3]

I think, ideally, there should be one main, fully Julia initial-value-problem package (the equivalent of Optim.jl in the JuliaOpt universe). The obvious choice is ODE.jl. However, @ChrisRackauckas is right that this PR has become too big for an easy merge. I would favor a split into LegacyODE.jl and ODE.jl as I think ODE.jl is a good name (if somewhat misleading) and that would make it clear on which package to focus development. IVP.jl would be a more correct alternative but the abbreviation is not as well known. What were @tkelman's objections to the name ODE.jl?

Concerning the points:

1. The idea is that the old syntax would still be working. This has been relaxed a little but 99% of old code should be running with this PR. I just checked the three examples and they run (except that ode45 is not exported thus needs ODE.ode45).

2. True. But if the old API is still valid then the old solver should just be copied and the ported at a later stage.

3. Yep, this is the main issue. It would be necessary to have community support for this change to be worth-while. At the moment this basically means you, @ChrisRackauckas. I think this maybe better discussed in a roadmap issue, I'll put one together in https://github.com/JuliaDiffEq/Roadmap.

4. This adds ForwardDiff.jl, which is ok for me. But I agree dependencies should be kept to a minimum and Julia-only packages.

5. Moving this to JuliaDiffEq would be good. We could just rename it to something that makes it clear that it is in development and postpone the decision about package names to when it becomes operational.

[mauro3]

Also, just having had a look through JuliaDiffEq: is there a point in having both https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl and ODE.jl (or a new incarnation thereof)?

Also @tkelman's reservation about ODE.jl is that acronyms should be avoided: SciML/DifferentialEquations.jl#59 (comment)

Edit: Just having read through SciML/DifferentialEquations.jl#59: I don't think that in the long run there is a point in having both OrdinaryDiffEq.jl+AlgebraicDiffEq.jl and ODE.jl. So, ideally we'd distill the best out of all of those packages and put it into one (or two). Whatever we do, it would be great/essential to have @ChrisRackauckas's backing.

[ChrisRackauckas]

@tkelman spells it out more here: SciML/DifferentialEquations.jl#59 (comment)

abbreviations are usually okay, but ODE and SDE are bad names that probably would be spelled out more if registered today. people not in the field will want to know what the package is for and possibly use it as well. longer names are clearer.

If that's the way all of Julia's naming is going, we should take that into consideration.

ideally, there should be one main, fully Julia initial-value-problem package (the equivalent of Optim.jl in the JuliaOpt universe)

There's a lot of signs pointing that this is the wrong direction. Indeed, even Optim.jl is splitting out functionality things like LsqFit.jl to get everything ease to document and add choices which reduce dependencies. The more recent approach is the individual components + metapackages approach exemplified by Plots.jl and JuliaML's Learn.jl. This is the direction that DifferentialEquations.jl has gone, following their lead. The reason is because it allows for an "ease of use" package which is fully documented, but then individual components for developers to pull as needed, and allows for dev docs/issues to be more targeted and thus have a lower barrier to entry.

I don't think that in the long run there is a point in having both OrdinaryDiffEq.jl+AlgebraicDiffEq.jl and ODE.jl

Depends. Sure, if you're just looking for one interface that does it all (and wraps everything else), has a bunch of features, and also offers the best possible performance in order to be easily usable for problems like tough PDEs or as a backend to DDE solvers, OrdinaryDiffEq.jl already does that. It also has a unified interface with SDEs, DAEs, PDEs, and in the near future, DDEs. If that's all that one's looking for (which is likely most users), then sure there's no point.

But I do see that there are other things which can be offered. For example, the OrdinaryDiffEq code has a lot of extra stuff in it because it provides adaptive stepsize (PI controlled) stabilization, specific dense outputs for different algorithms to match orders with minimal calculations, autodifferentiation using ForwardDiff and root solving using NLsolve. It also allows for other things like using the symbolically calculated inverse Jacobian from ParameterizedFunctions.jl to speed up stiff solvers. This adds complexity and dependencies. On the otherhand, ODE.jl is essentially dependency-free and very simple code. This PR has an interesting new iterator interface at the expense of much more complex code (at least to me, I find this much harder to follow). Of course, I am wrapping it all in OrdinaryDiffEq anyways so it is all accessible via one package and one function (currently only this PR is wrapped, but if they end up as separate packages with different implementations I'll extend the wrapper to both since that would take less than an hour), but it would offer the choice of using a separate component if necessary.

While I don't value the iterator interface myself, I think it's an interesting idea. But in the end does not get anything new (at least for diffeq solvers). It introduces a lot of code complexity in order to allow the interface, but I don't think that interface is very usable itself. While the premise is that you can "inject code", there's not enough available and so callback functions actually end up having more functionality for declaring things like events. If the algorithms have a lot of specific codes (specific dense output functions, differing amounts of caches, etc.), then using the iterator interface to something non-trivial like dynamically resizing the ODE is very difficult to do and is algorithm-specific (you change the integrator and then you have to change the code in the loop part?), so it's better to write a DSL and some macros anyways. So in the end I think I've shown that the callback + event DSL approach is more performant and easier to use. So I don't plan on adding an iterator interface to OrdinaryDiffEq, so it's perfectly fine to have another setup which has it available (and maybe someone will really show something that you can do using this interface which you can't do with a callback, proving me wrong). I will still wrap whatever other methods are out there, even if it doesn't use the full iterator format (and, if it's useful, can probably put a type of callback function in the iterator part, so the same functionality is available).

So, ideally we'd distill the best out of all of those packages and put it into one (or two)

I did pull everything I could find which was salvageable from ODE.jl already. For example, the Rosenbrock implementations served as good inspiration. However, other things like ode78 I don't see much use for (it's a Fehlberg method, and there are much more efficient methods these days). The benchmarks show that I pulled the best out because, well, if they didn't I'd keep working on them haha (they do show that OrdinaryDiffEq's native algorithms do have about a millisecond extra of startup time over the ODE.jl algorithms due to dynamic dispatch on the alg symbol. So if a problem is solved in less than milliseconds, you can measure ODE.jl as faster (sometimes, since there's a lot of jitter here), but other than that the ODE.jl algorithms are slower (this is measured against this PR)). That's why I personally am going the way of just working on OrdinaryDiffEq, modularizing, and adding dev docs to help pull other contributors in. I think it's much easier to get OrdinaryDiffEq "more public" than it is to get ODE/ODEIterators to have the functionality of OrdinaryDiffEq.

However, given that its structure can easily call outside packages, there's no limitation that others have to work directly inside of it as well. It's allows for a "work as you wish" strategy. I don't really see a problem with having some algorithms in a package which is separate from the main solver, as long as it's in JuliaDiffEq so that way it can be maintained (fix depwarns, change which versions of testing, etc., even if I'm not working on it directly, I like to keep it all up-to-date and working)