Merge pull request #459 from SciML/levy_area

Switch to LevyArea.jl

Author: ChrisRackauckas

.buildkite/pipeline.yml

@@ -41,6 +41,7 @@ steps:
     agents:
       os: "linux"
       queue: "juliaecosystem"
+      exclusive: true
     env:
       GROUP: 'SROCKC2WeakConvergence'
     timeout_in_minutes: 240

Project.toml

@@ -14,6 +14,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LevyArea = "2d8b4e74-eb68-11e8-0fb9-d5eb67b50637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
@@ -39,6 +40,7 @@ DocStringExtensions = "0.8"
 FillArrays = "0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12"
 FiniteDiff = "2"
 ForwardDiff = "0.10.3"
+LevyArea = "1.0.0"
 MuladdMacro = "0.2.1"
 NLsolve = "4"
 OrdinaryDiffEq = "6.4"

src/StochasticDiffEq.jl

@@ -63,6 +63,7 @@ using DocStringExtensions
 
   using SparseDiffTools: forwarddiff_color_jacobian!, ForwardColorJacCache
 
+  using LevyArea
 
   const CompiledFloats = Union{Float32,Float64}
 
@@ -164,7 +165,7 @@ using DocStringExtensions
 
   export RandomEM
 
-  export IteratedIntegralApprox, IICommutative, IIWiktorsson
+  export IteratedIntegralApprox, IICommutative, IILevyArea
 
   #General Functions
   export solve, init, solve!, step!

src/alg_utils.jl

@@ -123,7 +123,7 @@ alg_interpretation(alg::LambaEulerHeun) = :Stratonovich
 alg_interpretation(alg::KomBurSROCK2) = :Stratonovich
 alg_interpretation(alg::RKMil{interpretation}) where {interpretation} = interpretation
 alg_interpretation(alg::SROCK1{interpretation,E}) where {interpretation,E} = interpretation
-alg_interpretation(alg::RKMilCommute{interpretation}) where {interpretation} = interpretation
+alg_interpretation(alg::RKMilCommute) = alg.interpretation
 alg_interpretation(alg::RKMilGeneral) = alg.interpretation
 alg_interpretation(alg::ImplicitRKMil{CS,AD,F,P,FDT,ST,CJ,N,T2,Controller,interpretation}) where {CS,AD,F,P,FDT,ST,CJ,N,T2,Controller,interpretation} = interpretation
 

src/algorithms.jl

