Integrate LevyArea.jl into StochasticDiffEq.jl

[ChrisRackauckas]

https://github.com/stochastics-uni-luebeck/LevyArea.jl https://arxiv.org/abs/2201.08424

This paper and library is great @fkastner! As noted in the paper, LevyArea.jl has new and faster implementations of the Levy area calculations. We might as well use them in our higher order integrators. @frankschae could you take a look at how difficult that would be? I assume it wouldn't be too difficult as it should be a direct stand-in for StochasticDiffEq.get_iterated_I!(h,W,nothing,iip,p). Note that we should then update the citation materials in each of the higher order algorithms which use it to point to their paper.

(P.S. Note that you're citing most of the Julia libraries incorrectly. Julia has a standard of using Citation.bib files, so it's recommended you cite the associated papers for the library. Here for example it's: https://github.com/SciML/StochasticDiffEq.jl/blob/master/CITATION.bib. https://github.com/SebastianM-C/PkgCite.jl helps automate this.)

[fkastner]

Hi Chris, thanks for your interest in our package.

The integration should probably look something like this:

using LevyArea
import StochasticDiffEq: get_iterated_I

function get_iterated_I(dt, dW, ::Nothing, alg::LevyArea.AbstractIteratedIntegralAlgorithm, p=nothing, c=1, γ=1//1)
    if isnothing(p)
        ε = c*dt^(γ+1//2)
        p = terms_needed(length(dW), dt, ε, alg, MaxL2())
    end
    I = LevyArea.levyarea(dW/√dt, p, alg)
    I .= 0.5.*dW.*dW' .+ dt.*I
end

and then you need some way to choose between the different algorithms.

For a simple test I just overwrote get_WikJ:

using StochasticDiffEq

function StochasticDiffEq.get_WikJ(ΔW,prob,alg::RKMilGeneral)
    return MronRoe()
end

u₀=[0.0,0.0,0.0]
f(u,p,t) = [0.0,0.0,0.0]
g(u,p,t) = [1.0 0.0;0.0 1.0;0.0 u[1]]
dt = 1//2^10
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u₀,(0.0,1.0); noise_rate_prototype=zeros(3,2))
sol = solve(prob, RKMilGeneral(), dt=dt, adaptive=false);

BTW. is it correct that you currently never choose the Wiktorsson approximation even if I pass ii_approx=IIWiktorsson() to RKMilGeneral? Why not?

Re the citation: It is only StochasticDiffEq.jl where we forgot to cite the paper, right? Or is there another package we cited incorrectly?

[ChrisRackauckas]

and then you need some way to choose between the different algorithms.

Cool, I think that passing it along in the SDE algorithm struct would work well.

p = terms_needed(length(dW), dt, ε, alg, MaxL2())

Is that a cache that should be constructed once and updated?

BTW. is it correct that you currently never choose the Wiktorsson approximation even if I pass ii_approx=IIWiktorsson() to RKMilGeneral? Why not?

If that is true that's a bug @frankschae

Re the citation: It is only StochasticDiffEq.jl where we forgot to cite the paper, right? Or is there another package we cited incorrectly?

No, the citation in the paper:

[49] SciML. StochasticDiffEq.jl. Version 6.37.1. Aug. 11, 2021. url: https://github.com/SciML/StochasticDiffEq.jl.

Does not match the suggested citation of the repository, which is:

@article{DifferentialEquations.jl-2017,
 author = {Rackauckas, Christopher and Nie, Qing},
 doi = {10.5334/jors.151},
 journal = {The Journal of Open Research Software},
 keywords = {Applied Mathematics},
 note = {Exported from https://app.dimensions.ai on 2019/05/05},
 number = {1},
 pages = {},
 title = {DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia},
 url = {https://app.dimensions.ai/details/publication/pub.1085583166 and http://openresearchsoftware.metajnl.com/articles/10.5334/jors.151/galley/245/download/},
 volume = {5},
 year = {2017}
}

@article{rackauckas2017adaptive,
  title={Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory},
  author={Rackauckas, Christopher and Nie, Qing},
  journal={Discrete and continuous dynamical systems. Series B},
  volume={22},
  number={7},
  pages={2731},
  year={2017},
  publisher={NIH Public Access}
}

from the Citation.bib file. I think there's a few other cases where the citation is pointing to the github repo itself and not the paper for the github repo.

[frankschae]

Thanks for the code example @fkastner!

BTW. is it correct that you currently never choose the Wiktorsson approximation even if I pass ii_approx=IIWiktorsson() to RKMilGeneral? Why not?

That sounds bad. I'll take a look.

[fkastner]

p = terms_needed(length(dW), dt, ε, alg, MaxL2())

Is that a cache that should be constructed once and updated?

No, it's the number of terms in the approximation, just like in your implementation. I just wanted to provide something that follows your current interface as closely as possible.

[frankschae]

BTW. is it correct that you currently never choose the Wiktorsson approximation even if I pass ii_approx=IIWiktorsson() to RKMilGeneral? Why not?

This turned out to be poor naming. IIWiktorsson was only used to check if the Levy area has to be approximated (in contrast to IICommutative).

[fkastner]

BTW. is it correct that you currently never choose the Wiktorsson approximation even if I pass ii_approx=IIWiktorsson() to RKMilGeneral? Why not?

This turned out to be poor naming. IIWiktorsson was only used to check if the Levy area has to be approximated (in contrast to IICommutative).

The docs say, that you implemented Wiktorsson's approximation and that it can be chosen using this keyword. And it looks like you did implement it. But, if I understand correctly, you never return one of the WikJGeneral_* types from get_WikJ and thus don't ever choose Wiktorsson's approximation? And there is currently no way to explicitly provide the algorithm it should use for the Levy area approximation.

[frankschae]

Yes, the default choice is not the one due to Wiktorsson but the one due to Kloeden, Platen, Wright. This is an error in the docs. (We should correct it after #459) There is a way - the same way as you did it, i.e., overloading StochasticDiffEq.get_WikJ(ΔW,prob,alg).

[ChrisRackauckas]

Awesome, thanks @fkastner this is a major improvement!