Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Simulating jump problems require non-integer tspan values #359

Open
TorkelE opened this issue Oct 1, 2023 · 3 comments
Open

Simulating jump problems require non-integer tspan values #359

TorkelE opened this issue Oct 1, 2023 · 3 comments

Comments

@TorkelE
Copy link
Member

TorkelE commented Oct 1, 2023

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
Copy link
Member

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

@isaacsas
Copy link
Member

isaacsas commented Oct 2, 2023

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
Copy link
Member

I guess it needs a bit more modifications

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants