Merge branch 'master' into W2Ito1

SciML · Oct 29, 2024 · 402e10b · 402e10b
2 parents e4abed8 + 2930ba8
commit 402e10b
Show file tree

Hide file tree

Showing 17 changed files with 207 additions and 100 deletions.
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -21,8 +21,9 @@ jobs:
           - AlgConvergence3
           - WeakConvergence1
         version:
+          - 'lts'
           - '1'
-          - '1.6'
+          - 'pre'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "StochasticDiffEq"
 uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "6.65.1"
+version = "6.70.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -10,6 +10,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
@@ -34,9 +35,10 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Adapt = "3, 4"
 ArrayInterface = "6, 7"
 DataStructures = "0.18"
-DiffEqBase = "6.130.1"
+DiffEqBase = "6.154"
 DiffEqNoiseProcess = "5.13"
 DocStringExtensions = "0.8, 0.9"
+FastPower = "1"
 FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 JumpProcesses = "9"
@@ -45,18 +47,18 @@ LinearAlgebra = "1.6"
 Logging = "1.6"
 MuladdMacro = "0.2.1"
 NLsolve = "4"
-OrdinaryDiffEq = "6.52"
+OrdinaryDiffEq = "6.87"
 Random = "1.6"
 RandomNumbers = "1.5.3"
 RecursiveArrayTools = "2, 3"
 Reexport = "0.2, 1.0"
-SciMLBase = "2.0.6"
+SciMLBase = "2.51"
 SciMLOperators = "0.2.9, 0.3"
 SparseArrays = "1.6"
 SparseDiffTools = "2"
 StaticArrays = "0.11, 0.12, 1.0"
 UnPack = "0.1, 1.0"
-julia = "1.6"
+julia = "1.10"
 
 [extras]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"

diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl
@@ -14,7 +14,7 @@ using DocStringExtensions
          beta2_default, beta1_default, gamma_default,
          qmin_default, qmax_default, qsteady_min_default, qsteady_max_default,
          stepsize_controller!, accept_step_controller, step_accept_controller!,
-         step_reject_controller!, PIController, DummyController
+         step_reject_controller!, PIController, DummyController, issplit
 
   using UnPack, RecursiveArrayTools, DataStructures
   using DiffEqNoiseProcess, Random, ArrayInterface
@@ -35,6 +35,8 @@ using DocStringExtensions
   using LinearAlgebra, Random
 
   import ForwardDiff.Dual
+
+  import FastPower
 
   import DiffEqBase: step!, initialize!, DEAlgorithm,
                      AbstractSDEAlgorithm, AbstractRODEAlgorithm, DEIntegrator, AbstractDiffEqInterpolation,
@@ -47,7 +49,8 @@ using DocStringExtensions
                      resize_non_user_cache!,deleteat_non_user_cache!,addat_non_user_cache!,
                      terminate!,get_du, get_dt,get_proposed_dt,set_proposed_dt!,
                      u_modified!,savevalues!,add_tstop!,add_saveat!,set_reltol!,
-                     set_abstol!, postamble!, last_step_failed, has_Wfact, has_jac
+                     set_abstol!, postamble!, last_step_failed, has_Wfact, has_jac,
+                     get_tstops, get_tstops_array, get_tstops_max
 
   using DiffEqBase: check_error!, is_diagonal_noise, @..
 

