Mutation attempt in `ContinuousCallback` when solving `SDEProblem{StaticArraysCore.SVector{2, Float64}}`

Specific Domains Modelling & Simulations
, ,
1

I have an SDE where the drift function f is defined out-of-place and the state is a StaticArrays.SVector. I need to terminate! the integration using a ContinuousCallback. This raises an error, as the StochasticDiffEq package calls the mutating functionsde_interpolant!.

This error does not occur when solving the equivalent ODEProblem with the same callback or if I solve the problem in-place by using an StaticArrays.MVector

Before raising an issue, I wanted to double check that I was not doing something stupid. In the meantime, I will move to an in-place-formulation

Here a MWE.

import DifferentialEquations as DE
import StaticArrays as SA

function f(u::V, p, t) where V <: SA.StaticVector
    V(u[2], -p)
end
function g(::V, p, t) where V <: SA.StaticVector
    V(1., 0.)
end

function condition(u, t, integrator)
    u[1]
end

cb = DE.ContinuousCallback(condition, DE.terminate!)

u0 = SA.SVector(1., 0.)
tspan = (0.0, 1.)
p = 4.

## SDE out-of-place
fn = DE.SDEFunction{false}(f, g)
sdeprob = DE.SDEProblem(fn, u0, tspan, p)

DE.solve(sdeprob, DE.SOSRI()) # retcode: Success
DE.solve(sdeprob, DE.SOSRI(); callback = cb) # setindex!(::StaticArraysCore.SVector{2, Float64}, value, ::Int) is not defined


## ODE out-of-place
odeprob = DE.ODEProblem(f, u0, tspan, p)

DE.solve(odeprob, DE.Tsit5()) # retcode: Success
DE.solve(odeprob, DE.Tsit5(); callback = cb) # retcode: Terminated

## SDE in-place
function f!(du, u, p, t)
    du[1] = u[2]
    du[2] = -p
end
function g!(Σ, u, p, t)
    Σ[1] = 1.
    Σ[2] = 0.
end

u0ip = SA.MVector(u0)
fnip = DE.SDEFunction{true}(f!, g!)
sdeiipprob = DE.SDEProblem(fnip, u0ip, tspan, p)

DE.solve(sdeiipprob, DE.SOSRI()) # retcode: Success
DE.solve(sdeiipprob, DE.SOSRI(); callback = cb) # retcode: Terminated

The erroring call generates the following stacktrace.

ERROR: setindex!(::StaticArraysCore.SVector{2, Float64}, value, ::Int) is not defined.
 Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:44
  [2] setindex!(a::StaticArraysCore.SVector{2, Float64}, value::Float64, i::Int64)
    @ StaticArrays C:\Users\user\.julia\packages\StaticArrays\mI3IC\src\indexing.jl:3
  [...]
 [10] sde_interpolant!(out::StaticArraysCore.SVector{…}, Θ::Float64, dt::Float64, u0::StaticArraysCore.SVector{…}, u1::StaticArraysCore.SVector{…}, idxs::Nothing, deriv::Type{…})    
    @ StochasticDiffEq C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\dense.jl:46
 [11] sde_interpolant!(out::StaticArraysCore.SVector{…}, Θ::Float64, integrator::StochasticDiffEq.SDEIntegrator{…}, idxs::Nothing, deriv::Type)
    @ StochasticDiffEq C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\dense.jl:20
 [12] current_interpolant!
    @ C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\dense.jl:72 [inlined]
 [13] (::StochasticDiffEq.SDEIntegrator{…})(val::StaticArraysCore.SVector{…}, t::Float64, deriv::Type; idxs::Nothing)
    @ StochasticDiffEq C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\integrators\integrator_interface.jl:37
 [14] (::StochasticDiffEq.SDEIntegrator{…})(val::StaticArraysCore.SVector{…}, t::Float64, deriv::Type)
    @ StochasticDiffEq C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\integrators\integrator_interface.jl:33
 [15] change_t_via_interpolation!
    @ C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\integrators\integrator_interface.jl:12 [inlined]
 [16] apply_callback!(integrator::StochasticDiffEq.SDEIntegrator{…}, callback::SciMLBase.ContinuousCallback{…}, cb_time::Float64, prev_sign::Float64, event_idx::Float64)
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\5LeiG\src\callbacks.jl:462
 [17] handle_callbacks!
    @ C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\integrators\integrator_utils.jl:307 [inlined]
 [18] loopfooter!
    @ C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\integrators\integrator_utils.jl:197 [inlined]
 [19] solve!(integrator::StochasticDiffEq.SDEIntegrator{…})
    @ StochasticDiffEq C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\solve.jl:788
 [20] #__solve#74
    @ C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\solve.jl:7 [inlined]
 [21] __solve
    @ C:\Users\user\.julia\packages\StochasticDiffEq\rq8zz\src\solve.jl:1 [inlined]
 [22] #solve_call#22
    @ C:\Users\user\.julia\packages\DiffEqBase\5LeiG\src\solve.jl:172 [inlined]  
 [23] solve_call
    @ C:\Users\user\.julia\packages\DiffEqBase\5LeiG\src\solve.jl:137 [inlined]  
 [24] solve_up(prob::SciMLBase.SDEProblem{…}, sensealg::Nothing, u0::StaticArraysCore.SVector{…}, p::Float64, args::StochasticDiffEq.SOSRI; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\5LeiG\src\solve.jl:630 
 [25] solve_up
    @ C:\Users\user\.julia\packages\DiffEqBase\5LeiG\src\solve.jl:603 [inlined]
 [26] #solve#28
    @ C:\Users\user\.julia\packages\DiffEqBase\5LeiG\src\solve.jl:587 [inlined]