RFC: Duck typing and consistent return types for all solvers.

[acroy]

As discussed in #7 we want to allow more general types for y0. This PR removes the type declarations for the function arguments F, y0 and tspan. All solvers return two arrays of values Vector{eltype(tspan)} and Vector{typeof(y0)}. This will most certainly break existing codes.

There is one issue I would like to mention: the present version of oderosenbrock will IMO only work with y0::AbstractVector or y0::Number (mainly due to jacobi).

Here is a simple example for the new possibilities:

using ODE
H = [[0. 1.]; [1. 0.]];
rho0 = zeros(2,2);
rho0[1,1]=1.;
t,y = ode23((t,y)->(-im*H*y + im*y*H), [0., 2pi], complex(rho0));

using Winston
plot(t,  Float64[real(R[1,1]) for R in y], t, cos(t).^2, "bo")
oplot(t, Float64[real(R[2,2]) for R in y], "--", t, sin(t).^2, "rs")

This should work with all solvers except oderosenbrock.

[tomasaschan]

This looks good to me, although I'm leaning towards preferring consistent naming of arguments inside the methods as well - i.e. letting all methods take F, tspan, y0 rather than x0, and return yout rather than xout. I know it doesn't matter for client code, but it will make the code in ODE.jl slightly cleaner, IMO.

