Plan for modularization

[musm]

There's a lot of great functionality here...but there is almost too much functionality in one package 😄

One thing I think would make sense as a first step is to separate out PDE solvers from ODE solvers. This could be done alongside deprecating ODE.jl. IMO there should be just "one" ODE package that is "the" package for ODEs in Julia. If there are two ODE packages the tradeoffs between the two should be made clear in the README. I.e. why should I use this ODE package vs the other ODE package.

A lot people will likely want to use the ODE solvers without using the PDE solvers. From experience most of the time I have to write a specialized PDE solver for my use case, which will take advantage of an ODE solver.

Things to think about:

What exactly to separate out?
Should the SDE solvers be in their own package too (currently I think this would also be a good idea, but there might be an argument made against that ..., I suppose)

The ODE package should contain: (please fill in/update)

• ODE solvers (stiff/non-stiff/geometric integrators etc.)
• exponential integrators.

[ChrisRackauckas]

A few things to unravel here.

One thing is, it's not my position to be able to call ODE.jl officially deprecated. That would have to be something the community decides on. There is work in ODE.jl in the form of SciML/ODE.jl#49 to essentially re-do ODE.jl as iterators. who knows when that will actually "finish" but it is better than the current ODE.jl (which isn't saying much because the current ODE.jl isn't always type-stable...). So while it's not for me to to unilaterally decide what is best for the community, you can see that all of the benchmarks I have are against this new iterator update, and they show that the DifferentialEquations.jl ODE solvers are both faster and more accurate than any of those of ODE.jl. I can objectively say that there DifferentialEquations.jl has more algorithms, which are both faster and more accurate, and have more features (higher order dense output, plot recipes, multithreading, etc), and leave it to the community to make the decision of what's "official", whether there should be one, and what the criteria for choosing that one is. In the meantime, I'll just keep implementing/improving more algorithms.

That said, I agree that this package should be split, the real question is how. I like the experience that DifferentialEquations.jl offers: you do using DifferentialEquations and anything in that documentation should work. One way to continue to offer that is to use conditional dependencies, i.e. have DifferentialEquations install/import the companion packages on demand. However, that would have to wait because conditional dependencies are very fragile in Julia right now. So if I was to break it apart now, what I could do is make DifferentialEquations.jl be the common package of abstractions (like a DiffEqBase), and then name things like DiffEqODE, DiffEqSDE, etc.. Then to use the ODE solvers, you'd just need to do using DiffEqODE, and the SDE solvers would require using DiffEqSDE.

Even this leads to issues without conditional dependencies, since one of the ways I am adding for solving an SDE is to solve the related Kolmogorov PDE, and the PDE solvers really need the ODE solvers, and then allowing you to build SDE problems by adding noise to an ODE problems means that they all require each other and so we're at the same point where using one of them would use all of them (this is ONLY if they aren't conditional dependencies. Otherwise you could have them require each other just when needed...). This is the current dilemma which is why I haven't made a move on splitting it up yet (the actual splitting would be pretty easy: take no more than a few hours).

To really cut down on the size of the package, I should probably split the documentation, examples, and benchmarks out. I have it set so that way Github doesn't calculate those in the line total, but if I didn't have them ignored, you'd see that the majority of DifferentialEquations.jl is actually Jupyter notebooks (when counting lines / size since those things are so bloated). Would it be okay to have the user install another package for the examples? I guess I could just add that to the lines where I tell people how to open the example notebooks.

[ChrisRackauckas]

Another thing that heavily contributes to the size of the package is the ODE tableaus. I have a crap ton of them. I think that most people should stick to using the optimized ODE algorithms and the tableaus themselves are more of a research tool (it's like 100 tableaus or so, and they plot stability diagrams and you can benchmark them all against each other). I use some of them in the tests, though they are only used as options along with the alg=:ExplicitRK algorithm. I could spawn that out to a separate package as well. That would be easy as a conditional dependency since you pass in the tableau objects to the ODE solver, so it would work without DiffEqODE (or whatever happens) ever having to call using ODETableaus.

[musm]

To really cut down on the size of the package, I should probably split the documentation, examples, and benchmarks out. I have it set so that way Github doesn't calculate those in the line total, but if I didn't have them ignored, you'd see that the majority of DifferentialEquations.jl is actually Jupyter notebooks (when counting lines / size since those things are so bloated). Would it be okay to have the user install another package for the examples? I guess I could just add that to the lines where I tell people how to open the example notebooks.

How large are all the files? If it's a few MB I'd say it's not worth separating into a new package, which will be less likely to be installed

That said, I agree that this package should be split, the real question is how. I like the experience that DifferentialEquations.jl offers: you do using DifferentialEquations and anything in that documentation should work. One way to continue to offer that is to use conditional dependencies, i.e. have DifferentialEquations install/import the companion packages on demand. However, that would have to wait because conditional dependencies are very fragile in Julia right now. So if I was to break it apart now, what I could do is make DifferentialEquations.jl be the common package of abstractions (like a DiffEqBase), and then name things like DiffEqODE, DiffEqSDE, etc.. Then to use the ODE solvers, you'd just need to do using DiffEqODE, and the SDE solvers would require using DiffEqSDE.

Even this leads to issues without conditional dependencies, since one of the ways I am adding for solving an SDE is to solve the related Kolmogorov PDE, and the PDE solvers really need the ODE solvers, and then allowing you to build SDE problems by adding noise to an ODE problems means that they all require each other and so we're at the same point where using one of them would use all of them (this is ONLY if they aren't conditional dependencies. Otherwise you could have them require each other just when needed...). This is the current dilemma which is why I haven't made a move on splitting it up yet (the actual splitting would be pretty easy: take no more than a few hours).

