change names, delete path convergence test

Author: frankschae

src/algorithms.jl

@@ -91,6 +91,7 @@ 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{T} <: StochasticDiffEqAdaptiveAlgorithm 
   interpretation::Symbol
@@ -107,6 +108,9 @@ 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, TruncationType} <: StochasticDiffEqAdaptiveAlgorithm
   interpretation::Symbol

src/caches/basic_method_caches.jl

@@ -115,44 +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 = get_WikJ(ΔW,dt,prob,alg)
-  RKMilCommuteConstantCache{typeof(WikJ)}(WikJ)
+  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,dt,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
@@ -162,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,dt,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}
@@ -178,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,dt,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

src/iterated_integrals.jl

@@ -5,53 +5,50 @@ by the J.
 
 abstract type AbstractJ end
 
-# Wiktorssson Approximation (for diagonal and commutative noise processes where the Levy area approximation is not required.)
+# Iterated stochastic integral (for diagonal and commutative noise processes where the Levy area approximation is not required.)
+abstract type AbstractJDiagonal <: AbstractJ end
+abstract type AbstractJCommute <: AbstractJ end
 
-abstract type AbstractWikJ <: AbstractJ end
-abstract type AbstractWikJDiagonal <: AbstractWikJ end
-abstract type AbstractWikJCommute <: AbstractWikJ end
+struct JDiagonal_oop <: AbstractJDiagonal end
 
-struct WikJDiagonal_oop <: AbstractWikJDiagonal end
-
-mutable struct WikJDiagonal_iip{WikJType} <: AbstractWikJDiagonal
-  WikJ::WikJType
-  WikJDiagonal_iip(ΔW) = new{typeof(ΔW)}(false .* ΔW .* ΔW)
+mutable struct JDiagonal_iip{JType} <: AbstractJDiagonal
+  J::JType
+  JDiagonal_iip(ΔW) = new{typeof(ΔW)}(false .* ΔW .* ΔW)
 end
 
-struct WikJCommute_oop <: AbstractWikJCommute end
+struct JCommute_oop <: AbstractJCommute end
 
-mutable struct WikJCommute_iip{WikJType} <: AbstractWikJCommute
-  WikJ::WikJType
-  function WikJCommute_iip(ΔW)
-    WikJ = false .* vec(ΔW) .* vec(ΔW)'
-    new{typeof(WikJ)}(WikJ)
+mutable struct JCommute_iip{JType} <: AbstractJCommute
+  J::JType
+  function JCommute_iip(ΔW)
+    J = false .* vec(ΔW) .* vec(ΔW)'
+    new{typeof(J)}(J)
   end
 end
 