@@ -19,7 +19,7 @@ abstract type StochasticDiffEqJumpNewtonDiffusionAdaptiveAlgorithm{CS,AD,FDT,ST,
 
 abstract type IteratedIntegralApprox end
 struct IICommutative <:  IteratedIntegralApprox end
-struct IIWiktorsson <:  IteratedIntegralApprox end
+struct IILevyArea <:  IteratedIntegralApprox end
 
 ################################################################################
 
@@ -91,31 +91,40 @@ RKMilCommute: Nonstiff Method
 An explicit Runge-Kutta discretization of the strong order 1.0 Milstein method for commutative noise problems.
 Defaults to solving the Ito problem, but RKMilCommute(interpretation=:Stratonovich) makes it solve the Stratonovich problem.
 Uses a 1.5/2.0 error estimate for adaptive time stepping.
+Default: ii_approx=IICommutative() does not approximate the Levy area.
 """
-struct RKMilCommute{interpretation} <: StochasticDiffEqAdaptiveAlgorithm end
-RKMilCommute(;interpretation=:Ito) = RKMilCommute{interpretation}()
+struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm 
+  interpretation::Symbol
+  ii_approx::T
+end
+RKMilCommute(;interpretation=:Ito, ii_approx=IICommutative()) = RKMilCommute(interpretation,ii_approx)
 
 """
 Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
 Springer. Berlin Heidelberg (2011)
 
 RKMilGeneral: Nonstiff Method
-RKMilGeneral(;interpretation=:Ito, ii_approx=IIWiktorsson()
+RKMilGeneral(;interpretation=:Ito, ii_approx=IILevyArea()
 An explicit Runge-Kutta discretization of the strong order 1.0 Milstein method for general non-commutative noise problems.
 Allows for a choice of interpretation between :Ito and :Stratonovich.
 Allows for a choice of iterated integral approximation.
+Default: ii_approx=IILevyArea() uses LevyArea.jl to choose optimal algorithm. See
+Kastner, F. and Rößler, A., arXiv: 2201.08424
+Kastner, F. and Rößler, A., LevyArea.jl, 10.5281/ZENODO.5883748, https://github.com/stochastics-uni-luebeck/LevyArea.jl
 """
-struct RKMilGeneral{T<:IteratedIntegralApprox, TruncationType} <: StochasticDiffEqAdaptiveAlgorithm
+struct RKMilGeneral{T, TruncationType} <: StochasticDiffEqAdaptiveAlgorithm
   interpretation::Symbol
   ii_approx::T
   c::Int
   p::TruncationType
 end
-function RKMilGeneral(;interpretation=:Ito, ii_approx=IIWiktorsson(), c=1, p=nothing, dt=nothing)
+
+function RKMilGeneral(;interpretation=:Ito,ii_approx=IILevyArea(), c=1, p=nothing, dt=nothing)
   γ = 1//1
   p==true && (p = Int(floor(c*dt^(1//1-2//1*γ)) + 1))
-  RKMilGeneral(interpretation, ii_approx, c, p)
+  RKMilGeneral{typeof(ii_approx), typeof(p)}(interpretation, ii_approx, c, p)
 end
+
 """
 WangLi3SMil_A: Nonstiff Method
 Fixed step-size explicit 3-stage Milstein methods for Ito problem with strong and weak order 1.0

src/caches/basic_method_caches.jl

@@ -115,43 +115,44 @@ function alg_cache(alg::RKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototy
   RKMilCache(u,uprev,du1,du2,K,tmp,L)
 end
 
-struct RKMilCommuteConstantCache{WikType} <: StochasticDiffEqConstantCache
-  WikJ::WikType
+struct RKMilCommuteConstantCache{JalgType} <: StochasticDiffEqConstantCache
+  Jalg::JalgType
 end
-@cache struct RKMilCommuteCache{uType,rateType,rateNoiseType,WikType} <: StochasticDiffEqMutableCache
+@cache struct RKMilCommuteCache{uType,rateType,rateNoiseType,JalgType} <: StochasticDiffEqMutableCache
   u::uType
   uprev::uType
   du1::rateType
   du2::rateType
   K::rateType
   gtmp::rateNoiseType
   L::rateNoiseType
-  WikJ::WikType
+  Jalg::JalgType
   mil_correction::rateType
   Kj::uType
   Dgj::rateNoiseType
   tmp::uType
 end
 
-function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} WikJ = WikJ = get_WikJ(ΔW,prob,alg)
-  RKMilCommuteConstantCache{typeof(WikJ)}(WikJ)
+function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} 
+  Jalg = get_Jalg(ΔW,dt,prob,alg)
+  RKMilCommuteConstantCache{typeof(Jalg)}(Jalg)
 end
 
 function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
   du1 = zero(rate_prototype); du2 = zero(rate_prototype)
   K = zero(rate_prototype); gtmp = zero(noise_rate_prototype);
   L = zero(noise_rate_prototype); tmp = zero(rate_prototype)
-  WikJ = get_WikJ(ΔW,prob,alg)
+  Jalg = get_Jalg(ΔW,dt,prob,alg)
   mil_correction = zero(rate_prototype)
   Kj = zero(u); Dgj = zero(noise_rate_prototype)
-  RKMilCommuteCache(u,uprev,du1,du2,K,gtmp,L,WikJ,mil_correction,Kj,Dgj,tmp)
+  RKMilCommuteCache(u,uprev,du1,du2,K,gtmp,L,Jalg,mil_correction,Kj,Dgj,tmp)
 end
 
-struct RKMilGeneralConstantCache{WikType} <: StochasticDiffEqConstantCache
-  WikJ::WikType
+struct RKMilGeneralConstantCache{JalgType} <: StochasticDiffEqConstantCache
+  Jalg::JalgType
 end
 
-@cache struct RKMilGeneralCache{uType, rateType, rateNoiseType, WikType} <: StochasticDiffEqMutableCache
+@cache struct RKMilGeneralCache{uType, rateType, rateNoiseType, JalgType} <: StochasticDiffEqMutableCache
   u::uType
   uprev::uType
   tmp::uType
@@ -161,12 +162,12 @@ end
   L::rateNoiseType
   mil_correction::uType
   ggprime::rateNoiseType
-  WikJ::WikType
+  Jalg::JalgType
 end
 
-function alg_cache(alg::RKMilGeneral,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
-  WikJ = get_WikJ(ΔW,prob,alg)
-  RKMilGeneralConstantCache{typeof(WikJ)}(WikJ)
+function alg_cache(alg::RKMilGeneral, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype, jump_rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  Jalg = get_Jalg(ΔW, dt, prob, alg)
+  RKMilGeneralConstantCache{typeof(Jalg)}(Jalg)
 end
 
 function alg_cache(alg::RKMilGeneral,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
@@ -177,6 +178,6 @@ function alg_cache(alg::RKMilGeneral,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_
   L = zero(noise_rate_prototype)
   mil_correction = zero(u)
   ggprime = zero(noise_rate_prototype)
-  WikJ = get_WikJ(ΔW,prob,alg)
-  RKMilGeneralCache{typeof(u), typeof(rate_prototype), typeof(noise_rate_prototype), typeof(WikJ)}(u, uprev, tmp, du₁, du₂, K, L, mil_correction, ggprime, WikJ)
+  Jalg = get_Jalg(ΔW,dt,prob,alg)
+  RKMilGeneralCache{typeof(u), typeof(rate_prototype), typeof(noise_rate_prototype), typeof(Jalg)}(u, uprev, tmp, du₁, du₂, K, L, mil_correction, ggprime, Jalg)
 end