It's unfortunate that conditional modules are lacking at the moment. However, it still seems to me that the backbone are all the ODE solvers, is it not possible to split just the ODE solver out and have another package with all the other solvers? If I understand correctly your concern is the following: ODE solvers are needed for the PDE solvers, the SDE and ODE solvers are needed by the SPDE solvers.

Would it not be possible to split the SDE/ODE solvers from PDE related solvers?

Right now it's quite difficult to wrap my head around the code, with how it's currently structured.

[ChrisRackauckas]

How large are all the files? If it's a few MB I'd say it's not worth separating into a new package, which will be less likely to be installed

Benchmarks + Docs + Examples are 92 MB right now. Probably a lot of the .git history too (182 MB). Those two together are >95% of the size.

Would it not be possible to split the SDE/ODE solvers from PDE related solvers?

SDEs and PDEs are essentially the same via the Forward Kolmogorov equations. I am not using that quite yet because I need to make Grids.jl for efficient iterators in order to have the best possible FDM methods (going to make them right the first time, not the MATLAB-style with full matrices). However, that's where that's heading in the very near future. But then again, only one way to solve SDEs is through parabolic equations, and only one way to solve parabolic equations is through SDEs, but the domains for which these methods are best are large enough that I want to fully exploit this relationship (this was the reason for making all of this together in the first place).

Right now the different solvers are all in different folders (src/ode, src/sde, etc.). The only things that are together are the types (the solution types are all together, the problem types are all together, etc) in the src/general directory. Maybe splitting the types by equations would make it easier, and the general folder should just have the functions on the abstract types?