-function get_iterated_I(dt, dW, dZ, Wik::WikJDiagonal_oop, p=nothing, C=1, γ=1//1)
-  WikJ = 1//2 .* dW .* dW
-  WikJ
+function get_iterated_I(dt, dW, dZ, alg::JDiagonal_oop, p=nothing, c=1, γ=1//1)
+  1//2 .* dW .* dW
 end
 
-function get_iterated_I!(dt, dW, dZ, Wik::WikJDiagonal_iip, p=nothing, C=1, γ=1//1)
-  @unpack WikJ = Wik
+function get_iterated_I!(dt, dW, dZ, alg::JDiagonal_iip, p=nothing, c=1, γ=1//1)
+  @unpack J = alg
   if typeof(dW) <: Number
-    Wik.WikJ = 1//2 .* dW .^ 2
+    alg.J = 1//2 .* dW .^ 2
   else
-    @.. WikJ = 1//2 * dW^2
+    @.. J = 1//2 * dW^2
   end
   return nothing
 end
 
-function get_iterated_I(dt, dW, dZ, Wik::WikJCommute_oop, p=nothing, C=1, γ=1//1)
-  WikJ = 1//2 .* vec(dW) .* vec(dW)'
-  WikJ
+function get_iterated_I(dt, dW, dZ, alg::JCommute_oop, p=nothing, c=1, γ=1//1)
+  J = 1//2 .* vec(dW) .* vec(dW)'
+  J
 end
 
-function get_iterated_I!(dt, dW, dZ, Wik::WikJCommute_iip, p=nothing, C=1, γ=1//1)
-  @unpack WikJ = Wik
-  mul!(WikJ, vec(dW), vec(dW)')
-  @.. WikJ *= 1//2
+function get_iterated_I!(dt, dW, dZ, alg::JCommute_iip, p=nothing, c=1, γ=1//1)
+  @unpack J = alg
+  mul!(J, vec(dW), vec(dW)')
+  @.. J *= 1//2
   return nothing
 end
 
@@ -62,59 +59,59 @@ function get_iterated_I(dt, dW, dZ, alg::LevyArea.AbstractIteratedIntegralAlgori
     ε = c * dt^(γ + 1//2)
     p = terms_needed(length(dW), dt, ε, alg, MaxL2())
   end
-  I = LevyArea.levyarea(dW/√dt, p, alg)
+  I = LevyArea.levyarea(dW / √dt, p, alg)
   I .= 1//2 * dW .* dW' .+ dt .* I
 end
 
-mutable struct IteratedIntegralAlgorithm_iip{WikJType,algType} <: LevyArea.AbstractIteratedIntegralAlgorithm
-  WikJ::WikJType
-  alg::algType
-  function IteratedIntegralAlgorithm_iip(ΔW, alg)
-    WikJ = false .* vec(ΔW) .* vec(ΔW)'
-    new{typeof(WikJ),typeof(alg)}(WikJ, alg)
+mutable struct IteratedIntegralAlgorithm_iip{JType,levyalgType} <: LevyArea.AbstractIteratedIntegralAlgorithm
+  J::JType
+  levyalg::levyalgType
+  function IteratedIntegralAlgorithm_iip(ΔW, levyalg)
+    J = false .* vec(ΔW) .* vec(ΔW)'
+    new{typeof(J),typeof(levyalg)}(J, levyalg)
   end
 end
 
-function get_iterated_I!(dt, dW, dZ, Wik::IteratedIntegralAlgorithm_iip, p=nothing, c=1, γ=1//1)
-  @unpack WikJ, alg = Wik
+function get_iterated_I!(dt, dW, dZ, alg::IteratedIntegralAlgorithm_iip, p=nothing, c=1, γ=1//1)
+  @unpack J, levyalg = alg
   if isnothing(p)
     ε = c * dt^(γ + 1//2)
-    p = terms_needed(length(dW), dt, ε, alg, MaxL2())
+    p = terms_needed(length(dW), dt, ε, levyalg, MaxL2())
   end
-  WikJ .= LevyArea.levyarea(dW/√dt, p, alg)
-  WikJ .= 1//2 * dW .* dW' .+ dt .* WikJ
+  J .= LevyArea.levyarea(dW / √dt, p, levyalg)
+  J .= 1//2 * dW .* dW' .+ dt .* J
   return nothing
 end
 
 # Default algorithms, keep KPWJ_oop() to have a non-mutating version
-function get_WikJ(ΔW, dt, prob, alg)
+function get_Jalg(ΔW, dt, prob, alg)
   if alg.ii_approx isa IILevyArea
     if isinplace(prob)
       if typeof(ΔW) <: Number || is_diagonal_noise(prob)
-        return WikJDiagonal_iip(ΔW)
+        return JDiagonal_iip(ΔW)
       else
         # optimal_algorithm(dim, stepsize, eps=stepsize^(3/2), norm=MaxL2())
         return IteratedIntegralAlgorithm_iip(ΔW, LevyArea.optimal_algorithm(length(ΔW), dt))
       end
     else
       if typeof(ΔW) <: Number || is_diagonal_noise(prob)
-        return WikJDiagonal_oop()
+        return JDiagonal_oop()
       else
         return LevyArea.optimal_algorithm(length(ΔW), dt)
       end
     end
   elseif alg.ii_approx isa IICommutative
     if isinplace(prob)
       if typeof(ΔW) <: Number || is_diagonal_noise(prob)
-        return WikJDiagonal_iip(ΔW)
+        return JDiagonal_iip(ΔW)
       else
-        return WikJCommute_iip(ΔW)
+        return JCommute_iip(ΔW)
       end
     else
       if typeof(ΔW) <: Number || is_diagonal_noise(prob)
-        return WikJDiagonal_oop()
+        return JDiagonal_oop()
       else
-        return WikJCommute_oop()
+        return JCommute_oop()
       end
     end
   else
@@ -123,6 +120,6 @@ function get_WikJ(ΔW, dt, prob, alg)
 end
 
 # Specific Levy area alg for an SDE solver
-# function StochasticDiffEq.get_WikJ(ΔW,prob,alg::SOLVER)
+# function StochasticDiffEq.get_Jalg(ΔW,prob,alg::SOLVER)
 #  return MronRoe()
 # end