Proposal for GSoC 2016

[PerezHz]

Hi, everyone!

My name is Jorge Perez, I have been learning Julia since last year's summer, and would like to submit a proposal to work on the ODE.jl package, for this year's GSoC!

I'm currently a physics PhD student at UNAM, Mexico, under supervision of Luis Benet and David Sanders, authors of JuliaDiff/TaylorSeries.jl.

I have already talked a bit with @mauro3 on julia-users mailing list, but wanted also to contact you guys in order to introduce myself and see if there's anyone else interested in mentoring or co-mentoring this proposal.

The main idea

I have been looking at your work and it's really awesome! It would be an honor to be able to work with you guys! The main idea for my proposal is this: to implement automatic differentiation techniques to ODE.jl, in particular the Taylor method. As you know, AD techniques are essentially based on exploiting Taylor series expansions for a given initial value problem (IVP). The advantages from using these techniques include:

• Implementing high-(variable)-order, variable timestep, ODE solvers which use AD techniques to generate high-order (20th, 30th, 40th order) Taylor methods. Taylor methods offer the ability to reach machine-epsilon precision, with very competitive speeds.
• Derivative estimation (Jacobians, etc): instead of depending on finite-difference methods to calculate derivatives (which are error-prone as well as numerically unstable), AD techniques can calculate, for example, 1st and 2nd time-derivatives, or even Jacobians, automatically, to machine-epsilon precision. Looking at ODE.jl source code, AD could replace finite difference for example here, here, and here.
• Variable order and timestep: when integrating near singularities (i.e., regions where the "force" increases substantially) it's cheaper to increase the order than to reduce the timestep (see, for example section 3 of Jorba and Zou, 2004).
• Root-finding: Taylor methods allow to efficiently perform root finding (e.g., Newton-Raphson), without having to explicitly calculate derivatives, via high-order interpolation.
• Sympletic methods: Since Taylor methods are so precise, as a bonus they do preserve symplectic structure without asking them to preserve it: for energy-conserving problems, relative energy error is essentially a random walk around zero (each step of this random walk is of the order of eps(Float64)), due solely to rounding errors (can't ask for anything better, numerically!).

For the last two points I think AD could help ODE.jl a lot, as @pwl was discussing in #49 for root finding, and @tlycken, @stevengj, @mauro3 and @acroy were discussing in #7 and #10 for symplectic integrators.

A 100% Julia example

Right now I have an example on which I've been working (:100:% pure Julia!). There, I calculate times of lunar and solar eclipse occurrence, and compare those times with the ones obtained by NASA's Goddard Space Flight Center (GSFC) and Jet Propulsion Laboratory (JPL). For this example, I used a ~25th order, fixed timestep, Taylor method which integrated numerically a model of the Solar System including the Sun, the Moon and all the planets (Mercury through Neptune), and implemented a root-finding subroutine (which actually exploits TaylorSeries.jl) to calculate the time of each eclipse occurrence (ie., lunar conjunction or opposition wrt the Earth's orbit around the Sun) to an absolute precision of 1.0E-20. In a time span from September 2015, through December 2020, my code was able to reproduce each eclipse, with a difference of 25 seconds (max!) wrt NASA's results. This is ongoing work, but if you look carefully at my results, there seems to be a systematic linear drift from one eclipse to the next. Currently, I think this drift appears to be physical (it appears that polar flattening of both the Earth and Moon is really important!) rather than numerical (the relative error in the energy is less than 40*eps(Float64)).

I'm looking forward to receiving your comments, suggestions and ideas for this proposal!

Best regards,
Jorge

[tom26]

Just want to say I love the idea of using AD for numerical ODE algorithms :) Good luck sir.

[PerezHz]

Thanks, @tom26 !

[acroy]

It's great that your are interested in doing a GSoC project for ODE.jl! I really like your proposal to add AD ability (the finite difference scheme we use now always bothered me). Could you outline to what extent you would modify/improve the existing methods in ODE.jl and which new methods you would add?

[PerezHz]

Thank you for your reply @acroy ! I'm working on it, just give a couple of days and I'll post here my full GSoC proposal, but basically the idea would be to modify existing finite difference estimations (e.g., jacobians) to AD estimations.

I would also add a new ODE solver which would work using TaylorSeries.jl to implement a variable-order, variable-timestep, Taylor method (conforming to API structure), for non-stiff, explicit problems of the form dy/dt=F(y,t).

For the latter type of problems, I would also add a Newton-Raphson method for finding roots of the equation G(y)=0 (y is the state vector), which would return the value of the state vector y_, the value of the independent value t_, and the error at the event (=G(y*)).

I would of course document everything and give some cool, useful examples!

Anything you think would be important to add, please let me know

[obiajulu]

Hello,

My name is Joseph Obiajulu, and I'm currently a junior study mathematics and computer science at Princeton University who is interested in contributing to ODE.jl this summer as part of GSoC.

Reading through issues from over the last few years, I have seen that most of the development of ODE.jl has focused on implementing new solvers to and testing them. I was curious is there are any additional ode solvers that the Julia community would like to see implemented, and if there may be another need for a consorted testing effort (updating testing suite, etc).

Hopefully there is a project that I might be able to help progress!

Also, @jorgeaperezhz, your project idea sounds absolutely genius. Whether or not I join the ODE.jl team for the summer, I will be very curious to see how this project develops. Keep up the great work!

[ufechner7]

I might have mentioned this before, but having an interface to the RADAU solver in Julia would be very useful for solving very stiff DAE's fast. Assimulo wraps Radau5DAE and Radau5ODE for Python. This could be used as reference (see: http://www.jmodelica.org/assimulo ). The original Fortran code is available here: https://www.unige.ch/~hairer/software.html

[tomasaschan]

The purpose of ODE.jl is to provide a nice suite of pure-Julia implementations of ODE solvers; interfacing and calling out to third-party (non-Julia) libraries is not less valuable, but better located in library-specific projects such as https://github.com/JuliaLang/Sundials.jl.

[ufechner7]

Well, the algorithm of the RADAU solver could be implemented in Julia. Assimulo also has a Python implementation of the algorithm. It would be interesting to see, how fast the Julia version would be in comparison to the Fortran version. Furthermore the Julia version could be extended to use automatically calculated Jacobians, which would be awesome.

[tomasaschan]

Yes, for a pure-Julia implementation, this is definitely the place!

[acroy]

@obiajulu Thanks for your interest in contributing to ODE.jl. I think there is still need for a few more solvers. Apart from the already mentioned RADAU, there is also a PR (#72) for implementing RODAS and maybe @mauro3 would appreciate some help for that. Otherwise, it would be great to have a proper testing suite - especially when there will be more methods implemented.

[mauro3]

Thanks @PerezHz and @obiajulu for being interested in doing a GSoC project on ODE.jl!

@PerezHz: I like the idea of including Taylor Methods in ODE.jl and generally adding AD-tools. Concerning AD (and finite-diff) evaluation of Jacobians, I would really like to see matrix coloring methods implemented. I did some preliminary work on this (https://github.com/mauro3/MatrixColorings.jl). (Also, I suspect that even with Taylor Methods implemented there is a case for having symplectic solvers, as they might be more appropriate for some cases.) I think adding test/benchmark cases, such as your example, would make a proposal more rounded. Root-finding should be implemented in a general fashion such that all solvers supporting dense output can use it (although solvers/dense-outputters could use special root-finding methods if appropriate).

@obiajulu: yes we need more solvers! I would really like to see solvers which can be applied to space-discretised PDEs (e.g. my PR #72). For instance, IMEX methods could be interesting, but also fully implicit methods. Also, structure preserving/symplectic solvers would be good to have. I think working on implementing several solvers and related test-cases would be a great and very valuable GSoC project.

Last, I recently put some thoughts into how to take ODE.jl forward and wrote it down: https://github.com/mauro3/ODE.jl/blob/m3/future-plans/doc/ideas-for-future.org. It is a bit rough at the moment, suffers from a egocentric writing-style and has not received community feedback (except from @stevengj). So it might not reflect what the community thinks. However, I think it could be a value for GSoC proposals/ideas, thus I thought I better post it now than wait until it's too late. Please @PerezHz and @obiajulu, have a look at it. (@ the others, if you think this is valuable, then I can open an issue/PR to discuss it further)

[obiajulu]

Hello All,

I should start by thanking everyone for giving such great feedback as to what they think would make a good GSoC project over the last three days. I have researching and working on different project proposals, and wanted to update you all with a couple of GSoC ideas. I am especially interested to hear what you all think of these potential directions for the summer and see if there are any who would like to be a mentor on one of these projects. Also, I'm eager to help the Julia community, and some of these proposal may be a bit ambitious; I would greatly appreciate your help in predicting what could be done in one summer and what I must save for later contributions:

Project ideas:

1. New ODE.jl Solvers for Stiff ODEs and Improved Testing Suite
   1. Potentially Solvers to be implemented; I would work out with my mentor which ones to focus on of the ones below, or others which may be more in demand
      1. Native Implementation of LSODA or LSODA like method (Big, multistep project)
         1. LSODA automatically switches between Adam’s method and BDF method, and is a widely used package
         2. Implement Implicit Runge-Kutta methods
         3. Implement Adam-Bashforth methods
            1. Could use RK for initial values
         4. Implement BDF methods
      2. Realistically, the above project would probably be a full summer’s worth. But, here are other options for new solvers:
         1. Finish RODAS solver
            1. Work on this solver was started by, and I would work to finish PR#72
            2. Work on Some Sympletic Solvers
               1. 2nd order: Verlet method (work to finish PR WIP: symplectic solver based on velocity verlet #10)
               2. 1st order: Sympletic Euler Method
               3. 4th order: PRK (Partitioned Runge–Kutta)
2. Implementation of Numerical Methods Useful for Solving Time Dependent schrodinger equations (More of a project for JuliaQuantum project, but related to ODE.jl). Would focus on one or more of the following, depending on which the Julia and JuliaQuantum community most are interested in
   1. Implementation of symplectic splitting methods for
      1. PRK (Partitioned Runge–Kutta)
         1. http://personales.upv.es/serblaza/2002JCAM.pdf
         2. http://sanzserna.org/pdf/645_JMSS-A.%20Portillowavepacket.pdf
   2. Implementation of Lanczos method time propagation
      1. Motive: The Lanczos method “For time-dependent Hamiltonian operators the Lanczos method is the best choice, because it provides good accuracy and efficiency for a wide class of potentials” (source: http://arxiv.org/pdf/1407.7370.pdf)
      2. Reach goal: Other Krysolv space methods, like BICGSTAB(I) and/or other polynomial approximation method Chebseyv, potentially building off work done in ApproxFun.jl
   3. Implementation of Floquet Change of Basis Method
      1. For systems with periodic time-depedent schrodienger
      2. Motive: Issue Renamed tests.jl to runtests.jl #34 in QuDynamics.jl
   4. Other General Resources:
      1. http://lup.lub.lu.se/luur/download?func=downloadFile&recordOId=3363166&fileOId=3363176
      2. https://na.uni-tuebingen.de/ex/num4_ss10/lecture/chap3.pdf
      3. http://www.sam.math.ethz.ch/~mhg/pub/biksm.pdf
      4. http://qutip.org/docs/3.1.0/guide/dynamics/dynamics-floquet.html#
      5. https://en.wikipedia.org/wiki/Lanczos_algorithm
      6. https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method
3. Enhancing Testing Suite
   1. Motive: This project does seem to be in demand (see, for instance, Issue Improve the test suite #56 and https://github.com/mauro3/ODE.jl/blob/m3/future-plans/doc/ideas-for-future.org)
   2. Work would focus on improving IVPTestSuite.jl
4. Update Native Julia implementation of SDE Solvers
   1. Seems Euler–Maruyama Method already implemented in SDE.jl, but there are still some SDE methods that have not, namely Milstein, Strong Convergence methods and Weak Convergence methods referenced in http://www.caam.rice.edu/~cox/stoch/dhigham.pdf and http://math.gmu.edu/~tsauer/pre/wires.pdf

I am and will continue to work on these four proposals in parallel, though paying more attention to the first two, as they seem to be in higher demand. Also, I figure having multiple potential projects is always a good thing, especially if one finds themselves finishing one earlier than anticipated.

Also, I did accumulate a couple of comments and questions from various people’s posts and previous work:

@acroy: I think there is already a package SDE.jl. I don’t know how functional it is with the most recent release of Julia, or how well maintained it is. Do you think working on this SDE.jl would be a worthwhile GSoC?

@mschauer: I see that you were very much actively working on SDE.jl a few years ago, but seemed to be halted by technical problems. It seems, though, that this issue has been resolved with fixed-size arrays being implemented (PR #5857 in JuliaLang/julia). I was wondering what you think of me picking up where you left off and potentially have you mentor me this project of enhancing SDE.jl? Also, what are your thoughts on the SDE methods I suggested on implementing?

@mauro: What are your thoughts on some of the new solvers I proposed? Would you be interested in being a mentor for the implementation of any of them. Also, what do you think about native LSODA or LSODA like solver? Also, I had trouble finding the licensing details on LSODA, so I’m not sure if this could be a potential solver to directly implement in Julia, but I wondered if you might have more luck.

Once again, I would love to hear your feedback and be curious to see who might be interested in mentoring for one of these projects.

[PerezHz]

@acroy answering your questions about what existing methods I would modify, those would be mainly to replace finite-difference (FD) estimations by automatic-diff (AD):

hinit: ODE.jl, line 82

    d2 = norm(f1 - f0, Inf)/(tau*h0)

fdjacobian: ODE.jl, line 205

    dFdx = (F(t,x+dx)-ftx)./dx

fdjacobian: ODE.jl, line 218

        dFdx[:,j] = (F(t,x+dx)-ftx)./dx[j]

ode23s: ODE.jl, line 293

        T = h*d*(F(t + h/100, y) - F0)/(h/100)

I think that covers most FD estimations on ODE.jl. That would have to be done very closely with @mauro3, to ensure we include the MatrixColorings.jl approach along with AD techniques. About what new methods I would add to ODE.jl, those would be essentially

• A Taylor integrator for explicit ODEs
• A root finder/event finder based on AD
• A Taylor integrator for fully implicit, non-linear DAEs (e.g., this paper for equivalent work done on C++)
• A validated Taylor integrator for explicit ODEs

About the last point: a validated integrator means one which uses interval arithmetic (ValidatedNumerics.jl) to integrate ODEs. A validated integrator paves the way for computer-assisted proofs, which have been applied, e.g., to mathematically demonstrate the existence of the Lorenz attractor. Really powerful, but also very ambitious.

I hope that answers your questions about what my proposal is specifically; otherwise please let me know!

[PerezHz]

@mauro3 thank you for your reply!

I'll be happy to combine the MatrixColorings.jl approach along with AD techniques! About symplectic solvers, I definitely agree with you, having a Taylor integrator doesn't override the need for symplectic solvers; documentation-wise I do think it's important for the user to know that a Taylor integrator can be a decent alternative for a symplectic integrator, for high-accuracy problems. Thank you for your feedback and I will definitely include your suggestions about testing and benchmarking.

[PerezHz]

Thanks @obiajulu, you have some really great ideas and I hope we can both contribute to ODE.jl!