Merge pull request #496 from rmsrosa/noise_hygiene

deepcopy prob.noise

Author: ChrisRackauckas

Project.toml

@@ -34,7 +34,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Adapt = "3"
 ArrayInterface = "2.4, 3.0, 4, 5, 6"
 DataStructures = "0.18"
-DiffEqBase = "6.19"
+DiffEqBase = "6.104"
 JumpProcesses = "9"
 DiffEqNoiseProcess = "5.13"
 DocStringExtensions = "0.8, 0.9"

README.md

@@ -11,7 +11,7 @@ functionality should check out [DifferentialEquations.jl](https://github.com/Sci
 
 ## API
 
-StochasticDiffEq.jl is part of the JuliaDiffEq common interface, but can be used independently of DifferentialEquations.jl. The only requirement is that the user passes an StochasticDiffEq.jl algorithm to `solve`. For example, we can solve the [SDE tutorial from the docs](https://diffeq.sciml.ai/stable/tutorials/sde_example/) using the `SRIW1()` algorithm:
+StochasticDiffEq.jl is part of the JuliaDiffEq common interface, but can be used independently of [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). The only requirement is that the user passes an StochasticDiffEq.jl algorithm to `solve`. For example, we can solve the [SDE tutorial from the docs](https://diffeq.sciml.ai/stable/tutorials/sde_example/) using the `SRIW1()` algorithm:
 
 ```julia
 using StochasticDiffEq
@@ -66,7 +66,7 @@ end
 prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))
 ```
 
-Colored noise can be set using [an `AbstractNoiseProcess`](https://diffeq.sciml.ai/stable/features/noise_process/). For example, we can set the underlying noise process to a `GeometricBrownian` via:
+Colored noise can be set using [an `AbstractNoiseProcess`](https://diffeq.sciml.ai/stable/features/noise_process/). For example, we can set the underlying noise process to a `GeometricBrownianMotionProcess` via:
 
 ```julia
 μ = 1.0

src/solve.jl

@@ -5,6 +5,12 @@ function DiffEqBase.__solve(prob::DiffEqBase.AbstractRODEProblem,
                             kwargs...) where recompile_flag
   integrator = DiffEqBase.__init(prob,alg,timeseries,ts,recompile;kwargs...)
   solve!(integrator)
+  if typeof(prob) <: DiffEqBase.AbstractRODEProblem && typeof(prob.noise) == typeof(integrator.sol.W) && (!haskey(kwargs, :alias_noise) || kwargs[:alias_noise] === true)
+    # would be better to make the following a function `noise_deepcopy!(W::T, Z::T) where {T <: AbstractNoiseProcess}` in `DiffEqNoiseProcess.jl` or a proper `copy` overload, but this should do it for the moment
+    for x in fieldnames(typeof(prob.noise))
+        setfield!(prob.noise, x, deepcopy(getfield(integrator.sol.W, x)))
+    end
+  end
   integrator.sol
 end
 
@@ -412,7 +418,7 @@ function DiffEqBase.__init(
       =#
     end
   elseif typeof(prob) <: DiffEqBase.AbstractRODEProblem
-    W = prob.noise
+    W = (!haskey(kwargs, :alias_noise) || kwargs[:alias_noise] === true) ? deepcopy(prob.noise) : prob.noise
     if W.reset
       # Reseed
       if typeof(W) <: Union{NoiseProcess, NoiseTransport} && W.reseed

test/noise_type_test.jl

@@ -37,6 +37,26 @@ sol = solve(prob,SRA())
 
 @test length(sol.W[1]) == 4
 
+f(du,u,p,t) = (du.=1.01u)
+g(du,u,p,t) = (du.=0.1)
+Z = WienerProcess(0.0, [0.0])
+prob = SDEProblem(f,g,[1.0],(0.0,1.0),noise=Z)
+
+sol = solve(prob,EM(),dt=1/100)
+
+@test sol.W == prob.noise
+@test objectid(prob.noise) != objectid(sol.W)
+
+sol = solve(prob,EM(),dt=1/1000,alias_noise=false)
+
+@test sol.W == prob.noise
+@test objectid(prob.noise) == objectid(sol.W)
+
+sol = solve(prob,EM(),dt=1/1000, alias_noise=true)
+
+@test sol.W == prob.noise
+@test objectid(prob.noise) != objectid(sol.W)
+
 function g(du,u,p,t)
   @test typeof(du) <: SparseMatrixCSC
   du[1,1] = 0.3u[1]
@@ -72,3 +92,18 @@ tspan = (0.0,2.0)
 prob = SDEProblem(drift,vol,u0,tspan, noise=W)
 sol = solve(prob,EM(),dt=0.01)
 @test sol.W.curt ≈ last(tspan)
+
+@test typeof(sol.W) == typeof(prob.noise) <: NoiseFunction
+@test objectid(prob.noise) != objectid(sol.W)
+
+sol = solve(prob,EM(),dt=0.01,alias_noise=true)
+@test sol.W.curt ≈ last(tspan)
+
+@test typeof(sol.W) == typeof(prob.noise) <: NoiseFunction
+@test objectid(prob.noise) != objectid(sol.W)
+
+sol = solve(prob,EM(),dt=0.01,alias_noise=false)
+@test sol.W.curt ≈ last(tspan)
+
+@test typeof(sol.W) == typeof(prob.noise) <: NoiseFunction
+@test objectid(prob.noise) == objectid(sol.W)

test/reversal_tests.jl

@@ -85,6 +85,7 @@ W_reverse = reverse(W_forward)
 
     @test sol_forward(ts).u ≈ sol_reverse(ts).u rtol = 1e-3
     @test length(sol_forward.t) == length(sol_reverse.t)
+    GC.gc()
   end
 end
 
@@ -103,6 +104,7 @@ end
 
     @test sol_forward(ts).u ≈ sol_reverse(ts).u rtol = 1e-2
     @test length(sol_forward.t) == length(sol_reverse.t)
+    GC.gc()
   end
 end
 
@@ -138,5 +140,6 @@ end
       @test sol_forward(ts).u ≈ sol_reverse(ts).u rtol = 1e-2
     end
     @test length(sol_forward.t) == length(sol_reverse.t)
+    GC.gc()
   end
 end