Using Tau-leaping

[augustinas1]

It seems that tau-leaping methods are currently broken? I might be defining something incorrectly, though the errors are a bit cryptic... Using Catalyst v16.1.1 and JumpProcesses v9.26.0 in Julia v1.12.6.

using JumpProcesses
using Catalyst
using StochasticDiffEq

sir_model = @reaction_network begin
    β, S + I --> 2I
    ν, I --> R
end

p = (:β => 0.1 / 1000, :ν => 0.01)
u₀ = [:S => 990, :I => 10, :R => 0]
tspan = (0.0, 250.0)
jump_prob = JumpProblem(sir_model, u₀, tspan, p, aggregator=PureLeaping())

sol = solve(jump_prob, TauLeaping())

ERROR: MethodError: no method matching zero(::Nothing)
The function `zero` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  zero(::Type{Union{}}, Any...)
   @ Base number.jl:315
  zero(::Type{Missing})
   @ Base missing.jl:104
  zero(::Type{Dates.Time})
   @ Dates ~/.julia/juliaup/julia-1.12.6+0.aarch64.apple.darwin14/Julia-1.12.app/Contents/Resources/julia/share/julia/stdlib/v1.12/Dates/src/types.jl:460
  ...

Stacktrace:
 [1] alg_cache(alg::TauLeaping, prob::DiscreteProblem{…}, u::Vector{…}, ΔW::Nothing, ΔZ::Nothing, p::MTKParameters{…}, rate_prototype::Vector{…}, noise_rate_prototype::Nothing, jump_rate_prototype::Nothing, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Vector{…}, f::DiscreteFunction{…}, t::Float64, dt::Float64, ::Type{…}, verbose::DEVerbosity{…})
   @ StochasticDiffEqLeaping ~/.julia/packages/StochasticDiffEqLeaping/YRM93/src/caches/tau_caches.jl:27
 [2] _sde_init(_prob::JumpProblem{…}, alg::TauLeaping; 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, save_discretes::Bool, initialize_integrator::Bool, seed::UInt64, rng::Nothing, alias::Nothing, initializealg::DefaultInit, kwargs::@Kwargs{})
   @ StochasticDiffEqCore ~/.julia/packages/StochasticDiffEqCore/qCsWr/src/solve.jl:521
 [3] __init(_prob::JumpProblem{…}, alg::TauLeaping; kwargs::@Kwargs{…})
   @ StochasticDiffEqCore ~/.julia/packages/StochasticDiffEqCore/qCsWr/src/solve.jl:97
 [4] __solve(prob::JumpProblem{…}, alg::TauLeaping; merge_callbacks::Bool, kwargs::@Kwargs{})
   @ StochasticDiffEqCore ~/.julia/packages/StochasticDiffEqCore/qCsWr/src/solve.jl:33
 [5] __solve(prob::JumpProblem{…}, alg::TauLeaping)
   @ StochasticDiffEqCore ~/.julia/packages/StochasticDiffEqCore/qCsWr/src/solve.jl:27
 [6] solve(prob::JumpProblem{…}, args::TauLeaping; kwargs::@Kwargs{})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:656
 [7] solve(prob::JumpProblem{…}, args::TauLeaping)
   @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:655
 [8] top-level scope
   @ ~/Projects/draft/src/tau_leaping_test.jl:19

I get a different error trying to use CaoTauLeaping (which should still be compatible?):

sol = solve(jump_prob, CaoTauLeaping())

ERROR: FieldError: type Nothing has no field `dW`; Nothing has no fields at all.
Stacktrace:
 [1] getproperty
   @ ./Base_compiler.jl:54 [inlined]
 [2] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::StochasticDiffEqLeaping.TauLeapingCache{…})
   @ StochasticDiffEqLeaping ~/.julia/packages/StochasticDiffEqLeaping/YRM93/src/perform_step/tau_leaping.jl:25
 [3] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
   @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/rnOL4/src/solve.jl:889
 [4] __solve(prob::JumpProblem{…}, alg::CaoTauLeaping; merge_callbacks::Bool, kwargs::@Kwargs{})
   @ StochasticDiffEqCore ~/.julia/packages/StochasticDiffEqCore/qCsWr/src/solve.jl:34
 [5] __solve
   @ ~/.julia/packages/StochasticDiffEqCore/qCsWr/src/solve.jl:27 [inlined]
 [6] #solve#31
   @ ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:656 [inlined]
 [7] solve(prob::JumpProblem{…}, args::CaoTauLeaping)
   @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:655
 [8] top-level scope
   @ ~/Projects/draft/src/tau_leaping_test.jl:20

On a different note, trying to use SSAStepper also leads to an error, which makes sense given the PureLeaping aggregator, but isn't clear from the error message:

sol = solve(jump_prob, SSAStepper())

ERROR: BoundsError: attempt to access Tuple{} at index [0]
Stacktrace:
  [1] getindex(t::Tuple, i::Int64)
    @ Base ./tuple.jl:33
  [2] __init(jump_prob::JumpProblem{…}, alg::SSAStepper; save_start::Bool, save_end::Bool, seed::Nothing, alias_jump::Bool, saveat::Nothing, callback::Nothing, tstops::Nothing, numsteps_hint::Int64)
    @ JumpProcesses ~/.julia/packages/JumpProcesses/sLhTC/src/SSA_stepper.jl:0
  [3] init_call(_prob::JumpProblem{…}, args::SSAStepper; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:87
  [4] init_call(_prob::JumpProblem{…}, args::SSAStepper)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:67
  [5] init(prob::JumpProblem{…}, args::SSAStepper; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:106
  [6] init(prob::JumpProblem{…}, args::SSAStepper)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:105
  [7] __solve(jump_prob::JumpProblem{…}, alg::SSAStepper; kwargs::@Kwargs{})
    @ JumpProcesses ~/.julia/packages/JumpProcesses/sLhTC/src/SSA_stepper.jl:131
  [8] __solve(jump_prob::JumpProblem{…}, alg::SSAStepper)
    @ JumpProcesses ~/.julia/packages/JumpProcesses/sLhTC/src/SSA_stepper.jl:129
  [9] solve(prob::JumpProblem{…}, args::SSAStepper; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:656
 [10] solve(prob::JumpProblem{…}, args::SSAStepper)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/CJTzr/src/solve.jl:655
 [11] top-level scope
    @ ~/Projects/draft/src/tau_leaping_test.jl:21

[ChrisRackauckas]

Building for TauLeaping needs to build to regular jumps. I'm not sure Catalyst targets that, it would need to be an option in the JumpProblem construction since it fundamentally needs a different representation of the jumps as a function.

[isaacsas]

This would need to be supported at the ModelingToolkit level as we need codegen for RegularJumps.

Tau-leaping is on the TODO list for the summer, but realistically there is a lot of work still needed on the tau-leaping solvers before I would really recommend them for day to day use (i.e. good adaptive time-stepping, correctness validation, and some attention to non-negativity preservation).

[augustinas1]

I see, thanks both! I guess I'll leave this open as a feature request then

[ChrisRackauckas]

Yes please do. I'm hoping it does get that summer lovin'