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
         4. Could use RK for initial values
         5. 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
         2. Work on this solver was started by, and I would work to finish PR#72
         3. Work on Some Sympletic Solvers
            1. 2nd order: Verlet method (work to finish PR #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 #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 #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!

[acroy]

@perezhz : Your ideas sound very promising and it would be great to see those plans realized! The reason I was asking for details about improvements vs. methods is, that pointing this out clearly in your proposal can be helpful for reviewers. As @mauro3 mentioned, it would be good to also add "test/benchmark cases" to the proposal.

[acroy]

@obiajulu : Thanks for sharing your ideas. I think @mauro3 probably has the best overview to comment on the new solvers you suggested (in terms of priority and feasibility). Regarding your other points under 2 and 4:

• "Lanczos method time propagation": Note that there are (at least) two packages in this direction ExpmV.jl and Expokit.jl.
• "Floquet method": Personally I would find it very interesting to have a package, which implements the Floquet method (or the t,t' method) for solving ODEs/PDEs with time-periodic terms. However this view might not exactly reflect the needs/opinion of the community :-)
• "SDE Solvers": I know about SDE.jl and there are a few others. However, IMO the scope of those packages is different from a package a la ODE.jl for stochastic differential equations. As you say, one would need higher order methods, maybe even adaptive methods, a unified interface and (convergence) tests. I think such a package would be very useful for many applications (including JuliaQuantum related projects). A while ago there was also a brief discussion in julia-users on this.

Generally, to increase your chances of being selected I would recommend to do a little demo project, which shows your (Julia) coding fu.

[mauro3]

@obiajulu, yes I would be interested in mentoring or co-mentoring your projects-ideas 1 and/or 3. Thanks for giving a good outline of possible projects!

Concerning 1: solver-wise I would most like to see some implicit Runge-Kutta methods (RADAU and/or SDIRK). For a BDF-like method, MEBDFI (http://www.dm.uniba.it/~testset/solvers/mebdfi.php, BSD licence) could be a good one, and improves on BDF. (In particular as we got a BDF with DASSL.jl already). Both the RADAU and MEBDFI do well in the test cases run here: http://www.dm.uniba.it/~testset/report/testset.pdf, so presumably this means that they are often good. I have no clear picture of LSODA (BTW in public domain), but could be good too but maybe a bit dated? However, as warm-up an explicit solver might be good, for instance your proposed Adam-Bashforth or maybe a SSP R-K method (https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/TS/TSSSP.html).

Concerning 2: the above-mentioned test-suite https://www.dm.uniba.it/~testset/testsetivpsolvers/ would be a good model for expanding IVPTestSuite.jl. But I think just doing this would be a bit boring. Thus I suggest to maybe mix and match 1 and 2 a bit: improve 2 such that it can be used to validate the work on 1.

I second @acroy, on that you should do a demo Julia-project.

Also, time is getting short, I'm registered as a GSoC mentor now, if you upload a proposal I will comment on it. Note however, I will not be available on Friday!

[mschauer]

Hi @obiajulu Yes, I picked up the work on SDE.jl again, https://github.com/mschauer/Bridge.jl is going to replace SDE.jl. You have it right, Bridge.jl did diverge from SDE.jl because of #5857-related issues. I can certainly assist you on contributions on that area. @acroy has a point that this is a bit of a different beast than ODE.jl

[obiajulu]

@acroy Thanks for your feedback and sorry for the late reply. I personally like the idea of implementing the floquent method as well, but I'm not sure I could find someone to mentor the project, which may reflect the fact that it is not much in demand. Alas! I may still try to work on it on the side.

@mauro3 Thanks for your feedback and your interest in potentially mentoring one of my proposed GSoC projects. I'm currently doing the final work and research (I'm looking into MEBDFI) needed to finish my proposal and hope to upload it by Wednesday or Thursday.

@mschauer That's great. I'll definitely look into Bridge.jl Thanks for the lead.

[mauro3]

@obiajulu If you upload a proposal-draft, I can give you feedback.

[mauro3]

@obiajulu and @perezhz Are you aware of the "proof of enrollment" which needs to be submitted by the deadline tomorrow too? Here https://developers.google.com/open-source/gsoc/help/proof-of-enrollment

[PerezHz]

@mauro3 Thanks, yes, I'm aware of it and have already uploaded it!

[acroy]

@perezhz, @obiajulu : Good luck with your proposals! Looking forward to have you onboard.

@obiajulu : I don't know how much the Floquet method is used outside physics, but there it is certainly an important tool. Let's see how things work out, maybe we can come back to that idea after summer.

@mschauer : Great to hear that SDE.jl is continued. Would you nevertheless be interested in starting a dedicated stochastic differential equation solver package? We could start by collecting a few methods and finding a nice interface. Then we will see ...

[obiajulu]

@mauro3 Thanks for the heaps up! I already uploaded the proof of enrollment as well.

[mschauer]

@acroy, @obiajulu You could especially have a look at https://github.com/mschauer/Bridge.jl/blob/master/example/tutorial.jl, though a bit outdated, to see "what is already there"

[obiajulu]

@mauro3 I am more or less finished with my written proposal, and I uploaded a draft to GSoC. I know that you may not have time to give any feedback, but I figured you might be interested in at least taking a glance.

@mschauer Definitely will take a look

[mauro3]

@obiajulu Thanks for submitting a proposal!

@perezhz Maybe you can close this issue now?

[PerezHz]

Sure @mauro3 ! Thanks everyone for your support and feedback! I'm closing this issue

[mschauer]

@obiajulu Congrats, success with your project.