Merge pull request #2833 from SciML/sde_throws

Better error throwing on ill-formed SDESystems

Author: ChrisRackauckas

.github/workflows/Downstream.yml

@@ -32,7 +32,7 @@ jobs:
           - {user: SciML, repo: NeuralPDE.jl, group: NNPDE}
           - {user: SciML, repo: DataDrivenDiffEq.jl, group: Downstream}
           - {user: SciML, repo: StructuralIdentifiability.jl, group: All}
-          - {user: SciML, repo: ModelingToolkitStandardLibrary.jl}
+          - {user: SciML, repo: ModelingToolkitStandardLibrary.jl, group: Core}
           - {user: SciML, repo: ModelOrderReduction.jl, group: All}
           - {user: SciML, repo: MethodOfLines.jl, group: Interface}
           - {user: SciML, repo: MethodOfLines.jl, group: 2D_Diffusion}

src/systems/diffeqs/sdesystem.jl

@@ -140,6 +140,10 @@ struct SDESystem <: AbstractODESystem
             check_variables(dvs, iv)
             check_parameters(ps, iv)
             check_equations(deqs, iv)
+            check_equations(neqs, dvs)
+            if size(neqs, 1) != length(deqs)
+                throw(ArgumentError("Noise equations ill-formed. Number of rows must match number of drift equations. size(neqs,1) = $(size(neqs,1)) != length(deqs) = $(length(deqs))"))
+            end
             check_equations(equations(cevents), iv)
         end
         if checks == true || (checks & CheckUnits) > 0

test/sdesystem.jl

@@ -462,9 +462,11 @@ fdif!(du, u0, p, t)
     eqs_short = [D(x) ~ σ * (y - x),
         D(y) ~ x * (ρ - z) - y
     ]
-    sys1 = SDESystem(eqs_short, noiseeqs, t, [x, y, z], [σ, ρ, β], name = :sys1)
-    sys2 = SDESystem(eqs_short, noiseeqs, t, [x, y, z], [σ, ρ, β], name = :sys1)
-    @test_throws ArgumentError SDESystem([sys2.y ~ sys1.z], [], t, [], [],
+    noise_eqs = [y - x
+                 x - y]
+    sys1 = SDESystem(eqs_short, noise_eqs, t, [x, y, z], [σ, ρ, β], name = :sys1)
+    sys2 = SDESystem(eqs_short, noise_eqs, t, [x, y, z], [σ, ρ, β], name = :sys1)
+    @test_throws ArgumentError SDESystem([sys2.y ~ sys1.z], [sys2.y], t, [], [],
         systems = [sys1, sys2], name = :foo)
 end
 
@@ -615,6 +617,17 @@ prob = SDEProblem(sys1, sts .=> [1.0, 0.0, 0.0],
     (0.0, 100.0), ps .=> (10.0, 26.0))
 solve(prob, LambaEulerHeun(), seed = 1)
 
+# Test ill-formed due to more equations than states in noise equations
+
+@parameters p d
+@variables t X(t)
+eqs = [D(X) ~ p - d * X]
+noise_eqs = [sqrt(p), -sqrt(d * X)]
+@test_throws ArgumentError SDESystem(eqs, noise_eqs, t, [X], [p, d]; name = :ssys)
+
+noise_eqs = reshape([sqrt(p), -sqrt(d * X)], 1, 2)
+ssys = SDESystem(eqs, noise_eqs, t, [X], [p, d]; name = :ssys)
+
 # SDEProblem construction with StaticArrays
 # Issue#2814
 @parameters p d

test/symbolic_events.jl

@@ -458,7 +458,7 @@ let
 
     ∂ₜ = D
     eqs = [∂ₜ(A) ~ -k * A]
-    @named ssys = SDESystem(eqs, Equation[], t, [A], [k, t1, t2],
+    @named ssys = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
         discrete_events = [cb1, cb2])
     u0 = [A => 1.0]
     p = [k => 0.0, t1 => 1.0, t2 => 2.0]
@@ -468,7 +468,7 @@ let
     cond1a = (t == t1)
     affect1a = [A ~ A + 1, B ~ A]
     cb1a = cond1a => affect1a
-    @named ssys1 = SDESystem(eqs, Equation[], t, [A, B], [k, t1, t2],
+    @named ssys1 = SDESystem(eqs, [0.0], t, [A, B], [k, t1, t2],
         discrete_events = [cb1a, cb2])
     u0′ = [A => 1.0, B => 0.0]
     sol = testsol(
@@ -478,11 +478,11 @@ let
     # same as above - but with set-time event syntax
     cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once
     cb2‵ = [2.0] => affect2
-    @named ssys‵ = SDESystem(eqs, Equation[], t, [A], [k], discrete_events = [cb1‵, cb2‵])
+    @named ssys‵ = SDESystem(eqs, [0.0], t, [A], [k], discrete_events = [cb1‵, cb2‵])
     testsol(ssys‵, u0, p, tspan; paramtotest = k)
 
     # mixing discrete affects
-    @named ssys3 = SDESystem(eqs, Equation[], t, [A], [k, t1, t2],
+    @named ssys3 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
         discrete_events = [cb1, cb2‵])
     testsol(ssys3, u0, p, tspan; tstops = [1.0], paramtotest = k)
 
@@ -492,24 +492,24 @@ let
         nothing
     end
     cb2‵‵ = [2.0] => (affect!, [], [k], [k], nothing)
-    @named ssys4 = SDESystem(eqs, Equation[], t, [A], [k, t1],
+    @named ssys4 = SDESystem(eqs, [0.0], t, [A], [k, t1],
         discrete_events = [cb1, cb2‵‵])
     testsol(ssys4, u0, p, tspan; tstops = [1.0], paramtotest = k)
 
     # mixing with symbolic condition in the func affect
     cb2‵‵‵ = (t == t2) => (affect!, [], [k], [k], nothing)
-    @named ssys5 = SDESystem(eqs, Equation[], t, [A], [k, t1, t2],
+    @named ssys5 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
         discrete_events = [cb1, cb2‵‵‵])
     testsol(ssys5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)
-    @named ssys6 = SDESystem(eqs, Equation[], t, [A], [k, t1, t2],
+    @named ssys6 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
         discrete_events = [cb2‵‵‵, cb1])
     testsol(ssys6, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)
 
     # mix a continuous event too
     cond3 = A ~ 0.1
     affect3 = [k ~ 0.0]
     cb3 = cond3 => affect3
-    @named ssys7 = SDESystem(eqs, Equation[], t, [A], [k, t1, t2],
+    @named ssys7 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
         discrete_events = [cb1, cb2‵‵‵],
         continuous_events = [cb3])
     sol = testsol(ssys7, u0, p, (0.0, 10.0); tstops = [1.0, 2.0])