Also, we should probably provide some helper method to extract the different parts of the solution, because that is way too clumsy atm. Basically just a wrapper around the list comprehensions, so we can do getsolution(y, 1) instead of Float64[real(R[1,1] for R in y].

I think this is the way to go, even though it will break stuff - but I have no idea how much this package is used, and thereby no idea how much stuff it will break. But ODE.jl is in an early stage of development, so I don't think a lot of people are using it heavily yet...

[acroy]

@tlycken: I agree, we should have consistent naming throughout the package. I considered to rename variables already with this PR, but I think we should have a more thorough cleanup as a separate PR. Just to mention a few other points:

• Adjust the function descriptions to contain only method specific information (the interface will be the same for all solvers)
• Do we need Polynomial?
• Rename ode23 to ode23_bs (Bogacki and Shampine) and let const ode23 = ode_bs23
• Split ODE.jl into several files, maybe stiff.jl and nonstiff.jl ?
• Consistent error handling in ode23 and oderkf for h <= hmin

Regarding the conversion: I am not sure we should do anything to this effect. If y::T has some general type one might argue that the user (or inventor of T) knows best how to access it, e.g. by providing a custom getindex for Vector{T}. I actually tried a variant where yout::Vector{Vector{T}} was automatically converted to Matrix{T} (via hcat(yout...).'), any other yout::Vector{T} was not converted. In the former case this allows you to access yout as before (like yout[:,1] etc). It is convenient, but might be too confusing, so I took it out.

I have no idea how this is typically done, but we could copy the current master into a branch v0.1 to conserve the "old" behavior (just in case someone cares) and continue in master?

[tomasaschan]

@acroy Yeah, you're probably right - cleanup could (should) be separated into one or more different PR's, so that can wait.

Re. Polynomial.jl: stuff from that package is only used in ode_ms (at the very bottom) which is using an algorithm I'm not familiar with. It's possible we don't need the package to do the work, or that we can rationalize away the entire method - I don't know.

I think the "git way" of doing it would be to tag the current master as v0.1 and keep merging PR's into master. Next time we're in a state we like, we tag it v0.2, etc. Perhaps we should also update METADATA.jl to point to v0.1. (That said, I'm not a maintainer, so I can't do half of this stuff...)

[acroy]

Since we are breaking compatibility anyways, I would actually like to follow @stevengj's suggestion and change the order of the function arguments to F, y0, tspan with this PR.

[tomasaschan]

@acroy Coming from MATLAB, I find that argument order somewhat counterintuitive, but of course people with other backgrounds might find F, tspan, y0 counterintuitive. I tried to read back thruogh #7 to see what the reasoning behind changing the order was, but I didn't really understand if - and if so, why - it's better.

[acroy]

@tlycken : I think NumPy uses F, y0, tspan. One argument was, that F, y0, tspan matches the interface of quadgk, which allows you to use varargs for tspan; something like ode(F, y0, t0, t1, tspan...; keywords...).

Edit: See @stevengj's comment in #7.

[tomasaschan]

@acroy OK, if NumPy does it, it's just as standard as whatever else. I'm still not sold on t0, t1, tspan..., since that would place the times out of order (or am I still misunderstanding this?)

[acroy]

@tlycken : I think I didn't explain it well :-/ In the end the actual solver (oderkf) will use tspan=[t0 t1 tspan...], but the "frontend" ode(F, y0, t0, t1, tspan...; keywords...) allows you to write ode(F,y0,0.,12.) or ode(F,y0,0.,6.,12.) or ode(F,y0,(0.:1.:12.)...).

[tomasaschan]

@acroy OK. I still had to think it through a couple of times before I got it, but basically, we're saying that after we do tspan_impl = [t0 t1 tspan...] we'll have to handle the following cases:

if len(tspan_impl) == 2
    # given just t0 and tfinal
else
    # also given intermediate points
end

and since we require both t0 and t1 in the method signature, we know that tspan_impl will have at least two elements.

I guess I'm just worried that people might look at the method signature (i.e. using methods(ode)) and think they need to give the input as ode(F, y0, 1, 100, (2:99)...) rather than ode(F, y0, (1:100)...).

[acroy]

While implementing the tspan varargs I found that this doesn't work with our current definition of ode4s_s and ode4s_kr, because both can have an additional argument G which denotes the Jacobian of F. I would suggest to use a keyword instead (jac or DF).

[tomasaschan]

@acroy alternatively, place the jacobian right after F. I guess it depends on if we want to force the user to supply it, or just use a finite difference if they don't.

[acroy]

@tlycken : Well, that is how it is done now. The problem is that ode4s_s(F, G, x0, t0, t1, ts...) is not distinguishable from ode4s_s(F, x0, t0, t1, ts...) (and that it doesn't follow our "API specification").

[ivarne]

@acroy If G is a Function and x0 (or preferably y0) is something else, we can actually distinguish the two methods with Julia multiple dispatch. (Sad to raise an argument for the solution I do not like).

If we can require F to be a generic function, we can have:

type derivative{order:<Int} end
F(::derivative {1}, t, y) = ...
F(::derivative {2}, t, y, y_d1) = G(t, y, y_d1)

We also have 4 different automatic methods (symbolic, complex, dual numbers/power series, and finite difference). We can safely pick by trial an error, or keyword argument.

[tomasaschan]

Since we have alternative ways to find the Jacobian, a keyword argument jacobian is probably the best way to go.

[acroy]

@ivarne : Probably I am misunderstanding your suggestion, but the problem here is not that F (the rhs of the ODE) may depend on the Jacobian G, but that some methods need G itself to calculate the new solution. The easiest way out would IMO be to pass G via a keyword, then we would have ode( F, y0, t0, t1, ts...; jac=(t, y)->jacobian(F, t, y), other_kws...) and jacobian could be one of the 4 automatic methods.

[tomasaschan]

I know there have been other names for it, but I think clarity is nice and jacobian is completely unambiguous while still not too long. We could also have a few valid symbol arguments to let the user choose one of the built-in methods, e.g.

odes(F, y0, tspan...; jacobian = :finitedifference)
odes(F, y0, tspan...; jacobian = :symbolic)
odes(F, y0, tspan...; jacobian = :dualnumbers)

If typeof(jacobian) == Function we just call jacobian(F, y, t) (argument order same as for ode) to evaluate it, while for the other cases we can use our own implementations of the various approximation methods.

[ivarne]

@acroy There is some misunderstanding at least. I wrote my last answer on the phone, so I might have better success explaining my suggestion now.

My suggestion was to look at the function F. If F is a generic function, you can actually define multiple methods on F, and part of our specified interface will be which method we invoke on F. My thought was that the normal method invocation of F would be F(t, y). Then if the user wants, he can also define F(t,y,y_d1) to give the jacobian (maybe combined with setting jacobian = :function). In my first proposal I also suggested that we might have a type parameter for documentation, but that is not necessary.

F(t,y) = sin(y.*t)
F(t,y,y_der) = cos(y.*t).*y
t_values, y_values = ode(F, [0. 2.pi], 0., jacobian = :function)

Not sure if this is a good interface (it is definitely not intuitive for many new users), but it might be worth considering.

I like @tlycken's jacobian keyword. It gives a clear and somewhat obvious interface that allows the user a simple way to override our default choices.

[acroy]

I have changed the argument order and ode4s_s and ode4s_kr now use the keyword jacobian. I like the idea of using symbols to switch between different ways to calculate the Jacobian, but since we have only one method now, I just check if the user has provided a Function and otherwise use the (built-in) finite-difference method.

[acroy]

Mhmm, the Travis failure seems to be unrelated?!

[ivarne]

I restarted the build on the travis site. If it is transient, we might have a bug to report to core Julia

[acroy]

@ivarne : It is still failing at the same point with

using ode23
ERROR: error compiling __ode23#1__: unsupported or misplaced expression ... in function __ode23#1__
 in anonymous at no file:22
 in include at boot.jl:238
 in include_from_node1 at loading.jl:114
 in process_options at client.jl:303
 in _start at client.jl:389
at /home/travis/build/JuliaLang/ODE.jl/test/tests.jl:45

Which version does JULIAVERSION="juliareleases" actually give?

[ivarne]

JULIAVERSION="juliareleases" is currently 0.2.1. I tested locally with that version, and it seems to be failing for me also. This is probably a bug in Julia that has been fixed in master.

[tomasaschan]

That still means that if we merge this, we fail to support for Julia 0.2 in the next release of ODE.jl, right? I'm not a fan...

Clarification: I am a fan of getting the functionality from this PR, but not of dropping support for Julia 0.2...

[acroy]

Seems that Julia 0.2 doesn't like [t0 t1 ts...]. I changed it to hcat(t0, t1, ts...).

[tomasaschan]

@acroy Nice job finding that! =) hcat seems to do what we want for sensible inputs - I tried the following quick, and everything did what I expected:

f(t0, t1, ts...) = hcat(t0, t1, ts...)

f(0, 1)
f(0, 1, 2)
f(0, 1, 2, 3)
f((0:4)...)

For the last one, I first tried f(0:4) which of course didn't work. But I'm thinking, maybe it would be a good idea to provide an overload ode(F, y0, trange::Range1; kwargs...) = ode(F, y0, trange...; kwargs...) for convenience? The only difference that makes, is that you don't have to manually wrap ranges in ( )... to use them as tspans. (Of course, doesn't have to be part of this PR!)

[ivarne]

@tlycken That last part is just going to be confusing. Why should not ode(F, y0, 1:2:7) or the new float ranges ode(F, y0, 0:0.1:1) work, and when should we stop adding more convenience features that actually look like ode(F, y0, tspan)?

[tomasaschan]

@ivarne Sorry, consider the 1 in Range1 as a typo. My point was to add something that would work on all (1-dimensional and discrete) ranges, regardless of their type. But if there isn't a base class that catches all that behavior, I agree it's not worth it.

[ivarne]

Sorry for picking on a detail. We can definitely try to collect range types, but that does not solve the consistency problem (which I failed to illustrate).

Why should ode(F, y0, 1:2:7) work, but ode(F, y0, [1:2:7]) fail? If we accept 1 kind of iterable, we should accept any iterable!

[tomasaschan]

@ivarne It's good that you pick on these details - the two possible outcomes is that either a) we find a way to do this that doesn't cause problems later, or b) we find that it's not possible, and avoid causing problems now by simply not implementing it. So I'm all for the nit-picking! =)

