Merge pull request #49 from SciML/nondiag

fix wrapper handling of non-diagonal sparse noise

Author: ChrisRackauckas

src/integrators/interface.jl

@@ -107,7 +107,11 @@ function DiffEqBase.auto_dt_reset!(integrator::SDDEIntegrator)
   end
 
   # determine initial time step
-  sde_prob = SDEProblem(f,g, prob.u0, prob.tspan, prob.p)
+  sde_prob = SDEProblem(f,g, prob.u0, prob.tspan, prob.p;
+                        noise_rate_prototype = prob.noise_rate_prototype,
+                        noise = prob.noise,
+                        seed = prob.seed,
+                        prob.kwargs...)
   integrator.dt = StochasticDiffEq.sde_determine_initdt(integrator.u, integrator.t,
       integrator.tdir, integrator.opts.dtmax, integrator.opts.abstol, integrator.opts.reltol,
       integrator.opts.internalnorm, sde_prob, StochasticDiffEq.get_current_alg_order(getalg(integrator.alg), integrator.cache), integrator)
@@ -244,7 +248,7 @@ function DiffEqBase.u_modified!(integrator::SDDEIntegrator, bool::Bool)
 end
 
 get_proposed_dt(integrator::SDDEIntegrator) = integrator.dtpropose
-set_proposed_dt!(integrator::SDDEIntegrator,dt::Number) = (integrator.dtpropose = dt)
+set_proposed_dt!(integrator::SDDEIntegrator,dt::Number) = (integrator.dtpropose = dt; integrator.dtcache = dt)
 
 function set_proposed_dt!(integrator::SDDEIntegrator,integrator2::SDDEIntegrator)
   integrator.dtpropose = integrator2.dtpropose
@@ -337,3 +341,15 @@ end
   !isnothing(integrator.W) && DiffEqNoiseProcess.save_noise!(integrator.W)
   !isnothing(integrator.P) && DiffEqNoiseProcess.save_noise!(integrator.P)
 end
+
+@inline DiffEqBase.get_tmp_cache(integrator::SDDEIntegrator) =
+  get_tmp_cache(integrator, integrator.alg, integrator.cache)
+# avoid method ambiguity
+for typ in (StochasticDiffEq.StochasticDiffEqAlgorithm,StochasticDiffEq.StochasticDiffEqNewtonAdaptiveAlgorithm)
+  @eval @inline DiffEqBase.get_tmp_cache(integrator::SDDEIntegrator, alg::$typ, cache::StochasticDiffEq.StochasticDiffEqConstantCache) = nothing
+end
+@inline DiffEqBase.get_tmp_cache(integrator::SDDEIntegrator, alg, cache) = (cache.tmp,)
+@inline DiffEqBase.get_tmp_cache(integrator::SDDEIntegrator, alg::StochasticDiffEq.StochasticDiffEqNewtonAdaptiveAlgorithm, cache) =
+    (cache.nlsolver.tmp, cache.nlsolver.ztmp)
+@inline DiffEqBase.get_tmp_cache(integrator::SDDEIntegrator, alg::StochasticDiffEq.StochasticCompositeAlgorithm, cache) =
+    get_tmp_cache(integrator, alg.algs[1], cache.caches[1])

test/nondiagonal_sparse_noise.jl

@@ -0,0 +1,66 @@
+using StochasticDelayDiffEq
+using Random
+using SparseArrays
+
+function sir_dde!(du,u,h,p,t)
+    (S,I,R) = u
+    (β,c,τ) = p
+    N = S+I+R
+    infection = β*c*I/N*S
+    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
+    Nd = Sd+Id+Rd
+    recovery = β*c*Id/Nd*Sd
+    @inbounds begin
+        du[1] = -infection
+        du[2] = infection - recovery
+        du[3] = recovery
+    end
+    nothing
+end;
+
+# Define a sparse matrix by making a dense matrix and setting some values as not zero
+A = zeros(3,2)
+A[1,1] = 1
+A[2,1] = 1
+A[2,2] = 1
+A[3,2] = 1
+A = SparseArrays.sparse(A);
+
+# Make `g` write the sparse matrix values
+function sir_delayed_noise!(du,u,h,p,t)
+    (S,I,R) = u
+    (β,c,τ) = p
+    N = S+I+R
+    infection = β*c*I/N*S
+    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
+    Nd = Sd+Id+Rd
+    recovery = β*c*Id/Nd*Sd
+    du[1,1] = -sqrt(infection)
+    du[2,1] = sqrt(infection)
+    du[2,2] = -sqrt(recovery)
+    du[3,2] = sqrt(recovery)
+end;
+
+function condition(u,t,integrator) # Event when event_f(u,t) == 0
+  u[2]
+end;
+function affect!(integrator)
+  integrator.u[2] = 0.0
+end;
+cb = ContinuousCallback(condition,affect!);
+
+δt = 0.1
+tmax = 40.0
+tspan = (0.0,tmax)
+t = 0.0:δt:tmax;
+u0 = [990.0,10.0,0.0]; # S,I,R
+
+function sir_history(p, t)
+    [1000.0, 0.0, 0.0]
+end;
+
+p = [0.05,10.0,4.0]; # β,c,τ
+Random.seed!(1234);
+
+prob_sdde = SDDEProblem(sir_dde!,sir_delayed_noise!,u0,sir_history,tspan,p;noise_rate_prototype=A);
+sol_sdde = solve(prob_sdde,LambaEM(),callback=cb);

test/runtests.jl

@@ -3,3 +3,4 @@ using SafeTestsets
 @safetestset "SDDEProblem, solve" begin include("test_prob_sol.jl") end
 @safetestset "Analyticless Convergence Tests" begin include("analyticless_convergence_tests.jl") end
 @safetestset "Event handling" begin include("events.jl") end
+@safetestset "Non-Diagonal Sparse Noise" begin include("nondiagonal_sparse_noise.jl") end