[The thing that's probably hard to understand is the special casing: each optimized algorithm has its own interpolation scheme, and there's different default settings for widely used algorithms, etc. In fact, I think the ODE algorithms are the most in depth in this respect (ode_solve is just an outer wrapper, but it's able to handle so many different cases that it ends up being a lot of code). What I should really do is extend the contributor documentation more: detailing what each file/folder is doing.]

[musm]

Wow yeah....my .julia folder is ballooning out of hand in size, I suspect this isn't just a problem with your package.

For the SDEs issue, I'm thinking of lightweight SDE solvers such as Euler-maruyama and variants. What if the ODE package contained, stiff/regular/geometric/sde solvers (light weight)/?

[ChrisRackauckas]

The .julia folder is ballooning since each package is a Git repository, and Git repos retain history. Tom just was talking with the METADATA maintainers about pruning Plots.jl for this reason. I pruned the repo once before actually (I had to remove sensitive information, i.e. unpublished methods for SDE adaptivity and higher order SDE methods, both waiting on publication). This is an unsafe operation since it means older releases will not work. This also means it doesn't get smaller when you remove things from the repository... it's setup for small packages, but this fact, plus the lack of conditional dependencies, are causing clear problems in larger packages.

I agree that the FEM set should be its own package though. FEM usually will only be PDEs. The interactions between ODEs/SDEs -> PDEs is only for things like square domains which are better suited for FDM methods. The FDM methods should be pretty lightweight: just building an ODE to solve.

[musm]

Sounds good, this issue was inspired by the fact that 95% of the time myself and those I work with only use at best 3 different ode solvers, and 90% of the time it's just ode45

[ChrisRackauckas]

That's the problem: best is dependent on the problem. Dormand-Prince 4/5 solvers (ode45) are only "best" in cases where you need mild tolerances (1e-3 to 1e-5 accuracy?) and a quick interpolant (though this isn't really true, it's based on old wisdom. See how the newer Tsit5 algorithms tends to do better in this domain. Indeed, it was made to simply be a better 4/5 algorithm). If you need more accuracy but still in the floating point range, then Vern6 or Vern7 seem to be your best bet (they have order 6/7 interpolants respectively). If you need more than floating point accuracy and an accurate interpolation, then Vern9 is will do the best, but if you don't need the interpolation but need crazy amounts of accuracy at the specific steps (1e-40), then Feagin14 is best. That's just for nonstiff equations: for stiff equations it can be more complicated (though you can do auto-switching like LSODA). I think I should just make the defaults "smarter" and then the true power of having this framework will come through (i.e. choose the default algorithm based on user chosen tolerances, and do auto-switching to stiff solvers). Part of this would be pretty easy to implement when I find the time.

As you probably know, the story gets even more complicated for PDEs. SDEs are oddly different: I am finding that some of the higher order Runge-Kutta methods I am using with the special adaptivity just tend to work on tough equations better than other methods (since there really isn't a "stiff solver" do to the non-finite moments of the inverse normal distribution). But that's for a later time (after publication). Even then, it's faster to solve a PDE to get the distribution for an SDE than it is to do a Monte Carlo experiment, but only for low-dimensional SDEs. Etc., so many different cases.

---

That said, I think the discussion on separating pieces of a package like this should come down to two factors: size and how easy it breaks (i.e. does that section cause precompilation errors?). The code base is surprisingly small in size since the code itself is entirely Julia (it's about 30,000 lines of Julia, but since it's all text it's small). And since it's all Julia code (no C/Fortran dependencies), it doesn't cause build issues. That's why I haven't seen this as very pertinent.

The crucial things to have as conditional dependencies are things like Sundials and ODEInterface since, because they have to build binary dependencies, they can cause breakage (rather easily). That's why those must be optional. But the FEM methods? It's a 76 KB folder that doesn't cause compilation issues.

[musm]

That's the problem: best is dependent on the problem. Dormand-Prince 4/5 solvers (ode45) are only "best" in cases where you need mild tolerances (1e-3 to 1e-5 accuracy?) and a quick interpolant (though this isn't really true, it's based on old wisdom. See how the newer Tsit5 algorithms tends to do better in this domain. Indeed, it was made to simply be a better 4/5 algorithm). If you need more accuracy but still in the floating point range, then Vern6 or Vern7 seem to be your best bet (they have order 6/7 interpolants respectively). If you need more than floating point accuracy and an accurate interpolation, then Vern9 is will do the best, but if you don't need the interpolation but need crazy amounts of accuracy at the specific steps (1e-40), then Feagin14 is best. That's just for nonstiff equations: for stiff equations it can be more complicated (though you can do auto-switching like LSODA). I think I should just make the defaults "smarter" and then the true power of having this framework will come through (i.e. choose the default algorithm based on user chosen tolerances, and do auto-switching to stiff solvers). Part of this would be pretty easy to implement when I find the time.

Thanks for this. It might be interesting to think about an interface that is agnostic to the user's expertise.

For example: solve( ..... , :identifier)
where :identifier could be :stiff, :fast, :accurate, :general, :geometric
And based on the library writer (i.e. you) the best algorithm is selected. Automatic detection would also be an interesting long term idea. In other words the user could basically call plan_ode( my_ode_system) and this would spit out the best solver to call with solve. This second proposal would likely entail some heuristics to select the best algorithm for the problem at hand.

As you probably know, the story gets even more complicated for PDEs. SDEs are oddly different: I am finding that some of the higher order Runge-Kutta methods I am using with the special adaptivity just tend to work on tough equations better than other methods (since there really isn't a "stiff solver" do to the non-finite moments of the inverse normal distribution). But that's for a later time (after publication). Even then, it's faster to solve a PDE to get the distribution for an SDE than it is to do a Monte Carlo experiment, but only for low-dimensional SDEs. Etc., so many different cases.

That said, I think the discussion on separating pieces of a package like this should come down to two factors: size and how easy it breaks (i.e. does that section cause precompilation errors?). The code base is surprisingly small in size since the code itself is entirely Julia (it's about 30,000 lines of Julia, but since it's all text it's small). And since it's all Julia code (no C/Fortran dependencies), it doesn't cause build issues. That's why I haven't seen this as very pertinent.

The crucial things to have as conditional dependencies are things like Sundials and ODEInterface since, because they have to build binary dependencies, they can cause breakage (rather easily). That's why those must be optional. But the FEM methods? It's a 76 KB folder that doesn't cause compilation issues.

I'm not sure if "size" should account for whether separation makes sense. One could argue that a package should contain only those features that are essential for it's main function and that the smaller this is the more preferable. Whether this is relevant to DifferentialEquations.jl is another story, since we probably won't see any of this happening until Julia's conditional module's features are more mature.

[ChrisRackauckas]

Conditional imports looks like it's coming sooner rather than later: JuliaLang/julia#6195 (comment)

I'll hold out for that.

[ChrisRackauckas]

What names should I use? DiffEqODE, DiffEqSDE, etc.? Or DifferentialEquationsODE, DifferentialEquationsSDE, or should I ask Tony what he'd prefer?

I am thinking that instead I should just start doing the split, and for now just set it up with them requiring each other. So it's more of a semantic breakup until conditional modules come.

[musm]

Some additional thoughts:

ODESolve
SDESolve
PDESolve

Or:
ODESolvers
SDESolvers
PDESolvers

[ChrisRackauckas]

@tkelman, could you chime in? Is there a naming scheme you would suggest? Any you'd reject?

I'm also sending the Julia Praxis squad here for opinions.

[JeffreySarnoff]

Consider the appropriate ways to categorize your application with respect to the subfields it addresses. One way that I have found useful is to analyze the collective from the view that it is a whole, and use the elements that come of analysis as grist for resynthesis .. a remaking of the sense in which it is a whole; then use that as the basis for subselection.