# Oval2 Timings

##### Chris Rackauckas
using StochasticDiffEq, DiffEqProblemLibrary, Random, Base.Threads
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
prob = DiffEqProblemLibrary.SDEProblemLibrary.oval2ModelExample(largeFluctuations=true,useBigs=false)
prob_func(prob,i,repeat) = remake(prob,seed=i)
prob = EnsembleProblem(remake(prob,tspan=(0.0,1.0)),prob_func=prob_func)
js = 16:21
dts = 1.0 ./ 2.0 .^ (js)
trajectories = 1000
fails = Array{Int}(undef,length(dts),3)
times = Array{Float64}(undef,length(dts),3)

6×3 Array{Float64,2}:
6.95333e-310  6.95333e-310  6.95333e-310
6.95333e-310  6.95333e-310  6.95333e-310
6.95333e-310  6.95333e-310  6.95333e-310
6.95333e-310  6.95333e-310  6.95333e-310
6.95333e-310  6.95333e-310  6.95333e-310
6.95333e-310  6.95333e-310  6.95333e-310


## Timing Runs

sol = solve(prob,SRIW1(),EnsembleThreads(),abstol=2.0^(-13),reltol=2.0^(-7),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 4.593863072

sol = solve(prob,SRI(error_terms=2),EnsembleThreads(),abstol=2.0^(-13),reltol=2.0^(-7),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 7.535223829

sol = solve(prob,SRI(),EnsembleThreads(),abstol=2.0^(-14),reltol=2.0^(-18),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 82.958798375

sol = solve(prob,SRI(tableau=StochasticDiffEq.constructSRIOpt1()),EnsembleThreads(),abstol=2.0^(-7),reltol=2.0^(-4),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 0.848733018

sol = solve(prob,SOSRI(),EnsembleThreads(),abstol=2.0^(-7),reltol=2.0^(-4),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 0.478400502

sol = solve(prob,SOSRI(),EnsembleThreads(),abstol=2.0^(-7),reltol=2.0^(-6),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 0.525046647

sol = solve(prob,SOSRI(),EnsembleThreads(),abstol=2.0^(-12),reltol=2.0^(-15),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 12.166897742

sol = solve(prob,SOSRI(),EnsembleThreads(),abstol=2.0^(-13),reltol=2.0^(-7),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 1.22488415

sol = solve(prob,SOSRI(),EnsembleThreads(),abstol=2.0^(-12),reltol=2.0^(-15),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 12.088284559

sol = solve(prob,SOSRI2(),EnsembleThreads(),abstol=2.0^(-12),reltol=2.0^(-15),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 12.041467194

sol = solve(prob,SOSRI2(),EnsembleThreads(),abstol=2.0^(-13),reltol=2.0^(-11),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 6.167647134

sol = solve(prob,SOSRI2(),EnsembleThreads(),abstol=2.0^(-13),reltol=2.0^(-11),maxiters=Int(1e11),qmax=1.125,save_everystep=false,trajectories=Threads.nthreads())
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Adaptive Fails is $numfails. Elapsed time was$adaptive_time")

The number of Adaptive Fails is 0. Elapsed time was 5.762501974

for j in eachindex(js)
println("j = $j") sol =solve(prob,EM(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=Threads.nthreads()) t1 = @elapsed sol = solve(prob,EM(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=trajectories) numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories]) println("The number of Euler-Maruyama Fails is$numfails. Elapsed time was $t1") fails[j,1] = numfails times[j,1] = t1 end  j = 1 The number of Euler-Maruyama Fails is 20. Elapsed time was 13.253878568 j = 2 The number of Euler-Maruyama Fails is 5. Elapsed time was 26.154132367 j = 3 The number of Euler-Maruyama Fails is 1. Elapsed time was 51.547233009 j = 4 The number of Euler-Maruyama Fails is 0. Elapsed time was 101.9276721 j = 5 The number of Euler-Maruyama Fails is 0. Elapsed time was 202.383144237 j = 6 The number of Euler-Maruyama Fails is 0. Elapsed time was 402.002865982  for j in 1:4 println("j =$j")
t1 = @elapsed sol = solve(prob,SRIW1(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=trajectories)
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of SRIW1 Fails is $numfails. Elapsed time was$t1")
fails[j,3] = numfails
times[j,3] = t1
end

j = 1
The number of SRIW1 Fails is 971. Elapsed time was 0.498085818
j = 2
The number of SRIW1 Fails is 977. Elapsed time was 0.568782522
j = 3
The number of SRIW1 Fails is 978. Elapsed time was 0.531061161
j = 4
The number of SRIW1 Fails is 981. Elapsed time was 0.512728842

js = 17:21
dts = 1.0 ./2.0 .^ (js)
for j in 1:6
println("j = $j") sol =solve(prob,ImplicitEM(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=Threads.nthreads()) t1 = @elapsed sol = solve(prob,ImplicitEM(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=trajectories) numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories]) println("The number of Implicit-EM Fails is$numfails. Elapsed time was $t1") end  j = 1 The number of Implicit-EM Fails is 35. Elapsed time was 0.24586918 j = 2 The number of Implicit-EM Fails is 31. Elapsed time was 0.332139545 j = 3 The number of Implicit-EM Fails is 41. Elapsed time was 0.254886108 j = 4 The number of Implicit-EM Fails is 33. Elapsed time was 0.261064087 j = 5 The number of Implicit-EM Fails is 31. Elapsed time was 0.38440035 j = 6  ERROR: BoundsError: attempt to access 5-element Array{Float64,1} at index   js = 17:21 dts = 1.0 ./ 2.0 .^(js) for j in 1:6 println("j =$j")
t1 = @elapsed sol = solve(prob,ImplicitRKMil(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=trajectories)
numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories])
println("The number of Implicit-RKMil Fails is $numfails. Elapsed time was$t1")
end

j = 1
The number of Implicit-RKMil Fails is 28. Elapsed time was 4.387731263
j = 2
The number of Implicit-RKMil Fails is 31. Elapsed time was 4.509489773
j = 3
The number of Implicit-RKMil Fails is 31. Elapsed time was 4.509685423
j = 4
The number of Implicit-RKMil Fails is 16. Elapsed time was 4.532243258
j = 5
The number of Implicit-RKMil Fails is 23. Elapsed time was 4.402556784
j = 6

ERROR: BoundsError: attempt to access 5-element Array{Float64,1} at index 

for j in 1:6
println("j = $j") sol =solve(prob,RKMil(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=Threads.nthreads()) t1 = @elapsed sol = solve(prob,RKMil(),EnsembleThreads(),dt=dts[j],maxiters=Int(1e11),save_everystep=false,verbose=false,trajectories=trajectories) numfails = sum([Int(any(isnan,sol[i]) || sol[i].t[end] != 1) for i in 1:trajectories]) println("The number of RKMil Fails is$numfails. Elapsed time was $t1") fails[j,2] = numfails times[j,2] = t1 end  j = 1 The number of RKMil Fails is 3. Elapsed time was 2.131215569 j = 2 The number of RKMil Fails is 5. Elapsed time was 2.089111541 j = 3 The number of RKMil Fails is 4. Elapsed time was 2.080083787 j = 4 The number of RKMil Fails is 5. Elapsed time was 2.071570627 j = 5 The number of RKMil Fails is 4. Elapsed time was 2.043506636 j = 6  ERROR: BoundsError: attempt to access 5-element Array{Float64,1} at index   using Plots lw = 3 p2 = plot(dts,times,xscale=:log2,yscale=:log2,guidefont=font(16),tickfont=font(14),yguide="Elapsed Time (s)",xguide=L"Chosen$\Delta t\$",top_margin=50px,linewidth=lw,lab=["Euler-Maruyama" "RK-Mil" "RosslerSRI"],legendfont=font(14))

ERROR: LoadError: UndefVarError: @L_str not defined
in expression starting at none:1

plot!(dts,repmat([best_adaptive_time],11),linewidth=lw,line=:dash,lab="ESRK+RSwM3",left_margin=75px)

ERROR: UndefVarError: repmat not defined

scatter!([2.0^(-20);2.0^(-20);2.0^(-18)],[times[5,1];times[5,2];times[3,3]],markersize=20,c=:red,lab="")
plot(p2,size=(800,800))

ERROR: UndefVarError: p2 not defined

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


## Appendix

These benchmarks are a part of the DiffEqBenchmarks.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl

To locally run this tutorial, do the following commands:

using DiffEqBenchmarks
DiffEqBenchmarks.weave_file("StiffSDE","Oval2Timings.jmd")

Computer Information:

Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
JULIA_DEPOT_PATH = /builds/JuliaGPU/DiffEqBenchmarks.jl/.julia
JULIA_CUDA_MEMORY_LIMIT = 2147483648
JULIA_PROJECT = @.

Status: /builds/JuliaGPU/DiffEqBenchmarks.jl/benchmarks/StiffSDE/Project.toml
[10745b16-79ce-11e8-11f9-7d13ad32a3b2] Statistics