Solver problem with not-in-place non-diagonal noise SDEProblem

[alvestad10]

I ran into a problem today with the use of non-diagonal noise and not-in-place execution of the drift term and noise-coefficient. The problem seems to be on a subset of solvers including ImplicitEM(), SKenCarp() and SOSRA(), but working on solvers like EM(), EulerHeun() and RKMilGeneral().

A simplified code containing the problem

using StochasticDiffEq

function a_func(u,p,t)
    du1 = -(p[1] * u[1] - p[2] * u[2])
    du2 = -(p[1] * u[2] + p[2] * u[1])
    return [du1,du2]
end

function b_func(u,p,t)
    du11 = sqrt(2)*p[1] ; du12 = 0
    du21 = sqrt(2)*p[2] ; du22 = 0
    return [du11 du12
            du21 du22]
end

Ks = [1.,0.]

prob = SDEProblem(a_func,b_func,[0.,0.],(0.0,0.02),Ks,
                    noise_rate_prototype=similar(Ks,2,2))
sol = solve(prob,ImplicitEM(theta=0.5), #EM()
            dt=1e-2, saveat=0.01)

with stacktrace

ERROR: LoadError: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(2), Base.OneTo(2)), must have singleton at dim 2")
Stacktrace:
  [1] promote_shape
    @ ./indices.jl:183 [inlined]
  [2] promote_shape(a::Matrix{Float64}, b::Vector{Float64})
    @ Base ./indices.jl:169
  [3] +(A::Matrix{Float64}, Bs::Vector{Float64})
    @ Base ./arraymath.jl:45
  [4] muladd
    @ ./math.jl:1152 [inlined]
  [5] perform_step!(integrator::StochasticDiffEq.SDEIntegrator{ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, Val{:central}, Float64, :Predictive}, false, Vector{Float64}, Float64, Float64, Float64, Vector{Float64}, Float64, Float64, Float64, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, Nothing, SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(b_func), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Matrix{Float64}}, ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, Val{:central}, Float64, :Predictive}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats}, StochasticDiffEq.ImplicitEMConstantCache{OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, false, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonConstantCache{Float64, Float64, Matrix{Float64}, LU{Float64, Matrix{Float64}}, SciMLBase.UDerivativeWrapper{SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Float64}}}}}, SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(b_func), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, OrdinaryDiffEq.PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Nothing, 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, Int64, Float64, Float64, Float64, Tuple{}, Float64, Tuple{}}, Nothing, Float64, Nothing, Nothing}, cache::StochasticDiffEq.ImplicitEMConstantCache{OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, false, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonConstantCache{Float64, Float64, Matrix{Float64}, LU{Float64, Matrix{Float64}}, SciMLBase.UDerivativeWrapper{SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Float64}}}}})
    @ StochasticDiffEq ~/Developer/Stochastic/test_can_be_deleted/dev/StochasticDiffEq/src/perform_step/sdirk.jl:30
  [6] solve!(integrator::StochasticDiffEq.SDEIntegrator{ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, Val{:central}, Float64, :Predictive}, false, Vector{Float64}, Float64, Float64, Float64, Vector{Float64}, Float64, Float64, Float64, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.WHITE_NOISE_BRIDGE), false, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, false}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, Nothing, SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(b_func), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Matrix{Float64}}, ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, Val{:central}, Float64, :Predictive}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats}, StochasticDiffEq.ImplicitEMConstantCache{OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, false, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonConstantCache{Float64, Float64, Matrix{Float64}, LU{Float64, Matrix{Float64}}, SciMLBase.UDerivativeWrapper{SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Float64, Vector{Float64}}}}}, SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(b_func), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, OrdinaryDiffEq.PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Nothing, 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, Int64, Float64, Float64, Float64, Tuple{}, Float64, Tuple{}}, Nothing, Float64, Nothing, Nothing})
    @ StochasticDiffEq ~/Developer/Stochastic/test_can_be_deleted/dev/StochasticDiffEq/src/solve.jl:607
  [7] #__solve#100
    @ ~/Developer/Stochastic/test_can_be_deleted/dev/StochasticDiffEq/src/solve.jl:7 [inlined]
  [8] #solve_call#42
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:61 [inlined]
  [9] solve_up(prob::SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, Nothing, SDEFunction{false, typeof(a_func), typeof(b_func), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(b_func), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Matrix{Float64}}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::ImplicitEM{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, Val{:central}, Float64, :Predictive}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:dt, :saveat), Tuple{Float64, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:87
 [10] #solve#43
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:73 [inlined]
 [11] top-level scope
    @ ~/Developer/Stochastic/test_can_be_deleted/test.jl:26

[ChrisRackauckas]

Thanks for the report. Yeah these are missing a handling for the non-diagonal derivative approximation.