automated noise rate prototype

Author: ChrisRackauckas

src/systems/diffeqs/sdesystem.jl

@@ -161,6 +161,7 @@ function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,p=parammap;
                                     version = nothing, tgrad=false,
                                     jac = false, Wfact = false,
                                     checkbounds = false, sparse = false,
+                                    sparsenoise = sparse,
                                     linenumbers = true, parallel=SerialForm(),
                                     kwargs...)
 ```
@@ -172,6 +173,7 @@ function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,parammap=DiffEqBa
                                     version = nothing, tgrad=false,
                                     jac = false, Wfact = false,
                                     checkbounds = false, sparse = false,
+                                    sparsenoise = sparse,
                                     linenumbers = true, parallel=SerialForm(),
                                     kwargs...) where iip
 
@@ -180,7 +182,16 @@ function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,parammap=DiffEqBa
                         sparse=sparse)
     u0 = varmap_to_vars(u0map,states(sys))
     p = varmap_to_vars(parammap,parameters(sys))
-    SDEProblem(f,f.g,u0,tspan,p;kwargs...)
+    if typeof(sys.noiseeqs) <: AbstractVector
+        noise_rate_prototype = nothing
+    elseif sparsenoise
+        I,J,V = findnz(SparseArrays.sparse(sys.noiseeqs))
+        noise_rate_prototype = SparseArrays.sparse(I,J,zero(eltype(u0)))
+    else
+        noise_rate_prototype = zeros(eltype(u0),size(sys.noiseeqs))
+    end
+
+    SDEProblem(f,f.g,u0,tspan,p;noise_rate_prototype=noise_rate_prototype,kwargs...)
 end
 
 function DiffEqBase.SDEProblem(sys::SDESystem, args...; kwargs...)

test/sdesystem.jl

@@ -1,5 +1,5 @@
 using ModelingToolkit, StaticArrays, LinearAlgebra
-using DiffEqBase, StochasticDiffEq
+using DiffEqBase, StochasticDiffEq, SparseArrays
 using Test
 
 # Define some variables
@@ -43,3 +43,25 @@ f = SDEFunction(de)
 prob = SDEProblem(SDEFunction(de),f.g,[1.0,0.0,0.0],(0.0,100.0),(10.0,26.0,2.33),
                   noise_rate_prototype = zeros(3,3))
 sol = solve(prob,EM(),dt=0.001)
+
+u0map = [
+    x => 1.0,
+    y => 0.0,
+    z => 0.0
+]
+
+parammap = [
+    σ => 10.0,
+    β => 26.0,
+    ρ => 2.33
+]
+
+prob = SDEProblem(de,u0map,(0.0,100.0),parammap)
+@test size(prob.noise_rate_prototype) == (3,3)
+@test prob.noise_rate_prototype isa Matrix
+sol = solve(prob,EM(),dt=0.001)
+
+prob = SDEProblem(de,u0map,(0.0,100.0),parammap,sparsenoise=true)
+@test size(prob.noise_rate_prototype) == (3,3)
+@test prob.noise_rate_prototype isa SparseMatrixCSC
+sol = solve(prob,EM(),dt=0.001)