Incorrect(?) dynamics of SDE systems in latest version

[MasonProtter]

Sorry this isn't the most minimal possible MWE, but I got it down as small as I could. Here's a distillation of a test case from Neuroblox that passed on v9.75, and now fails on v9.76:

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

function test_mtk(t_dur = 1.0)
    function make_sys(; name, θ=1.0, ϕ=0.1)
        sts = @variables x(t)=0.0 [output=true] y(t)=0.0 jcn(t) [input=true]
        p = @parameters θ=θ ϕ=ϕ
        @brownian ξ
        eqs = [
            D(x) ~ y,
            D(y) ~ θ*(1-x^2)*y - x + ϕ*ξ + jcn
        ]
        System(eqs, t, sts, p; name)
    end
    @named sys1 = make_sys()
    @named sys2 = make_sys()

    connection_eqs = [
        sys1.jcn ~ sys2.x 
        sys2.jcn ~ sys1.x
    ]
    @named conn_sys = System(connection_eqs, t, [], [])
    sys = compose(conn_sys, [sys1, sys2])
    
    sys_final = structural_simplify(sys)
    prob = SDEProblem(sys_final, [], (0.0, t_dur), []; seed=1234)
    sol = solve(prob, LambaEM())
    sol(t_dur, idxs=[sys1.x, sys1.y, sys2.x, sys2.y])
end

function test_normal(t_dur = 1.0)
    function f!(du, u, p, t)
        (; θ1, ϕ1, θ2, ϕ2) = p
        (x1, y1, x2, y2) = u
        jcn1 = x2
        jcn2 = x1
        #=D(x1)=# du[1] = y1
        #=D(y1)=# du[2] = θ1*(1-x1^2)*y1 - x1 + jcn1
        #=D(x2)=# du[3] = y2
        #=D(y2)=# du[4] = θ2*(1-x2^2)*y2 - x2 + jcn2
        nothing
    end
    function noise!(du, u, p, t)
        (; ϕ1, ϕ2) = p
        du[2] = ϕ1
        du[4] = ϕ2
        nothing
    end
    let
        p = (;θ1=1.0, ϕ1=0.1, θ2=1.0, ϕ2=0.1)
        u0 = zeros(4)
        prob = SDEProblem(f!, noise!, u0, (0.0, t_dur), p; seed=1234, noise_rate_prototype=nothing)
        sol = solve(prob, LambaEM())
        sol(t_dur)
    end
end

On v9.75 I get:

julia> tmtk = test_mtk()
4-element Vector{Float64}:
  0.08387167746864717
  0.0024813763968508764
 -0.048014755423273825
 -0.17317419636309844

julia> tn = test_normal()
4-element Vector{Float64}:
  0.08387167746864717
  0.0024813763968508764
 -0.048014755423273825
 -0.17317419636309844

julia> tmtk ≈ tn
true

and on v9.76 I get

julia> tmtk = test_mtk()
4-element Vector{Float64}:
 -0.06679234625972856
 -0.0907440827828222
  0.07981500480635247
 -0.01403600107650567

julia> tn = test_normal()
4-element Vector{Float64}:
  0.08387167746864717
  0.0024813763968508764
 -0.048014755423273825
 -0.17317419636309844

julia> tmtk ≈ tn
false

This only occurs for me with stochastic systems, so I suspect that something went wrong in the handling of stochastic noise, but it could also potentially have to do with stuff like the handling of random seeds or something, I'm not yet sure.

Here's the project info from the 9.76 test:

(jl_8lIddP) pkg> st
Status `/tmp/jl_8lIddP/Project.toml`
  [961ee093] ModelingToolkit v9.76.0
  [789caeaf] StochasticDiffEq v6.76.0

[ChrisRackauckas]

The noise equations need to permute the same way.

[AayushSabharwal]

It seems the noise equations are permuted correctly. Here's a side by side comparison with MTK@9.76 on the left and MTK@9.75 on the right

I'm still figuring out what the issue is

[AayushSabharwal]

This seems to be the fact that reordering equations changes how the noise is applied in the SDE:

function test_normal(t_dur = 1.0)
    function f!(du, u, p, t)
        (; θ1, ϕ1, θ2, ϕ2) = p
        (x1, y1, x2, y2) = u
        jcn1 = x2
        jcn2 = x1
        #=D(x1)=# du[1] = y1
        #=D(y1)=# du[2] = θ1*(1-x1^2)*y1 - x1 + jcn1
        #=D(x2)=# du[3] = y2
        #=D(y2)=# du[4] = θ2*(1-x2^2)*y2 - x2 + jcn2
        nothing
    end
    function noise!(du, u, p, t)
        (; ϕ1, ϕ2) = p
        du[2] = ϕ1
        du[4] = ϕ2
        nothing
    end
    let
        p = (;θ1=1.0, ϕ1=0.1, θ2=1.0, ϕ2=0.1)
        u0 = zeros(4)
        prob = SDEProblem(f!, noise!, u0, (0.0, t_dur), p; seed=1234, noise_rate_prototype=nothing)
        sol = solve(prob, LambaEM())
        sol(t_dur)
    end
end

This returns

4-element Vector{Float64}:
  0.08387167746864717
  0.0024813763968508764
 -0.048014755423273825
 -0.17317419636309844

As in the MWE. But if I reorder equations:

function test_normal2(t_dur = 1.0)
    function f!(du, u, p, t)
        (; θ1, ϕ1, θ2, ϕ2) = p
        (x1, x2, y1, y2) = u
        jcn1 = x2
        jcn2 = x1
        #=D(x1)=# du[1] = y1
        #=D(y1)=# du[3] = θ1*(1-x1^2)*y1 - x1 + jcn1
        #=D(x2)=# du[2] = y2
        #=D(y2)=# du[4] = θ2*(1-x2^2)*y2 - x2 + jcn2
        nothing
    end
    function noise!(du, u, p, t)
        (; ϕ1, ϕ2) = p
        du[3] = ϕ1
        du[4] = ϕ2
        nothing
    end
    let
        p = (;θ1=1.0, ϕ1=0.1, θ2=1.0, ϕ2=0.1)
        u0 = zeros(4)
        prob = SDEProblem(f!, noise!, u0, (0.0, t_dur), p; seed=1234, noise_rate_prototype=nothing)
        sol = solve(prob, LambaEM())
        sol(t_dur)
    end
end

This returns

4-element Vector{Float64}:
 -0.016854339441662165
 -0.055214384218212176
 -0.1174609017866221
 -0.21335087334447642

[AayushSabharwal]

These are the plots of test_mtk in MTK@9.76 and MTK@9.75 respectively. The shapes are similar, but the labeling is different which is confusing.

[ChrisRackauckas]

Oh, yes 😅 . If you reorder the equations, you'll get a different Brownian motion because the order that the random numbers are applied to the equations is different 😅. So it's no bug, but definitely a tricky case to actually test well. I'd recommend just turning off the sorting for that test? Or directly building a NoiseProcessGrid and then permuting its indices.

[MasonProtter]

I see, thanks for the info @AayushSabharwal!