diff --git a/src/alg_utils.jl b/src/alg_utils.jl
@@ -134,22 +134,22 @@ beta1_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm
 
 isdtchangeable(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true
 
-alg_interpretation(alg::StochasticDiffEqAlgorithm) = :Ito
-alg_interpretation(alg::EulerHeun) = :Stratonovich
-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) = 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
-
-alg_interpretation(alg::RS1) = :Stratonovich
-alg_interpretation(alg::RS2) = :Stratonovich
-
-alg_interpretation(alg::NON) = :Stratonovich
-alg_interpretation(alg::COM) = :Stratonovich
-alg_interpretation(alg::NON2) = :Stratonovich
+SciMLBase.alg_interpretation(alg::StochasticDiffEqAlgorithm) = SciMLBase.AlgorithmInterpretation.Ito
+SciMLBase.alg_interpretation(alg::EulerHeun) = SciMLBase.AlgorithmInterpretation.Stratonovich
+SciMLBase.alg_interpretation(alg::LambaEulerHeun) = SciMLBase.AlgorithmInterpretation.Stratonovich
+SciMLBase.alg_interpretation(alg::KomBurSROCK2) = SciMLBase.AlgorithmInterpretation.Stratonovich
+SciMLBase.alg_interpretation(alg::RKMil{interpretation}) where {interpretation} = interpretation
+SciMLBase.alg_interpretation(alg::SROCK1{interpretation,E}) where {interpretation,E} = interpretation
+SciMLBase.alg_interpretation(alg::RKMilCommute) = alg.interpretation
+SciMLBase.alg_interpretation(alg::RKMilGeneral) = alg.interpretation
+SciMLBase.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
+
+SciMLBase.alg_interpretation(alg::RS1) = SciMLBase.AlgorithmInterpretation.Stratonovich
+SciMLBase.alg_interpretation(alg::RS2) = SciMLBase.AlgorithmInterpretation.Stratonovich
+
+SciMLBase.alg_interpretation(alg::NON) = SciMLBase.AlgorithmInterpretation.Stratonovich
+SciMLBase.alg_interpretation(alg::COM) = SciMLBase.AlgorithmInterpretation.Stratonovich
+SciMLBase.alg_interpretation(alg::NON2) = SciMLBase.AlgorithmInterpretation.Stratonovich
 
 alg_compatible(prob, alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true
 alg_compatible(prob, alg::StochasticDiffEqAlgorithm) = false

diff --git a/src/algorithms.jl b/src/algorithms.jl
@@ -77,49 +77,49 @@ Springer. Berlin Heidelberg (2011)
 
 RKMil: Nonstiff Method
 An explicit Runge-Kutta discretization of the strong order 1.0 Milstein method.
-Defaults to solving the Ito problem, but RKMil(interpretation=:Stratonovich) makes it solve the Stratonovich problem.
+Defaults to solving the Ito problem, but RKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich) makes it solve the Stratonovich problem.
 Only handles scalar and diagonal noise.
 """
 struct RKMil{interpretation} <: StochasticDiffEqAdaptiveAlgorithm end
-RKMil(;interpretation=:Ito) = RKMil{interpretation}()
+RKMil(;interpretation=SciMLBase.AlgorithmInterpretation.Ito) = RKMil{interpretation}()
 
 """
 Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations.
 Springer. Berlin Heidelberg (2011)
 
 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.
