Hi Chris, thanks for elaborating (maybe you should copy your message into some vision.md file in your repo). Your plan sounds cool, good luck with it! I would definitely use a plug and play ODE/PDE solver infrastructure.
I agree that at the fast pace you're moving at the moment it is best to just have your own package(s). However, if you code up more ODE/DAE solvers, it would be great if those are usable without all the machinery of DifferentialEquations.jl and if they eventually make into ODE.jl (or some other "standard" ODE package). A few more comments in-line below. Cheers! Mauro On Wed, 2016-06-08 at 20:40, Chris Rackauckas <rackd...@gmail.com> wrote: > Thanks for chiming in! There's a lot to say about this topic, and I think it's > great! I felt like Julia's differential equation solvers were behind where > they > should be, which was crazy to me since there is so many state-of-the-art > mathematical packages in other areas (optimization, stats, etc.) I think we are missing someone who actually does research in this area. As far as I can tell, all of us ODE.jl-ers are just "users", but not researching new methods. > There are a few reasons why I didn't go about helping ODE.jl. I think the API > is the biggest one. My muse here was JuMP. I used JuMP for a few really > difficult optimization problems, and just swapped around the solvers by > switching a keyword argument, and it handled the rest. I thought it was > beautiful and the idea was simple: put the information for a problem into a > type, and then just wrap solvers. Once I had that idea and gave it a try on > the > FEM "scratch work" I had been playing with, I just thought that this is why > Julia will be the future. So my idea was to structure problem specifying and > solving as follows: defining a problem is defining a type, from there you > dispatch to a solver which is preamble, solver algorithm, put solution to > object, and then have helpers make the solver object easy to use. The key here > is that the solver algorithm is just the middle: you can just swap out the > middle part in seconds, not just as a user, but as a coder. If you look at > sdeSolvers.jl, you can see what my final plan is going towards: define all of > the solvers as drop in code, and choose one in a conditional (though with the > loop as part of the "def" code, right now that's not the case, and there > should > be easy loop pre-ambles post-ambles for saving outputs etc). Once you have > that > structure, adding new solvers is easy! So I added all of the solvers in > odeSolvers (i.e. all of the standard solvers which I think ODE.jl has as well, > plus a few) in one morning, just adding a new "way to handle the middle". > Having the special output type also made convergence testing and plotting free > (plotting is then super nice because the solution types have recipes to Tom's > Plots.jl!), and from there things were just in motion. I think one of the key features here is that whilst JuMP is all fancy, its components are pretty down to earth. If you just want to optimize a function without using JuMP's DSL then just use, e.g., NLopt.jl. Maybe something to keep in mind when designing DifferentialEquations. > Once I had that structure, I wanted it all to be similar so that way they > could > interact well. Even though they do currently use it, there will be (S)PDE > solvers which allow for using method of lines discretizations which then just > plug into (S)ODE solvers. And what about all of those C/Fortran > solvers? (Do you know https://github.com/luchr/ODEInterface.jl and TS of https://github.com/JuliaParallel/PETSc.jl ?) > They can also just be added as part of the routine, and to switch over > to it the user just swaps a keyword argument. This freedom of solvers > has both practical, academic, and educational advantages. The > practical advantage is that with users able to change 10 characters to > switch the solvers, inevitably adding new solvers becomes "free new > features". Also, these features can be existing solvers. For example, > if you guys are doing well, and option may be to use ODE.jl solvers > (if not just to benchmark between them, our solvers will most likely > converge over time to about the same speed). But academically, this > structure let me (and anyone else) just experiment with different > solvers. Everything else is handled for you (plotting, convergence > tests, setting up the internals, saving everything), so you literally > just write the loop. That's why I am using it to develop SDE > algorithms. And lastly, education. For example, it will only take > about 10 lines to implement a Simpson's Rule solver for ODEs. This is > a quiz problem I have to my undergrads: what order is it, and what's > the region of stability? It's always unstable, so you'd never want to > use it on a real problem. I plan to just add that as an option in the > solver because it won't slow anything else down: it'll just be a bad > option. The documentation will say it's a bad method, but this lets > the curious students see what happens. > > So once I had this JuMP-inspired idea, I was scared of ODE.jl's pace > of development. I didn't want to wait to talk about API decisions and > all that: So I made a prototype. If you have special API requirements, then it would still be good to hear them. > Note that at first I wasn't going to step on ODE.jl's toes so I put in > the README that it doesn't handle ODEs, but then I needed one to > generate some silly figure, and so for fun I used to SDE solver > template to then make the ODEProblem, coded RK solvers from memory, > opened up Hairer's book and implemented Rosenbrock solvers, and tada I > had an ODE solver suite some hours later. Will the whole infrastructure for SDE solvers be usable for ODE/DAE solvers too? With infrastructure I mean: adaptive step size and order selection, dense output, root finding/events, etc. > I don't think then that DifferentialEquations.jl and ODE,jl are diametrically > opposed in any way. DifferentialEquations.jl is more this high-level interface > to how to specify differential equation problems, interface with solvers, and > then give a structured solution which is "easy to play with" (I think making > an > animation of a solution using animate(sol) is just too cool). It has a bunch > of > solvers because I was testing it all out and because that's what my PhD > research is (stochastic solvers specifically, so of course I implement > deterministic solvers all the time though!). ODE.jl is a set of ODE solvers, > and can actually be used by DifferentialEquations.jl (though it hasn't yet). > > How this should all be organized? I don't really know. I mentioned it on the > newly formed JuliaMath organization. But I don't think anyone really cared, > but > that's also because I didn't "have" anything. Now that I have all of this > implemented, the question is coming back (as it should). My take on it is that > in this form, "differential equation problems" are all pretty much the same: > they are a problem specification via some functions, some kind of mesh (either > just time, or space and time), an associated solvers, and some helper > functions > for doing standard analyses of the solution. These pieces will then play with > each other: solving an SDE can be solving a PDE if you want to do it via the > Fokker-Plank / Forward Kolmogorov / Feynman-Kac way, solving a PDE might just > discretize to an ODE problem. or you can solve some PDEs probabilistically by > solving a Monte Carlo simulation of SDEs. That's why I kept them together. But > to be honest,I don't care about organization as much as long as it's easy to > just implement every possible way of solving things, > > There's a lot left to do. I've only had this package for less than a month, > and > was taking 6 courses + 2 seminars per quarter this year, so I didn't get as > much done as I hoped. However, I am turning in that last final in a few days, > so I will be freed up to start working more on this. Also, I have a lot of > private codes which are being submitted soon and once those are published we > will really have some state-of-the-art functionality. Next year I have courses > and plan to take this head on. So I am just going to keep doing as much as > possible, and we as a community can discuss what we want to do with it.
