search

SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
141 stars 35 forks source link

Simulating jump problems require non-integer tspan values #359

Open TorkelE opened 1 year ago

TorkelE commented 1 year ago

Simple example:

using Catalyst, JumpProcesses

rn = @reaction_network begin
    (p,d), 0 <--> X
end
u0 = [:X => 1.0]
tspan = (0.0, 10.0)
p = [:p => 1.0, :d => 0.1]

dprob = DiscreteProblem(rn, u0, tspan, p)
jprob = JumpProblem(rn, dprob, Direct())
sol = solve(jprob, SSAStepper())

If I change to tspan = (0, 10):

using Catalyst, JumpProcesses

rn = @reaction_network begin
    (p,d), 0 <--> X
end
u0 = [:X => 1.0]
tspan = (0, 10)
p = [:p => 1.0, :d => 0.1]

dprob = DiscreteProblem(rn, u0, tspan, p)
jprob = JumpProblem(rn, dprob, Direct())
sol = solve(jprob, SSAStepper())

I instead get an error:

ERROR: InexactError: Int64(1.1)
Stacktrace:

Int64 at [float.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

convert at [number.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

setindex! at [array.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

time_to_next_jump(p::JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, u::Vector{Float64}, params::Vector{Float64}, t::Int64) at [direct.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

generate_jumps!(p::JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, integrator::JumpProcesses.SSAIntegrator{DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#213#214", Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#685"{Bool, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}, Dict{Any, Any}}, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}}, Vector{Float64}, Int64, Int64, Vector{Float64}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Int64}, Nothing, DiscreteProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#213#214", Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#685"{Bool, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}, Dict{Any, Any}}, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, SSAStepper, SciMLBase.ConstantInterpolation{Vector{Int64}, Vector{Vector{Float64}}}, DiffEqBase.Stats, Nothing}, DiscreteCallback{JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, typeof(SciMLBase.FINALIZE_DEFAULT)}, Nothing, NamedTuple{(:callback,), Tuple{CallbackSet{Tuple{}, Tuple{}}}}, Vector{Int64}}, u::Vector{Float64}, params::Vector{Float64}, t::Int64) at [direct.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

initialize! at [direct.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

AbstractSSAJumpAggregator at [ssajump.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

__init(jump_prob::JumpProblem{true, DiscreteProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#213#214", Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#685"{Bool, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}, Dict{Any, Any}}, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Direct, CallbackSet{Tuple{}, Tuple{DiscreteCallback{JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, typeof(SciMLBase.FINALIZE_DEFAULT)}}}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, Tuple{}, Nothing, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Random.TaskLocalRNG, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}}, alg::SSAStepper; save_start::Bool, save_end::Bool, seed::Nothing, alias_jump::Bool, saveat::Nothing, callback::Nothing, tstops::Vector{Int64}, numsteps_hint::Int64) at [SSA_stepper.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

__init at [SSA_stepper.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

#init_call#30 at [solve.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

init_call at [solve.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

#init#32 at [solve.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

init at [solve.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

#__solve#221 at [SSA_stepper.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

__solve at [SSA_stepper.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

#solve#47 at [solve.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

solve(prob::JumpProblem{true, DiscreteProblem{Vector{Float64}, Tuple{Int64, Int64}, true, Vector{Float64}, DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#213#214", Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#685"{Bool, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}, Dict{Any, Any}}, JumpSystem{RecursiveArrayTools.ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Direct, CallbackSet{Tuple{}, Tuple{DiscreteCallback{JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, typeof(SciMLBase.FINALIZE_DEFAULT)}}}, JumpProcesses.DirectJumpAggregation{Int64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Tuple{}, Tuple{}, Random.TaskLocalRNG}, Tuple{}, Nothing, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Float64}}, Random.TaskLocalRNG, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}}, args::SSAStepper) at [solve.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

top-level scope at [playground.jl](vscode-file://vscode-app/snap/code/140/usr/share/code/resources/app/out/vs/code/electron-sandbox/workbench/workbench.html)

The problem also happens when Catalyst is not used, e.g. this also give a similar error

begin
    using JumpProcesses
    # here we order S = 1, I = 2, and R = 3
    # substrate stoichiometry:
    substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
        [2 => 1]]            # 1*I
    # net change by each jump type
    netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
        [2 => -1, 3 => 1]]   # I -> I-1, R -> R+1
    # rate constants for each jump
    p = (0.1 / 1000, 0.01)

    # p[1] is rate for S+I --> 2I, p[2] for I --> R
    pidxs = [1, 2]

    maj = MassActionJump(substoich, netstoich; param_idxs = pidxs)

    u₀ = [999, 1, 0]       #[S(0),I(0),R(0)]
    tspan = (0, 250)
    dprob = DiscreteProblem(u₀, tspan, p)

    # use the Direct method to simulate
    jprob = JumpProblem(dprob, Direct(), maj)

    # solve as a pure jump process, i.e. using SSAStepper
    sol = solve(jprob, SSAStepper())
end

This works fine with e.g. ODEs:

rn = @reaction_network begin
    (p,d), 0 <--> X
end
u0 = [:X => 1.0]
tspan = (0, 10)
p = [:p => 1.0, :d => 0.1]

oprob = ODEProblem(rn, u0, tspan, p)
sol = solve(oprob, Tsit5())
ChrisRackauckas commented 1 year ago

Yes JumpProcesses.jl should adapt the conversions done in DiffEqBase.

isaacsas commented 1 year ago

What conversion should we use? On quick glance it looks like the type of dt is used in DiffEqBase, but there is no dt here.

ChrisRackauckas commented 1 year ago

I guess it needs a bit more modifications