GSOC 2014

[fsaporito]

Hi,

I am Francesco Saporito, and I was intested in applying for the ODE Solver project. I have been poking around the code and reading the Julia manual, and I am getting toghether some ideas for the project.
I wanted to know what the community had in mind regarding the ODE Solver and what the GSOC student it's expected to do to the project. In the project idea it was hinted that a major rewrite of the code was involved, but I wanted to know what you believe that has to be done for the solver.

Thank you for the suggestions.

Saporito Francesco

[acroy]

Glad to hear that you are interested in ODE.jl! As you have probably seen we are currently trying to revise the package to some extent and would be most happy to get some help. I don't know if there is already a definite GSOC ODE project, but I am sure we can come up with something that suits your interests.

I hope that the others will join the discussion with suggestions, but some potential topics that come to my mind would be

• Implementation of an adaptive solver for stiff ODEs. Such methods typically need the Jacobian matrix of the right-hand side and it might be interesting to see if/how one could make use of DualNumbers.jl, which provides a different approach to numerical differentiation than finite differences with some epsilon.
• "ODE solver facility". In WIP: Api changes #7 @tlycken suggested to introduce some kind of "solver object", which allows you to specify the problem independently of the method.
• Test & performance suite for ODE solvers. @tlycken started some work in WIP: Test framework #3, which (I believe) is based on DETEST. In general it would be nice to have a proper test suite for our solvers, which allows it to check correctness, but also gives some performance benchmarks.

But you mentioned that you already have some ideas?

[fsaporito]

Thank you for the reply.
The project is among the ones on the julia's website list (http://julialang.org/gsoc/2014/).
The description says:

"Project: Native Julia solvers for ordinary differential equations

Julia should have a full set of ODE solvers, as they are vital for numeric programming. This project should build on the work of ODE.jl. The best way forward will probably be to implement the solvers in Julia. See #75.

Expected Results: A set of production-quality ODE solvers."

This is what I have thought so far:

I would start by implementing some basic ODE solvers, like implicit and explicit Euler, Crank-Nicolson and Heun. This is actually more useful for myself than to the library, because it's lets me enter in the julia's development frame of mind and can be used to assert that I am suitable for the project. (This will come before the actual project and before the GSOC, let's call it a warm up :) ).

The first thing to do in my opinion is to improve the test suite and the documentation, to allow a much smoother development of the library.

The library should then expand by implementing some other methods, with maybe some redundance in problem solving (different methods can solve the same thing) but with the possibility for the user to chose from a group of solvers.
Maybe we can even create a comparison between the methods for the specific problem posed by the user, so that it's possible to see which method is better for the task (where better can be intended as speed, precision, etc).

For the new methods to be implemented, I don't know much on the subject yet, but I think that maybe finite volume and finite element can be added.

What do you think?

[tomasaschan]

Fun! ODE.jl could really benefit from growing its number of active developers, so it's nice to see that there is interest! =)

As you'll see if you read through the open issues, much of the current discussion is around API design and basic functionality. Many good ideas have been raised, but at least some of them might require a fair amount of work to actually be implemented - however, I don't know if maybe work on many of those things will already be started by the time GSoC starts. A good way for you to familiarize yourself with Julia and with ODE.jl would be to get involved in the discussion ;)

@acroy's suggestions for topics are great - in particular, the library is really in need for some rigorous testing, both for accuracy and performance, so if I got to pick, that would be my first priority for a GSoC project. However, there are some other things we could definitely use as well:

• Even without the dual numbers approach for differentiation, stiff solvers are really a requirement for this package to be "feature-complete". The usual way of handling them is to have the user provide an analytical jacobian of the system as an extra function argument, i.e. instead of ode(f, tspan, y0), one has to call odes(f, tspan, y0, J).
• Some of the current implementations are taken from places where it's a little unclear what licensing requirements are. Therefore, a pure white-paper implementation of e.g. a general RK method (similar to oderkf in its current form) would be a welcome addition. This is perhaps a too small a project for GSoC (I don't know) but it for example be combined with stiff solvers above.
• Documentation is also needed, but since the API is changing a lot right now, maybe we should make all the decisions before we document them ;)
• All Runge-Kutta methods are finite-difference based - it might be nice to add some finite volume or finite element methods as well.
• We currently lack BVP solvers completely.