Wouldn't it be possible to duck-type it, though?

f(t0, t1, ts...) = hcat(t0,t1,ts...)
f(ts) = length(ts) > 2 ? f(ts[1], t[2], t[3:end]...) : f(ts[1], ts[2])

f(1,5)
#  1x5 Array{Int64,2}:
#   1  2  3  4  5
f(1,2,(3:5)...)
#  1x5 Array{Int64,2}:
#   1  2  3  4  5
f((1:5)...)
# 1x5 Array{Int64,2}:
#  1  2  3  4  5
f(1:5)
# 1x5 Array{Int64,2}:
#  1  2  3  4  5
f([1:5])
# 1x5 Array{Int64,2}:
#  1  2  3  4  5

If f is called with a non-expanded ts as the only argument, and ts is something that can't be expanded into at least two parts (ts[1] and ts[2]), we get a BoundsError (which we could of course wrap with something more instructive, but at least it doesn't try to hcat). Multiple dispatch works because in the case of a range, array or whatever iterable we pass for timespan, it is one argument less than the default f(t0, t1, ts...).

I'm still not saying this is a need-to-have, just that it's a nice-to-have that we shouldn't dismiss without good reason. Feel free to nit-pick once again =)

[ivarne]

@tlycken Thanks, I think you managed to answer my nitpick concern. An iterable is only one argument, but the range is start, [middle points], end.

