Merge pull request #75 from frankschae/ensemble_work_precision

Work Precision Set for EnsembleProblem

Author: ChrisRackauckas

.travis.yml

@@ -11,10 +11,7 @@ julia:
 #    - julia: nightly
 notifications:
   email: false
-# uncomment the following lines to override the default test script
-#script:
-#  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
-#  - julia -e 'Pkg.clone(pwd()); Pkg.build("DiffEqDevTools"); Pkg.test("DiffEqDevTools"; coverage=true)'
+script: julia -e 'using Pkg; Pkg.build(); Pkg.test(coverage=false)'
 after_success:
   # push coverage results to Coveralls
   - julia -e 'cd(Pkg.dir("DiffEqDevTools")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'

src/benchmark.jl

@@ -394,6 +394,122 @@ function WorkPrecisionSet(prob::AbstractRODEProblem,abstols,reltols,setups,test_
   WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,error_estimate,numruns_error)
 end
 
+
+function WorkPrecisionSet(prob::AbstractEnsembleProblem,abstols,reltols,setups,test_dt=nothing;
+                          numruns=5, trajectories=1000,
+                          print_names=false,names=nothing,appxsol_setup=nothing,expected_value=nothing,
+                          error_estimate=:weak_final,ensemblealg=EnsembleThreads(),
+                          kwargs...)
+
+  @assert names === nothing || length(setups) == length(names)
+
+  weak_timeseries_errors = error_estimate ∈ WEAK_TIMESERIES_ERRORS
+  weak_dense_errors = error_estimate ∈ WEAK_DENSE_ERRORS
+
+  N = length(setups); M = length(abstols)
+  times = Array{Float64}(undef,M,N)
+  solutions = Array{Any}(undef,M,N)
+  if names === nothing
+    names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
+  end
+  time_tmp = Vector{Float64}(undef,numruns)
+
+  # First calculate all of the errors
+  for k in 1:N
+    for j in 1:M
+      if !haskey(setups[1],:dts)
+        sol = solve(prob,setups[k][:alg],ensemblealg;
+          setups[k]...,
+          abstol=abstols[j],
+          reltol=reltols[j],
+          timeseries_errors=false,
+          dense_errors = false,
+          trajectories=Int(trajectories),kwargs...)
+      else
+        sol = solve(prob,setups[k][:alg],ensemblealg;
+          setups[k]...,
+          abstol=abstols[j],
+          reltol=reltols[j],
+          dt=setups[k][:dts][j],
+          timeseries_errors=false,
+          dense_errors = false,
+          trajectories=Int(trajectories),kwargs...)
+      end
+      solutions[j,k] = sol
+    end
+    @info "$(setups[k][:alg]) ($k/$N)"
+  end
+
+  if error_estimate ∈ WEAK_ERRORS
+    if expected_value != nothing
+      errors = [[LinearAlgebra.norm(Statistics.mean(solutions[i,j].u .- expected_value))
+        for i in 1:M] for j in 1:N]
+    else
+      sol = solve(prob,appxsol_setup[:alg],ensemblealg;kwargs...,appxsol_setup...,
+          timeseries_errors=false,dense_errors = false,trajectories=Int(trajectories))
+      errors = [[LinearAlgebra.norm(Statistics.mean(solutions[i,j].u .- sol.u))
+        for i in 1:M] for j in 1:N]
+    end
+  else
+    error("use RODEProblem instead of EnsembleProblem for strong errors.")
+  end
+
+  local _sol
+
+  # Now time it
+  for k in 1:N
+    # precompile
+    GC.gc()
+    if !haskey(setups[1],:dts)
+      _sol = solve(prob,setups[k][:alg],ensemblealg;
+        setups[k]...,
+        abstol=abstols[1],
+        reltol=reltols[1],
+        timeseries_errors=false,
+        dense_errors = false,
+        trajectories=Int(trajectories),kwargs...)
+    else
+      _sol = solve(prob,setups[k][:alg],ensemblealg;
+        setups[k]...,
+        abstol=abstols[1],
+        reltol=reltols[1],
+        dt=setups[k][:dts][1],
+        timeseries_errors=false,
+        dense_errors = false,
+        trajectories=Int(trajectories),kwargs...)
+    end
+    #x = isempty(_sol.t) ? 0 : round(Int,mean(_sol.t) - sum(_sol.t)/length(_sol.t))
+    GC.gc()
+    for j in 1:M
+      for i in 1:numruns
+        time_tmp[i] = @elapsed if !haskey(setups[k],:dts)
+          sol = solve(prob,setups[k][:alg],ensemblealg;
+              setups[k]...,
+              abstol=abstols[j],
+              reltol=reltols[j],
+              timeseries_errors=false,
+              dense_errors = false,
+              trajectories=Int(trajectories),kwargs...)
+        else
+          sol = solve(prob,setups[k][:alg],ensemblealg;
+              setups[k]...,
+              abstol=abstols[j],
+              reltol=reltols[j],
+              dt=setups[k][:dts][j],
+              timeseries_errors=false,
+              dense_errors = false,
+              trajectories=Int(trajectories),kwargs...)
+        end
+      end
+      times[j,k] = mean(time_tmp) #+ x
+      GC.gc()
+    end
+  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,error_estimate,Int(trajectories))
+end
+
 function get_sample_errors(prob::AbstractRODEProblem,setup,test_dt=nothing;
                           appxsol_setup=nothing,
                           numruns,error_estimate=:final,

test/analyticless_convergence_tests.jl

@@ -40,6 +40,44 @@ sim2 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=100,
 @test abs(sim2.𝒪est[:final]-1.5) < 0.3
 @show sim2.𝒪est[:final]
 
+
+# EnsembleProblem
+
+function prob_func(prob, i, repeat)
+    remake(prob,seed=seeds[i])
+end
+
+u₀ = [1.0,1.0]
+function f2!(du,u,p,t)
+  du[1] = -273//512*u[1]
+  du[2] = -1//160*u[1]-(-785//512+sqrt(2)/8)*u[2]
+end
+function g2!(du,u,p,t)
+  du[1,1] = 1//4*u[1]
+  du[1,2] = 1//16*u[1]
+  du[2,1] = (1-2*sqrt(2))/4*u[1]
+  du[2,2] = 1//10*u[1]+1//16*u[2]
+end
+dts = 1 .//2 .^(3:-1:0)
+tspan = (0.0,3.0)
+
+h2(z) = z^2 # but apply it only to u[1]
+
+prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2))
+
+numtraj = Int(1e5)
+seed = 100
+Random.seed!(seed)
+seeds = rand(UInt, numtraj)
+ensemble_prob = EnsembleProblem(prob;
+        output_func = (sol,i) -> (h2(sol[end][1]),false),
+        prob_func = prob_func
+        )
+sim = test_convergence(dts,ensemble_prob,DRI1(),save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,weak_timeseries_errors=false,weak_dense_errors=false,expected_value=exp(-3.0))
+
+@test abs(sim.𝒪est[:weak_final]-2.0) < 0.3
+@show sim.𝒪est[:weak_final]
+
 ### SDDE
 
 function hayes_modelf(du,u,h,p,t)

