Merge remote-tracking branch 'origin/master'

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "DiffEqDevTools"
 uuid = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "2.12.0"
+version = "2.15.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

src/benchmark.jl

@@ -36,35 +36,46 @@ function Shootout(prob,setups;appxsol=nothing,names=nothing,error_estimate=:fina
   effratios = Matrix{Float64}(undef,N,N)
   timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
   dense_errors = error_estimate ∈ DENSE_ERRORS
-  if names == nothing
+  if names === nothing
     names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
   end
   for i in eachindex(setups)
     sol = solve(prob,setups[i][:alg];timeseries_errors=timeseries_errors,
-    dense_errors = dense_errors,kwargs...,setups[i]...) # Compile and get result
+    dense_errors = dense_errors,kwargs...,setups[i]...)
 
-    if appxsol != nothing
-      errsol = appxtrue(sol,appxsol)
+    if :prob_choice ∈ keys(setups[i])
+      cur_appxsol = appxsol[setups[i][:prob_choice]]
+    else
+      cur_appxsol = appxsol
+    end
+
+    if cur_appxsol != cur_appxsol
+      errsol = appxtrue(sol,cur_appxsol)
       errors[i] = errsol.errors[error_estimate]
       solutions[i] = errsol
     else
       errors[i] = sol.errors[error_estimate]
       solutions[i] = sol
     end
 
-    benchmark_f = let prob=prob,alg=setups[i][:alg],sol=sol,kwargs=kwargs
-      function benchmark_f()
-        @elapsed solve(prob,alg,(sol.u),(sol.t),(sol.k);
-              timeseries_errors = false,
-              dense_errors = false, kwargs...)
-      end
+    if haskey(setups[i], :prob_choice)
+      _prob = prob[setups[i][:prob_choice]]
+    else
+      _prob = prob
+    end
+
+    benchmark_f = let _prob=_prob,alg=setups[i][:alg],sol=sol,kwargs=kwargs
+      () -> @elapsed solve(_prob, alg, sol.u, sol.t, sol.k;
+                           timeseries_errors = false,
+                           dense_errors = false, kwargs...)
     end
+    benchmark_f() # pre-compile
 
-    b_t =  benchmark_f()
+    b_t = benchmark_f()
     if b_t > seconds
       times[i] = b_t
     else
-      times[i] = minimum([b_t;map(i->benchmark_f(),2:numruns)])
+      times[i] = mapreduce(i -> benchmark_f(), min, 2:numruns; init = b_t)
     end
 
     effs[i] = 1/(errors[i]*times[i])
