[WIP] Port of DOP853

[helgee]

This is an almost complete port of DOP853 (the dense output functionality is still missing) with quite good performance (https://groups.google.com/forum/#!searchin/julia-users/dop853/julia-users/41B3qULjF-0/WqGa-Lc04AkJ), approximately 2-3x slower than Fortran (without compiler optimizations) for longer runs. Since the 5-th order version DOPRI5 looks rather similar implementing it as well should not be a huge amount of work

See also: #25

To Do:

• API adjustment
• Dense output
• DOPRI5
• Write tests
• Documentation
• What else?
• Is it a good idea to separate different solver families in different source files?

[stevengj]

In general, I would prefer something a little cleaner than a line-for-line port of Fortran 77 code, although that is a good baseline starting point for comparison to later efforts.

[ivarne]

• Documentation

We do not have a style of documentation yet, but it would be very nice if you added a file (eg. /doc/DOP853.md) where you explain the method and what kind of systems it performs well on. Links to the original fortran source, and other citations will also be useful.

We should also have a /doc/index.md file with links to the other files, but that is a separate task.

[acroy]

Regarding the tests, you can have a look at test/tests.jl. It should be relatively easy to adapt that code for dop853. If you add the dop853-tests to that file, we will also get feedback via Travis when you update the PR.

[stevengj]

There is an interesting efficiency question here: the problem with k=f(x,y,params) is that it forces Julia to allocate k each time you call f, rather than returning it.

I think that probably the best thing would be to have a low-level interface in which you pass a subtype of abstract ODEProblem with methods like F!(p::ODEProblem, y, x, t) that act in-place on x. The high-level interface could use a simple wrapper class like:

type ODEProblemFunction <: ODEProblem
     f::Function
end
F!(p::ODEProblemFunction, y, x, t) = copy!(y, p.f(x,t))
odeXX(f::Function, ....) = odeXX(ODEProblemFunction(f), ...)

The ODEProblem could also have methods to describe more complicated DAE and implicit problems like those discussed in #30, methods to compute the Jacobian (or both the Jacobian and f at the same time), etc.

That way, if the user needs greater efficiency, they could define their own ODEProblem subtype and pass it directly.

[ivarne]

@stevengj Good to see that you are considering the techniques I suggested in #31. I'm not sure your approach will give as good type inference as the one I suggested, but it looks prettier so I'd love to be proven wrong.

[stevengj]

@ivarne, a generic ODEProblemFunction will have to do runtime type inference to figure out the return type of f, but if the user defines their own ODEProblemFoo type with their own F! method I don't see why the inference couldn't be done at compile time.

It is just nice to have an interface where the user can pass a function for simplicity; I view forcing the user to define their own type as a low-level interface that users can ignore unless it is needed.

[ivarne]

@stevengj Sorry, I read your comment again and you are completely right. We had very similar ideas, but I focused too much on the implementation detail for allowing anonymous functions in my comment.

[tomasaschan]

There seems to be a lot of code here that just defines a million variables for various parameters. Is there some logical way to group them to make the code more readable?

I'm not familiar enough with these methods to be able to look at it at a glance and determine, but the current effort for oderkf is to use one single "core" method that does the actual work, and then a few wrapper methods that pass the relevant coefficients as parameters - for example, as of #16, that lets the same code be both 2nd and 4th order adaptive RK solver. Adding new solvers to the RK-family should be as easy as adding a new table of coefficients and a wrapper method that passes the correct table. Is the same approach possible with these solvers? I.e., would it be possible to (sometime) refactor this so that we could define

function odedop(APIvars, solverparams)
   # main implementation of DOP solvers like DOPRI5 and DOP853
end

dopri5_p = [] # parameters for DOPRI5
dop853_p = [] # parameters for DOP853

dop853(APIvars) = odedop(APIvars, dop853_p)
dopri5(APIvars) = odedop(APIvars, dopri5_p)

If that is possible, it would make for much easier addition of new solvers to the family. It will also make it easier to "hide" all the ugly parameter definition stuff from the implementation code and increase legibility, since we can just define all the parameter objects last.

[helgee]

First of all many thanks for the valuable feedback everybody!

@tlycken Sounds like a good idea to me.

@stevengj & @acroy Did I understand correctly that the current approach to passing parameters to the callback requires that the callback function is defined in the same scope as the parameters? I think this might be impractical if a library contains "reusable" callbacks and only the parameters are read from input files. Are there any reasons why we should not use a params keyword? When defining your own ODEProblem this is of course not necessary since you can just add the parameters to the data structure.

[ivarne]

Reasons to not use a params keyword:

1. Dictionaries are slow, and keyword arguments give poor type inference. While type assertions can partially solve the inference problem, they makes the code look ugly.
2. You will have a runtime check to see if you should pass the params to the callback. Alternatively you can use macros to generate code specifically for both cases.

@helgee I think I understand why you ask, but I don't see the issue with scope and ODEProblem. Both ODEProblem and the F! function is defined in the ODE module, and can be independently extended from any module/context. You will just have to import them, then extend them.

[helgee]

@ivarne You misunderstood me there. When one uses the low-level API with ODEProblem there is no issue, this only concerns the high-level wrapper. So you would suggest to always use the low-level API in a library context?

[ivarne]

I would really prefer to have only one interface, but unfortunately the most intuitive interface is slow, and we need better options for performance critical cases. The question is how advanced cases we need to support for the simple interface, and when we want to say, "use this faster interface instead for complex problems."

[stevengj]

@helgee, we don't need a params parameter because Julia has anonymous functions and lexical scoping. i.e. you can always pass (t, y) -> f(t, y, params). Or you can declare an new ODEProblem type to carry around extra state if you want to use the (proposed) low-level interface.

My suggestion is that all of our library implementations use the low-level (object-passing) interface, and that we provide (one-line) wrappers to support the high-level (function-passing) interface(s).

[berceanu]

What is the current status of DOP853? Is anyone still actively developing it? What remains to be done (dense output, etc)?

[tomasaschan]

@berceanu I don't know if @helgee is still doing any work on this, and I haven't taken a close look at the code to see what has already been done and what has not. I suspect you could help out e.g. by creating pull requests against his fork - when he merges them into his branch, the commits will automatically be added to this pull request.

Take a look at doc/api.md for the latest discussion about what the API should look like.

[berceanu]

@tlycken I'm quite new to github. Could you please explain the "pull requests against his fork" suggestion in more detail? I already have a new version of his code, so how can I contribute?

[ivarne]

@berceanu It is hard to know what level of detail you want the explanation on. Please ask (or search) if you need more specific directions.

• Make a fork of ODE.jl on github. (click fork)
• Add the aproperiate remotes to the ODE.jl git repository in .julia

cd ~/.julia/ODE.jl
git remote add helgee [email protected]:helgee/ODE.jl.git
git remote add fork [email protected]:berceanu/ODE.jl.git
git checkout -b DOP853 --track helgee/master

• Do your changes
• Commit your work

git commit

• Push your commits to your fork

git push fork DOP853

• Create a pull request on @helgee's repository. And post a link to it here. (I think this should work)

[tomasaschan]

@berceanu I'll make an attempt to explain what I mean =)

Basically, the most common way to organize work in open-source github projects is to have a main repo (the origin), which the core developers work on. If someone outside the team (say, for example, you :) ) wants to contribute to the project, they create a fork from origin - let's refer to it as fork1->origin to be clear with where it's forked from. After pushing a few commits to fork1->origin, you can create a pull request in which (s)he basically notifies the maintainers that "hey, I did some cool stuff - why don't you pull that into origin?"