Edit: now I saw that you've posted another comment. Yes, your ideas sound good - as I said above, a good test suite would really add a lot of value to the package. If you'd like to work on that, that'd be really great! As for warm-up, why don't you start going through the issues and submit pull requests already now? ;)

[jiahao]

If you can show us reasonable implementations of explicit and implicit Euler in your GSoC application, that would go a very long way toward convincing us that we should have you on board.

cc: @StefanKarpinski

[fsaporito]

@tlycken I thought of implementing something related to the project by starting from the basic level, to allow a sort of educational pathway in the project course and to start getting the hang of things. Starting with actual issues, would mean a more step learning curve with the julia language and a more difficult project start in my opinion.

[acroy]

@sapo93 Then @jiahao's suggestion sounds like the perfect starting point!

[fsaporito]

Give me a couple of days to sort out the material and I'll start working with Euler :)

[jiahao]

@sapo93 we look forward to your application.

[stevengj]

I want to second @tlycken's call for testing, and let me expand on that by saying that we really need rigorous testing and comparison of our solvers to competing solvers (in Matlab, Mathematica, Sundials, etcetera).

Without rigorous comparisons (on a variety of test cases exhibiting different characteristics, including both "easy" and "hard" cases), it is difficult to know whether our ODE solvers are any good.

[mauro3]

I just commented over in issue #20 that PETSc has nice ODE/DAE solver setup (EDIT: now the separate issue #30). Both from an API point of view and from what solvers it contains. Also it has a liberal licence, so code could probably be directly based off that.

[fsaporito]