@@ -86,10 +97,10 @@ function ShootoutSet(probs,setups;probaux=nothing,
   N = length(probs)
   shootouts = Vector{Shootout}(undef,N)
   winners = Vector{String}(undef,N)
-  if names == nothing
+  if names === nothing
     names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
   end
-  if probaux == nothing
+  if probaux === nothing
     probaux = Vector{Dict{Symbol,Any}}(undef,N)
     for i in 1:N
       probaux[i] = Dict{Symbol,Any}()
@@ -143,7 +154,6 @@ mutable struct WorkPrecisionSet
   prob
   setups
   names
-  sample_error
   error_estimate
   numruns
 end
@@ -153,68 +163,80 @@ function WorkPrecision(prob,alg,abstols,reltols,dts=nothing;
   N = length(abstols)
   errors = Vector{Float64}(undef,N)
   times = Vector{Float64}(undef,N)
-  if name == nothing
+  if name === nothing
     name = "WP-Alg"
   end
-  timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
-  dense_errors = error_estimate ∈ DENSE_ERRORS
-  for i in 1:N
-    # Calculate errors and precompile
-    if dts == nothing
-      sol = solve(prob,alg;kwargs...,abstol=abstols[i],
-      reltol=reltols[i],timeseries_errors=timeseries_errors,
-      dense_errors = dense_errors) # Compile and get result
-    else
-      sol = solve(prob,alg;kwargs...,abstol=abstols[i],
-      reltol=reltols[i],dt=dts[i],timeseries_errors=timeseries_errors,
-      dense_errors = dense_errors) # Compile and get result
-    end
 
-    if appxsol != nothing
-      errsol = appxtrue(sol,appxsol)
-      errors[i] = mean(errsol.errors[error_estimate])
-    else
-      errors[i] = mean(sol.errors[error_estimate])
-    end
+  if haskey(kwargs, :prob_choice)
+    _prob = prob[kwargs[:prob_choice]]
+  else
+    _prob = prob
+  end
 
-    benchmark_f = let dts=dts,prob=prob,alg=alg,sol=sol,abstols=abstols,reltols=reltols,kwargs=kwargs
-      function benchmark_f()
-        if dts == nothing
-          @elapsed solve(prob,alg,(sol.u),(sol.t),(sol.k);
-                abstol=(abstols[i]),
-                reltol=(reltols[i]),
-                timeseries_errors = false,
-                dense_errors = false, kwargs...)
+  let _prob = _prob
+    timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
+    dense_errors = error_estimate ∈ DENSE_ERRORS
+    for i in 1:N
+      if dts === nothing
+        sol = solve(_prob,alg;kwargs...,abstol=abstols[i],
+        reltol=reltols[i],timeseries_errors=timeseries_errors,
+        dense_errors = dense_errors)
+      else
+        sol = solve(_prob,alg;kwargs...,abstol=abstols[i],
+        reltol=reltols[i],dt=dts[i],timeseries_errors=timeseries_errors,
+        dense_errors = dense_errors)
+      end
+
+      if haskey(kwargs, :prob_choice)
+        cur_appxsol = appxsol[kwargs[:prob_choice]]
+      else
+        cur_appxsol = appxsol
+      end
+
+      if cur_appxsol !== nothing
+        errsol = appxtrue(sol,cur_appxsol)
+        errors[i] = mean(errsol.errors[error_estimate])
+      else
+        errors[i] = mean(sol.errors[error_estimate])
+      end
+
+      benchmark_f = let dts=dts,_prob=_prob,alg=alg,sol=sol,abstols=abstols,reltols=reltols,kwargs=kwargs
+        if dts === nothing
+          () -> @elapsed solve(_prob, alg, sol.u, sol.t, sol.k;
+                               abstol = abstols[i],
+                               reltol = reltols[i],
+                               timeseries_errors = false,
+                               dense_errors = false, kwargs...)
         else
-          @elapsed solve(prob,alg,(sol.u),(sol.t),(sol.k);
-                abstol=(abstols[i]),
-                reltol=(reltols[i]),
-                dt=(dts[i]),
-                timeseries_errors = false,
-                dense_errors = false, kwargs...)
+          () -> @elapsed solve(_prob, alg, sol.u, sol.t, sol.k;
+                               abstol = abstols[i],
+                               reltol = reltols[i],
+                               dt = dts[i],
+                               timeseries_errors = false,
+                               dense_errors = false, kwargs...)
         end
       end
-    end
+      benchmark_f() # pre-compile
 
-    b_t =  benchmark_f()
-    if b_t > seconds
-      times[i] = b_t
-    else
-      times[i] = minimum([b_t;map(i->benchmark_f(),2:numruns)])
+      b_t = benchmark_f()
+      if b_t > seconds
+        times[i] = b_t
+      else
+        times[i] = mapreduce(i -> benchmark_f(), min, 2:numruns; init = b_t)
+      end
     end
   end
   return WorkPrecision(prob,abstols,reltols,errors,times,name,N)
 end
 
-function WorkPrecisionSet(prob::Union{AbstractODEProblem,AbstractDDEProblem,
-                                      AbstractDAEProblem},
+function WorkPrecisionSet(prob,
                           abstols,reltols,setups;
                           print_names=false,names=nothing,appxsol=nothing,
                           error_estimate=:final,
                           test_dt=nothing,kwargs...)
   N = length(setups)
   wps = Vector{WorkPrecision}(undef,N)
-  if names == nothing
+  if names === nothing
     names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
   end
   for i in 1:N
@@ -231,7 +253,7 @@ function WorkPrecisionSet(prob::Union{AbstractODEProblem,AbstractDDEProblem,
                                  name=names[i],kwargs...,setups[i]...)
     end
   end
-  return WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,nothing,error_estimate,nothing)
+  return WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,error_estimate,nothing)
 end
 
 @def error_calculation begin
@@ -285,7 +307,8 @@ end
 function WorkPrecisionSet(prob::AbstractRODEProblem,abstols,reltols,setups,test_dt=nothing;
                           numruns=20,numruns_error = 20,
                           print_names=false,names=nothing,appxsol_setup=nothing,
-                          error_estimate=:final,parallel_type = :none,kwargs...)
+                          error_estimate=:final,parallel_type = :none,
+                          kwargs...)
 
   timeseries_errors = DiffEqBase.has_analytic(prob.f) && error_estimate ∈ TIMESERIES_ERRORS
   weak_timeseries_errors = error_estimate ∈ WEAK_TIMESERIES_ERRORS
@@ -294,7 +317,7 @@ function WorkPrecisionSet(prob::AbstractRODEProblem,abstols,reltols,setups,test_
   N = length(setups); M = length(abstols)
   times = Array{Float64}(undef,M,N)
   tmp_solutions = Array{Any}(undef,numruns_error,M,N)
-  if names == nothing
+  if names === nothing
     names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
   end
   time_tmp = Vector{Float64}(undef,numruns)
@@ -310,8 +333,7 @@ function WorkPrecisionSet(prob::AbstractRODEProblem,abstols,reltols,setups,test_
       @error_calculation
     end
   end
-  analytical_solution_ends = [tmp_solutions[i,1,1].u_analytic[end] for i in 1:numruns_error]
-  sample_error = 1.96std(norm.(analytical_solution_ends))/sqrt(numruns_error)
+
   _solutions_k = [[EnsembleSolution(tmp_solutions[:,j,k],0.0,true) for j in 1:M] for k in 1:N]
   solutions = [[DiffEqBase.calculate_ensemble_errors(sim;weak_timeseries_errors=weak_timeseries_errors,weak_dense_errors=weak_dense_errors) for sim in sol_k] for sol_k in _solutions_k]
   if error_estimate ∈ WEAK_ERRORS
@@ -365,46 +387,62 @@ function WorkPrecisionSet(prob::AbstractRODEProblem,abstols,reltols,setups,test_
   end
 
   wps = [WorkPrecision(prob,abstols,reltols,errors[i],times[:,i],names[i],N) for i in 1:N]
-  WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,sample_error,error_estimate,numruns_error)
+  WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,error_estimate,numruns_error)
 end
 
-@def sample_errors begin
-  if !DiffEqBase.has_analytic(prob.f)
-    true_sol = solve(prob,appxsol_setup[:alg];kwargs...,appxsol_setup...,
-                     save_everystep=false)
-    analytical_solution_ends[i] = norm(true_sol.u[end])
-  else
-    _dt = prob.tspan[2] - prob.tspan[1]
-    if typeof(prob.u0) <: Number
-      W = sqrt(_dt)*randn()
-    else
-      W = sqrt(_dt)*randn(size(prob.u0))
+function get_sample_errors(prob::AbstractRODEProblem,setup,test_dt=nothing;
+                          appxsol_setup=nothing,
+                          numruns,error_estimate=:final,
+                          sample_error_runs = Int(1e7),
+                          solution_runs,
+                          parallel_type = :none,kwargs...)
+
+  maxnumruns = findmax(numruns)[1]
+
+  tmp_solutions_full = map(1:solution_runs) do i
+    @info "Solution Run: $i"
+    # Use the WorkPrecision stuff to calculate the errors
+    tmp_solutions = Array{Any}(undef,maxnumruns,1,1)
+    setups = [setup]
+    abstols = [1e-2] # Standard default
+    reltols = [1e-2] # Standard default
+    M = 1; N = 1
+    timeseries_errors = false; dense_errors = false
+    if parallel_type == :threads
+      Threads.@threads for i in 1:maxnumruns
+        @error_calculation
+      end
+    elseif parallel_type == :none
+      for i in 1:maxnumruns
+        @error_calculation
+      end
     end
-    analytical_solution_ends[i] = norm(prob.f.analytic(prob.u0,prob.p,prob.tspan[2],W))
+    tmp_solutions = vec(tmp_solutions)
   end
-end
 
-function get_sample_errors(prob::AbstractRODEProblem,test_dt=nothing;
-                          appxsol_setup=nothing,
-                          numruns=20,std_estimation_runs = maximum(numruns),
-                          error_estimate=:final,parallel_type = :none,kwargs...)
-  _std_estimation_runs = Int(std_estimation_runs)
-  analytical_solution_ends = Vector{typeof(norm(prob.u0))}(undef,_std_estimation_runs)
-  if parallel_type == :threads
-    Threads.@threads for i in 1:_std_estimation_runs
-      @sample_errors
-    end
-  elseif parallel_type == :none
-    for i in 1:_std_estimation_runs
-      @info "Standard deviation estimation: $i/$_std_estimation_runs"
-      @sample_errors
+  if DiffEqBase.has_analytic(prob.f)
+    analytical_mean_end = mean(1:sample_error_runs) do i
+      _dt = prob.tspan[2] - prob.tspan[1]
+      if typeof(prob.u0) <: Number
+        W = sqrt(_dt)*randn()
+      else
+        W = sqrt(_dt)*randn(size(prob.u0))
+      end
+      prob.f.analytic(prob.u0,prob.p,prob.tspan[2],W)
     end
+  else
+    # Use the mean of the means as the analytical mean
+    analytical_mean_end = mean(mean(tmp_solutions[i].u[end] for i in 1:length(tmp_solutions)) for tmp_solutions in tmp_solutions_full)
   end
-  est_std = std(analytical_solution_ends)
-  if typeof(numruns) <: Number
-    return 1.96est_std/sqrt(numruns)
+
+  if numruns isa Number
+    mean_solution_ends = [mean([tmp_solutions[i].u[end] for i in 1:maxnumruns]) for tmp_solutions in tmp_solutions_full]
+    return sample_error = 1.96std(norm(mean_sol_end - analytical_mean_end) for mean_sol_end in mean_solution_ends)/sqrt(numruns)
   else
-    return [1.96est_std/sqrt(_numruns) for _numruns in numruns]
+    map(1:length(numruns)) do i
+      mean_solution_ends = [mean([tmp_solutions[i].u[end] for i in 1:numruns[i]]) for tmp_solutions in tmp_solutions_full]
+      sample_error = 1.96std(norm(mean_sol_end - analytical_mean_end) for mean_sol_end in mean_solution_ends)/sqrt(numruns[i])
+    end
   end
 end
 

src/convergence.jl

@@ -36,7 +36,7 @@ function ConvergenceSimulation(solutions,convergence_axis;
 end
 
 function test_convergence(dts::AbstractArray,prob::Union{AbstractRODEProblem,AbstractSDEProblem},
-                          alg;trajectories=10000,save_everystep=true,timeseries_steps=1,
+                          alg;trajectories,save_everystep=true,timeseries_steps=1,
                           timeseries_errors=save_everystep,adaptive=false,
                           weak_timeseries_errors=false,weak_dense_errors=false,kwargs...)
   N = length(dts)
@@ -66,8 +66,13 @@ function analyticless_test_convergence(dts::AbstractArray,
   for j in 1:trajectories
     @info "Monte Carlo iteration: $j/$trajectories"
     t = prob.tspan[1]:test_dt:prob.tspan[2]
-    brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
-    brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
+    if prob.noise_rate_prototype === nothing
+      brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
+      brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
+    else
+      brownian_values = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
+      brownian_values2 = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
+    end
     np = NoiseGrid(t,brownian_values,brownian_values2)
     _prob = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,
                        noise=np,