From 96f6ef0a2c82aa4dac416f199189ec9c06cc8e48 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 30 Jun 2024 03:20:25 +0800 Subject: [PATCH 01/36] Fix trait usage downstream --- Project.toml | 2 +- src/StochasticDiffEq.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index f1f8ba5e..1c1980cc 100644 --- a/Project.toml +++ b/Project.toml @@ -45,7 +45,7 @@ LinearAlgebra = "1.6" Logging = "1.6" MuladdMacro = "0.2.1" NLsolve = "4" -OrdinaryDiffEq = "6.52" +OrdinaryDiffEq = "6.85" Random = "1.6" RandomNumbers = "1.5.3" RecursiveArrayTools = "2, 3" diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index 61cff965..d9cc79d3 100644 --- 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 From c2c34b44b624f6f5a9790291c5388cadd6d0578f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 30 Jun 2024 03:28:33 +0800 Subject: [PATCH 02/36] bump Julia version --- .github/workflows/CI.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ce50fc33..d1ff8473 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,7 +22,6 @@ jobs: - WeakConvergence1 version: - '1' - - '1.6' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 From 35fb560f7953f7d9c89a0a1b05429c056c06f0c0 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 30 Jun 2024 03:28:56 +0800 Subject: [PATCH 03/36] bump julia version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1c1980cc..845c7009 100644 --- a/Project.toml +++ b/Project.toml @@ -56,7 +56,7 @@ 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" From 43f1ef9fe4a92ba56af713768763530208da0542 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 30 Jun 2024 05:43:46 +0800 Subject: [PATCH 04/36] fix utility tests --- test/utility_tests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index 2311ff57..f52561df 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -60,7 +60,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator println(Alg) Random.seed!(0); sol1 = solve(prob1, Alg(theta=1); adaptive=false, dt=0.01) Random.seed!(0); sol2 = solve(prob2, Alg(theta=1); adaptive=false, dt=0.01) - @test sol1(1.0) ≈ sol2(1.0) rtol=1e-4 + @test sol1(1.0) ≈ sol2(1.0) rtol=1e-2 Random.seed!(0); sol1_ip = solve(prob1_ip, Alg(theta=1); adaptive=false, dt=0.01) Random.seed!(0); sol2_ip = solve(prob2_ip, Alg(theta=1); adaptive=false, dt=0.01) @test sol1_ip(1.0) ≈ sol2_ip(1.0) rtol=1e-4 @@ -76,7 +76,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator println(SKenCarp) Random.seed!(0); sol1 = solve(prob1, SKenCarp(); adaptive=false, dt=0.01) Random.seed!(0); sol2 = solve(prob2, SKenCarp(); adaptive=false, dt=0.01) - @test sol1(1.0) ≈ sol2(1.0) rtol=1e-4 + @test sol1(1.0) ≈ sol2(1.0) rtol=1e-2 Random.seed!(0); sol1_ip = solve(prob1_ip, SKenCarp(); adaptive=false, dt=0.01) Random.seed!(0); sol2_ip = solve(prob2_ip, SKenCarp(); adaptive=false, dt=0.01) @test sol1_ip(1.0) ≈ sol2_ip(1.0) rtol=1e-3 From 4af4ab8382c257bf2b802aafca4f15baa20cf08b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 30 Jun 2024 06:05:43 +0800 Subject: [PATCH 05/36] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 845c7009..672c48da 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StochasticDiffEq" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" authors = ["Chris Rackauckas "] -version = "6.65.1" +version = "6.66.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 8129a7193886521390417c44e41561308679eb0e Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 25 Jul 2024 12:15:52 -0400 Subject: [PATCH 06/36] don't try to use `calc_W!` on out of place SDE. How did this ever work correctly? --- test/utility_tests.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index f52561df..77985858 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -17,8 +17,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator jac=(u,p,t) -> A) prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) - W = integrator.cache.nlsolver.cache.W - calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false) + W = calc_W(integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false) @test convert(AbstractMatrix, W) ≈ concrete_W @test W \ u0 ≈ concrete_W \ u0 From 1443a30b0a8f822251f07917cae92d227bb24346 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 25 Jul 2024 12:34:42 -0400 Subject: [PATCH 07/36] fix direct push I accidentally pushed to main, and forgot to import `calc_W`, so CI is now failing. This PR should fix the CI. --- test/utility_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index 77985858..8bb4e3bc 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -1,5 +1,5 @@ using StochasticDiffEq, LinearAlgebra, SparseArrays, Random, LinearSolve, Test -using StochasticDiffEq.OrdinaryDiffEq: WOperator, calc_W! +using StochasticDiffEq.OrdinaryDiffEq: WOperator, calc_W!, calc_W using StochasticDiffEq.SciMLOperators: MatrixOperator @testset "Derivative Utilities" begin From 06ecdf2ae8f7ce8c9c811f17fafafabe21f0715f Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 25 Jul 2024 13:56:52 -0400 Subject: [PATCH 08/36] fix --- test/utility_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index 8bb4e3bc..51b32aa1 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -17,7 +17,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator jac=(u,p,t) -> A) prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) - W = calc_W(integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false) + W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, false) @test convert(AbstractMatrix, W) ≈ concrete_W @test W \ u0 ≈ concrete_W \ u0 From f4b8a313ecfec516d6f7aaf26a614bbf17ba8b9f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 26 Jul 2024 09:21:18 -0400 Subject: [PATCH 09/36] Update Project.toml --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 672c48da..ffca010d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StochasticDiffEq" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" authors = ["Chris Rackauckas "] -version = "6.66.0" +version = "6.67.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -45,7 +45,7 @@ LinearAlgebra = "1.6" Logging = "1.6" MuladdMacro = "0.2.1" NLsolve = "4" -OrdinaryDiffEq = "6.85" +OrdinaryDiffEq = "6.87" Random = "1.6" RandomNumbers = "1.5.3" RecursiveArrayTools = "2, 3" From 8c5826b954120b9a0480ba16260db3fdc1b2bcb8 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 2 Aug 2024 12:09:49 -0400 Subject: [PATCH 10/36] test W-transform the solvers already use the W-transform so we should be testing with it (especially because I'm about to delete the non-W-transformed version in https://github.com/SciML/OrdinaryDiffEq.jl/pull/2297) --- test/utility_tests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index 51b32aa1..acc8b78b 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -8,7 +8,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator mm = [2.0 0.0; 0.0 1.0] u0 = [1.0, 1.0]; tmp = zeros(2) tspan = (0.0,1.0); dt = 0.01; dtgamma = 0.5dt - concrete_W = -mm + dtgamma * A + concrete_W = -mm/dtgamma + A # Out-of-place _f = (u,p,t) -> A*u; _g = (u,p,t) -> σ*u @@ -17,7 +17,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator jac=(u,p,t) -> A) prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) - W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, false) + W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, true) @test convert(AbstractMatrix, W) ≈ concrete_W @test W \ u0 ≈ concrete_W \ u0 @@ -29,7 +29,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) W = integrator.cache.nlsolver.cache.W - calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false) + calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, true) # Did not update because it's an array operator # We don't want to build Jacobians when we have operators! @@ -37,7 +37,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp != concrete_W \ u0 # But jacobian2W! will update the cache - StochasticDiffEq.OrdinaryDiffEq.jacobian2W!(integrator.cache.nlsolver.cache.W._concrete_form, mm, dtgamma, integrator.cache.nlsolver.cache.W.J.A, false) + StochasticDiffEq.OrdinaryDiffEq.jacobian2W!(integrator.cache.nlsolver.cache.W._concrete_form, mm, dtgamma, integrator.cache.nlsolver.cache.W.J.A, true) @test convert(AbstractMatrix, integrator.cache.nlsolver.cache.W) == concrete_W ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp == concrete_W \ u0 end From 73c7840db212006344319bf2475094efe0c2f46b Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 2 Aug 2024 14:51:52 -0400 Subject: [PATCH 11/36] fix calc_W! test with Rosenbrock transform --- test/utility_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index acc8b78b..c16a05a4 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -29,7 +29,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) W = integrator.cache.nlsolver.cache.W - calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, true) + calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, #=repeat_step=#false, #=W_transform=#true) # Did not update because it's an array operator # We don't want to build Jacobians when we have operators! From ef4a16c5f4169e39c9dbb763503aad24fbfd6b33 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 2 Aug 2024 15:13:35 -0400 Subject: [PATCH 12/36] fix oop also --- test/utility_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index c16a05a4..532a936b 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -17,7 +17,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator jac=(u,p,t) -> A) prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) - W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, true) + W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, #=repeat_step=#false, #=W_transform=#true) @test convert(AbstractMatrix, W) ≈ concrete_W @test W \ u0 ≈ concrete_W \ u0 From 0098ddc79d6139cb0d8c835bba476a81a78f1f28 Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 8 Aug 2024 15:45:36 -0400 Subject: [PATCH 13/36] broadcast Base.:(/) for non-AbstractArray --- src/solve.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solve.jl b/src/solve.jl index d06354c8..98ac1031 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -283,7 +283,7 @@ function DiffEqBase.__init( uprev = recursivecopy(u) if !(uType <: AbstractArray) - rand_prototype = zero(u/u) # Strip units and type info + rand_prototype = zero(u ./ u) # Strip units and type info randType = typeof(rand_prototype) else randElType = uBottomEltypeNoUnits # Strip units and type info From c5750221fc7099cf947d104df45b97d2f49156de Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 8 Aug 2024 15:52:33 -0400 Subject: [PATCH 14/36] add test --- test/noindex_tests.jl | 88 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 87 insertions(+), 1 deletion(-) diff --git a/test/noindex_tests.jl b/test/noindex_tests.jl index b98512a2..31223ddb 100644 --- a/test/noindex_tests.jl +++ b/test/noindex_tests.jl @@ -1,4 +1,5 @@ -using StochasticDiffEq, Test, Random, DiffEqNoiseProcess +using StochasticDiffEq, Test, Random, DiffEqNoiseProcess, + RecursiveArrayTools, LinearAlgebra Random.seed!(100) struct NoIndexArray{T, N} <: AbstractArray{T, N} @@ -50,3 +51,88 @@ for alg in algs @test_nowarn sol(0.1) @test_nowarn sol(similar(prob.u0), 0.1) end + + +struct CustomArray{T, N} + x::Array{T, N} +end +Base.size(x::CustomArray) = size(x.x) +Base.axes(x::CustomArray) = axes(x.x) +Base.ndims(x::CustomArray) = ndims(x.x) +Base.ndims(::Type{<:CustomArray{T,N}}) where {T,N} = N +Base.zero(x::CustomArray) = CustomArray(zero(x.x)) +Base.zero(::Type{<:CustomArray{T,N}}) where {T,N} = CustomArray(zero(Array{T,N})) +Base.similar(x::CustomArray, dims::Union{Integer, AbstractUnitRange}...) = CustomArray(similar(x.x, dims...)) +Base.copyto!(x::CustomArray, y::CustomArray) = CustomArray(copyto!(x.x, y.x)) +Base.copy(x::CustomArray) = CustomArray(copy(x.x)) +Base.length(x::CustomArray) = length(x.x) +Base.isempty(x::CustomArray) = isempty(x.x) +Base.eltype(x::CustomArray) = eltype(x.x) +Base.zero(x::CustomArray) = CustomArray(zero(x.x)) +Base.fill!(x::CustomArray, y) = CustomArray(fill!(x.x, y)) +Base.getindex(x::CustomArray, i) = getindex(x.x, i) +Base.setindex!(x::CustomArray, v, idx) = setindex!(x.x, v, idx) +Base.mapreduce(f, op, x::CustomArray; kwargs...) = mapreduce(f, op, x.x; kwargs...) +Base.any(f::Function, x::CustomArray; kwargs...) = any(f, x.x; kwargs...) +Base.all(f::Function, x::CustomArray; kwargs...) = all(f, x.x; kwargs...) +Base.similar(x::CustomArray, t) = CustomArray(similar(x.x, t)) +Base.:(==)(x::CustomArray, y::CustomArray) = x.x == y.x +Base.:(*)(x::Number, y::CustomArray) = CustomArray(x*y.x) +Base.:(/)(x::CustomArray, y::Number) = CustomArray(x.x/y) +LinearAlgebra.norm(x::CustomArray) = norm(x.x) + +struct CustomStyle{N} <: Broadcast.BroadcastStyle where {N} end +CustomStyle(::Val{N}) where N = CustomStyle{N}() +CustomStyle{M}(::Val{N}) where {N,M} = NoIndexStyle{N}() +Base.BroadcastStyle(::Type{<:CustomArray{T,N}}) where {T,N} = CustomStyle{N}() +Broadcast.BroadcastStyle(::CustomStyle{N}, ::Broadcast.DefaultArrayStyle{0}) where {N} = CustomStyle{N}() +Base.similar(bc::Base.Broadcast.Broadcasted{CustomStyle{N}}, ::Type{ElType}) where {N, ElType} = CustomArray(similar(Array{ElType, N}, axes(bc))) +Base.Broadcast._broadcast_getindex(x::CustomArray, i) = x.x[i] +Base.Broadcast.extrude(x::CustomArray) = x +Base.Broadcast.broadcastable(x::CustomArray) = x + +@inline function Base.copyto!(dest::CustomArray, bc::Base.Broadcast.Broadcasted{<:CustomStyle}) + axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc)) + bc′ = Base.Broadcast.preprocess(dest, bc) + dest′ = dest.x + @simd for I in 1:length(dest′) + @inbounds dest′[I] = bc′[I] + end + return dest +end +@inline function Base.copy(bc::Base.Broadcast.Broadcasted{<:CustomStyle}) + bcf = Broadcast.flatten(bc) + x = find_x(bcf) + data = zeros(eltype(x), size(x)) + @inbounds @simd for I in 1:length(x) + data[I] = bcf[I] + end + return CustomArray(data) +end +find_x(bc::Broadcast.Broadcasted) = find_x(bc.args) +find_x(args::Tuple) = find_x(find_x(args[1]), Base.tail(args)) +find_x(x) = x +find_x(::Any, rest) = find_x(rest) +find_x(x::CustomArray, rest) = x.x + +RecursiveArrayTools.recursive_unitless_bottom_eltype(x::CustomArray) = eltype(x) +RecursiveArrayTools.recursivecopy!(dest::CustomArray, src::CustomArray) = copyto!(dest, src) +RecursiveArrayTools.recursivecopy(x::CustomArray) = copy(x) +RecursiveArrayTools.recursivefill!(x::CustomArray, a) = fill!(x, a) +DiffEqNoiseProcess.wiener_randn!(rng::AbstractRNG,rand_vec::CustomArray) = randn!(rng,rand_vec.x) + +Base.show_vector(io::IO, x::CustomArray) = Base.show_vector(io, x.x) + +Base.show(io::IO, x::CustomArray) = (print(io, "CustomArray");show(io, x.x)) +function Base.show(io::IO, ::MIME"text/plain", x::CustomArray) + println(io, Base.summary(x), ":") + Base.print_array(io, x.x) +end + +prob = SDEProblem((du, u, p, t)->copyto!(du, u),(du, u, p, t)->copyto!(du, u), CustomArray(ones(10)), (0.0,1.0)) + +for alg in algs + sol_ca = @test_nowarn solve(prob, alg) + @test_nowarn sol_ca(0.1) + @test_nowarn sol_ca(similar(prob.u0), 0.1) + end \ No newline at end of file From 4d72b60dcec9aacb85635612dfec96f4204b2d09 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 13 Aug 2024 08:36:00 -0400 Subject: [PATCH 15/36] add tstops api functions --- src/StochasticDiffEq.jl | 3 ++- src/integrators/integrator_interface.jl | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index d9cc79d3..af7b4f54 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -47,7 +47,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/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl index 6ff668f0..d2b41808 100644 --- 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)) \ No newline at end of file From fe083b7096eb8ea2a4cc770429fdf1727c5d0cd4 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 13 Aug 2024 08:37:10 -0400 Subject: [PATCH 16/36] eof --- src/integrators/integrator_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl index d2b41808..aed5d584 100644 --- a/src/integrators/integrator_interface.jl +++ b/src/integrators/integrator_interface.jl @@ -393,4 +393,4 @@ 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)) \ No newline at end of file +DiffEqBase.get_tstops_max(integ::SDEIntegrator) = maximum(get_tstops_array(integ)) From eae6f16e0a5b3954b891a6cb3ce6898cc4e40e32 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 18 Aug 2024 11:28:53 -0400 Subject: [PATCH 17/36] update DiffEqBase --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ffca010d..9ae11da1 100644 --- a/Project.toml +++ b/Project.toml @@ -34,7 +34,7 @@ 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" FiniteDiff = "2" From e02a78b70e124d0a8549d48103c33760125e0fa5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 18 Aug 2024 13:10:32 -0400 Subject: [PATCH 18/36] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9ae11da1..1c3c5541 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StochasticDiffEq" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" authors = ["Chris Rackauckas "] -version = "6.67.0" +version = "6.68.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 6eefb4e4af7e5174d0a26a18b0b48a8095f36d80 Mon Sep 17 00:00:00 2001 From: Hossein Pourbozorg Date: Fri, 30 Aug 2024 14:06:46 +0330 Subject: [PATCH 19/36] `alg_interpretation` from SciMLBase --- src/alg_utils.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 0d49bd48..589f938c 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -133,22 +133,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) = :Ito +SciMLBase.alg_interpretation(alg::EulerHeun) = :Stratonovich +SciMLBase.alg_interpretation(alg::LambaEulerHeun) = :Stratonovich +SciMLBase.alg_interpretation(alg::KomBurSROCK2) = :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) = :Stratonovich +SciMLBase.alg_interpretation(alg::RS2) = :Stratonovich + +SciMLBase.alg_interpretation(alg::NON) = :Stratonovich +SciMLBase.alg_interpretation(alg::COM) = :Stratonovich +SciMLBase.alg_interpretation(alg::NON2) = :Stratonovich alg_compatible(prob, alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true alg_compatible(prob, alg::StochasticDiffEqAlgorithm) = false From dbfbb77d4b084fbe7eed8b1e43961a1bed8f004e Mon Sep 17 00:00:00 2001 From: Hossein Pourbozorg Date: Fri, 30 Aug 2024 18:13:08 +0330 Subject: [PATCH 20/36] full qualified calls --- src/perform_step/SROCK_perform_step.jl | 12 ++++++------ src/perform_step/low_order.jl | 24 ++++++++++++------------ src/perform_step/sdirk.jl | 8 ++++---- 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/perform_step/SROCK_perform_step.jl b/src/perform_step/SROCK_perform_step.jl index 046cfc7b..5e9eba74 100644 --- 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) == :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) == :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) == :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) == :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) == :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) == :Ito if W.dW isa Number || is_diagonal_noise(integrator.sol.prob) integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁) @.. uᵢ₋₂ = uᵢ₋₁ + sqrt(abs(dt))*gₘ₋₂ diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index e5c85818..054ee5f5 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -208,11 +208,11 @@ end K = @.. uprev + dt * du1 L = integrator.g(uprev,p,t) mil_correction = zero(u) - if alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == :Ito utilde = K + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(W.dW.^2 .- abs(dt))./2 - elseif alg_interpretation(integrator.alg) == :Stratonovich + elseif SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich utilde = uprev + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(W.dW.^2)./2 @@ -241,11 +241,11 @@ end integrator.g(L,uprev,p,t) @.. K = uprev + dt * du1 @.. du2 = zero(eltype(u)) # This makes it safe to re-use the array - if alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == :Ito @.. tmp = K + integrator.sqdt * L integrator.g(du2,tmp,p,t) @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2 - abs(dt)) - elseif alg_interpretation(integrator.alg) == :Stratonovich + elseif SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich @.. tmp = uprev + integrator.sqdt * L integrator.g(du2,tmp,p,t) @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2) @@ -275,7 +275,7 @@ end J = get_iterated_I(dt, dW, W.dZ, Jalg) mil_correction = zero(u) - if alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == :Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) J = J .- 1//2 .* abs(dt) else @@ -289,7 +289,7 @@ end K = uprev + dt*du1 if is_diagonal_noise(integrator.sol.prob) - tmp = (alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + tmp = (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L gtmp = integrator.g(tmp,p,t) Dgj = (gtmp - L)/sqdt ggprime_norm = integrator.opts.internalnorm(Dgj,t) @@ -343,7 +343,7 @@ end J = Jalg.J @.. mil_correction = zero(u) - if alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == :Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. J -= 1 // 2 * abs(dt) else @@ -357,7 +357,7 @@ end @.. K = uprev + dt*du1 if is_diagonal_noise(integrator.sol.prob) - tmp .= (alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + tmp .= (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L integrator.g(gtmp,tmp,p,t) @.. Dgj = (gtmp - L)/sqdt ggprime_norm = integrator.opts.internalnorm(Dgj,t) @@ -397,7 +397,7 @@ end J = get_iterated_I(dt, dW, W.dZ, Jalg, integrator.alg.p, integrator.alg.c, alg_order(integrator.alg)) - if alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == :Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) J = J .- 1//2 .* abs(dt) else @@ -413,7 +413,7 @@ end if dW isa Number || is_diagonal_noise(integrator.sol.prob) K = @.. uprev + dt*du₁ - utilde = (alg_interpretation(integrator.alg) == :Ito ? K : uprev) + L*integrator.sqdt + utilde = (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t) .- L) ./ (integrator.sqdt) mil_correction = ggprime .* J u = K + L .* dW + mil_correction @@ -467,7 +467,7 @@ end @.. mil_correction = zero(eltype(u)) ggprime_norm = zero(eltype(ggprime)) - if alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == :Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. J -= 1 // 2 * abs(dt) else @@ -478,7 +478,7 @@ end if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. K = uprev + dt*du₁ @.. du₂ = zero(eltype(u)) - tmp .= (alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + tmp .= (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L integrator.g(du₂,tmp,p,t) @.. ggprime = (du₂ - L)/sqdt ggprime_norm = integrator.opts.internalnorm(ggprime,t) diff --git a/src/perform_step/sdirk.jl b/src/perform_step/sdirk.jl index 00f346b2..7a325568 100644 --- a/src/perform_step/sdirk.jl +++ b/src/perform_step/sdirk.jl @@ -24,14 +24,14 @@ end if cache isa ImplicitRKMilConstantCache || integrator.opts.adaptive == true - if alg_interpretation(alg) == :Ito || + if SciMLBase.alg_interpretation(alg) == :Ito || cache isa ImplicitEMConstantCache K = @.. uprev + dt * ftmp utilde = K + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime .* (integrator.W.dW.^2 .- abs(dt))./2 gtmp += mil_correction - elseif alg_interpretation(alg) == :Stratonovich || + elseif SciMLBase.alg_interpretation(alg) == :Stratonovich || cache isa ImplicitEulerHeunConstantCache utilde = uprev + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) @@ -154,12 +154,12 @@ end if cache isa ImplicitRKMilCache gtmp3 = cache.gtmp3 - if alg_interpretation(alg) == :Ito + if SciMLBase.alg_interpretation(alg) == :Ito @.. z = uprev + dt * tmp + integrator.sqdt * gtmp integrator.g(gtmp3,z,p,t) @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation @.. gtmp2 += gtmp3*(dW.^2 - abs(dt))/2 - elseif alg_interpretation(alg) == :Stratonovich + elseif SciMLBase.alg_interpretation(alg) == :Stratonovich @.. z = uprev + integrator.sqdt * gtmp integrator.g(gtmp3,z,p,t) @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation From ac6c27fd1db0130eb15c0194635fc4412515db22 Mon Sep 17 00:00:00 2001 From: Hossein Pourbozorg Date: Fri, 30 Aug 2024 20:26:02 +0330 Subject: [PATCH 21/36] use `EnumX` --- Project.toml | 2 ++ src/StochasticDiffEq.jl | 1 + src/alg_utils.jl | 20 +++++++------ src/algorithms.jl | 26 ++++++++--------- src/perform_step/SROCK_perform_step.jl | 12 ++++---- src/perform_step/low_order.jl | 28 +++++++++---------- src/perform_step/sdirk.jl | 12 ++++---- .../sde_complex_adaptive_mean_test.jl | 12 ++++---- test/reversal_tests.jl | 10 +++---- test/stratonovich_convergence_tests.jl | 20 ++++++------- test/zerod_noise_test.jl | 6 ++-- 11 files changed, 77 insertions(+), 72 deletions(-) diff --git a/Project.toml b/Project.toml index 1c3c5541..4162b434 100644 --- a/Project.toml +++ b/Project.toml @@ -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" +EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" @@ -37,6 +38,7 @@ DataStructures = "0.18" DiffEqBase = "6.154" DiffEqNoiseProcess = "5.13" DocStringExtensions = "0.8, 0.9" +EnumX = "1" FiniteDiff = "2" ForwardDiff = "0.10.3" JumpProcesses = "9" diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index af7b4f54..f7a1350a 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -73,6 +73,7 @@ else end import SciMLBase + import EnumX using SparseDiffTools: forwarddiff_color_jacobian!, ForwardColorJacCache diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 589f938c..c5bc8199 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -133,22 +133,24 @@ beta1_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm isdtchangeable(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true -SciMLBase.alg_interpretation(alg::StochasticDiffEqAlgorithm) = :Ito -SciMLBase.alg_interpretation(alg::EulerHeun) = :Stratonovich -SciMLBase.alg_interpretation(alg::LambaEulerHeun) = :Stratonovich -SciMLBase.alg_interpretation(alg::KomBurSROCK2) = :Stratonovich +EnumX.@enumx AlgorithmInterpretation Ito Stratonovich + +SciMLBase.alg_interpretation(alg::StochasticDiffEqAlgorithm) = AlgorithmInterpretation.Ito +SciMLBase.alg_interpretation(alg::EulerHeun) = AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::LambaEulerHeun) = AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::KomBurSROCK2) = 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) = :Stratonovich -SciMLBase.alg_interpretation(alg::RS2) = :Stratonovich +SciMLBase.alg_interpretation(alg::RS1) = AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::RS2) = AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::NON) = :Stratonovich -SciMLBase.alg_interpretation(alg::COM) = :Stratonovich -SciMLBase.alg_interpretation(alg::NON2) = :Stratonovich +SciMLBase.alg_interpretation(alg::NON) = AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::COM) = AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::NON2) = 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 index b3880059..ceaca812 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -77,11 +77,11 @@ 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=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=AlgorithmInterpretation.Ito) = RKMil{interpretation}() """ Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations. @@ -89,37 +89,37 @@ 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=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 + interpretation::AlgorithmInterpretation.T ii_approx::T end -RKMilCommute(;interpretation=:Ito, ii_approx=IICommutative()) = RKMilCommute(interpretation,ii_approx) +RKMilCommute(;interpretation=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=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 AlgorithmInterpretation.Ito and 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::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=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=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=AlgorithmInterpretation.Ito,eigen_est=nothing) = SROCK1{interpretation,typeof(eigen_est)}(eigen_est) # Weak Order 2 for Alg in [:SROCK2, :KomBurSROCK2, :SROCKC2] @@ -701,7 +701,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=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. """ @@ -721,7 +721,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=AlgorithmInterpretation.Ito) = ImplicitRKMil{chunk_size,autodiff, typeof(linsolve),typeof(precs),diff_type, OrdinaryDiffEq._unwrap_val(standardtag), diff --git a/src/perform_step/SROCK_perform_step.jl b/src/perform_step/SROCK_perform_step.jl index 5e9eba74..285a90a3 100644 --- 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 SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich + if SciMLBase.alg_interpretation(integrator.alg) == 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) && SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich + if (i > mdeg - 2) && SciMLBase.alg_interpretation(integrator.alg) == 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) && SciMLBase.alg_interpretation(integrator.alg) == :Ito + elseif (i == mdeg) && SciMLBase.alg_interpretation(integrator.alg) == 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 SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich + if SciMLBase.alg_interpretation(integrator.alg) == 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) && SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich + if (i > mdeg - 2) && SciMLBase.alg_interpretation(integrator.alg) == 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) && SciMLBase.alg_interpretation(integrator.alg) == :Ito + elseif (i == mdeg) && SciMLBase.alg_interpretation(integrator.alg) == 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ₘ₋₂ diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index 054ee5f5..297cb7e4 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -208,16 +208,16 @@ end K = @.. uprev + dt * du1 L = integrator.g(uprev,p,t) mil_correction = zero(u) - if SciMLBase.alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito utilde = K + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(W.dW.^2 .- abs(dt))./2 - elseif SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich + elseif SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Stratonovich utilde = uprev + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(W.dW.^2)./2 else - error("Alg interpretation invalid. Use either :Ito or :Stratonovich") + error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") end u = K+L.*W.dW+mil_correction @@ -241,16 +241,16 @@ end integrator.g(L,uprev,p,t) @.. K = uprev + dt * du1 @.. du2 = zero(eltype(u)) # This makes it safe to re-use the array - if SciMLBase.alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito @.. tmp = K + integrator.sqdt * L integrator.g(du2,tmp,p,t) @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2 - abs(dt)) - elseif SciMLBase.alg_interpretation(integrator.alg) == :Stratonovich + elseif SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Stratonovich @.. tmp = uprev + integrator.sqdt * L integrator.g(du2,tmp,p,t) @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2) else - error("Alg interpretation invalid. Use either :Ito or :Stratonovich") + error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") end @.. u = K+L*W.dW + tmp if integrator.opts.adaptive @@ -275,7 +275,7 @@ end J = get_iterated_I(dt, dW, W.dZ, Jalg) mil_correction = zero(u) - if SciMLBase.alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) J = J .- 1//2 .* abs(dt) else @@ -289,7 +289,7 @@ end K = uprev + dt*du1 if is_diagonal_noise(integrator.sol.prob) - tmp = (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + tmp = (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L gtmp = integrator.g(tmp,p,t) Dgj = (gtmp - L)/sqdt ggprime_norm = integrator.opts.internalnorm(Dgj,t) @@ -343,7 +343,7 @@ end J = Jalg.J @.. mil_correction = zero(u) - if SciMLBase.alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. J -= 1 // 2 * abs(dt) else @@ -357,7 +357,7 @@ end @.. K = uprev + dt*du1 if is_diagonal_noise(integrator.sol.prob) - tmp .= (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + tmp .= (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L integrator.g(gtmp,tmp,p,t) @.. Dgj = (gtmp - L)/sqdt ggprime_norm = integrator.opts.internalnorm(Dgj,t) @@ -397,7 +397,7 @@ end J = get_iterated_I(dt, dW, W.dZ, Jalg, integrator.alg.p, integrator.alg.c, alg_order(integrator.alg)) - if SciMLBase.alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) J = J .- 1//2 .* abs(dt) else @@ -413,7 +413,7 @@ end if dW isa Number || is_diagonal_noise(integrator.sol.prob) K = @.. uprev + dt*du₁ - utilde = (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) + L*integrator.sqdt + utilde = (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t) .- L) ./ (integrator.sqdt) mil_correction = ggprime .* J u = K + L .* dW + mil_correction @@ -467,7 +467,7 @@ end @.. mil_correction = zero(eltype(u)) ggprime_norm = zero(eltype(ggprime)) - if SciMLBase.alg_interpretation(integrator.alg) == :Ito + if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. J -= 1 // 2 * abs(dt) else @@ -478,7 +478,7 @@ end if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. K = uprev + dt*du₁ @.. du₂ = zero(eltype(u)) - tmp .= (SciMLBase.alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + tmp .= (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L integrator.g(du₂,tmp,p,t) @.. ggprime = (du₂ - L)/sqdt ggprime_norm = integrator.opts.internalnorm(ggprime,t) diff --git a/src/perform_step/sdirk.jl b/src/perform_step/sdirk.jl index 7a325568..7e39cc6a 100644 --- a/src/perform_step/sdirk.jl +++ b/src/perform_step/sdirk.jl @@ -24,21 +24,21 @@ end if cache isa ImplicitRKMilConstantCache || integrator.opts.adaptive == true - if SciMLBase.alg_interpretation(alg) == :Ito || + if SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Ito || cache isa ImplicitEMConstantCache K = @.. uprev + dt * ftmp utilde = K + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime .* (integrator.W.dW.^2 .- abs(dt))./2 gtmp += mil_correction - elseif SciMLBase.alg_interpretation(alg) == :Stratonovich || + elseif SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Stratonovich || cache isa ImplicitEulerHeunConstantCache utilde = uprev + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(integrator.W.dW.^2)./2 gtmp += mil_correction else - error("Alg interpretation invalid. Use either :Ito or :Stratonovich") + error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") end end @@ -154,18 +154,18 @@ end if cache isa ImplicitRKMilCache gtmp3 = cache.gtmp3 - if SciMLBase.alg_interpretation(alg) == :Ito + if SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Ito @.. z = uprev + dt * tmp + integrator.sqdt * gtmp integrator.g(gtmp3,z,p,t) @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation @.. gtmp2 += gtmp3*(dW.^2 - abs(dt))/2 - elseif SciMLBase.alg_interpretation(alg) == :Stratonovich + elseif SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Stratonovich @.. z = uprev + integrator.sqdt * gtmp integrator.g(gtmp3,z,p,t) @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation @.. gtmp2 += gtmp3*(dW.^2)/2 else - error("Alg interpretation invalid. Use either :Ito or :Stratonovich") + error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") end end diff --git a/test/adaptive/sde_complex_adaptive_mean_test.jl b/test/adaptive/sde_complex_adaptive_mean_test.jl index 542fb669..5edfd51f 100644 --- a/test/adaptive/sde_complex_adaptive_mean_test.jl +++ b/test/adaptive/sde_complex_adaptive_mean_test.jl @@ -44,20 +44,20 @@ let out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMil{:Stratonovich}(), dt=1e-4, callback=scb1, seed = i, adaptive = false) + solve(prob1, RKMil{AlgorithmInterpretation.Stratonovich}(), dt=1e-4, callback=scb1, seed = i, adaptive = false) avg1 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMil{:Stratonovich}(), tstops = T, callback=scb1, + solve(prob1, RKMil{AlgorithmInterpretation.Stratonovich}(), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg2 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob2, RKMil{:Stratonovich}(), tstops = T, callback=scb1, + solve(prob2, RKMil{AlgorithmInterpretation.Stratonovich}(), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg3 .+= out1.saveval ./ Ntraj end @@ -88,20 +88,20 @@ let out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMilGeneral(interpretation = :Stratonovich), dt=1e-4, callback=scb1, seed = i, adaptive = false) + solve(prob1, RKMilGeneral(interpretation = AlgorithmInterpretation.Stratonovich), dt=1e-4, callback=scb1, seed = i, adaptive = false) avg1 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMilGeneral(interpretation = :Stratonovich), tstops = T, callback=scb1, + solve(prob1, RKMilGeneral(interpretation = AlgorithmInterpretation.Stratonovich), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg2 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob2, RKMilGeneral(interpretation = :Stratonovich), tstops = T, callback=scb1, + solve(prob2, RKMilGeneral(interpretation = AlgorithmInterpretation.Stratonovich), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg3 .+= out1.saveval ./ Ntraj end diff --git a/test/reversal_tests.jl b/test/reversal_tests.jl index 06116faf..9268101b 100644 --- a/test/reversal_tests.jl +++ b/test/reversal_tests.jl @@ -19,14 +19,14 @@ Stratonovich_solver = [ # non-stiff methods EulerHeun(), LambaEulerHeun(), - RKMil(interpretation=:Stratonovich), - RKMilCommute(interpretation=:Stratonovich), - RKMilGeneral(interpretation=:Stratonovich), + RKMil(interpretation=AlgorithmInterpretation.Stratonovich), + RKMilCommute(interpretation=AlgorithmInterpretation.Stratonovich), + RKMilGeneral(interpretation=AlgorithmInterpretation.Stratonovich), # S-Rock methods - SROCK1(interpretation=:Stratonovich), + SROCK1(interpretation=AlgorithmInterpretation.Stratonovich), # stiff methods ImplicitEulerHeun(), - ImplicitRKMil(interpretation=:Stratonovich), + ImplicitRKMil(interpretation=AlgorithmInterpretation.Stratonovich), ISSEulerHeun(), ] diff --git a/test/stratonovich_convergence_tests.jl b/test/stratonovich_convergence_tests.jl index 4b8c14e3..37cfb1c8 100644 --- a/test/stratonovich_convergence_tests.jl +++ b/test/stratonovich_convergence_tests.jl @@ -22,16 +22,16 @@ sim = test_convergence(dts,prob,ImplicitEulerHeun(theta=1),trajectories=Int(5e2 sim = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),trajectories=Int(5e2)) +sim = test_convergence(dts,prob,RKMil(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,RKMilGeneral(interpretation=:Stratonovich),trajectories=Int(5e2)) +sim = test_convergence(dts,prob,RKMilGeneral(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=:Stratonovich),trajectories=Int(5e2)) +sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,SROCK1(interpretation=:Stratonovich),trajectories=Int(2e2)) +sim = test_convergence(dts,prob,SROCK1(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(2e2)) @test abs(sim.𝒪est[:l2]-1) < 0.15 sim = test_convergence(dts,prob,KomBurSROCK2(),trajectories=Int(2e2)) @@ -61,25 +61,25 @@ sim = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),trajectories println("RKMils") -sim = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),trajectories=Int(1e2)) +sim = test_convergence(dts,prob,RKMil(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,RKMilGeneral(interpretation=:Stratonovich),trajectories=Int(1e2)) +sim = test_convergence(dts,prob,RKMilGeneral(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=:Stratonovich), +sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=AlgorithmInterpretation.Stratonovich), trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,ImplicitRKMil(theta=1,interpretation=:Stratonovich), +sim = test_convergence(dts,prob,ImplicitRKMil(theta=1,interpretation=AlgorithmInterpretation.Stratonovich), trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,interpretation=:Stratonovich), +sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,interpretation=AlgorithmInterpretation.Stratonovich), trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,SROCK1(interpretation=:Stratonovich),trajectories=Int(1e2)) +sim = test_convergence(dts,prob,SROCK1(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 sim = test_convergence(dts,prob,KomBurSROCK2(),trajectories=Int(2e2)) diff --git a/test/zerod_noise_test.jl b/test/zerod_noise_test.jl index 9a51d7f9..20ac4d17 100644 --- a/test/zerod_noise_test.jl +++ b/test/zerod_noise_test.jl @@ -10,16 +10,16 @@ end u0 = [1.0] prob = SDEProblem{true}(f, g, u0, (0.0, 0.1)) -sol_ito = solve(prob, RKMil{:Ito}()) +sol_ito = solve(prob, RKMil{AlgorithmInterpretation.Ito}()) @test length(sol_ito) < 100 -sol_strato = solve(prob, RKMil{:Stratonovich}(); dt=1e-2) +sol_strato = solve(prob, RKMil{AlgorithmInterpretation.Stratonovich}(); dt=1e-2) @test length(sol_strato) < 100 sol_ito = solve(prob, RKMil()) @test length(sol_ito) < 100 -sol_strato = solve(prob, RKMil(interpretation=:Stratonovich); dt=1e-2) +sol_strato = solve(prob, RKMil(interpretation=AlgorithmInterpretation.Stratonovich); dt=1e-2) @test length(sol_strato) < 100 sol_leh = solve(prob, LambaEulerHeun()) From ef189b824ec7d76cdcb55fe950d8dcc15e301c20 Mon Sep 17 00:00:00 2001 From: Hossein Pourbozorg Date: Fri, 30 Aug 2024 20:37:31 +0330 Subject: [PATCH 22/36] move temporary --- src/StochasticDiffEq.jl | 2 ++ src/alg_utils.jl | 2 -- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index f7a1350a..7ca7f147 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -91,6 +91,8 @@ end seed_multiplier() = 1 end + EnumX.@enumx AlgorithmInterpretation Ito Stratonovich + include("misc_utils.jl") include("algorithms.jl") include("options_type.jl") diff --git a/src/alg_utils.jl b/src/alg_utils.jl index c5bc8199..e3bbde00 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -133,8 +133,6 @@ beta1_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm isdtchangeable(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true -EnumX.@enumx AlgorithmInterpretation Ito Stratonovich - SciMLBase.alg_interpretation(alg::StochasticDiffEqAlgorithm) = AlgorithmInterpretation.Ito SciMLBase.alg_interpretation(alg::EulerHeun) = AlgorithmInterpretation.Stratonovich SciMLBase.alg_interpretation(alg::LambaEulerHeun) = AlgorithmInterpretation.Stratonovich From c1f29d66bf00049fb8f8e49cb0a9c6409edccf3a Mon Sep 17 00:00:00 2001 From: Hossein Pourbozorg Date: Sat, 31 Aug 2024 15:51:09 +0330 Subject: [PATCH 23/36] move the type to SciMLBase --- src/StochasticDiffEq.jl | 2 -- src/alg_utils.jl | 18 ++++++------ src/algorithms.jl | 26 ++++++++--------- src/perform_step/SROCK_perform_step.jl | 12 ++++---- src/perform_step/low_order.jl | 28 +++++++++---------- src/perform_step/sdirk.jl | 12 ++++---- .../sde_complex_adaptive_mean_test.jl | 12 ++++---- test/reversal_tests.jl | 10 +++---- test/stratonovich_convergence_tests.jl | 20 ++++++------- test/zerod_noise_test.jl | 6 ++-- 10 files changed, 72 insertions(+), 74 deletions(-) diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index 7ca7f147..f7a1350a 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -91,8 +91,6 @@ end seed_multiplier() = 1 end - EnumX.@enumx AlgorithmInterpretation Ito Stratonovich - include("misc_utils.jl") include("algorithms.jl") include("options_type.jl") diff --git a/src/alg_utils.jl b/src/alg_utils.jl index e3bbde00..b7b15256 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -133,22 +133,22 @@ beta1_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm isdtchangeable(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true -SciMLBase.alg_interpretation(alg::StochasticDiffEqAlgorithm) = AlgorithmInterpretation.Ito -SciMLBase.alg_interpretation(alg::EulerHeun) = AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::LambaEulerHeun) = AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::KomBurSROCK2) = AlgorithmInterpretation.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) = AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::RS2) = AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::RS1) = SciMLBase.AlgorithmInterpretation.Stratonovich +SciMLBase.alg_interpretation(alg::RS2) = SciMLBase.AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::NON) = AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::COM) = AlgorithmInterpretation.Stratonovich -SciMLBase.alg_interpretation(alg::NON2) = 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 index ceaca812..c6f12490 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -77,11 +77,11 @@ 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=AlgorithmInterpretation.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=AlgorithmInterpretation.Ito) = RKMil{interpretation}() +RKMil(;interpretation=SciMLBase.AlgorithmInterpretation.Ito) = RKMil{interpretation}() """ Kloeden, P.E., Platen, E., Numerical Solution of Stochastic Differential Equations. @@ -89,37 +89,37 @@ 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=AlgorithmInterpretation.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::AlgorithmInterpretation.T + interpretation::SciMLBase.AlgorithmInterpretation.T ii_approx::T end -RKMilCommute(;interpretation=AlgorithmInterpretation.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=AlgorithmInterpretation.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 AlgorithmInterpretation.Ito and AlgorithmInterpretation.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::AlgorithmInterpretation.T + interpretation::SciMLBase.AlgorithmInterpretation.T ii_approx::T c::Int p::TruncationType end -function RKMilGeneral(;interpretation=AlgorithmInterpretation.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=AlgorithmInterpretation.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=AlgorithmInterpretation.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] @@ -701,7 +701,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=AlgorithmInterpretation.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. """ @@ -721,7 +721,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=AlgorithmInterpretation.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/perform_step/SROCK_perform_step.jl b/src/perform_step/SROCK_perform_step.jl index 285a90a3..9825d8d5 100644 --- 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 SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.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) && SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.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) && SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.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 SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.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) && SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.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) && SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.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ₘ₋₂ diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index 297cb7e4..72b06620 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -208,16 +208,16 @@ end K = @.. uprev + dt * du1 L = integrator.g(uprev,p,t) mil_correction = zero(u) - if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito utilde = K + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(W.dW.^2 .- abs(dt))./2 - elseif SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Stratonovich + elseif SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Stratonovich utilde = uprev + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(W.dW.^2)./2 else - error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") + error("Algorithm interpretation invalid. Use either SciMLBase.AlgorithmInterpretation.Ito or SciMLBase.AlgorithmInterpretation.Stratonovich") end u = K+L.*W.dW+mil_correction @@ -241,16 +241,16 @@ end integrator.g(L,uprev,p,t) @.. K = uprev + dt * du1 @.. du2 = zero(eltype(u)) # This makes it safe to re-use the array - if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito @.. tmp = K + integrator.sqdt * L integrator.g(du2,tmp,p,t) @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2 - abs(dt)) - elseif SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Stratonovich + elseif SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Stratonovich @.. tmp = uprev + integrator.sqdt * L integrator.g(du2,tmp,p,t) @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2) else - error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") + error("Algorithm interpretation invalid. Use either SciMLBase.AlgorithmInterpretation.Ito or SciMLBase.AlgorithmInterpretation.Stratonovich") end @.. u = K+L*W.dW + tmp if integrator.opts.adaptive @@ -275,7 +275,7 @@ end J = get_iterated_I(dt, dW, W.dZ, Jalg) mil_correction = zero(u) - if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) J = J .- 1//2 .* abs(dt) else @@ -289,7 +289,7 @@ end K = uprev + dt*du1 if is_diagonal_noise(integrator.sol.prob) - tmp = (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L + tmp = (SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L gtmp = integrator.g(tmp,p,t) Dgj = (gtmp - L)/sqdt ggprime_norm = integrator.opts.internalnorm(Dgj,t) @@ -343,7 +343,7 @@ end J = Jalg.J @.. mil_correction = zero(u) - if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. J -= 1 // 2 * abs(dt) else @@ -357,7 +357,7 @@ end @.. K = uprev + dt*du1 if is_diagonal_noise(integrator.sol.prob) - tmp .= (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L + tmp .= (SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L integrator.g(gtmp,tmp,p,t) @.. Dgj = (gtmp - L)/sqdt ggprime_norm = integrator.opts.internalnorm(Dgj,t) @@ -397,7 +397,7 @@ end J = get_iterated_I(dt, dW, W.dZ, Jalg, integrator.alg.p, integrator.alg.c, alg_order(integrator.alg)) - if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) J = J .- 1//2 .* abs(dt) else @@ -413,7 +413,7 @@ end if dW isa Number || is_diagonal_noise(integrator.sol.prob) K = @.. uprev + dt*du₁ - utilde = (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) + L*integrator.sqdt + utilde = (SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito ? K : uprev) + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t) .- L) ./ (integrator.sqdt) mil_correction = ggprime .* J u = K + L .* dW + mil_correction @@ -467,7 +467,7 @@ end @.. mil_correction = zero(eltype(u)) ggprime_norm = zero(eltype(ggprime)) - if SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. J -= 1 // 2 * abs(dt) else @@ -478,7 +478,7 @@ end if dW isa Number || is_diagonal_noise(integrator.sol.prob) @.. K = uprev + dt*du₁ @.. du₂ = zero(eltype(u)) - tmp .= (SciMLBase.alg_interpretation(integrator.alg) == AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L + tmp .= (SciMLBase.alg_interpretation(integrator.alg) == SciMLBase.AlgorithmInterpretation.Ito ? K : uprev) .+ integrator.sqdt .* L integrator.g(du₂,tmp,p,t) @.. ggprime = (du₂ - L)/sqdt ggprime_norm = integrator.opts.internalnorm(ggprime,t) diff --git a/src/perform_step/sdirk.jl b/src/perform_step/sdirk.jl index 7e39cc6a..d475a009 100644 --- a/src/perform_step/sdirk.jl +++ b/src/perform_step/sdirk.jl @@ -24,21 +24,21 @@ end if cache isa ImplicitRKMilConstantCache || integrator.opts.adaptive == true - if SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Ito || + if SciMLBase.alg_interpretation(alg) == SciMLBase.AlgorithmInterpretation.Ito || cache isa ImplicitEMConstantCache K = @.. uprev + dt * ftmp utilde = K + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime .* (integrator.W.dW.^2 .- abs(dt))./2 gtmp += mil_correction - elseif SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Stratonovich || + elseif SciMLBase.alg_interpretation(alg) == SciMLBase.AlgorithmInterpretation.Stratonovich || cache isa ImplicitEulerHeunConstantCache utilde = uprev + L*integrator.sqdt ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt) mil_correction = ggprime.*(integrator.W.dW.^2)./2 gtmp += mil_correction else - error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") + error("Algorithm interpretation invalid. Use either SciMLBase.AlgorithmInterpretation.Ito or SciMLBase.AlgorithmInterpretation.Stratonovich") end end @@ -154,18 +154,18 @@ end if cache isa ImplicitRKMilCache gtmp3 = cache.gtmp3 - if SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Ito + if SciMLBase.alg_interpretation(alg) == SciMLBase.AlgorithmInterpretation.Ito @.. z = uprev + dt * tmp + integrator.sqdt * gtmp integrator.g(gtmp3,z,p,t) @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation @.. gtmp2 += gtmp3*(dW.^2 - abs(dt))/2 - elseif SciMLBase.alg_interpretation(alg) == AlgorithmInterpretation.Stratonovich + elseif SciMLBase.alg_interpretation(alg) == SciMLBase.AlgorithmInterpretation.Stratonovich @.. z = uprev + integrator.sqdt * gtmp integrator.g(gtmp3,z,p,t) @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation @.. gtmp2 += gtmp3*(dW.^2)/2 else - error("Algorithm interpretation invalid. Use either AlgorithmInterpretation.Ito or AlgorithmInterpretation.Stratonovich") + error("Algorithm interpretation invalid. Use either SciMLBase.AlgorithmInterpretation.Ito or SciMLBase.AlgorithmInterpretation.Stratonovich") end end diff --git a/test/adaptive/sde_complex_adaptive_mean_test.jl b/test/adaptive/sde_complex_adaptive_mean_test.jl index 5edfd51f..4b1b5650 100644 --- a/test/adaptive/sde_complex_adaptive_mean_test.jl +++ b/test/adaptive/sde_complex_adaptive_mean_test.jl @@ -44,20 +44,20 @@ let out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMil{AlgorithmInterpretation.Stratonovich}(), dt=1e-4, callback=scb1, seed = i, adaptive = false) + solve(prob1, RKMil{SciMLBase.AlgorithmInterpretation.Stratonovich}(), dt=1e-4, callback=scb1, seed = i, adaptive = false) avg1 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMil{AlgorithmInterpretation.Stratonovich}(), tstops = T, callback=scb1, + solve(prob1, RKMil{SciMLBase.AlgorithmInterpretation.Stratonovich}(), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg2 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob2, RKMil{AlgorithmInterpretation.Stratonovich}(), tstops = T, callback=scb1, + solve(prob2, RKMil{SciMLBase.AlgorithmInterpretation.Stratonovich}(), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg3 .+= out1.saveval ./ Ntraj end @@ -88,20 +88,20 @@ let out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMilGeneral(interpretation = AlgorithmInterpretation.Stratonovich), dt=1e-4, callback=scb1, seed = i, adaptive = false) + solve(prob1, RKMilGeneral(interpretation = SciMLBase.AlgorithmInterpretation.Stratonovich), dt=1e-4, callback=scb1, seed = i, adaptive = false) avg1 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob1, RKMilGeneral(interpretation = AlgorithmInterpretation.Stratonovich), tstops = T, callback=scb1, + solve(prob1, RKMilGeneral(interpretation = SciMLBase.AlgorithmInterpretation.Stratonovich), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg2 .+= out1.saveval ./ Ntraj out1 = SavedValues(Float64,ComplexF64) scb1 = SavingCallback(fout, out1, saveat=T, save_everystep=false, save_start=false) - solve(prob2, RKMilGeneral(interpretation = AlgorithmInterpretation.Stratonovich), tstops = T, callback=scb1, + solve(prob2, RKMilGeneral(interpretation = SciMLBase.AlgorithmInterpretation.Stratonovich), tstops = T, callback=scb1, save_everystep=false, save_start=false) avg3 .+= out1.saveval ./ Ntraj end diff --git a/test/reversal_tests.jl b/test/reversal_tests.jl index 9268101b..778c7b56 100644 --- a/test/reversal_tests.jl +++ b/test/reversal_tests.jl @@ -19,14 +19,14 @@ Stratonovich_solver = [ # non-stiff methods EulerHeun(), LambaEulerHeun(), - RKMil(interpretation=AlgorithmInterpretation.Stratonovich), - RKMilCommute(interpretation=AlgorithmInterpretation.Stratonovich), - RKMilGeneral(interpretation=AlgorithmInterpretation.Stratonovich), + RKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), + RKMilCommute(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), + RKMilGeneral(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), # S-Rock methods - SROCK1(interpretation=AlgorithmInterpretation.Stratonovich), + SROCK1(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), # stiff methods ImplicitEulerHeun(), - ImplicitRKMil(interpretation=AlgorithmInterpretation.Stratonovich), + ImplicitRKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), ISSEulerHeun(), ] diff --git a/test/stratonovich_convergence_tests.jl b/test/stratonovich_convergence_tests.jl index 37cfb1c8..1ac25ed0 100644 --- a/test/stratonovich_convergence_tests.jl +++ b/test/stratonovich_convergence_tests.jl @@ -22,16 +22,16 @@ sim = test_convergence(dts,prob,ImplicitEulerHeun(theta=1),trajectories=Int(5e2 sim = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,RKMil(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) +sim = test_convergence(dts,prob,RKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,RKMilGeneral(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) +sim = test_convergence(dts,prob,RKMilGeneral(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) +sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(5e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,SROCK1(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(2e2)) +sim = test_convergence(dts,prob,SROCK1(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(2e2)) @test abs(sim.𝒪est[:l2]-1) < 0.15 sim = test_convergence(dts,prob,KomBurSROCK2(),trajectories=Int(2e2)) @@ -61,25 +61,25 @@ sim = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),trajectories println("RKMils") -sim = test_convergence(dts,prob,RKMil(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) +sim = test_convergence(dts,prob,RKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,RKMilGeneral(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) +sim = test_convergence(dts,prob,RKMilGeneral(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 -sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=AlgorithmInterpretation.Stratonovich), +sim = test_convergence(dts,prob,ImplicitRKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,ImplicitRKMil(theta=1,interpretation=AlgorithmInterpretation.Stratonovich), +sim = test_convergence(dts,prob,ImplicitRKMil(theta=1,interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,interpretation=AlgorithmInterpretation.Stratonovich), +sim = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich), trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 -sim = test_convergence(dts,prob,SROCK1(interpretation=AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) +sim = test_convergence(dts,prob,SROCK1(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.1 sim = test_convergence(dts,prob,KomBurSROCK2(),trajectories=Int(2e2)) diff --git a/test/zerod_noise_test.jl b/test/zerod_noise_test.jl index 20ac4d17..80e104ec 100644 --- a/test/zerod_noise_test.jl +++ b/test/zerod_noise_test.jl @@ -10,16 +10,16 @@ end u0 = [1.0] prob = SDEProblem{true}(f, g, u0, (0.0, 0.1)) -sol_ito = solve(prob, RKMil{AlgorithmInterpretation.Ito}()) +sol_ito = solve(prob, RKMil{SciMLBase.AlgorithmInterpretation.Ito}()) @test length(sol_ito) < 100 -sol_strato = solve(prob, RKMil{AlgorithmInterpretation.Stratonovich}(); dt=1e-2) +sol_strato = solve(prob, RKMil{SciMLBase.AlgorithmInterpretation.Stratonovich}(); dt=1e-2) @test length(sol_strato) < 100 sol_ito = solve(prob, RKMil()) @test length(sol_ito) < 100 -sol_strato = solve(prob, RKMil(interpretation=AlgorithmInterpretation.Stratonovich); dt=1e-2) +sol_strato = solve(prob, RKMil(interpretation=SciMLBase.AlgorithmInterpretation.Stratonovich); dt=1e-2) @test length(sol_strato) < 100 sol_leh = solve(prob, LambaEulerHeun()) From bcd441f0df443598e8d287992968295271c4f0b5 Mon Sep 17 00:00:00 2001 From: Hossein Pourbozorg Date: Sat, 31 Aug 2024 15:56:14 +0330 Subject: [PATCH 24/36] add `SciMLBase` to tests --- test/adaptive/sde_complex_adaptive_mean_test.jl | 1 + test/reversal_tests.jl | 1 + test/stratonovich_convergence_tests.jl | 2 ++ test/zerod_noise_test.jl | 1 + 4 files changed, 5 insertions(+) diff --git a/test/adaptive/sde_complex_adaptive_mean_test.jl b/test/adaptive/sde_complex_adaptive_mean_test.jl index 4b1b5650..a427961c 100644 --- a/test/adaptive/sde_complex_adaptive_mean_test.jl +++ b/test/adaptive/sde_complex_adaptive_mean_test.jl @@ -1,4 +1,5 @@ using StochasticDiffEq, DiffEqCallbacks, DiffEqNoiseProcess, Test, Random +import SciMLBase Random.seed!(100) # Definitions according to vector/matrix representations of operators from QuantumOptics diff --git a/test/reversal_tests.jl b/test/reversal_tests.jl index 778c7b56..9164d99d 100644 --- a/test/reversal_tests.jl +++ b/test/reversal_tests.jl @@ -3,6 +3,7 @@ using Random using SDEProblemLibrary # automatically construct SDE transformation for Ito reversal using ModelingToolkit +import SciMLBase # tested solvers additive_noise_solver = [ diff --git a/test/stratonovich_convergence_tests.jl b/test/stratonovich_convergence_tests.jl index 1ac25ed0..2e457dfb 100644 --- a/test/stratonovich_convergence_tests.jl +++ b/test/stratonovich_convergence_tests.jl @@ -1,5 +1,7 @@ using StochasticDiffEq, Test, Random, DiffEqDevTools using SDEProblemLibrary: prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich +import SciMLBase + Random.seed!(100) dts = 1 ./2 .^(10:-1:2) #14->7 good plot diff --git a/test/zerod_noise_test.jl b/test/zerod_noise_test.jl index 80e104ec..65271678 100644 --- a/test/zerod_noise_test.jl +++ b/test/zerod_noise_test.jl @@ -1,4 +1,5 @@ using StochasticDiffEq, Test +import SciMLBase function f(du, u, p, t) du[1] = u[1] From bd686bef312db6d5a71b12aa315217a270dfb66d Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 1 Sep 2024 09:56:20 -0400 Subject: [PATCH 25/36] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4162b434..8566e3f0 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ 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" From bb193ccb420a03c93aa741e634152cb1518cb914 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 1 Sep 2024 17:23:09 -0400 Subject: [PATCH 26/36] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8566e3f0..9938cb99 100644 --- a/Project.toml +++ b/Project.toml @@ -38,7 +38,6 @@ DataStructures = "0.18" DiffEqBase = "6.154" DiffEqNoiseProcess = "5.13" DocStringExtensions = "0.8, 0.9" -EnumX = "1" FiniteDiff = "2" ForwardDiff = "0.10.3" JumpProcesses = "9" From fd814f79617f41a1debb67ef71b014032e1ad711 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 1 Sep 2024 17:23:13 -0400 Subject: [PATCH 27/36] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9938cb99..c166a205 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" From 0995b681c4392bbbdca205cc4df818a3d81141dc Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 1 Sep 2024 17:23:25 -0400 Subject: [PATCH 28/36] Update src/StochasticDiffEq.jl --- src/StochasticDiffEq.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index f7a1350a..af7b4f54 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -73,7 +73,6 @@ else end import SciMLBase - import EnumX using SparseDiffTools: forwarddiff_color_jacobian!, ForwardColorJacCache From 870b882b62429003198581faa4b124120565302f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 1 Sep 2024 22:07:58 -0400 Subject: [PATCH 29/36] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c166a205..d9d97c03 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StochasticDiffEq" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" authors = ["Chris Rackauckas "] -version = "6.68.0" +version = "6.69.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 4c0109f21d567c13aba2cbda1531d6b718b5a3bd Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 19 Sep 2024 12:01:40 -0400 Subject: [PATCH 30/36] Fix `calc_W` for W_transform refactor CI didn't catch this for some reason. --- test/utility_tests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index 532a936b..50815468 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -17,7 +17,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator jac=(u,p,t) -> A) prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) - W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, #=repeat_step=#false, #=W_transform=#true) + W = calc_W(integrator, integrator.cache.nlsolver, dtgamma, #=repeat_step=#false) @test convert(AbstractMatrix, W) ≈ concrete_W @test W \ u0 ≈ concrete_W \ u0 @@ -29,7 +29,7 @@ using StochasticDiffEq.SciMLOperators: MatrixOperator prob = SDEProblem(fun, u0, tspan) integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt) W = integrator.cache.nlsolver.cache.W - calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, #=repeat_step=#false, #=W_transform=#true) + calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, #=repeat_step=#false) # Did not update because it's an array operator # We don't want to build Jacobians when we have operators! From 2b3ec9cca8c3663bea498d0d21b10a05a0d46120 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 19 Sep 2024 13:18:35 -0400 Subject: [PATCH 31/36] horrible nasty hack --- test/utility_tests.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/utility_tests.jl b/test/utility_tests.jl index 50815468..f55a48c3 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -1,7 +1,14 @@ using StochasticDiffEq, LinearAlgebra, SparseArrays, Random, LinearSolve, Test using StochasticDiffEq.OrdinaryDiffEq: WOperator, calc_W!, calc_W using StochasticDiffEq.SciMLOperators: MatrixOperator +using OrdinaryDiffEq +#horid nasty hack to deal with temporary calc_W refactor +# if there is a method that takes a W_transform argument, define the version that doesn't to set W_transform to true +if hasmethod(calc_W, (Any, Any, Any, Any, Any)) + OrdinaryDiffEq.calc_W(integ, nlsolver, dgamma, repeat_step::Bool) = OrdinaryDiffEq.calc_W(integ, nlsolver, dgamma, repeat_step, true) + OrdinaryDiffEq.calc_W!(integ, nlsolver, cache, dgamma, repeat_step::Bool) = OrdinaryDiffEq.calc_W(integ, nlsolver, dgamma, cacherepeat_step, true) +end @testset "Derivative Utilities" begin @testset "calc_W!" begin A = [-1.0 0.0; 0.0 -0.5]; σ = [0.9 0.0; 0.0 0.8] From f10deb09baadae8f81e44fee6f507fdf17dd0528 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 20 Sep 2024 08:49:45 +0200 Subject: [PATCH 32/36] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d9d97c03..b9ad4cf1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StochasticDiffEq" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" authors = ["Chris Rackauckas "] -version = "6.69.0" +version = "6.69.1" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From d7b8d0cad67ac9a2c6ad8bd77bca14c290a50f7c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 20 Oct 2024 07:30:25 -0400 Subject: [PATCH 33/36] Switch to FastPower.jl --- Project.toml | 1 + src/StochasticDiffEq.jl | 2 ++ src/integrators/stepsize_controllers.jl | 4 ++-- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index b9ad4cf1..f659d4b7 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index af7b4f54..2db86dc3 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -35,6 +35,8 @@ using DocStringExtensions using LinearAlgebra, Random import ForwardDiff.Dual + + import FastPower import DiffEqBase: step!, initialize!, DEAlgorithm, AbstractSDEAlgorithm, AbstractRODEAlgorithm, DEIntegrator, AbstractDiffEqInterpolation, diff --git a/src/integrators/stepsize_controllers.jl b/src/integrators/stepsize_controllers.jl index c1ee9ece..01a0206b 100644 --- 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 From f90a9bc80d4d34679f2f1b783ab1da277e10c8b6 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 20 Oct 2024 07:32:07 -0400 Subject: [PATCH 34/36] add compat --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index f659d4b7..6122e412 100644 --- a/Project.toml +++ b/Project.toml @@ -38,6 +38,7 @@ DataStructures = "0.18" DiffEqBase = "6.154" DiffEqNoiseProcess = "5.13" DocStringExtensions = "0.8, 0.9" +FastPower = "1" FiniteDiff = "2" ForwardDiff = "0.10.3" JumpProcesses = "9" From 2b8dea4f58a7eaf338fdc389a35b78b47329cc1a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 20 Oct 2024 08:49:53 -0400 Subject: [PATCH 35/36] Update CI.yml --- .github/workflows/CI.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index d1ff8473..ab7e13ad 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,7 +21,9 @@ jobs: - AlgConvergence3 - WeakConvergence1 version: + - 'lts' - '1' + - 'pre' steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v2 From 2930ba8c7bc113aa57c5b8c41c3dcc6dc2c17dd9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 20 Oct 2024 09:11:22 -0400 Subject: [PATCH 36/36] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6122e412..f47223a1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StochasticDiffEq" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" authors = ["Chris Rackauckas "] -version = "6.69.1" +version = "6.70.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"