Out-of-place CoupledSDEs erroring with trajectory

[reykboerner]

Describe the bug
Using trajectory on a CoupledSDEs with specified covariance works for in-place deterministic functions but not for out-of-place. The code below throws the following error:

DimensionMismatch: Sizes (Size(2,), Size(2, 2)) of input arrays do not match

Minimal Working Example

using DynamicalSystemsBase
using StochasticDiffEq

ρ = 0.3
Σ = [1 ρ; ρ 1]
f!(du, u, p, t) = du .= 1.01u
f(u, p, t) = 1.01u
sde_iip = CoupledSDEs(f!, zeros(2); covariance=Σ)
sde_oop = CoupledSDEs(f, zeros(2); covariance=Σ)

trajectory(sde_iip, 1) # works
trajectory(sde_oop, 1) # throws error

Package versions

DynamicalSystemsBase v3.11.2
StochasticDiffEq v6.70.0
Julia v1.10.0

[Datseris]

always paste the error please.

Regardless, this should have been caught by the test suite. Tests for SDEs need to be increased.

[reykboerner]

Full error:

ERROR: LoadError: DimensionMismatch: Sizes (Size(2,), Size(2, 2)) of input arrays do not match
Stacktrace:
  [1] _throw_size_mismatch(::SVector{2, Float64}, ::Vararg{Any})
    @ StaticArrays ~/.julia/packages/StaticArrays/9Yt0H/src/traits.jl:117
  [2] same_size
    @ ~/.julia/packages/StaticArrays/9Yt0H/src/traits.jl:111 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/StaticArrays/9Yt0H/src/mapreduce.jl:76 [inlined]
  [4] _map
    @ ~/.julia/packages/StaticArrays/9Yt0H/src/mapreduce.jl:42 [inlined]
  [5] map
    @ ~/.julia/packages/StaticArrays/9Yt0H/src/mapreduce.jl:39 [inlined]
  [6] +(a::SVector{2, Float64}, b::StaticArraysCore.SizedMatrix{2, 2, Float64, 2, Matrix{Float64}})
    @ StaticArrays ~/.julia/packages/StaticArrays/9Yt0H/src/linalg.jl:12
  [7] perform_step!(integrator::StochasticDiffEq.SDEIntegrator{SOSRA, false, SVector{2, Float64}, Float64, Float64, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, SVector{2, Float64}, RODESolution{Float64, 2, Vector{SVector{2, Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{SVector{2, Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, @Kwargs{}, Matrix{Float64}}, SOSRA, StochasticDiffEq.LinearInterpolationData{Vector{SVector{2, Float64}}, Vector{Float64}}, SciMLBase.DEStats, Nothing, Nothing}, StochasticDiffEq.ThreeStageSRAConstantCache{Float64, Float64}, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, OrdinaryDiffEqCore.PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Bool, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, Nothing, Float64, Nothing, Nothing}, cache::StochasticDiffEq.ThreeStageSRAConstantCache{Float64, Float64})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/5lNtE/src/perform_step/sra.jl:181
  [8] step!
    @ ~/.julia/packages/StochasticDiffEq/5lNtE/src/iterator_interface.jl:11 [inlined]
  [9] step!(::CoupledSDEs{false, 2, StochasticDiffEq.SDEIntegrator{SOSRA, false, SVector{2, Float64}, Float64, Float64, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, SVector{2, Float64}, RODESolution{Float64, 2, Vector{SVector{2, Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{SVector{2, Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, @Kwargs{}, Matrix{Float64}}, SOSRA, StochasticDiffEq.LinearInterpolationData{Vector{SVector{2, Float64}}, Vector{Float64}}, SciMLBase.DEStats, Nothing, Nothing}, StochasticDiffEq.ThreeStageSRAConstantCache{Float64, Float64}, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, OrdinaryDiffEqCore.PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Bool, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, Nothing, Float64, Nothing, Nothing}, SciMLBase.NullParameters})
    @ StochasticSystemsBase ~/.julia/packages/DynamicalSystemsBase/r9Ln5/ext/src/CoupledSDEs.jl:171
 [10] trajectory_continuous(ds::CoupledSDEs{false, 2, StochasticDiffEq.SDEIntegrator{SOSRA, false, SVector{2, Float64}, Float64, Float64, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, SVector{2, Float64}, RODESolution{Float64, 2, Vector{SVector{2, Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{SVector{2, Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, @Kwargs{}, Matrix{Float64}}, SOSRA, StochasticDiffEq.LinearInterpolationData{Vector{SVector{2, Float64}}, Vector{Float64}}, SciMLBase.DEStats, Nothing, Nothing}, StochasticDiffEq.ThreeStageSRAConstantCache{Float64, Float64}, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, OrdinaryDiffEqCore.PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Bool, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, Nothing, Float64, Nothing, Nothing}, SciMLBase.NullParameters}, T::Int64; Dt::Float64, Δt::Float64, Ttr::Float64, accessor::Nothing, kw::@Kwargs{})
    @ DynamicalSystemsBase ~/.julia/packages/DynamicalSystemsBase/r9Ln5/src/core/trajectory.jl:85
 [11] trajectory_continuous
    @ ~/.julia/packages/DynamicalSystemsBase/r9Ln5/src/core/trajectory.jl:75 [inlined]
 [12] #trajectory#5
    @ ~/.julia/packages/DynamicalSystemsBase/r9Ln5/src/core/trajectory.jl:44 [inlined]
 [13] trajectory
    @ ~/.julia/packages/DynamicalSystemsBase/r9Ln5/src/core/trajectory.jl:36 [inlined]
 [14] trajectory(ds::CoupledSDEs{false, 2, StochasticDiffEq.SDEIntegrator{SOSRA, false, SVector{2, Float64}, Float64, Float64, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, SVector{2, Float64}, RODESolution{Float64, 2, Vector{SVector{2, Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), Nothing, false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Vector{Float64}}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{SVector{2, Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, Nothing, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, @Kwargs{}, Matrix{Float64}}, SOSRA, StochasticDiffEq.LinearInterpolationData{Vector{SVector{2, Float64}}, Vector{Float64}}, SciMLBase.DEStats, Nothing, Nothing}, StochasticDiffEq.ThreeStageSRAConstantCache{Float64, Float64}, SDEFunction{false, SciMLBase.FullSpecialize, typeof(f), StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, StochasticSystemsBase.var"#8#12"{Float64, Matrix{Float64}}, Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, OrdinaryDiffEqCore.PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Bool, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, Nothing, Float64, Nothing, Nothing}, SciMLBase.NullParameters}, T::Int64)
    @ DynamicalSystemsBase ~/.julia/packages/DynamicalSystemsBase/r9Ln5/src/core/trajectory.jl:36
 [15] macro expansion
    @ show.jl:1229 [inlined]
 [16] top-level scope
    @ ~/Library/CloudStorage/OneDrive-UniversityofReading/Documents/phd/research/software/bugtest.jl:17

[reykboerner]

The issue disappears when changing the SDE solver to LambdaEM or EM. The following change to the MWE above runs as expected:

diffeq = (alg = LambaEM(), dt=0.1)
sde_oop = CoupledSDEs(f, zeros(2); covariance=Σ, diffeq=diffeq)

So it seems related to the default SOSRA solver - @oameye do you have an idea?

[oameye]

SOSRA() only works for diagonal noise. See SciML/StochasticDiffEq.jl#578. I was planning to make a warning on SciML side of things. I started that with this PR. However, I still need to add the traits in StochasticDiffEq.jl.

[reykboerner]

Okay thanks! Should we add a note in the docs or a warning in CoupledSDEs until that is resolved?

[Datseris]

I vote no: all efforts should go so that the warning occurs only once and only at the correct source.

[oameye]

It should be solved now with SciML/StochasticDiffEq.jl#588.