Detect diagonal noise in SDESystem

[MasonProtter]

Fixes #2490

see #2882 (comment) for current state of PR

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

The idea here is that without this change, we were basically not communicating to the solver when the noise equations are either a diagonal matrix (using N independent @brownians once), or a Nx1 matrix (only one @brownian variable which also means diagonal noise rather confusingly which means scalar noise)

Example:

using ModelingToolkit, StochasticDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

# Define some variables
@parameters σ ρ β
@variables x(t) y(t) z(t)
@brownian a
eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
    D(y) ~ x * (ρ - z) - y + 0.1a * y,
    D(z) ~ x * y - β * z + 0.1a * z]

@mtkbuild de = System(eqs, t)

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)
sol = solve(prob, SOSRI())

Result before:

ERROR: The algorithm is not compatible with the chosen noise type. Please see the documentation on the solver methods
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] __init(_prob::SDEProblem{…}, alg::SOSRI, timeseries_init::Vector{…}, ts_init::Vector{…}, ks_init::Type, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_noise::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Rational{…}, qsteady_min::Int64, qsteady_max::Int64, beta2::Nothing, beta1::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, delta::Rational{…}, maxiters::Int64, dtmax::Float64, dtmin::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, force_dtmin::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, initialize_integrator::Bool, seed::UInt64, alias_u0::Bool, alias_jumps::Bool, kwargs::@Kwargs{})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/tdY4b/src/solve.jl:111
  [3] __solve(prob::SDEProblem{…}, alg::SOSRI, timeseries::Vector{…}, ts::Vector{…}, ks::Nothing, recompile::Type{…}; kwargs::@Kwargs{…})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/tdY4b/src/solve.jl:6
  [4] solve_call(_prob::SDEProblem{…}, args::SOSRI; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Sueu7/src/solve.jl:612
  [5] solve_call(_prob::SDEProblem{…}, args::SOSRI)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Sueu7/src/solve.jl:569
  [6] solve_up(prob::SDEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::SOSRI; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Sueu7/src/solve.jl:1062
  [7] solve_up(prob::SDEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::SOSRI)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Sueu7/src/solve.jl:1048
  [8] solve(prob::SDEProblem{…}, args::SOSRI; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Sueu7/src/solve.jl:985
  [9] solve(prob::SDEProblem{…}, args::SOSRI)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/Sueu7/src/solve.jl:975
 [10] top-level scope
    @ REPL[152]:1

Result after:

retcode: Success
Interpolation: 1st order linear
t: 18416-element Vector{Float64}:
   0.0
   0.0007146315987816352
   0.0008575579185379623
   0.0010183500282638302
   0.0011992411517054316
   0.0014027436655772332
   0.00163168399368301
   0.001888330123096511
   0.002169732452494152
   0.0024747179844340924
   0.002801083673034135
   0.0031467337392412443
   0.0035101593109016778
   0.0038893633340266673
   0.004282801251770871
   0.004689082950137011
   0.005106359961984013
   0.005534011073886903
   0.005970909411562779
   0.006416219161068274
   0.0068685852392637164
   ⋮
  99.86081178121842
  99.86691462888227
  99.87284456492715
  99.87822415992166
  99.88366393532804
  99.88978368266024
  99.89589654984118
  99.90277352541975
  99.90884369645995
  99.91567263888018
  99.92335519910294
  99.93199807935355
  99.93899925091941
  99.94687556893099
  99.95573642669403
  99.96333445262874
  99.9718822318053
  99.98149848337891
  99.9917836663466
 100.0
u: 18416-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0]
 [0.9935723575521516, 0.0016623161802513916, 5.887476766870236e-7]
 [0.9931413985835341, 0.001996060215547187, 8.452725403365441e-7]
 [0.9897372191696431, 0.002364209074699878, 1.1910128901146396e-6]
 [0.9863754156964627, 0.0027870349754163344, 1.6442415456946722e-6]
 [0.9816453017746493, 0.0032485982077592327, 2.2360388646766377e-6]
 [0.9798823115339813, 0.00376419538431509, 3.002680985394169e-6]
 [0.9769627498611689, 0.004348290343345211, 3.994742748676605e-6]
 [0.9746752673052053, 0.005000608975232769, 5.2452143832247885e-6]
 [0.975561937347947, 0.005704683532084442, 6.780230831634797e-6]
 [0.9718800868550801, 0.006431432457643715, 8.644914220182036e-6]
 [0.9689080631202299, 0.007244604580622524, 1.0807603471826052e-5]
 [0.9652030007703591, 0.008079161111140092, 1.3371366562575522e-5]
 [0.9590883231049572, 0.008941237325611936, 1.6297379837288528e-5]
 [0.9587393808617176, 0.009811626838090969, 1.9656634680658225e-5]
 [0.9531404195509355, 0.010691604729587828, 2.342847672417148e-5]
 [0.9466340429697403, 0.011613839787154832, 2.7603782378153965e-5]
 [0.9387485145643154, 0.012547710361083157, 3.218172151733701e-5]
 [0.9366750081050457, 0.01348022317055927, 3.703196748463813e-5]
 [0.9342421952139918, 0.014485371063599746, 4.240227183615855e-5]
 [0.9292640616155787, 0.015474981880609914, 4.8194527448965316e-5]
 ⋮
 [5.897812975515105, 6.076645906038707, 1.3286115780193821]
 [5.880039741850383, 6.128572486802333, 1.3365183956845863]
 [5.903266999793311, 6.148150568928711, 1.3399252122262055]
 [5.891448398400042, 6.050028253586368, 1.347903994703874]
 [5.957592323941897, 6.084432178115382, 1.3654956169934376]
 [5.937188474950169, 6.0530838520784505, 1.3667764881070934]
 [5.974121750739623, 6.0357310672631845, 1.3496066140028748]
 [5.959039281116, 6.091080956514798, 1.3704623941006386]
 [6.002353668596953, 6.05426843580676, 1.3823165588516568]
 [6.032884184201906, 5.978308329753016, 1.3839087768390836]
 [6.054295488965005, 5.944376129570052, 1.398591650635967]
 [6.106997470914939, 5.940932286509221, 1.391178399422503]
 [6.045950238598028, 5.949308391126071, 1.3770320989414266]
 [5.997594101075954, 5.994909616470001, 1.355085339471976]
 [5.922569829217446, 6.061347469973805, 1.3627802914472535]
 [5.946840201546448, 6.023691890475072, 1.3700820457176575]
 [5.94821208242336, 5.964936071946218, 1.3800562379600143]
 [5.906584837349537, 5.962649759958491, 1.3818425027449064]
 [5.91918916372324, 6.058807477925927, 1.370347328607035]
 [6.060812868419076, 5.93461436780659, 1.3770640471338815]

[ChrisRackauckas]

or a Nx1 matrix (only one @brownian variable which also means diagonal noise rather confusingly)

For reference, the reason for the noise structure is that the Levy area, computed as the correlation between the Brownian motions, is the thing that is difficult and thus requires splitting SDE solvers into non-commutative, commutative, and "diagonal" methods, where diagonal can specialize the Levy area calculation to effectively be trivial and thus achieve higher order. If you have one Brownian motion, then the Levy area is also trivial, and so the diagonal noise solvers also are higher order on scalar noise with a small change to the form (which the solvers do), and so diagonal noise and scalar noise hit the same higher order solver code path than the other noise types.

Additionally there's additive noise which can hit different solver types, but that does not have signaling at the problem level and is just a solver choice. I would like MTK to signal it somehow, since it's somewhat common and the best solvers are those for additive noise (since that is the most trivial noise case, so you get much higher order).

[ChrisRackauckas]

Downstream failure is unrelated. @baggepinnen do you know what that is?

[MasonProtter]

Before you merge, I just realized one thing I should do to make it a bit more robust

[MasonProtter]

Okay, that should be better. The change I just pushed can now detect cases like when the noise equations would be

[a 0
 b 0
 0 c]

which is also actually diagonal noise.

[ChrisRackauckas]

Why would that be diagonal? 1 and 2 share the same brownian motion in that case?

[MasonProtter]

Wait, so if the noise equation is

[a
 b
 c]

then all three are independent, but if its

[a 0
 b 0
 0 c]

then a and b are shared?

[ChrisRackauckas]

Yes, see https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/#Example-4:-Systems-of-SDEs-with-Non-Diagonal-Noise for how it is all defined.

[MasonProtter]

to be concrete, what I mean is that both of these should describe diagonal noise, right?

    @parameters σ ρ β
    @variables x(t) y(t) z(t)
    @brownian a
    eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
           D(y) ~ x * (ρ - z) - y + 0.1a * y,
           D(z) ~ x * y - β * z + 0.1a * z]

    @mtkbuild de = System(eqs, t)

and

    @parameters σ ρ β
    @variables x(t) y(t) z(t)
    @brownian a b
    eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
            D(y) ~ x * (ρ - z) - y + 0.1a * y,
            D(z) ~ x * y - β * z + 0.1b * z]

    @mtkbuild de = System(eqs, t)