Efficiency is also a concern. I tried to benchmark your f(ts) against y(ts) = ts' on ts = [1:5000] and the performance difference was not insignificant. It will probably not matter much when we combine it with a ODE solver, but I do not understand the value of deviating from matlab's array tspan with t..., if we are going to support any iterable anyway.

[acroy]

I don't see a problem in providing ode(F, y0, tspan; kws...) and ode(F, y0, t0, t1, ts...; kws...). The "backend" functions (like oderkf) anyhow use tspan.

[acroy]

I have added variants odeX(F, y0, tspan; kws...) for the "frontends". Should we (try to) address #26 with this PR or do it separately?

[ivarne]

I think we should try to push things forward and try to change things so that things become more like the situation we want. We have tagged a release and the next one is announced to be breaking. I opened #26 in a separate because I think it deserves a separate discussion and solution.

[tomasaschan]

I agree. There are lots of things we want to fix and change, but there's no reason to try to do them all in one PR. Better merge this, so work on other things can be based on this - especially since this PR makes changes to the signatures of almost all our methods.

[acroy]

Just noticed that with the nightlies Travis takes 1:20 min longer ... This might be connected to #29.

[ivarne]

@acroy That is a lot! It might also be connected to just Julia taking longer to install. The new sys.so file that contains precompilled parts of the standard library needs to be generated on install.

Edit: Also, I do not think Travis runs only one test on each machine at once, so the timings might be quite random.

[acroy]

@ivarne: I also see a (dramatic) slowdown on my machine compared to some earlier version of julia. I suspect it is connected to the issue you referenced in #29, which effects generic matmul (that is my guess at least).

[acroy]

I could get rid of some more deprecation warnings. However, we still suffer from JuliaLang/julia#6118 in ode45, ode4s and ode4ms, but I would like to address that separately (at least ode45 in #16, but probably also the others).

[ivarne]

Is this ready to be merged? It looks good to me, and definitely brings us closer to the recent API discussions.

[acroy]

Maybe I can squash the commits into two. Otherwise I think the PR does what we wanted. The side-effects have to be handled separately.

[ivarne]

Now all code using ODE.jl will break when they update.

[acroy]

Yup. Maybe we should take some part of api.md and put it into README.md so that people can see how they have to change their code?

[ivarne]

Yea. We should at least put in a WARNING and a link.