update to seeded RNGs

Author: ChrisRackauckas

REQUIRE

@@ -1,6 +1,6 @@
-julia 0.6.0-pre
+julia 0.6.0
 Parameters 0.5.0
-DiffEqBase 1.10.0
+DiffEqBase 1.13.0
 RecursiveArrayTools 0.8.0
 DataStructures 0.4.6
 Juno 0.2.5

src/StochasticDiffEq.jl

@@ -6,6 +6,8 @@ module StochasticDiffEq
 
   import Base: start, next, done, eltype
 
+  import RandomNumbers: Xorshifts
+
   using Reexport
   @reexport using DiffEqBase
 

src/integrators/integrator_interface.jl

@@ -91,9 +91,9 @@ end
 
 @inline function fill_new_noise_caches!(integrator,c,scaling_factor,idxs)
   if isinplace(integrator.W)
-    integrator.W.dist(@view(c[2][idxs]),integrator.W,scaling_factor)
+    integrator.W.dist(@view(c[2][idxs]),integrator.W,scaling_factor,integrator.W.rng)
     if alg_needs_extra_process(integrator.alg)
-      integrator.W.dist(@view(c[3][idxs]),integrator.W,scaling_factor)
+      integrator.W.dist(@view(c[3][idxs]),integrator.W,scaling_factor,integrator.W.rng)
     end
   else
     c[2][idxs] .= integrator.noise(length(idxs),integrator,scaling_factor)

src/solve.jl

@@ -271,26 +271,32 @@ function init{uType,tType,isinplace,algType<:Union{AbstractRODEAlgorithm,Abstrac
     randType = typeof(rand_prototype) # Strip units and type info
   end
 
+  prob.seed == 0 ? seed = rand(UInt64) : seed = prob.seed
+
   if typeof(prob.noise) <: Void
     if isinplace
       if alg_needs_extra_process(alg)
         W = WienerProcess!(t,rand_prototype,rand_prototype,
                            save_everystep=save_everystep,
-                           timeseries_steps=timeseries_steps)
+                           timeseries_steps=timeseries_steps,
+                           rng = Xorshifts.Xoroshiro128Plus(seed))
       else
         W = WienerProcess!(t,rand_prototype,
                            save_everystep=save_everystep,
-                           timeseries_steps=timeseries_steps)
+                           timeseries_steps=timeseries_steps,
+                           rng = Xorshifts.Xoroshiro128Plus(seed))
       end
     else
       if alg_needs_extra_process(alg)
         W = WienerProcess(t,rand_prototype,rand_prototype,
                            save_everystep=save_everystep,
-                           timeseries_steps=timeseries_steps)
+                           timeseries_steps=timeseries_steps,
+                           rng = Xorshifts.Xoroshiro128Plus(seed))
       else
         W = WienerProcess(t,rand_prototype,
                            save_everystep=save_everystep,
-                           timeseries_steps=timeseries_steps)
+                           timeseries_steps=timeseries_steps,
+                           rng = Xorshifts.Xoroshiro128Plus(seed))
       end
     end
   else
@@ -318,7 +324,7 @@ function init{uType,tType,isinplace,algType<:Union{AbstractRODEAlgorithm,Abstrac
 
   sol = build_solution(prob,alg,ts,timeseries,W=W,
                 calculate_error = false,
-                interp = id, dense = dense)
+                interp = id, dense = dense, seed = seed)
 
   if recompile_flag == true
     FType = typeof(f)

test/iif_methods.jl

@@ -23,16 +23,16 @@ sol2 = solve(prob2,EM(),dt=1/10)
 srand(100)
 dts = 1./2.^(7:-1:4) #14->7 good plot
 println("IIF scalar")
-sim  = test_convergence(dts,prob,IIF1M(),numMonte=Int(2e1))
-@test abs(sim.𝒪est[:l2]-1) < 0.2 # closer to 1 at this part
+sim  = test_convergence(dts,prob,IIF1M(),numMonte=Int(1e2))
+@test abs(sim.𝒪est[:l2]-0.5) < 0.2 # closer to 1 at this part
 sim  = test_convergence(dts,no_noise_prob,IIF1M(),numMonte=Int(2e1))
-@test abs(sim.𝒪est[:l2]-1) < 0.2 # closer to 1 at this part
+@test abs(sim.𝒪est[:l2]-1.0) < 0.2 # closer to 1 at this part
 
 dts = 1./2.^(7:-1:4) #14->7 good plot
 println("IIF no noise scalar")
 srand(100)
-sim  = test_convergence(dts,prob,IIF2M(),numMonte=Int(2e1))
-@test abs(sim.𝒪est[:l2]-1) < 0.2 # closer to 1 at this part
+sim  = test_convergence(dts,prob,IIF2M(),numMonte=Int(1e2))
+@test abs(sim.𝒪est[:l2]-0.5) < 0.2 # closer to 1 at this part
 sim  = test_convergence(dts,no_noise_prob,IIF2M(),numMonte=Int(1e1))
 @test abs(sim.𝒪est[:l2]-2) < 0.2 # closer to 1 at this part