RFC: new adaptive (and fixed step) Runge Kutta solver

[mauro3]

As promised in #61, here a RK implementation sans license restrictions. It is a bit longer than the current code but has some more features as well. There is probably a bug lurking somewhere as it is less accurate than the current solvers.

Extra features:

• Tableaus are a datatype. This would allow dispatch on different methods. Also, this "solves" what was tried in PR WIP: Make coefficients Rational for better generics #60: the method will run using the type of number used for the coefficients.
• dense output (using 3rd order polynomials)

Improvement in ODE.jl to have hinit return tdir as well.

TODO:

• figure out why less accurate than current solvers.
• run through IVPtestsuite.jl
• doc update (there is no doc yet)
• fix performance issues
• update to conform to internal of other ODE solvers
• remove old solvers
• squash commits

This touches on issue #9 and fixes #25 (although that was fixed already, no?).

Extra features, maybe in this PR, maybe later:

• Tableaus in terms of Rational` and automatic conversion to suitable type
• iterator protocol
• maybe implement dense output in terms of Polynomials.jl, root finding.

[coveralls]

Coverage decreased (-41.03%) to 55.72% when pulling 8de3ff7 on mauro3:m3/rk into 162c881 on JuliaLang:master.

[coveralls]

Coverage decreased (-6.76%) to 90.0% when pulling a971878 on mauro3:m3/rk into 162c881 on JuliaLang:master.

[tomasaschan]

This will break on 0.4 with the newly introduced tuple changes; at the moment, it should probably be order::(@compat(Tuple{Vararg{Int}})) but I'm unsure on how this works in this context; the outer parenthesis might throw this off.

It will work on 0.3 as is, though :)

[mauro3]

It sure failed on Travis!

[tomasaschan]

Actually, this doesn't solve the problem that #60 addresses, since you've effectively converted all the tables to Float64 anyway in the constructors of TableauRKExplicit (when saving the arrays). It's a good way on its way, though, and the approach is better than what I had in #60, so I'd gladly yield that for this! :)

[mauro3]

Ah, yes, you're right. Another step is needed. But my idea was that the parameter T of the tableau would determine what type is used in the solver. As it doesn't make any sense to have the tabelau in a different type to the numerics. So, we'd need to write it:

bt_rk45_coef = (          [  0           0           0            0         0     0
                             1//4        0           0            0         0     0
                             3//32       9//32       0            0         0     0
                          1932//2197 -7200//2197  7296//2197      0         0     0
                           439//216     -8        3680//513    -845//4104   0     0
                            -8//27       2       -3544//2565   1859//4104 -11//40 0 ],
                           [25//216      0        1408//2565   2197//4104  -1//5  0
                            16//135      0        6656//12825 28561//56430 -9//50 2//55],
                            [0,          1//4,       3//8,       12//13,    1,    1//2])

bt_rk45 = TableauRKExplicit(:fehlberg,(4,5),Float64, bt_rk45_coef...)
bt_rk45_f32 = TableauRKExplicit(:fehlberg,(4,5),Float32, bt_rk45_coef...)

[coveralls]

Coverage decreased (-6.76%) to 90.0% when pulling a971878 on mauro3:m3/rk into 162c881 on JuliaLang:master.

[acroy]

The problem is, that the type of the solution is determined by the type of the time-step (most likely float), the type of the tableau coefficients and the type of F (rhs of the ODE). All three are somehow multiplied to give the new solution. So I am not sure it makes much sense to store the coefficients as Rational ...

Nice work BTW!

[coveralls]

Coverage decreased (-7.06%) to 89.7% when pulling 5123853 on mauro3:m3/rk into 162c881 on JuliaLang:master.

[mauro3]

My idea was that 99% of users will work with Float64. For the other 1%, they can instantiate something like bt_rk45_f32 where T==Float32. As my solver is written, then initial conditions and time vector are, if necessary, converted to that T (actually this should be tested). It would be up to the user to make sure that there is no conversion or no loss in the conversion. So, this is not fully automatic but can be extended if needed. Note that inside bt_rk45, which is used by default, the coefficients are stored as Float64.

And thanks. But, I just had a look at performance with IVPtestSuite and it's atrocious, so still lots to be done...

[stevengj]

const?

[tomasaschan]

The reason to store the coefficients as rational is that they are the least invasive when promoting with other number types. Then, we can do something conceptually equivalent to

TTime = eltype(tspan)
TRHS = typeof(F(0, #= whatever args we need =#))
TTableau = typeof(one(TTime) * one(TRHS) * one(Rational{Int})

and then (possibly) convert the coefficients.

This way, TTableau will automatically be Float32 for a tspan/F combo which returns Float32, remain Rational{Int} for such inputs etc. This will mean one conversion per call to the solver, but we'll avoid all kinds of promotion slowdown in the inner loops.

Interpolations.jl does very similar things to handle typing of coefficients gracefully.

[pwl]

Sorry for nitpicking but I am in the 1% not using Float64. Is this new solver generic in the sense that it would be possible to do computations using BigFloat without loosing precision on the way (via e.g. multiplications by explicit Float64 constants)? So far it kind of looks like Float64 is hard-wired to the coefficient tables now. I could definitely live with initial conversions at the expense slower solver startup.

[tomasaschan]

@mauro3 The reason this won't work is that it's impossible for us to prepare in this way for the case where someone has a self-created number type which doesn't promote cleanly with Float64 - consider for example a type HumongousFloat that has 1024 bits of precision instead of BigFloat's 256. We can't create versions of these tables for such types beforehand, but we can convert to these types on invocation. Storing the coefficients as rationals ensures that it's possible to convert to any concievable number type without loss of precision.

@pwl I don't think the current version of this PR works like that - and it might turn out to be out of scope for this PR to make it work like that. However, I do think it's valuable enough that we should solve this at some point. We'll see where this leads to and take it from there, I guess.

However important generics and sensible promotion rules are, I think the main merit of this PR is to provide a white-room Runge-Kutta implementation, allowing us to greatly simplify the license of this package. That's no small thing, indeed! :)

[pwl]

I just wanted to ensure that the 1% is well represented. :) And I already have my own working fork so there is no hurry in implementing this feature. Don't get me wrong, I see this PR as a remarkable piece of work.

[mauro3]

After fixing a bug in the step control, this now works. Here the performance tests from IVPTestSuites.jl:

The new solvers (*_v2) about a factor 2 faster, memory footprint is also down by a factor 2-5 (not shown). However the slope of the old 78 solver is slightly shallower.

All solvers new solvers a factor 2-5 faster.

And for the fixed step solver:

The solver is significantly better than the old ones. This is a bit surprising as the ode4 and ode4_v2 should be the same alogrithm?!

[mauro3]

@tlycken, yes I agree, that was the idea behind the little code example above, in bt_rk45_coef are the rational coefficients. But maybe we should discuss this in another PR.

[acroy]

Very nice :-)

The solver is significantly better than the old ones. This is a bit surprising as the ode4 and ode4_v2 should be the same alogrithm?!

Well, in the old code many copies of the solution are made, maybe this explains the difference? Partly, this was necessary to make the code work with user-defined types. At the moment, y doesn't have to be <:AbstractArray for example. Is this possible with your code as well?

[acroy]

Its probably better to use transpose here. But what happens if y0 is not an array?

[mauro3]

Changed to .' (I think .' is sugar for transpose, right?)

[mauro3]

@acroy: we should add tests for non-abstractarrays if we want to support them. But what is the use case?

[mauro3]

@acroy: and no, at the moment it doesn't work quite work for those general cases. It also stores the ys in a Matrix internally. I can update that.

Related, I find the current interface returning a Vector of something quite annoying. I can see that for the general case that is needed. But I always have to hcat(y...). Maybe return a Matrix for the normal case and only the Vector of Vectors otherwise?

[mauro3]

@acroy: before updating the PR to allow for more general datatypes, a few questions:

What is the "interface" a typeof(y) has to support? Looking through the code, it looks like just arithmetic. However, considering the speed and memory improvements provided by looping and indexing, it seems a shame to let those improvements slip. So, could I use indexing? Presumably it would be ok to use arrays for intermediate storage?

[acroy]

I agree that the typical use case will involve scalars or dense vectors, but there might be good reasons to have a user-defined type, which doesn't support the full array interface. For example, in JuliaQuantum/QuBase.jl we have QuArrays which are not <:AbstractArray. Also AFAIK @Jutho's tensors are not abstract arrays. Nevertheless, we definitely want to use ODE.jl to solve some differential equations.

Related, I find the current interface returning a Vector of something quite annoying.

Yeah, we discussed this when we put together the API specification (probably in #20). I think one point was, that one wants to use yout[end] as new input or to evaluate F. We could handle the array-input case separately, but this would also break the API consistency.

[acroy]

@mauro3: I think it is just norm, +, -, multiplication with/division by scalars. At some point (#37) I came up with an example gist to demonstrate the basics. But I agree that looping and indexing is really helpful, so maybe we should revise the requirements?

[tomasaschan]

I haven't done any profiling, but if memory management is the bottleneck here, we could very probably use the iterator versions (#49) to let the user take care of it in cases where they want to:

• We keep tout, yout = ode45(F, y0, tspan) returning a Vector{typeof(y0)}
• If that is too slow for the user, they can pre-allocate space for their results in a matrix, and use the iterator interface to fill it.
• Internally, we use the iterator interface to fill the Vector, so we don't duplicate code

If we choose that path (and there is a clean-enough way to do it for the user, so we can show it off in examples), we can move that issue out of scope for this PR :)

[mauro3]

I'm trying to adapt the solvers to the preferred API and internals but I find it quite painful and I'm wondering whether this in not a case of over-engineering? There are a few types involved:

Tt = typeof(tspan)
Ty = typeof(y0)
Tf = typeof(F(tspan(1),y0))

Et = eltype(tspan)
Ey = eltype(y0)
Ef = eltype(Tf)

E*: There should be promotion of the eltypes (#26,): Tt should probably be some type of floating point number, for sure a Real (maybe FixedPoint, Rational could maybe work but as @stevengj said, it doesn't make sense), Ty should be some "continuous" Number, Tf also.

T*: The T* are either equal to the eltypes, if scalar, or some kind of container type. In the current solvers Tt is discarded, the outputs are kept in a Vector{Ty} (this is where I experience my pain) and it is assumed Tf==Ty (and thus Ef==Ey).

Anyway, I wonder whether there is any need to support other containers than AbstractArrays? And in particular whether internal, temporary storage arrays can be just Arrays or not. If not, what is the use case?

The reason I am wondering these things is that I would like to see an API developed, similar to Sundials, where all storage can be preallocated and functions do in-place updated. I don't see how this is compatible with above internals. I don't think an iterator method would work, at least not if it is not in-place.

Sorry for the rant. I'll ponder this a bit longer...

[pwl]

The assumption that Tf==Ty is sensible, we always need as many equations as there are variables so it should be always possible to map the return value Tf back into Ty. As for Tt I think it should be an (abstract) vector of Real type (interpreted as either [start,stop] or a sequence of predetermined output times). Can you see any other uses for Tt?

On the other hand the in-place updates are a must have feature if we want good performance in the iterator version and I don't see how this could be realized with custom types of Ty and Tf. Maybe by defining an indexing methods on Ty and preallocating instances of Ty as "working arrays"? Would that be feasible? Can we simply go with Arrays for this PR and think of this as a separate feature (as @tlycken suggests) or does it have to be a part of this PR for some technical reasons?

[pwl]

And just a comment, in implicit methods Ty can only be an (abstract) vector because we have to solve an implicit equation for the next step via Newton iterations, which require a Jacobian for the whole system to be a 2d array. In principle it might be possible to prepare a specialized array type for each call but I guess it won't be very efficient.

My opinion is that types other than vectors would be very impractical to work with inside ODE.jl for reasons already mentioned and I think we should refrain from implementing API this way. Anyway, with iterators in place one can easily produce wrappers to generate user defined types from the array returned at each step via imap for example, same goes for the wrapper around F (I am doing this in my code and it works well enough). If somebody wants to define a specialized type he can also create these wrappers without much effort.

[mauro3]

I removed ode54 and renamed it to ode45. So the exported API should be non-breaking.

@berceanu: I definitely cannot reproduce your error. I ran your topo-condensate code and it works for me (well the @test_approx_eq_eps tests fail). Let me know if you have any luck with a fresh clone of ODE.jl.

[coveralls]

Coverage decreased (-3.11%) to 93.65% when pulling 7785e41 on mauro3:m3/rk into 162c881 on JuliaLang:master.

[berceanu]

I did

rm -rf .julia/v0.3/ODE/
git clone [email protected]:mauro3/ODE.jl.git ODE
git checkout m3/rk

and using ODE still gives the same error. Am I missing something?

On 05/22/2015 10:53 AM, Mauro wrote:

I removed |ode54| and renamed it to |ode45|. So the exported API
should be non-breaking.

@berceanu https://github.com/berceanu: I definitely cannot reproduce
your error. I ran your topo-condensate code and it works for me (well
the |@test_approx_eq_eps| tests fail). Let me know if you have any
luck with a fresh clone of ODE.jl.

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

[mauro3]

I did those steps as well, just to be sure, and still no error here. I also tried with a build of the realease-0.3 branch with:

   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "help()" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.3.9-pre+27 (2015-05-20 15:14 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 8418a3a* (1 day old release-0.3)
|__/                   |  x86_64-unknown-linux-gnu

with no error.

What operating system are you on? You compiled from source with the release-0.3 branch, right? Can you try it with 0.3.8 or the current head of relase-0.3? Are the julia-tests all passing? Or from a julia-0.3 binary?

If you got a 0.4 version around, can you try that?

[berceanu]

OK so I did a Pkg.update(), which bumped Compat from 0.3.7 to 0.4.4
(among other things) and now using ODE no longer crashes with the code
in your branch.

On 05/22/2015 01:31 PM, Mauro wrote:

I did those steps as well, just to be sure, and still no error here. I
also tried with a build of the realease-0.3 branch with:

| _ _ ()_ | A fresh approach to technical computing
() | () () | Documentation: http://docs.julialang.org
_ _ | | __ _ | Type "help()" for help.
| | | | | | |/ ` | |
| | || | | | (| | | Version 0.3.9-pre+27 (2015-05-20 15:14 UTC)
/ |_'|||__'| | Commit 8418a3a* (1 day old release-0.3)
|__/ | x86_64-unknown-linux-gnu
|

with no error.

What operating system are you on? You compiled from source with the
|release-0.3| branch, right? Can you try it with 0.3.8 or the current
head of relase-0.3? Or from a julia-0.3 binary?

If you got a 0.4 version around, can you try that?

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

[mauro3]

Cool, thanks! I think that must have been the new Tuple types which were not included in the old compat. Are those assertion errors I got running your code anything to worry about?

[berceanu]

Well, it looks like with the new ode45, for my test case I get the same
results for aout[end] with accuracy of 1e-6. The test that fails
required 1e-15, so I guess it's not something to worry about.
The other difference seems to be the number of time steps taken, before
I had length(tout) = 11 and with the new solver I have 12. Again,
probably nothing to fret over.

On 05/22/2015 02:14 PM, Mauro wrote:

Cool, thanks! I think that must have been the new Tuple types which
were not included in the old compat. Are those assertion errors I got
running your code anything to worry about?

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

[kmsquire]

I think that must have been the new Tuple types which were not included in the old compat.

Perhaps add a minimum version to REQUIRES for Compat.jl?

[mauro3]

@berceanu how does the performance compare? Time and memory allocated?

@kmsquire: done. Also increased julia version to 0.3

(and changed one copy to deepcopy)

[coveralls]

Coverage decreased (-3.11%) to 93.65% when pulling 9eb745b on mauro3:m3/rk into 162c881 on JuliaLang:master.

[acroy]

Ok. I just tagged a new version for ODE.jl/master. From my side we can merge this. @pwl also gave +1 and there were no objections. Let's do it then.

[acroy]

Merged by popular demand :-)

[mauro3]

:-)

[pwl]

Great! Tanks @mauro3!

[acroy]

@mauro3 : I knew I forgot something: could you please update LICENSE.md?

[mauro3]

It's here #73. Let me know what you think.

[tomasaschan]

💯 👍 ❗

Great job landing this!