+Defaults to solving the Ito problem, but RKMilCommute(interpretation=SciMLBase.AlgorithmInterpretation.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
+struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm 
+  interpretation::SciMLBase.AlgorithmInterpretation.T
   ii_approx::T
 end
-RKMilCommute(;interpretation=:Ito, ii_approx=IICommutative()) = RKMilCommute(interpretation,ii_approx)
+RKMilCommute(;interpretation=SciMLBase.AlgorithmInterpretation.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=IILevyArea()
+RKMilGeneral(;interpretation=SciMLBase.AlgorithmInterpretation.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 interpretation between SciMLBase.AlgorithmInterpretation.Ito and SciMLBase.AlgorithmInterpretation.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
+  interpretation::SciMLBase.AlgorithmInterpretation.T
   ii_approx::T
   c::Int
   p::TruncationType
 end
 
-function RKMilGeneral(;interpretation=:Ito,ii_approx=IILevyArea(), c=1, p=nothing, dt=nothing)
+function RKMilGeneral(;interpretation=SciMLBase.AlgorithmInterpretation.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{typeof(ii_approx), typeof(p)}(interpretation, ii_approx, c, p)
@@ -160,13 +160,13 @@ struct WangLi3SMil_F <: StochasticDiffEqAlgorithm end
 """
 SROCK1: S-ROCK Method
 Is a fixed step size stabilized explicit method for stiff problems.
-Defaults to solving th Ito problem but SROCK1(interpretation=:Stratonovich) can make it solve the Stratonovich problem.
+Defaults to solving th Ito problem but SROCK1(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich) can make it solve the Stratonovich problem.
 Strong order of convergence is 0.5 and weak order 1, but is optimised to get order 1 in case os scalar/diagonal noise.
 """
 struct SROCK1{interpretation,E} <: StochasticDiffEqAlgorithm
   eigen_est::E
 end
-SROCK1(;interpretation=:Ito,eigen_est=nothing) = SROCK1{interpretation,typeof(eigen_est)}(eigen_est)
+SROCK1(;interpretation=SciMLBase.AlgorithmInterpretation.Ito,eigen_est=nothing) = SROCK1{interpretation,typeof(eigen_est)}(eigen_est)
 
 # Weak Order 2
 for Alg in [:SROCK2, :KomBurSROCK2, :SROCKC2]
@@ -714,7 +714,7 @@ ImplicitEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
 ImplicitRKMil: Stiff Method
 An order 1.0 drift-implicit method.
 This is a theta method which defaults to theta=1 or the Trapezoid method on the drift term.
-Defaults to solving the Ito problem, but ImplicitRKMil(interpretation=:Stratonovich) makes it solve the Stratonovich problem.
+Defaults to solving the Ito problem, but ImplicitRKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich) makes it solve the Stratonovich problem.
 This method defaults to symplectic=false, but when true and theta=1/2 this is the implicit Midpoint method on the drift term and is symplectic in distribution.
 Handles diagonal and scalar noise. Uses a 1.5/2.0 heuristic for adaptive time stepping.
 """
@@ -734,7 +734,7 @@ ImplicitRKMil(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                           extrapolant=:constant,
                           theta = 1,symplectic = false,
                           new_jac_conv_bound = 1e-3,
-                          controller = :Predictive,interpretation=:Ito) =
+                          controller = :Predictive,interpretation=SciMLBase.AlgorithmInterpretation.Ito) =
                           ImplicitRKMil{chunk_size,autodiff,
                           typeof(linsolve),typeof(precs),diff_type,
                           OrdinaryDiffEq._unwrap_val(standardtag),

diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl
@@ -390,3 +390,7 @@ function DiffEqBase.set_u!(integrator::SDEIntegrator, u)
   integrator.u = u
   u_modified!(integrator, true)
 end
+
+DiffEqBase.get_tstops(integ::SDEIntegrator) = integ.opts.tstops
+DiffEqBase.get_tstops_array(integ::SDEIntegrator) = get_tstops(integ).valtree
+DiffEqBase.get_tstops_max(integ::SDEIntegrator) = maximum(get_tstops_array(integ))
diff --git a/src/integrators/stepsize_controllers.jl b/src/integrators/stepsize_controllers.jl
@@ -1,7 +1,7 @@
 
 function stepsize_controller!(integrator::SDEIntegrator, controller::PIController, alg)
-    integrator.q11 = DiffEqBase.value(DiffEqBase.fastpow(integrator.EEst,controller.beta1))
-    integrator.q = DiffEqBase.value(integrator.q11/DiffEqBase.fastpow(integrator.qold,controller.beta2))
+    integrator.q11 = DiffEqBase.value(FastPower.fastpower(integrator.EEst,controller.beta1))
+    integrator.q = DiffEqBase.value(integrator.q11/FastPower.fastpower(integrator.qold,controller.beta2))
     @fastmath integrator.q = DiffEqBase.value(max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),integrator.q/integrator.opts.gamma)))
 end
 

diff --git a/src/perform_step/SROCK_perform_step.jl b/src/perform_step/SROCK_perform_step.jl
@@ -14,7 +14,7 @@
   cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
   ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))
 
-  if alg_interpretation(integrator.alg) == :Stratonovich
+  if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Stratonovich
     α  = cosh(mdeg*cosh_inv)/(2*ω₀*cosh((mdeg-1)*cosh_inv))
     γ  = 1/(2*α)
     β  = -γ
@@ -42,7 +42,7 @@
     k = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
 
     u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
-    if (i > mdeg - 2) && alg_interpretation(integrator.alg) == :Stratonovich
+    if (i > mdeg - 2) && SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Stratonovich
       if i == mdeg - 1
         gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
         if W.dW isa Number || !is_diagonal_noise(integrator.sol.prob)
@@ -58,7 +58,7 @@
           u .+= (β .* gₘ₋₂ .+ γ .* gₘ₋₁) .* W.dW
         end
       end
-    elseif (i == mdeg) && alg_interpretation(integrator.alg) == :Ito
+    elseif (i == mdeg) && SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito
       if W.dW isa Number
         gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
         uᵢ₋₂ = uᵢ₋₁ + sqrt(abs(dt))*gₘ₋₂
@@ -105,7 +105,7 @@ end
   cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
   ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))
 
-  if alg_interpretation(integrator.alg) == :Stratonovich
+  if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Stratonovich
     α  = cosh(mdeg*cosh_inv)/(2*ω₀*cosh((mdeg-1)*cosh_inv))
     γ  = 1/(2*α)
     β  = -γ
@@ -132,7 +132,7 @@ end
     κ = - Tᵢ₋₂/Tᵢ
     integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
     @.. u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
-    if (i > mdeg - 2) && alg_interpretation(integrator.alg) == :Stratonovich
+    if (i > mdeg - 2) && SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Stratonovich
       if i == mdeg - 1
         integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
         if W.dW isa Number || is_diagonal_noise(integrator.sol.prob)
@@ -152,7 +152,7 @@ end
           @.. u += γ*k
         end
       end
-    elseif (i == mdeg) && alg_interpretation(integrator.alg) == :Ito
+    elseif (i == mdeg) && SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito
       if W.dW isa Number || is_diagonal_noise(integrator.sol.prob)
         integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
         @.. uᵢ₋₂ = uᵢ₋₁ + sqrt(abs(dt))*gₘ₋₂