I looked over chapter six of the manual and the General linear methods explanation posted in #30 . It seams to me a worthwhile idea for the GSOC project.
If the community feels that this is the right track for the library, I will apply for implementing this and for adding a test suite (which I usually do as I code, only this time it wouldn't be a use and trash test suite)
.
P.s. In the admission draft, its suggested to list some people that the student feels should be mentoring him for the project. Is there already someone that took interest in mentoring this project? in my post in the mailing list (https://groups.google.com/forum/#!topic/julia-dev/7o4gMdB8PIc) alessandro adrioni (https://github.com/andrioni) wrote that he was able to serve as a mentor for the project, but I never saw him write on the discussions regarding the library. Someone that steps forward?

[jiahao]

I can mentor, but you'd hit the jackpot with @stevengj if he were available.

[acroy]

For GSOC I would concentrate on one specific solver and implement that in a correct and performant way (with support for the features we are discussing in doc/api.md. Personally, I would like to see an implicit solver (like RADAU or RODAS), but there might be other suggestions. That should be enough work for the summer and if there is still time we will certainly find some extra things to do.

Regarding the test suite, I think one should not underestimate that task. We are not talking about the kind of really simple tests we have now. As @stevengj pointed out, we need rigorous testing and comparison of our solvers for a variety of test cases (involving different ODEs, parameters, accuracies). The tricky part is probably to obtain a sensible and useful analysis of the test results. I would consider this a separate project in it's own right.

[edit ivarne: fixed link to api.md]

[mauro3]

The referenced RADAU & RODAS work for banded Jacobians only. It would be good to have an implementation which works with a general sparse Jacobian.

[Nicksol]

Hi everyone,

I am Niluka Piyasinghe. I am interested about ODE solver GSOC project.

I am familiar with ODE's over 1 year since its included in academic curriculum and know some solving method like D operator, undetermined coefficient, Laplace transformation and matrix method.

[edit ivarne: formatting]

[Nicksol]

Hi everyone ,

I am Niluka piyasinghe.I am interested about ODE solver GSOC project.

I know many ODE solving methods like D operator, Laplace transformation, undermined coefficient and matrix method. These days I am trying to choose suitable solving method for proposal. Any suggestions?

[edit ivarne: formatting]

[mauro3]

@Nicksol: I think the methods you mention are used in analytic solutions of ODEs and not usually in numerical solutions. As mentioned earlier in this thread, interest would be in having a solver which works for stiff ODEs, i.e. an implicit solver, with adaptive time step. For example BDFs, implicit Runge-Kutta (RADAU is of that type), and general linear methods (DIMSIMs, IRKS). Have a look at Wikipedia for a quick overview.

Edit: and such a solver should work for DAEs as well. But I think most implicit ones can do that.

[Nicksol]

Here are my ideas as GSoC project

As the discussions in JuliaLang/julia#75 ODE23 solver seems work correctly . But there are some issues with ODE45 solver.Also there are some concerns about the licensing.

So as suggest in #25 , i will work on implementing pure Julia solver to ODE45 based on http://www.unige.ch/~hairer/prog/nonstiff/dopri5.f (There a C and Matlab versions of these in http://www.unige.ch/~hairer/software.html website.Since i am coding with these languages over 2 years.It won't be hard understand it.)

I will make it according to API changes describe in https://github.com/JuliaLang/ODE.jl/blob/master/doc/api.md

I am looking forward for suggestions.

Thanks.

[Nicksol]

Or I could focus on @acroy suggestions and implement RADAU solver with features in https://github.com/JuliaLang/ODE.jl/blob/master/doc/api.md .Since these API changes are tentative, I can come with changes according to it in the implementation as the project goes on.

Before submit any proposal to GSOC ,I would like to know ideas of the comminuty and looking forward for a mentor since students can suggest some people as mentor to the project when submitting the proposal.

[tomasaschan]

It's fun that there is so much interest in this package for GSoC!

I would like to re-iterate the suggestion to work on a good test suite for the package. There has been some discussion in #16 scratching the surface of how to work with performance tests, and there might be some salvageable ideas in #3 even if the PR is so old that most of the code is probably unusable. There also a few comments on this subject higher up in this thread, and the DETEST test suite seems to be a more or less mandatory resource in this context.

Specifying the exact requirements on the test suite must of course be part of the project itself, but I can think of a few things that should be considered:

• The nature of the ODE's and solvers in question: a high-order solver should be expected to provide better accuracy than a low-order solver, and only stiff solvers should be expected to be able to solve stiff problems. Calling explicit solvers on implicit methods shouldn't be possible, etc. The test suite needs to take this into account, and pair problems and solvers in an adequate way.

• There are several reasons to test this package:

  • Correctness of implementation: does evaluating this solver on that equation with those parameters really give the correct result? These are tests that exist as safe-guards against introducing bugs when updating the current implementations, e.g. for adding new features or API choices.
  • Quality of algorithms: is this solver a good type of solver for that type of problem? These tests exist so that if someone finds and implements a new algorithm for a problem type for which we already have a solver, we can compare the two implementations.
  • Performance tests: does the current implementation of this algorithm perform well (enough)? These tests exist primarily so that we can spot regressions between versions, and do something about them. (Note that regression can happen without us changing a thing - as was evident when the cause of Tests failing on nightly #29 was identified.)

  All three types of tests provide value to the test suite, but not all of them need to run e.g. every time the package is installed. Therefore, part of the project would also be to identify which tests should be run in which context (on installation, on Travis, for http://status.julialang.org, possibly for http://speed.julialang.org etc).

It's difficult to know exactly how much time it would take to implement all this - and it will probably depend a lot on how much background the student has both in numerics and in testing. It is of course very possible to do only part of this within the frames of a GSoC project.

[ufechner7]

Hello,
I want to suggest to have a look at the Assimulo solver suite (Python) for inspiration. They have a RADAU5DAE solver, implemented in Python. Documentation: http://www.jmodelica.org/assimulo_home/DAE_Radau5DAE.html
Source code: http://www.jmodelica.org/assimulo_home/_modules/assimulo/solvers/radau5.html

It can calculate the Jacobians numerically, but can also work with user provided Jacobians.
I like the option to set a threshold on how often the Jacobians are recalculated. It works very fast for our simulation problem, a stiff kite-power system model.

The also have a nice implicit Euler solver:
https://trac.jmodelica.org/assimulo/browser/trunk/src/solvers/euler.pyx

Regards: Uwe Fechner, TU Delft

[acroy]

I suggested to implement a stiff solver, because (a) we don't have an adaptive solver for stiff problems at the moment, (b) it would allow us to explore some additional aspects of the API (like handling the Jacobian) and (c) gives some ideas about stiff test-problems. Radau IIA (which is implemented in RADAU5 or RADAU5DAE mentioned above (which is GPL)) seems to be a reasonable (and manageable) starting point.

(I haven't looked a the PETSc code, but I am sort of afraid that one has to understand the PETSc machinery before one can take out one part of it, but maybe this isn't really a problem.)

[Nicksol]

I agree with @acroy.Since Julia isn't equiped with a solver for stiff problems, it's a good idea to come with solver.I am going to submit a GSOC proposal based on RADAU

Thing is i am going to need some time get more familiar with these numerical ODE solving theorems. so i will learn more in up coming weeks.(Not this week,I have mid semester exams.)

I will implement Euler method in my GSOC application as @jiahao suggested just to show that i can complete the task.

[dcarrera]

I see a lot of people linking to the codes at http://www.unige.ch/~hairer/software.html

I cannot find a license for these codes. Are they under some permissive license?

[ufechner7]

Why not asking the author directly:

http://www.unige.ch/~hairer/welcome.html

Am 06.11.2014 um 21:07 schrieb Daniel Carrera:

I see a lot of people linking to the codes at http://www.unige.ch/~hairer/software.html
http://www.unige.ch/%7Ehairer/software.html

I cannot find a license for these codes. Are they under some permissive license?

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

[ufechner7]

Actually, the license is included on his web-page:

Copyright (c) 2004, Ernst Hairer

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

• Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer.
• Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions and the following disclaimer in the
  documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS
IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

[dcarrera]

Excellent. Thanks.

[espintos]

I was wondering what's the state of this project so far, is it still moving forward?
Is there any interest for someone to mentor in the context of a Julia Summer of Code project?

I would really like to apply, I'm of course no expert in the field but I have experience using propietary software (gProms) to solve large DAE systems (with stiff ODEs) for polymerization reactions.

Although not related to solving ODEs in Julia, my experience with the language is writing a kinetic MC simulation using the Gillespie Algorithm for the above DAE system and comparing the results between the two. This work in Julia, is drafted in this notebook, I know it's a little messy and confusing, but the simulation is working and giving good results. I'm working right now at uploading a better version of the notebook.

As a minor plus, I'm a TA for Numerical Methods in chemical engineering class.

[acroy]

@espintos : First of all thanks for your interest in working on ODEs in Julia! We are always happy about new contributors.

Regarding this particular issue, I think it is fair to say that it is a bit outdated since there have been quite some developments in the last year or so, for instance,

• @mauro3's PR RFC: new adaptive (and fixed step) Runge Kutta solver #68 will (soon) replace our current set of (explicit) RK methods by a more efficient implementation, which is also free of license problems.
• Moreover, he has set up a test suite for initial value problems.
• With ode23s we have an adaptive solver for stiff problems.

This doesn't mean that everything is done of course, but we need to reconsider the project suggestions put forward in this issue. Some of the problems we have to tackle next are performance for custom scalar-like types after #68 is merged, mutating rhs functions, support of implicit ODEs, ... but I am not sure if those are good JSoC projects.

(I know that there is a JSoC project suggestion for "adding RK45 to Base". Personally, I don't think that this is necessary, but others might think differently)

[ChrisRackauckas]

Closing this because it's not relevant anymore. We can have a thread discuss new GSoC ideas, but this has served its purpose.