I'd find it rather distressing if the first was considered diagonal, but the second was considered correlated...

[ChrisRackauckas]

The first is a scalar noise SDE, which uses the simplified function form of the diagonal noise https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/#Example-3:-Systems-of-SDEs-with-Scalar-Noise but with a scalar noise process.

The second is non-diagonal and more specifically non-commutative (with commutative defined in https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/#special_noise_forms).

[ChrisRackauckas]

Actually I see your row equation thing is missing the part where "but with a scalar noise process". MTK does not really handle that part of the construction right now, so maybe it's best to keep that to non-diagonal until that can be specialized more.

[MasonProtter]

Hm okay I think there's something I'm just not understanding here. Should I just revert to doing

        # Fix for https://github.com/SciML/ModelingToolkit.jl/issues/2490
        noise_eqs = isdiag(sorted_g_rows) # Does each row have only one non-zero entry?
            diag(sorted_g_rows)
        elseif sorted_g_rows isa AbstractMatrix && size(sorted_g_rows, 2) == 1
            sorted_g_rows[:, 1] # Take a vector slice so solver knows there's no mixing
        else
            sorted_g_rows
        end

?

[ChrisRackauckas]

only make the diagonal one separate for now. Scalar noise needs to change the WienerProcess that is passed and I'm not sure that's handled

[MasonProtter]

Okay, sounds good.

[MasonProtter]

Okay, update on this. I had previously misunderstood what was going on, and when a user uses just one @brownian variable in all equations, that is indeed meant to be scalar noise, not separate diagonal noise functions.

I was confused about this because of the differences in meaning between a Vector and a Matrix of noise equations. We could maybe clear this up in the future by replacing the role of Vector here with LinearAlgebra.Diagonal.

In order to correctly handle scalar noise, I've added another field to SDESystem that tells it if a Vector of noise equations should be treated a scalar noise or not, and then when that gets turned into a SDEProblem, it'll pass the SDEProblem a WienerProcess(0.0, 0.0, 0.0) if the noise is scalar.

Since this package doesn't depend on packages which define WienerProcess, I've added a package extension in order to conditionally define it (this'll happen if the user does using StochasticDiffEq), but maybe I should just add a regular dependency on DiffEqNoiseProcess.jl instead?

[MasonProtter]

Since this package doesn't depend on packages which define WienerProcess, I've added a package extension in order to conditionally define it (this'll happen if the user does using StochasticDiffEq), but maybe I should just add a regular dependency on DiffEqNoiseProcess.jl instead?

Ah, I just noticed that before you could actually create SDEProblems without having StochasticDiffEq loaded, but with this change you'd get an error if you defined one with scalar noise. I guess I have to add a dependency then.