test/analyticless_stochastic_wp.jl

@@ -37,3 +37,68 @@ se2 = get_sample_errors(prob2,setups[1],test_dt,appxsol_setup = appxsol_setup,
                         numruns=[5,10,25,50,100],solution_runs=20)
 
 @test all(se[1:5]-se2 .< 1e-1)
+
+
+# Ensemble Problem with non-commutative noise process
+
+function prob_func(prob, i, repeat)
+    remake(prob,seed=seeds[i])
+end
+
+u₀ = [1.0,1.0]
+function f2!(du,u,p,t)
+  du[1] = -273//512*u[1]
+  du[2] = -1//160*u[1]-(-785//512+sqrt(2)/8)*u[2]
+  return nothing
+end
+function g2!(du,u,p,t)
+  du[1,1] = 1//4*u[1]
+  du[1,2] = 1//16*u[1]
+  du[2,1] = (1-2*sqrt(2))/4*u[1]
+  du[2,2] = 1//10*u[1]+1//16*u[2]
+  return nothing
+end
+dts = 1 .//2 .^(3:-1:0)
+tspan = (0.0,3.0)
+
+h2(z) = z^2 # but apply it only to u[1]
+
+prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2))
+
+numtraj = Int(1e5)
+seed = 100
+Random.seed!(seed)
+seeds = rand(UInt, numtraj)
+ensemble_prob = EnsembleProblem(prob;
+        output_func = (sol,i) -> (h2(sol[end][1]),false),
+        prob_func = prob_func
+        )
+
+reltols = 1.0 ./ 4.0 .^ (1:4)
+abstols = reltols#[0.0 for i in eachindex(reltols)]
+setups = [
+        Dict(:alg=>EM(),:dts=>dts),
+        Dict(:alg=>SimplifiedEM(),:dts=>dts),
+        Dict(:alg=>DRI1(),:dts=>dts, :adaptive=>false)
+        ]
+test_dt = 1//1000
+appxsol_setup = Dict(:alg=>EM(), :dt=>test_dt)
+
+# without analytical expectation value
+wp1 = @time WorkPrecisionSet(ensemble_prob,abstols,reltols,setups,test_dt;
+                            maxiters = 1e7,
+                            verbose=false,save_everystep=false, save_start=false,
+                            appxsol_setup = appxsol_setup,
+                            trajectories=numtraj, error_estimate=:weak_final)
+
+# with analytical expectation value
+wp2 = @time WorkPrecisionSet(ensemble_prob,abstols,reltols,setups,test_dt;
+                            maxiters = 1e7,
+                            verbose=false,save_everystep=false, save_start=false,
+                            appxsol_setup = appxsol_setup,expected_value=exp(-3.0),
+                            trajectories=numtraj, error_estimate=:weak_final)
+
+err1 = [wp1.wps[i].errors for i=1:length(setups)]
+err2 = [wp2.wps[i].errors for i=1:length(setups)]
+
+@test isapprox(err1, err2, atol=1e-4)