SDE simulations broken

[TorkelE]

Seems SDE simulations are broken?

MWE:

using ModelingToolkit, StochasticDiffEq

# Create SDE system.
@parameters p d 
@variables t X(t)
D = Differential(t)
eqs = [D(X) ~ p - d*X]
noise_eqs = [sqrt(p), -sqrt(d*X)]
@named ssys = SDESystem(eqs, noise_eqs, t, [X], [p, d])
ssys = complete(ssys)

# Create and solves a SDEProblem.
u0 = [X => 10.0]
tspan = (0.0, 10.0)
ps = [p => 2.0, d => 0.5]
sprob = SDEProblem(ssys, u0, tspan, ps)
solve(sprob, ImplicitEM())

yields a

ERROR: MethodError: no method matching issplit(::ImplicitEM{…})

Closest candidates are:
  issplit(::OrdinaryDiffEq.OrdinaryDiffEqStabilizedIRK.IRKC)
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/lib/OrdinaryDiffEqStabilizedIRK/src/alg_utils.jl:5
  issplit(::Union{OrdinaryDiffEq.CFNLIRK3, OrdinaryDiffEq.CNAB2, OrdinaryDiffEq.CNLF2, OrdinaryDiffEq.KenCarp3, OrdinaryDiffEq.KenCarp4, OrdinaryDiffEq.KenCarp47, OrdinaryDiffEq.KenCarp5, OrdinaryDiffEq.KenCarp58, OrdinaryDiffEq.SBDF})
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/alg_utils.jl:98
  issplit(::Union{OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, OrdinaryDiffEq.DAEAlgorithm})
   @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/alg_utils.jl:97

