From 32ceae8e9d699f9fb55e3c3b49fbe8e3431887d2 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Wed, 3 Aug 2022 07:52:54 -0300 Subject: [PATCH 1/8] first steps new rode solver --- src/StochasticDiffEq.jl | 1 + src/algorithms.jl | 2 + src/caches/basic_method_caches.jl | 16 ++++ src/perform_step/low_order.jl | 15 ++++ test/rode_linear_tests.jl | 5 ++ test/runtests copy.jl | 133 ++++++++++++++++++++++++++++++ test/runtests.jl | 2 +- 7 files changed, 173 insertions(+), 1 deletion(-) create mode 100644 test/runtests copy.jl diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index 72be39028..78e9a1df8 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -166,6 +166,7 @@ using DocStringExtensions StochasticDiffEqRODECompositeAlgorithm export RandomEM + export RandomHeun export IteratedIntegralApprox, IICommutative, IILevyArea diff --git a/src/algorithms.jl b/src/algorithms.jl index 4e9b75084..6e12881f0 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -848,6 +848,8 @@ end struct RandomEM <: StochasticDiffEqRODEAlgorithm end +struct RandomHeun <: StochasticDiffEqRODEAlgorithm end + const SplitSDEAlgorithms = Union{IIF1M,IIF2M,IIF1Mil,SKenCarp,SplitEM} @doc raw""" diff --git a/src/caches/basic_method_caches.jl b/src/caches/basic_method_caches.jl index 1eebfa252..b120222d0 100644 --- a/src/caches/basic_method_caches.jl +++ b/src/caches/basic_method_caches.jl @@ -72,6 +72,22 @@ function alg_cache(alg::RandomEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prot RandomEMCache(u,uprev,tmp,rtmp) end +struct RandomHeunConstantCache <: StochasticDiffEqConstantCache end +@cache struct RandomHeunCache{uType,rateType} <: StochasticDiffEqMutableCache + u::uType + uprev::uType + tmp::uType + rtmp1::rateType + rtmp2::rateType +end + +alg_cache(alg::RandomHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = RandomHeunConstantCache() + +function alg_cache(alg::RandomHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} + tmp = zero(u); rtmp1 = zero(rate_prototype); rtmp2 = zero(rate_prototype) + RandomHeunCache(u,uprev,tmp,rtmp1,rtmp2) +end + struct SimplifiedEMConstantCache <: StochasticDiffEqConstantCache end @cache struct SimplifiedEMCache{randType,uType,rateType,rateNoiseType} <: StochasticDiffEqMutableCache u::uType diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index 46cc0a66d..378283191 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -121,6 +121,21 @@ end @.. u = uprev + dt * rtmp end +@muladd function perform_step!(integrator,cache::RandomHeunConstantCache) + @unpack t,dt,uprev,u,W,p,f = integrator + ftmp = integrator.f(uprev,p,t,W.curW) + tmp = @.. uprev + dt * ftmp + u = uprev .+ (dt / 2) .* (ftmp .+ integrator.f(tmp,p,t+dt, W.curW .+ W.dW)) # Need to check the last argument, if it should also be in terms of an intermediate step + integrator.u = u +end + +@muladd function perform_step!(integrator,cache::RandomHeunCache) + @unpack rtmp1, rtmp2 = cache + @unpack t,dt,uprev,u,W,p,f = integrator + integrator.f(rtmp1,uprev,p,t,W.curW) + @.. u = uprev + dt * rtmp1 +end + # weak approximation EM @muladd function perform_step!(integrator,cache::SimplifiedEMConstantCache) @unpack t,dt,uprev,u,W,p,f = integrator diff --git a/test/rode_linear_tests.jl b/test/rode_linear_tests.jl index 5ff1dcf99..23a6c6e99 100644 --- a/test/rode_linear_tests.jl +++ b/test/rode_linear_tests.jl @@ -5,17 +5,20 @@ u0 = 1.00 tspan = (0.0,1.0) prob = RODEProblem(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100) +sol2 = solve(prob,RandomHeun(),dt=1/100) f(u,p,t,W,du) = (du.=1.01u.+0.87u.*W) u0 = ones(4) prob = RODEProblem(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100) +sol2 = solve(prob,RandomHeun(),dt=1/100) f(u,p,t,W) = 2u*sin(W) u0 = 1.00 tspan = (0.0,5.0) prob = RODEProblem{false}(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100) +sol2 = solve(prob,RandomHeun(),dt=1/100) function f(du,u,p,t,W) du[1] = 2u[1]*sin(W[1] - W[2]) @@ -25,6 +28,7 @@ u0 = [1.00;1.00] tspan = (0.0,5.0) prob = RODEProblem(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100) +sol2 = solve(prob,RandomHeun(),dt=1/100) function f(du,u,p,t,W) du[1] = -2W[3]*u[1]*sin(W[1] - W[2]) @@ -34,3 +38,4 @@ u0 = [1.00;1.00] tspan = (0.0,5.0) prob = RODEProblem(f,u0,tspan,rand_prototype=zeros(3)) sol = solve(prob,RandomEM(),dt=1/100) +sol2 = solve(prob,RandomHeun(),dt=1/100) diff --git a/test/runtests copy.jl b/test/runtests copy.jl new file mode 100644 index 000000000..a8b505c3f --- /dev/null +++ b/test/runtests copy.jl @@ -0,0 +1,133 @@ +using SafeTestsets +using Pkg + +function activate_gpu_env() + Pkg.activate("gpu") + Pkg.develop(PackageSpec(path=dirname(@__DIR__))) + Pkg.instantiate() +end + +const LONGER_TESTS = false + +const GROUP = get(ENV, "GROUP", "All") + +const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") + +@time begin + if GROUP == "All" || GROUP == "Interface" + @time @safetestset "First Rand Tests" begin include("first_rand_test.jl") end + @time @safetestset "Inference Tests" begin include("inference_test.jl") end + @time @safetestset "Linear RODE Tests" begin include("rode_linear_tests.jl") end + @time @safetestset "Complex Number Tests" begin include("complex_tests.jl") end + @time @safetestset "Static Array Tests" begin include("static_array_tests.jl") end + @time @safetestset "Noise Type Tests" begin include("noise_type_test.jl") end + @time @safetestset "Mass matrix tests" begin include("mass_matrix_tests.jl") end + #@time @safetestset "Sparse Jacobian tests" begin include("sparsediff_tests.jl") end + @time @safetestset "Outofplace Arrays Tests" begin include("outofplace_arrays.jl") end + @time @safetestset "tdir Tests" begin include("tdir_tests.jl") end + @time @safetestset "tstops Tests" begin include("tstops_tests.jl") end + @time @safetestset "saveat Tests" begin include("saveat_tests.jl") end + @time @safetestset "Oval2" begin include("oval2_test.jl") end + end + + if GROUP == "All" || GROUP == "Interface2" + @time @safetestset "Basic Tau Leaping Tests" begin include("tau_leaping.jl") end + @time @safetestset "Linear SDE Tests" begin include("sde/sde_linear_tests.jl") end + @time @safetestset "Two-dimensional Linear SDE Tests" begin include("sde/sde_twodimlinear_tests.jl") end + @time @safetestset "Element-wise Tolerances Tests" begin include("tolerances_tests.jl") end + @time @safetestset "Zero'd Noise Tests" begin include("zerod_noise_test.jl") end + @time @safetestset "Scalar Tests" begin include("scalar_noise.jl") end + @time @safetestset "Stiffness Detection Test" begin include("stiffness_detection_test.jl") end + @time @safetestset "Adaptive SDE Linear Tests" begin include("adaptive/sde_linearadaptive_tests.jl") end + end + + if GROUP == "All" || GROUP == "Interface3" + @time @safetestset "Composite Tests" begin include("composite_algorithm_test.jl") end + @time @safetestset "Events Tests" begin include("events_test.jl") end + @time @safetestset "Cache Tests" begin include("cache_test.jl") end + @time @safetestset "Adaptive Complex Mean Test" begin include("adaptive/sde_complex_adaptive_mean_test.jl") end + @time @safetestset "Utility Tests" begin include("utility_tests.jl") end + @time @safetestset "Non-diagonal SDE Tests" begin include("nondiagonal_tests.jl") end + @time @safetestset "No Index Tests" begin include("noindex_tests.jl") end + @time @safetestset "Multiple Dimension Linear Adaptive Test" begin include("adaptive/sde_twodimlinearadaptive_tests.jl") end + @time @safetestset "Autostepsize Test" begin include("adaptive/sde_autostepsize_test.jl") end + @time @safetestset "Additive Lorenz Attractor Test" begin include("adaptive/sde_lorenzattractor_tests.jl") end + @time @safetestset "Stochastic iterated integrals" begin include("levy_areas.jl") end + end + + if !is_APPVEYOR && (GROUP == "All" || GROUP == "AlgConvergence") + @time @safetestset "Convergence Tests" begin include("sde/sde_convergence_tests.jl") end + @time @safetestset "Dynamical SDE Tests" begin include("sde/sde_dynamical.jl") end + end + + if !is_APPVEYOR && GROUP == "AlgConvergence2" + @time @safetestset "IIF Convergence Tests" begin include("iif_methods.jl") end + @time @safetestset "Commutative Noise Methods Tests" begin include("commutative_tests.jl") end + @time @safetestset "Multivariate Geometric Tests" begin include("multivariate_geometric.jl") end + end + + if !is_APPVEYOR && GROUP == "AlgConvergence3" + @time @safetestset "Rossler Order Tests" begin include("sde/sde_rosslerorder_tests.jl") end + @time @safetestset "ODE Convergence Regression Tests" begin include("ode_convergence_regression.jl") end + @time @safetestset "Additive SDE Tests" begin include("sde/sde_additive_tests.jl") end + @time @safetestset "Split Tests" begin include("split_tests.jl") end + @time @safetestset "Stratonovich Convergence Tests" begin include("stratonovich_convergence_tests.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakConvergence1" + @time @safetestset "Multidimensional IIP Weak Convergence Tests" begin include("weak_convergence/multidim_iip_weak.jl") end + @time @safetestset "Platen's PL1WM weak second order" begin include("weak_convergence/PL1WM.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakConvergence2" + @time @safetestset "Roessler weak SRK Tests" begin include("weak_convergence/srk_weak_final.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakConvergence3" + @time @safetestset "Roessler weak SRK (non-diagonal) Tests" begin include("weak_convergence/srk_weak_final_non_diagonal.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakConvergence4" + @time @safetestset "Weak Stratonovich (non-diagonal) Tests" begin include("weak_convergence/weak_strat_non_diagonal.jl") end + @time @safetestset "SIE SME weak Tests" begin include("weak_convergence/SIE_SME.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakConvergence5" + @time @safetestset "Weak Stratonovich Tests" begin include("weak_convergence/weak_strat.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakConvergence6" + @time @safetestset "Roessler weak SRK diagonal Tests" begin include("weak_convergence/srk_weak_diagonal_final.jl") end + end + + if !is_APPVEYOR && GROUP == "OOPWeakConvergence" + @time @safetestset "OOP Weak Convergence Tests" begin include("weak_convergence/oop_weak.jl") end + @time @safetestset "Additive Weak Convergence Tests" begin include("weak_convergence/additive_weak.jl") end + end + + if !is_APPVEYOR && GROUP == "IIPWeakConvergence" + #activate_gpu_env() + @time @safetestset "IIP Weak Convergence Tests" begin include("weak_convergence/iip_weak.jl") end + end + + if !is_APPVEYOR && GROUP == "SROCKC2WeakConvergence" + #activate_gpu_env() + @time @safetestset "SROCKC2 Weak Convergence Tests" begin include("weak_convergence/weak_srockc2.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakAdaptiveCPU" + @time @safetestset "CPU Weak adaptive step size Brusselator " begin include("adaptive/sde_weak_brusselator_adaptive.jl") end + @time @safetestset "CPU Weak adaptive" begin include("adaptive/sde_weak_adaptive.jl") end + end + + if !is_APPVEYOR && GROUP == "WeakAdaptiveGPU" + activate_gpu_env() + @time @safetestset "GPU Weak adaptive step size scalar noise SDE" begin include("gpu/sde_weak_scalar_adaptive_gpu.jl") end + @time @safetestset "GPU Weak adaptive" begin include("gpu/sde_weak_adaptive_gpu.jl") end + end + + if !is_APPVEYOR && GROUP == "Multithreaded" + @time @safetestset "Mulithreaded Jump Thread Safety Tests" begin include("multithreaded_jump_test.jl") end + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index a8b505c3f..aeb829376 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,7 @@ end const LONGER_TESTS = false -const GROUP = get(ENV, "GROUP", "All") +const GROUP = "Interface" # get(ENV, "GROUP", "All") const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") From 913719502ccd981e7dfda1059a662d27b8dbe30e Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Thu, 11 Aug 2022 12:27:47 -0300 Subject: [PATCH 2/8] alg --- src/alg_utils.jl | 1 + src/perform_step/low_order.jl | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 6e25c42b8..fc1910bb4 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -59,6 +59,7 @@ alg_order(alg::IIF1Mil) = 1 // 1 alg_order(alg::EulerHeun) = 1 // 2 alg_order(alg::LambaEulerHeun) = 1 // 2 alg_order(alg::RandomEM) = 1 // 2 +alg_order(alg::RandomHeun) = 3 // 2 alg_order(alg::SimplifiedEM) = 1 // 2 alg_order(alg::RKMil) = 1 // 1 alg_order(alg::RKMilCommute) = 1 // 1 diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index 378283191..df075f9a1 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -125,7 +125,7 @@ end @unpack t,dt,uprev,u,W,p,f = integrator ftmp = integrator.f(uprev,p,t,W.curW) tmp = @.. uprev + dt * ftmp - u = uprev .+ (dt / 2) .* (ftmp .+ integrator.f(tmp,p,t+dt, W.curW .+ W.dW)) # Need to check the last argument, if it should also be in terms of an intermediate step + u = uprev .+ (dt / 2) .* (ftmp .+ integrator.f(tmp,p,t+dt, W.curW .+ W.dW)) integrator.u = u end @@ -134,6 +134,7 @@ end @unpack t,dt,uprev,u,W,p,f = integrator integrator.f(rtmp1,uprev,p,t,W.curW) @.. u = uprev + dt * rtmp1 + integrator.f(rtmp1,u,p,t,W.curW) end # weak approximation EM From 4ed7a4362ce44b187b95171741e4d9cc3326ada1 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Thu, 11 Aug 2022 16:44:27 -0300 Subject: [PATCH 3/8] fix alg_order for RandomHeun --- src/alg_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index fc1910bb4..7b35e580e 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -59,7 +59,7 @@ alg_order(alg::IIF1Mil) = 1 // 1 alg_order(alg::EulerHeun) = 1 // 2 alg_order(alg::LambaEulerHeun) = 1 // 2 alg_order(alg::RandomEM) = 1 // 2 -alg_order(alg::RandomHeun) = 3 // 2 +alg_order(alg::RandomHeun) = 1 // 2 alg_order(alg::SimplifiedEM) = 1 // 2 alg_order(alg::RKMil) = 1 // 1 alg_order(alg::RKMilCommute) = 1 // 1 From 205d4c3215c2aa87676173074d7255f51bee6fa9 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Thu, 11 Aug 2022 18:13:56 -0300 Subject: [PATCH 4/8] implement step for RandomHeunCache --- src/perform_step/low_order.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index df075f9a1..2b3c0540f 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -125,16 +125,17 @@ end @unpack t,dt,uprev,u,W,p,f = integrator ftmp = integrator.f(uprev,p,t,W.curW) tmp = @.. uprev + dt * ftmp - u = uprev .+ (dt / 2) .* (ftmp .+ integrator.f(tmp,p,t+dt, W.curW .+ W.dW)) + u = uprev .+ (dt/2) .* (ftmp .+ integrator.f(tmp,p,t+dt, W.curW .+ W.dW)) integrator.u = u end @muladd function perform_step!(integrator,cache::RandomHeunCache) - @unpack rtmp1, rtmp2 = cache + @unpack tmp, rtmp1, rtmp2 = cache @unpack t,dt,uprev,u,W,p,f = integrator integrator.f(rtmp1,uprev,p,t,W.curW) - @.. u = uprev + dt * rtmp1 - integrator.f(rtmp1,u,p,t,W.curW) + @.. tmp = uprev + dt * rtmp1 + integrator.f(rtmp2,tmp,p,t+dt,W.curW+W.dW) + @.. u = uprev + (dt/2) * (rtmp1 + rtmp2) end # weak approximation EM From 4a96cd6c4ac76846ec3a90f14e092ee396a7a3e0 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Sat, 13 Aug 2022 08:57:50 -0300 Subject: [PATCH 5/8] add newline --- src/StochasticDiffEq.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index 78e9a1df8..4469d7d5e 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -166,6 +166,7 @@ using DocStringExtensions StochasticDiffEqRODECompositeAlgorithm export RandomEM + export RandomHeun export IteratedIntegralApprox, IICommutative, IILevyArea From 6d81d6471bfc6bc22cfd9e7f98352ce24b45e760 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Sat, 13 Aug 2022 09:03:22 -0300 Subject: [PATCH 6/8] back with original runtests.jl --- test/runtests copy.jl | 133 ------------------------------------------ test/runtests.jl | 2 +- 2 files changed, 1 insertion(+), 134 deletions(-) delete mode 100644 test/runtests copy.jl diff --git a/test/runtests copy.jl b/test/runtests copy.jl deleted file mode 100644 index a8b505c3f..000000000 --- a/test/runtests copy.jl +++ /dev/null @@ -1,133 +0,0 @@ -using SafeTestsets -using Pkg - -function activate_gpu_env() - Pkg.activate("gpu") - Pkg.develop(PackageSpec(path=dirname(@__DIR__))) - Pkg.instantiate() -end - -const LONGER_TESTS = false - -const GROUP = get(ENV, "GROUP", "All") - -const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") - -@time begin - if GROUP == "All" || GROUP == "Interface" - @time @safetestset "First Rand Tests" begin include("first_rand_test.jl") end - @time @safetestset "Inference Tests" begin include("inference_test.jl") end - @time @safetestset "Linear RODE Tests" begin include("rode_linear_tests.jl") end - @time @safetestset "Complex Number Tests" begin include("complex_tests.jl") end - @time @safetestset "Static Array Tests" begin include("static_array_tests.jl") end - @time @safetestset "Noise Type Tests" begin include("noise_type_test.jl") end - @time @safetestset "Mass matrix tests" begin include("mass_matrix_tests.jl") end - #@time @safetestset "Sparse Jacobian tests" begin include("sparsediff_tests.jl") end - @time @safetestset "Outofplace Arrays Tests" begin include("outofplace_arrays.jl") end - @time @safetestset "tdir Tests" begin include("tdir_tests.jl") end - @time @safetestset "tstops Tests" begin include("tstops_tests.jl") end - @time @safetestset "saveat Tests" begin include("saveat_tests.jl") end - @time @safetestset "Oval2" begin include("oval2_test.jl") end - end - - if GROUP == "All" || GROUP == "Interface2" - @time @safetestset "Basic Tau Leaping Tests" begin include("tau_leaping.jl") end - @time @safetestset "Linear SDE Tests" begin include("sde/sde_linear_tests.jl") end - @time @safetestset "Two-dimensional Linear SDE Tests" begin include("sde/sde_twodimlinear_tests.jl") end - @time @safetestset "Element-wise Tolerances Tests" begin include("tolerances_tests.jl") end - @time @safetestset "Zero'd Noise Tests" begin include("zerod_noise_test.jl") end - @time @safetestset "Scalar Tests" begin include("scalar_noise.jl") end - @time @safetestset "Stiffness Detection Test" begin include("stiffness_detection_test.jl") end - @time @safetestset "Adaptive SDE Linear Tests" begin include("adaptive/sde_linearadaptive_tests.jl") end - end - - if GROUP == "All" || GROUP == "Interface3" - @time @safetestset "Composite Tests" begin include("composite_algorithm_test.jl") end - @time @safetestset "Events Tests" begin include("events_test.jl") end - @time @safetestset "Cache Tests" begin include("cache_test.jl") end - @time @safetestset "Adaptive Complex Mean Test" begin include("adaptive/sde_complex_adaptive_mean_test.jl") end - @time @safetestset "Utility Tests" begin include("utility_tests.jl") end - @time @safetestset "Non-diagonal SDE Tests" begin include("nondiagonal_tests.jl") end - @time @safetestset "No Index Tests" begin include("noindex_tests.jl") end - @time @safetestset "Multiple Dimension Linear Adaptive Test" begin include("adaptive/sde_twodimlinearadaptive_tests.jl") end - @time @safetestset "Autostepsize Test" begin include("adaptive/sde_autostepsize_test.jl") end - @time @safetestset "Additive Lorenz Attractor Test" begin include("adaptive/sde_lorenzattractor_tests.jl") end - @time @safetestset "Stochastic iterated integrals" begin include("levy_areas.jl") end - end - - if !is_APPVEYOR && (GROUP == "All" || GROUP == "AlgConvergence") - @time @safetestset "Convergence Tests" begin include("sde/sde_convergence_tests.jl") end - @time @safetestset "Dynamical SDE Tests" begin include("sde/sde_dynamical.jl") end - end - - if !is_APPVEYOR && GROUP == "AlgConvergence2" - @time @safetestset "IIF Convergence Tests" begin include("iif_methods.jl") end - @time @safetestset "Commutative Noise Methods Tests" begin include("commutative_tests.jl") end - @time @safetestset "Multivariate Geometric Tests" begin include("multivariate_geometric.jl") end - end - - if !is_APPVEYOR && GROUP == "AlgConvergence3" - @time @safetestset "Rossler Order Tests" begin include("sde/sde_rosslerorder_tests.jl") end - @time @safetestset "ODE Convergence Regression Tests" begin include("ode_convergence_regression.jl") end - @time @safetestset "Additive SDE Tests" begin include("sde/sde_additive_tests.jl") end - @time @safetestset "Split Tests" begin include("split_tests.jl") end - @time @safetestset "Stratonovich Convergence Tests" begin include("stratonovich_convergence_tests.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakConvergence1" - @time @safetestset "Multidimensional IIP Weak Convergence Tests" begin include("weak_convergence/multidim_iip_weak.jl") end - @time @safetestset "Platen's PL1WM weak second order" begin include("weak_convergence/PL1WM.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakConvergence2" - @time @safetestset "Roessler weak SRK Tests" begin include("weak_convergence/srk_weak_final.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakConvergence3" - @time @safetestset "Roessler weak SRK (non-diagonal) Tests" begin include("weak_convergence/srk_weak_final_non_diagonal.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakConvergence4" - @time @safetestset "Weak Stratonovich (non-diagonal) Tests" begin include("weak_convergence/weak_strat_non_diagonal.jl") end - @time @safetestset "SIE SME weak Tests" begin include("weak_convergence/SIE_SME.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakConvergence5" - @time @safetestset "Weak Stratonovich Tests" begin include("weak_convergence/weak_strat.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakConvergence6" - @time @safetestset "Roessler weak SRK diagonal Tests" begin include("weak_convergence/srk_weak_diagonal_final.jl") end - end - - if !is_APPVEYOR && GROUP == "OOPWeakConvergence" - @time @safetestset "OOP Weak Convergence Tests" begin include("weak_convergence/oop_weak.jl") end - @time @safetestset "Additive Weak Convergence Tests" begin include("weak_convergence/additive_weak.jl") end - end - - if !is_APPVEYOR && GROUP == "IIPWeakConvergence" - #activate_gpu_env() - @time @safetestset "IIP Weak Convergence Tests" begin include("weak_convergence/iip_weak.jl") end - end - - if !is_APPVEYOR && GROUP == "SROCKC2WeakConvergence" - #activate_gpu_env() - @time @safetestset "SROCKC2 Weak Convergence Tests" begin include("weak_convergence/weak_srockc2.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakAdaptiveCPU" - @time @safetestset "CPU Weak adaptive step size Brusselator " begin include("adaptive/sde_weak_brusselator_adaptive.jl") end - @time @safetestset "CPU Weak adaptive" begin include("adaptive/sde_weak_adaptive.jl") end - end - - if !is_APPVEYOR && GROUP == "WeakAdaptiveGPU" - activate_gpu_env() - @time @safetestset "GPU Weak adaptive step size scalar noise SDE" begin include("gpu/sde_weak_scalar_adaptive_gpu.jl") end - @time @safetestset "GPU Weak adaptive" begin include("gpu/sde_weak_adaptive_gpu.jl") end - end - - if !is_APPVEYOR && GROUP == "Multithreaded" - @time @safetestset "Mulithreaded Jump Thread Safety Tests" begin include("multithreaded_jump_test.jl") end - end - -end diff --git a/test/runtests.jl b/test/runtests.jl index aeb829376..a8b505c3f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,7 @@ end const LONGER_TESTS = false -const GROUP = "Interface" # get(ENV, "GROUP", "All") +const GROUP = get(ENV, "GROUP", "All") const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") From cf5c12b0862fdea95056a05e3b4f298f312e4e37 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Sat, 13 Aug 2022 16:10:12 -0300 Subject: [PATCH 7/8] add basic comparison tests for RODE solvers --- test/rode_linear_tests.jl | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/test/rode_linear_tests.jl b/test/rode_linear_tests.jl index 23a6c6e99..98dce55f2 100644 --- a/test/rode_linear_tests.jl +++ b/test/rode_linear_tests.jl @@ -1,34 +1,43 @@ -using StochasticDiffEq +using StochasticDiffEq, DiffEqNoiseProcess, Random, Test +Random.seed!(100) f(u,p,t,W) = 1.01u.+0.87u.*W u0 = 1.00 tspan = (0.0,1.0) prob = RODEProblem(f,u0,tspan) -sol = solve(prob,RandomEM(),dt=1/100) +sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) +prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) +@test abs(sol[end] - sol2[end]) < 0.1 -f(u,p,t,W,du) = (du.=1.01u.+0.87u.*W) +f!(du,u,p,t,W) = (du.=1.01u.+0.87u.*W) u0 = ones(4) -prob = RODEProblem(f,u0,tspan) -sol = solve(prob,RandomEM(),dt=1/100) +prob = RODEProblem(f!,u0,tspan) +sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) +prob = RODEProblem(f!,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) +@test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end]) f(u,p,t,W) = 2u*sin(W) u0 = 1.00 tspan = (0.0,5.0) prob = RODEProblem{false}(f,u0,tspan) -sol = solve(prob,RandomEM(),dt=1/100) +sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) +prob = RODEProblem{false}(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) +@test abs(sol[end]-sol2[end]) < 0.1 function f(du,u,p,t,W) du[1] = 2u[1]*sin(W[1] - W[2]) du[2] = -2u[2]*cos(W[1] + W[2]) end u0 = [1.00;1.00] -tspan = (0.0,5.0) +tspan = (0.0,4.0) prob = RODEProblem(f,u0,tspan) -sol = solve(prob,RandomEM(),dt=1/100) +sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) +prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) +@test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end]) function f(du,u,p,t,W) du[1] = -2W[3]*u[1]*sin(W[1] - W[2]) @@ -37,5 +46,7 @@ end u0 = [1.00;1.00] tspan = (0.0,5.0) prob = RODEProblem(f,u0,tspan,rand_prototype=zeros(3)) -sol = solve(prob,RandomEM(),dt=1/100) +sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) +prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) +@test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end]) From d949ca16f323c6a94d951cb54cdd2c9554a39a66 Mon Sep 17 00:00:00 2001 From: Ricardo Rosa Date: Sun, 14 Aug 2022 08:02:25 -0300 Subject: [PATCH 8/8] remove alloc in perform_step for nonscalar noise --- src/caches/basic_method_caches.jl | 7 ++++--- src/perform_step/low_order.jl | 14 ++++++++++---- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/caches/basic_method_caches.jl b/src/caches/basic_method_caches.jl index b120222d0..fbdef1b39 100644 --- a/src/caches/basic_method_caches.jl +++ b/src/caches/basic_method_caches.jl @@ -73,19 +73,20 @@ function alg_cache(alg::RandomEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prot end struct RandomHeunConstantCache <: StochasticDiffEqConstantCache end -@cache struct RandomHeunCache{uType,rateType} <: StochasticDiffEqMutableCache +@cache struct RandomHeunCache{uType,rateType,randType} <: StochasticDiffEqMutableCache u::uType uprev::uType tmp::uType rtmp1::rateType rtmp2::rateType + wtmp::randType end alg_cache(alg::RandomHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = RandomHeunConstantCache() function alg_cache(alg::RandomHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} - tmp = zero(u); rtmp1 = zero(rate_prototype); rtmp2 = zero(rate_prototype) - RandomHeunCache(u,uprev,tmp,rtmp1,rtmp2) + tmp = zero(u); rtmp1 = zero(rate_prototype); rtmp2 = zero(rate_prototype); wtmp = zero(ΔW) + RandomHeunCache(u,uprev,tmp,rtmp1,rtmp2,wtmp) end struct SimplifiedEMConstantCache <: StochasticDiffEqConstantCache end diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index 2b3c0540f..e14524577 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -124,17 +124,23 @@ end @muladd function perform_step!(integrator,cache::RandomHeunConstantCache) @unpack t,dt,uprev,u,W,p,f = integrator ftmp = integrator.f(uprev,p,t,W.curW) - tmp = @.. uprev + dt * ftmp - u = uprev .+ (dt/2) .* (ftmp .+ integrator.f(tmp,p,t+dt, W.curW .+ W.dW)) + tmp = @.. uprev + dt * ftmp + wtmp = @.. W.curW + W.dW + u = uprev .+ (dt/2) .* (ftmp .+ integrator.f(tmp,p,t+dt, wtmp)) integrator.u = u end @muladd function perform_step!(integrator,cache::RandomHeunCache) - @unpack tmp, rtmp1, rtmp2 = cache + @unpack tmp, rtmp1, rtmp2, wtmp = cache @unpack t,dt,uprev,u,W,p,f = integrator integrator.f(rtmp1,uprev,p,t,W.curW) @.. tmp = uprev + dt * rtmp1 - integrator.f(rtmp2,tmp,p,t+dt,W.curW+W.dW) + if W.dW isa Number + wtmp = W.curW + W.dW + else + @.. wtmp = W.curW + W.dW + end + integrator.f(rtmp2,tmp,p,t+dt,wtmp) @.. u = uprev + (dt/2) * (rtmp1 + rtmp2) end