Develop and wrap a PODES Implementation

[ChrisRackauckas]

From #47, it would be nice to have PODES in DifferentialEquations.jl. Adding this functionality can look a lot like how PyDSTools -> Bifurcate -> DifferentialEquations. The workflow to implement this looks something like:

• Make a port of PODES from MATLAB to Julia (discuss licensing issues with the developers first!). Make this as a package, PODES.jl, which is a standalone "does PODES things".
• Wrap the PODES method in DiffEqUncertainty.jl, a package for uncertainty quantification analysis. Make it easy to use DiffEq problem definitions and ParameterizedFunctions with PODES here.
• When complete, document DiffEqUncertainty.jl in DifferentialEquations.jl and officially add it to the ecosystem.

This structure allows for other forms of uncertainty quantification, and lets PODES be a "a backend" for it (probably the first of many). I'm not sure I can get to this for awhile, but I'll build the framework so that way if someone ports PODES then it's already ready.

[finmod]

Good roadmap for uncertainty quantification, particularly the first step of porting Matlab PODES to PODES.jl. The standard problems in PODES should also be ported to DiffEq standard problems. Since the number of equations and the number of parameters to be estimated jumps up from just one or two to quite a few (between five and ten) the benchmark with Matlab and others will become interesting.

[finmod]

Chris, Can you have a look at this link: http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-research/girolami/probints . This is a significant development. Have tried to set up the code in python but dependencies provide the usual headaches. It may be easier in julia

[ChrisRackauckas]

It may be easier in julia

In fact, this is pretty trivial to do in Julia. To get a trajectory for this solution, just modify the standard callback:

probabilistic_callback = @ode_callback begin
  u .= +.(u,σ*dt^(2*alg.order+1)*randn(u))
  @ode_savevalues
end

and then just

sol = solve(prob,callback=probabilistic_callback)

gives one solution to one trajectory according to their numerical methods. This will only work with methods which implement the callbacks and define the order in Alg (all of the OrdinaryDiffEq methods, and hopefully Sundials soon).

To make it even better, we could implement a method which does this in parallel on a cluster using the same method as this:

https://github.com/JuliaDiffEq/StochasticDiffEq.jl/blob/master/src/stochastic_utils.jl#L27

and then on the resulting on Monte Carlo solution type put a plot recipe that makes it output the same plot at as in the paper (just plot all of the Monte Carlo solutions with transparency using the plot recipe for ODE solutions and stack them!). In total, this is under like 100 lines. I'll see when I get around to it but this is, but there's all of the guts for how to implement it.

[ChrisRackauckas]

What are the advantages and disadvantages of this method vs that of PODES? This is really easy to do, while the PODES would take more work. Is there a reason to have both?

[finmod]

The answer to this question is in the slide presentation, in particular slide 8 and the concluding slides on TODO.

This direction of research constitutes outstanding thesis material and you are well-positioned to contribute significantly with the JuliaDiffEq ecosystem which provides a level playing field in Julia to assess the efficiency of old (GBS-RK4 Ernst Hairer et al.) and new solvers. In addition, you are also incorporating the tools to assess the properties of the parameter estimates.

The probabilistic measure approach is an important development of the work of Martin Hairer and it is no coincidence that all the contributors are clustering in the new institute in Warwick. PODES may now appear to be redundant in part but it may still be necessary to validate the algorithm if only to dismiss it in favour of PROBINTS. It is difficult to pass judgment since PODES is in Matlab and PROBINTS rests on MIT Uncertainty Quantification (MUQ) library in Python/C++. Note that Sundials is also used there. There is also an uncertainty quantification library in Julia, Alea.jl, but with no docs.

[ChrisRackauckas]

The new event handling framework will be able to handle this and will be able to be uniformly implemented on "most" solvers, so I'm waiting on that to implement the PROBINTS algorithm in DiffEqSensitivity.jl

[ChrisRackauckas]

For an update on this issue, see SciML/Roadmap#10 (comment)

[ChrisRackauckas]

I put together the ProbInts functionality in DiffEqUncertainty.jl. It's documented here:

http://docs.juliadiffeq.org/latest/analysis/uncertainty_quantification.html

I am not sure that my noise scaling is right. Is it supposed to be in the exponent or not? Could you take a look at it and make sure it matches the paper? I'll open this as a new issue over in DiffEqUncertainty.jl since it's now a developer issue instead of a feature request.

[gerrygralton]

Chris, Can you have a look at this link: http://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-research/girolami/probints . This is a significant development. Have tried to set up the code in python but dependencies provide the usual headaches. It may be easier in julia

Hey, do either of you, @ChrisRackauckas or @finmod, still have a copy of the code that was meant to be at that website? I'm interested in implementing the PDE stuff (Section 5) but the website has been taken down and the paper in that section is light on detail. I asked Girolami but he doesn't have a copy and I can't get in touch with Conrad.
Cheers,
Gerry

[ChrisRackauckas]

I do not have it.

[finmod]

You could follow this link: https://www.asc.ohio-state.edu/chkrebtii.1/podes/Example_Lorenz_IVP.html which leads to: https://github.com/ochkrebtii/uqdes. This must be a consolidated version of the original PODES file published on the warwick.ac.uk site. All of this is in MATLAB. Finally, if you fail to find happiness there. I have the original PODES folder in Matlab but it is 650 MB! You would have to tell me where to upload it.

More generally, do give a look at DiffEqUncertainty.jl in Julia in the SciML ecosystem. You may never look back at anything back!