To facilitate code review and collaboration, Github was designed so that a pull request isn't static - if you push more commits to the branch in fork1->origin from which you created the pull request, they're automatically appended to the pull request so we can all discuss the new changes.

This opens a nice possibility to "chain" pull requests - if you see someone working on a PR, and you want to contribute something to their contribution, you can fork their fork and request that your changes are merged into there, instead of into the original source. In other words, it is entirely possible to create fork2->fork1 (which is really fork2->fork1->origin) and start a code review process there. When the owner of fork1->origin decides that the contributions in the PR from fork2 are good, (s)he can merge them into fork1 - and they are then automatically appended to the PR toward origin, since they add commits to the branch.

Finally, it's noteworthy that you don't really have to create fork2->fork1->origin in order to make a pull request into fork1 - you could just as well create a PR from fork2->origin to fork1->origin. The only thing a PR is, really, is a message saying "hey, I want to merge this branch of mine into that branch of yours". Doing it via forks just makes it extra easy, since the user interface is built around this workflow.

I've used this pattern a couple of times when I wanted to help finish the work on a PR that someone else started, without having to first wait for that PR to get merged.

Depending on whether @helgee still wants to invest anything in this PR, this is either a good idea or not - if Helge is still in, this approach helps the two of you to finish this functionality without having to wait for action from us ODE.jl people. However, if he's not, you'll get stuck waiting even longer for him to get back to you, than you'd have to wait for us. In that case, it might be better to merge his branch into your fork, do some more work, and create a new pull request against ODE.jl.

Is this making any sense to you? =)

[helgee]

Please excuse the long radio silence but my day job was keeping me busy.

I am still interested in getting this merged and got started with the dense output functionality some weeks ago. Porting the numerical bits is actually straightforward. The challenge lies in defining a sensible and suffciently julian API and I have not made much progress there.

[acroy]

Although I think dense output would/will be a really great feature, I just wanted to point out that DOP853 without dense output still requires some work to match our API draft. IMHO the first step should really be to hook DOP853 into our "test suite", either by writing some wrapper functions or by adapting the tests.

But of course I don't want to discourage anyone from contributing to the dense output!

