Merge pull request #26 from devmotion/abstol_reltol

Fix and test element-wise tolerances for SRI, SRIW1, SRA, SRA1, and RKMil

Author: ChrisRackauckas

src/initdt.jl

@@ -1,7 +1,7 @@
 function sde_determine_initdt{tType,uType}(u0::uType,t::tType,tdir,dtmax,abstol,reltol,internalnorm,prob,order)
   f = prob.f
   g = prob.g
-  d₀ = internalnorm(u0./(abstol+abs.(u0)*reltol))
+  d₀ = internalnorm(u0./(abstol.+abs.(u0).*reltol))
   if !isinplace(prob)
     f₀ = f(t,u0)
     g₀ = 3g(t,u0)

src/integrators/low_order.jl

@@ -159,8 +159,8 @@ end
     @inbounds u[i] = K[i]+L[i]*W.dW[i] + tmp[i]
   end
   if integrator.opts.adaptive
-    @tight_loop_macros for i in eachindex(u)
-      @inbounds tmp[i] = @muladd(tmp[i])/@muladd(integrator.opts.abstol + max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
+    @tight_loop_macros for (i,atol,rtol) in zip(eachindex(u),Iterators.cycle(integrator.opts.abstol),Iterators.cycle(integrator.opts.reltol))
+      @inbounds tmp[i] = @muladd(tmp[i])/@muladd(atol + max(abs(uprev[i]),abs(u[i]))*rtol)
     end
     integrator.EEst = integrator.opts.internalnorm(tmp)
   end

src/integrators/sra.jl

@@ -85,8 +85,9 @@ end
   end
 
   if integrator.opts.adaptive
-    for i in eachindex(u)
-      @inbounds tmp[i] = @muladd(integrator.opts.delta*E₁[i]+E₂[i])/@muladd(integrator.opts.abstol + max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
+    @tight_loop_macros for (i,atol,rtol,δ) in zip(eachindex(u),Iterators.cycle(integrator.opts.abstol),
+						  Iterators.cycle(integrator.opts.reltol),Iterators.cycle(integrator.opts.delta))
+      @inbounds tmp[i] = @muladd(δ*E₁[i]+E₂[i])/@muladd(atol + max(abs(uprev[i]),abs(u[i]))*rtol)
     end
     integrator.EEst = integrator.opts.internalnorm(tmp)
   end
@@ -186,8 +187,9 @@ end
   end
 
   if integrator.opts.adaptive
-    @tight_loop_macros for i in eachindex(u)
-      @inbounds tmp[i] = (integrator.opts.delta*E₁[i]+E₂[i])/(integrator.opts.abstol + max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
+    @tight_loop_macros for (i,atol,rtol,δ) in zip(eachindex(u),Iterators.cycle(integrator.opts.abstol),
+						  Iterators.cycle(integrator.opts.reltol),Iterators.cycle(integrator.opts.delta))
+      @inbounds tmp[i] = @muladd(δ*E₁[i]+E₂[i])/@muladd(atol + max(abs(uprev[i]),abs(u[i]))*rtol)
     end
     integrator.EEst = integrator.opts.internalnorm(tmp)
   end

src/integrators/sri.jl

@@ -113,8 +113,9 @@ end
   end
 
   if integrator.opts.adaptive
-    @tight_loop_macros for i in eachindex(u)
-      @inbounds tmp[i] = @muladd(integrator.opts.delta*E₁[i]+E₂[i])/@muladd(integrator.opts.abstol + max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
+    @tight_loop_macros for (i,atol,rtol,δ) in zip(eachindex(u),Iterators.cycle(integrator.opts.abstol),
+						  Iterators.cycle(integrator.opts.reltol),Iterators.cycle(integrator.opts.delta))
+      @inbounds tmp[i] = @muladd(δ*E₁[i]+E₂[i])/@muladd(atol + max(abs(uprev[i]),abs(u[i]))*rtol)
     end
     integrator.EEst = integrator.opts.internalnorm(tmp)
   end
@@ -208,8 +209,9 @@ end
   end
 
   if integrator.opts.adaptive
-    @tight_loop_macros for i in eachindex(u)
-      @inbounds tmp[i] = @muladd(integrator.opts.delta*E₁[i]+E₂[i])/@muladd(integrator.opts.abstol + max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
+    @tight_loop_macros for (i,atol,rtol,δ) in zip(eachindex(u),Iterators.cycle(integrator.opts.abstol),
+						  Iterators.cycle(integrator.opts.reltol),Iterators.cycle(integrator.opts.delta))
+      @inbounds tmp[i] = @muladd(δ*E₁[i]+E₂[i])/@muladd(atol + max(abs(uprev[i]),abs(u[i]))*rtol)
     end
     integrator.EEst = integrator.opts.internalnorm(tmp)
   end

src/options_type.jl

@@ -1,10 +1,10 @@
-type SDEOptions{uEltype,uEltypeNoUnits,tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsType,ECType,SType,MI}
+type SDEOptions{tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsType,ECType,SType,MI,A,R,D}
   maxiters::MI
   timeseries_steps::Int
   save_everystep::Bool
   adaptive::Bool
-  abstol::uEltype
-  reltol::uEltypeNoUnits
+  abstol::A
+  reltol::R
   gamma::tTypeNoUnits
   qmax::tTypeNoUnits
   qmin::tTypeNoUnits
@@ -24,7 +24,7 @@ type SDEOptions{uEltype,uEltypeNoUnits,tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsT
   dense_errors::Bool
   beta1::tTypeNoUnits
   beta2::tTypeNoUnits
-  delta::uEltypeNoUnits
+  delta::D
   qoldinit::tTypeNoUnits
   dense::Bool
   save_start::Bool

src/solve.jl

@@ -210,15 +210,15 @@ function init{uType,tType,isinplace,algType<:Union{AbstractRODEAlgorithm,Abstrac
 
   uEltype = recursive_eltype(u)
 
-  opts = SDEOptions(maxiters,timeseries_steps,save_everystep,adaptive,uEltype(uEltype(1)*abstol),
-    uEltypeNoUnits(reltol),tTypeNoUnits(gamma),tTypeNoUnits(qmax),tTypeNoUnits(qmin),
+  opts = SDEOptions(maxiters,timeseries_steps,save_everystep,adaptive,map(uEltype,abstol),
+    map(uEltypeNoUnits,reltol),tTypeNoUnits(gamma),tTypeNoUnits(qmax),tTypeNoUnits(qmin),
     dtmax,dtmin,internalnorm,save_idxs,
     tstops_internal,saveat_internal,d_discontinuities_internal,
     userdata,
     progress,progress_steps,
     progress_name,progress_message,
     timeseries_errors,dense_errors,
-    tTypeNoUnits(beta1),tTypeNoUnits(beta2),uEltypeNoUnits(delta),tTypeNoUnits(qoldinit),
+    tTypeNoUnits(beta1),tTypeNoUnits(beta2),map(uEltypeNoUnits,delta),tTypeNoUnits(qoldinit),
     dense,save_start,save_noise,
     callbacks_internal,isoutofdomain,unstable_check,verbose,calck,force_dtmin,
     advance_to_tstop,stop_at_next_tstop)

test/runtests.jl

@@ -26,6 +26,7 @@ LONGER_TESTS && @time @testset "Weak Convergence Tests" begin include("weak_conv
 @time @testset "Composite Tests" begin include("composite_algorithm_test.jl") end
 @time @testset "Events Tests" begin include("events_test.jl") end
 @time @testset "Cache Tests" begin include("cache_test.jl") end
+@time @testset "Element-wise Tolerances Tests" begin include("tolerances_tests.jl") end
 
 #Adaptive SDE
 @time @testset "Adaptive SDE Linear Tests" begin include("adaptive/sde_linearadaptive_tests.jl") end

test/tolerances_tests.jl

@@ -0,0 +1,38 @@
+using StochasticDiffEq, DiffEqProblemLibrary, Base.Test
+
+#=
+f = (t,u) -> begin
+  du = similar(u)
+  prob_sde_2Dlinear.f(t,u,du)
+  return du
+end
+σ = (t,u) -> begin
+  du = similar(u)
+  prob_sde_2Dlinear.g(t,u,du)
+  return du
+end
+
+probs = [prob_sde_2Dlinear,
+         SDEProblem(f,σ,prob_sde_2Dlinear.u0,prob_sde_2Dlinear.tspan)]
+=#
+
+probs = [prob_sde_2Dlinear]
+algs = [SRI(),SRIW1(),SRA1(),SRA(),RKMil()]
+
+for alg in algs, prob in probs
+  dt = typeof(alg)<:StochasticDiffEqAdaptiveAlgorithm ? 0.0 : 0.1
+  srand(100)
+  sol = solve(prob,alg;dt=dt)
+
+  # Vector of element-wise absolute tolerances
+  srand(100)
+  sol2 = solve(prob,alg;dt=dt,abstol=fill(1e-2,4,2))
+  
+  @test sol.t == sol2.t && sol.u == sol2.u
+
+  # Vector of element-wise relative tolerances
+  srand(100)
+  sol2 = solve(prob,alg;dt=dt,reltol=fill(1e-2,4,2))
+
+  @test sol.t == sol2.t && sol.u == sol2.u
+end