Stacktrace:
  [1] islinearfunction(f::Function, alg::ImplicitEM{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/derivative_utils.jl:437
  [2] build_J_W(alg::ImplicitEM{…}, u::Vector{…}, uprev::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64, dt::Float64, f::SDEFunction{…}, ::Type{…}, ::Val{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/derivative_utils.jl:868
  [3] build_nlsolver(alg::ImplicitEM{…}, nlalg::OrdinaryDiffEq.NLNewton{…}, u::Vector{…}, uprev::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64, dt::Float64, f::SDEFunction{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, γ::Float64, c::Float64, α::Int64, ::Val{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/nlsolve/utils.jl:178
  [4] build_nlsolver
    @ ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/nlsolve/utils.jl:146 [inlined]
  [5] build_nlsolver
    @ ~/.julia/packages/OrdinaryDiffEq/HQ92J/src/nlsolve/utils.jl:136 [inlined]
  [6] alg_cache(alg::ImplicitEM{…}, prob::SDEProblem{…}, u::Vector{…}, ΔW::Vector{…}, ΔZ::Nothing, p::ModelingToolkit.MTKParameters{…}, rate_prototype::Vector{…}, noise_rate_prototype::Vector{…}, jump_rate_prototype::Nothing, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Vector{…}, f::SDEFunction{…}, t::Float64, dt::Float64, ::Type{…})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/aR0ZE/src/caches/sdirk_caches.jl:23
  [7] __init(_prob::SDEProblem{…}, alg::ImplicitEM{…}, 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/aR0ZE/src/solve.jl:474
  [8] __init (repeats 2 times)
    @ ~/.julia/packages/StochasticDiffEq/aR0ZE/src/solve.jl:18 [inlined]
  [9] __solve(prob::SDEProblem{…}, alg::ImplicitEM{…}, timeseries::Vector{…}, ts::Vector{…}, ks::Nothing, recompile::Type{…}; kwargs::@Kwargs{…})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/aR0ZE/src/solve.jl:6
 [10] __solve (repeats 5 times)
    @ ~/.julia/packages/StochasticDiffEq/aR0ZE/src/solve.jl:1 [inlined]
 [11] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:612 [inlined]
 [12] solve_call
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:569 [inlined]
 [13] solve_up(prob::SDEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::ImplicitEM{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:1080
 [14] solve_up
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:1066 [inlined]
 [15] #solve#51
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:1003 [inlined]
 [16] solve(prob::SDEProblem{…}, args::ImplicitEM{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:993
 [17] top-level scope
    @ ~/Desktop/Julia Playground/Environment - Temporary/playground.jl:18
Some type information was truncated. Use `show(err)` to see complete types.

on

(Environment - Temporary) pkg> st
Status `~/Desktop/Julia Playground/Environment - Temporary/Project.toml`
  [961ee093] ModelingToolkit v9.21.0
  [789caeaf] StochasticDiffEq v6.65.1

[ChrisRackauckas]

I know, it requires a hotfix. I mentioned it to @oscardssmith we need to test green ASAP.

[ChrisRackauckas]

Nothing to do with MTK so closing.

[TorkelE]

If we can get this fixed our plan was to release Catalyst v14 today

[oscardssmith]

just to clarify, is all that needs to happen here for issplit(::ImplicitEM) = false to be defined? I'm very unsure what change caused this, but it should be easy to add...

[ChrisRackauckas]

it was one of the refactors in the OrdinaryDiffEq lib splits

[ChrisRackauckas]

StochasticDiffEqAlgorithm should have a whole thing that goes to false except on KenCarpS

[TorkelE]

Presumably the change would need to be general to all see silvers as well, not just defining something for ImplicitEM

[oscardssmith]

note that fixing this bug will be difficult since it seems to be running into an out of bounds error first (that is a segfault if you don't start julia with check-bounds=true

julia> solve(sprob, ImplicitEM())
ERROR: BoundsError: attempt to access 1-element Vector{Float64} at index [2]
Stacktrace:
  [1] setindex!
    @ ./array.jl:1021 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/SymbolicUtils/c0xQb/src/code.jl:418 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/Symbolics/Eas9m/src/build_function.jl:546 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/SymbolicUtils/c0xQb/src/code.jl:375 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:163 [inlined]
  [6] macro expansion
    @ ./none:0 [inlined]
  [7] generated_callfunc
    @ ./none:0 [inlined]
  [8] RuntimeGeneratedFunction
    @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150 [inlined]
  [9] g
    @ ~/.julia/packages/ModelingToolkit/Nsxrf/src/systems/diffeqs/sdesystem.jl:410 [inlined]
 [10] sde_determine_initdt(u0::Vector{…}, t::Float64, tdir::Float64, dtmax::Float64, abstol::Float64, reltol::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), prob::SDEProblem{…}, order::Rational{…}, integrator::StochasticDiffEq.SDEIntegrator{…})
    @ StochasticDiffEq ~/.julia/dev/StochasticDiffEq/src/initdt.jl:34
 [11] auto_dt_reset!
    @ ~/.julia/dev/StochasticDiffEq/src/integrators/integrator_interface.jl:355 [inlined]
 [12] handle_dt!(integrator::StochasticDiffEq.SDEIntegrator{…})
    @ StochasticDiffEq ~/.julia/dev/StochasticDiffEq/src/solve.jl:643
 [13] __init(_prob::SDEProblem{…}, alg::ImplicitEM{…}, 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/dev/StochasticDiffEq/src/solve.jl:596
 [14] __init (repeats 2 times)
    @ ~/.julia/dev/StochasticDiffEq/src/solve.jl:18 [inlined]
 [15] __solve(prob::SDEProblem{…}, alg::ImplicitEM{…}, timeseries::Vector{…}, ts::Vector{…}, ks::Nothing, recompile::Type{…}; kwargs::@Kwargs{…})
    @ StochasticDiffEq ~/.julia/dev/StochasticDiffEq/src/solve.jl:6
 [16] __solve (repeats 5 times)
    @ ~/.julia/dev/StochasticDiffEq/src/solve.jl:1 [inlined]
 [17] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:612 [inlined]
 [18] solve_call
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:569 [inlined]
 [19] solve_up(prob::SDEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::Vector{…}, args::ImplicitEM{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:1080
 [20] solve_up
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:1066 [inlined]
 [21] #solve#51
    @ ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:1003 [inlined]
 [22] solve(prob::SDEProblem{…}, args::ImplicitEM{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/c8MAQ/src/solve.jl:993
 [23] top-level scope
    @ REPL[20]:1
Some type information was truncated. Use `show(err)` to see complete types.

[TorkelE]

This did work relatively recently though, shouldn't it be possible to identify the change which first caused the issue to appear? I'm happy to help, but not sure which updates have been made in the last week or two

[ChrisRackauckas]

Your model is not correctly formed. It's supposed to be:

using ModelingToolkit, StochasticDiffEq

# Create SDE system.
@parameters p d
@variables t X(t)
D = Differential(t)
eqs = [D(X) ~ p - d*X]
noise_eqs = reshape([sqrt(p), -sqrt(d*X)],1,2)
@named ssys = SDESystem(eqs, noise_eqs, t, [X], [p, d])
ssys = complete(ssys)

# Create and solves a SDEProblem.
u0 = [X => 10.0]
tspan = (0.0, 10.0)
ps = [p => 2.0, d => 0.5]
sprob = SDEProblem(ssys, u0, tspan, ps)
solve(sprob, ImplicitEM())

[TorkelE]

It works now, Thanks a lot! Back working on the Catalyst release, hopefully out very soon.