[acroy]

Are you using the same tolerances? (oderkf has reltol=1e-5 and abstol=1e-8)

[helgee]

Yes, the tolerances are the same.

[coveralls]

Coverage decreased (-4.64%) when pulling 134ced5 on helgee:master into 913a3a5 on JuliaLang:master.

[helgee]

The Julia version begins to deviate from the Fortran version after 4 integration steps. The time steps are identical though. Could this be due to differences in the floating point model?
https://gist.github.com/helgee/f2a4d40298ae4cb6f089

@acroy The Fortran version is three orders of magnitude more precise for problem A3. I guess there is something wrong with the Julia code then.

[acroy]

AFAICT, dopri5.f uses EXPO1=0.2D0-BETA*0.75D0 here. Not sure it matters ...

[helgee]

I tracked down where the solutions start to deviate in the A3 problem and this either an upstream Julia problem or I am doing something awfully wrong.

At that point Julia and Fortran do the same calculations with numerically indistinguishable numbers and get slightly different results. The error then probably grows due to the numerical method. I have also tried the specific calculation in C and Python and they agree with the Fortran result.

Have a look at the code here: https://gist.github.com/helgee/f179b9610aa104e38966

Julia output:

1.2620843374762893
1.2620843374762893

C, Fortran, Python output:

1.2620843374762898
1.2620843374762898

What is going on?

[ivarne]

The difference is that C, Fortran and Python uses system libm implementation for cos, but Julia uses openlibm.

You can compile julia with USE_SYSTEM_LIBM=1 in your Make.user file, or ccall the system library directly, to get a fair comparison between Julia and C, not openlibm vs system libm:

f(t,y) = y*ccall(("cos","libc"),Float64,(Float64,),t)

The reason Julia uses openlibm is that people were reporting bugs in their system libm implementations and blamed Julia for being slow or wrong.

[helgee]

@ivarne Correct, I'll keep on searching then.

[acroy]

That might explain some deviations between Julia and C etc, but it doesn't explain why the deviation from the known solution in Julia is so large. Does dopri5 also give a large deviation?

[helgee]

@acroy No. See the updated results: https://gist.github.com/helgee/0eef2ca36afb4339710d

[coveralls]

Coverage decreased (-4.66%) when pulling 00f4a0f on helgee:master into 913a3a5 on JuliaLang:master.

[helgee]

Well, I was actually writing the wrong variable into the output array for points=:all.

Now the Fortran and Julia versions agree perfectly for A3 (with OSX's libm) and the maximum deviation from the known solution is at ~1e-6.
There are still differences for the gravity example and I will keep on digging.

[helgee]

Now both integrators deliver the exact same results as the Fortran versions. Embarrasingly there was also a bug in my Fortran code, which I used for checking the results. A lot of searching in the wrong place :-P

[coveralls]

Coverage decreased (-4.66%) when pulling c0dd3b5 on helgee:master into 913a3a5 on JuliaLang:master.

[acroy]

In any case it is great that you found the error!

On Wednesday, July 2, 2014, Helge Eichhorn [email protected] wrote:

Now both integrators deliver the exact same results as the Fortran
versions. Embarrasingly there was also a bug in my Fortran code, which I
used for checking the results. A lot of searching in the wrong place :-P

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

[acroy]

I think it is better to use

push!(tout, x)
push!(yout, y1)

in order to be type stable (otherwise the return type depends on the value of points). Note that ATM :last is not in the current API draft.

[helgee]

I apologize for letting this PR go stale. Should I continue with this or would it make more sense to just add the 8th-order Dormand-Prince coefficients to the new Runge-Kutta solvers?

In that case I would still need facilities for event detection (e.g. a user function that gets called after every integration step and dense output) for ODE.jl to be useful to me. But until that is available I will use my wrapper for the Fortran codes (https://github.com/helgee/Dopri.jl).

[tomasaschan]

Adding the 8th-order coefficients would definitely be a nice contribution regardless of other things.

I personally hope that event detection can be handled with the approach outlined in #49, with support across all solvers in ODE.jl, but that's quite a large undertaking. Also, I don't have the opportunity to use ODE.jl regularly anymore - and thus my ability to invest time in it is also limited - so my opinion on the matter carries very little weight nowadays :)

[acroy]

I think we should stick to the current implementation for anything related to (explicit) RK methods. If you could add the 8th order coefficients this would be great.

If you have ideas or needs for event functions, please post them in #11.

[helgee]

Alright. I will close this then.

[mauro3]

It is easy to add the method to the existing RK solvers, you just need to provide a tableau like:
https://github.com/JuliaLang/ODE.jl/blob/01620ad46321b2088f4b90be1099b903c43dd788/src/runge_kutta.jl#L135

At the moment dense output is only 3rd order. To change that would require some code changes. I hope to implement a general scheme at some point but cannot promise a time frame. Same goes for event location.