search

SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.38k stars 196 forks source link

SDE simulations broken #2829

Closed TorkelE closed 4 days ago

TorkelE commented 4 days ago

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 commented 4 days ago

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

ChrisRackauckas commented 4 days ago

Nothing to do with MTK so closing.

TorkelE commented 4 days ago

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

oscardssmith commented 4 days ago

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 commented 4 days ago

it was one of the refactors in the OrdinaryDiffEq lib splits

ChrisRackauckas commented 4 days ago

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

TorkelE commented 4 days ago

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

oscardssmith commented 2 days ago

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 commented 2 days ago

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 commented 2 days ago

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 commented 2 days ago

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