Catalyst.jl icon indicating copy to clipboard operation
Catalyst.jl copied to clipboard

Published 20 hours ago • verified publisher icon SciML

VariableRateJump Problems

Open arnold-c opened this issue 2 years ago • 12 comments

Copying this question/feature request from the discourse discussion here.

I'm trying to code up a JumpProcess of an SIR-style model that includes a VariableRateJump e.g. the transmission parameter varies throughout the year in a sinusoidal fashion. I have what I think is working code manually defining the rates and affects in JumpProcesses.jl directly, but I'm struggling to translate this to Catalyst.jl. Is this something that can be done in Catalyst.jl?

I am aware of this issue on the GitHub repo for MTK.jl, but couldn't translate this to my problem in either MTK or Catalyst.

JumpProcesses.jl Code 

#%%
using JumpProcesses, DifferentialEquations

#%%
u₀ = [999, 1, 0]
δt = 0.005
tlower = 0.0
tmax = 250.0
tspan = (tlower, tmax)
tlength = length(tlower:δt:tmax)

dur_inf = 8
R₀ = 10.0
γ = 1 / dur_inf
μ = 1 / (62 * 365)
β = 0.001250441891294741
# Adjust the scale of the seasonal variation in infectivity i.e. A_scale scales the amplitude of cosine function to 1/A_scale. The maximum infectivity is β
A_scale = 2.0
p = (β, γ, μ, A_scale)

#%%
birth_rate(u, p, t) = p[3] * (u[1] + u[2] + u[3])  # μ*N
birth_affect!(integrator) = integrator.u[1] += 1  # S -> S + 1
birth_jump = ConstantRateJump(birth_rate, birth_affect!)

recov_rate(u, p, t) = p[2] * u[2]         # γ*I
function recov_affect!(integrator)
    integrator.u[2] -= 1
    integrator.u[3] += 1
    return nothing
end
recov_jump = ConstantRateJump(recov_rate, recov_affect!)

S_death_rate(u, p, t) = p[3] * u[1]  # μ*S
S_death_affect!(integrator) = integrator.u[1] -= 1  # S -> S + 1
S_death_jump = ConstantRateJump(S_death_rate, S_death_affect!)

I_death_rate(u, p, t) = p[3] * u[2]  # μ*I
I_death_affect!(integrator) = integrator.u[2] -= 1  # S -> S + 1
I_death_jump = ConstantRateJump(I_death_rate, I_death_affect!)

R_death_rate(u, p, t) = p[3] * u[3]  # μ*R
R_death_affect!(integrator) = integrator.u[3] -= 1  # S -> S + 1
R_death_jump = ConstantRateJump(R_death_rate, R_death_affect!)

# Place infection at the end as it is a VariableRateJump, which is ordered after ConstantRateJumps in the dependency graph.
# Amplitude = 0.5 * (cos(2pi * t / 365) + 1)) * 1/scale + (scale - 1)/scale
function infec_rate(u, p, t)
    # β*S*I * 0.5*cos(2πt/365) + 0.5
    Amplitude = 0.5 * (cos(2π * t / 365) + 1) * 1 / p[4] + (p[4] - 1) / p[4]
    return p[1] * u[1] * u[2] * Amplitude
end

# Lower bound on infection rate
lrate(u, p, t) = 0.0
# Upper bound on infection rate
urate(u, p, t) = p[1] * u[1] * u[2]  # β*S*I
# Interval that the infection rate must be within the bounds for
rateinterval(u, p, t) = 1 / (2 * urate(u, p, t))

function infec_affect!(integrator)
    integrator.u[1] -= 1     # S -> S - 1
    integrator.u[2] += 1     # I -> I + 1
    return nothing
end
infec_jump = VariableRateJump(infec_rate, infec_affect!; lrate, urate, rateinterval)

jumps = JumpSet(
    birth_jump, recov_jump, S_death_jump, I_death_jump, R_death_jump, infec_jump
)

# ConstantRateJumps are ordered before VariableRateJumps in the dependency graph, otherwise within the same ordering presented to the JumpProblem, i.e., birth, recovery ...
# Each row of the dependency graph is a list of rates that must be recalculated when the row's jump is executed, as it affects the underlying states each of the rates depends on.
dep_graph = [
    [1, 3, 6],          # Birth, S death, infection
    [2, 1, 4, 5, 6],    # Recovery, birth, I death, R death, infection
    [3, 1, 6],          # S death, birth, infection
    [4, 1, 2, 6],       # I death, birth, recovery, infection
    [5, 1, 2, 6],       # R death, birth, recovery, infection
    [6, 1, 2, 3, 4],    # Infection, birth, recovery, S death, I death
]

#%%
season_infec_prob = DiscreteProblem(u₀, tspan, p)
# Bounded VariableJumpRate problems require the Coevolve() algorithm
season_infec_jump_prob = JumpProblem(season_infec_prob, Coevolve(), jumps; dep_graph)
season_infec_sol = solve(season_infec_jump_prob, SSAStepper())
Core Catalyst.jl Code 

#%%
using DifferentialEquations, Catalyst, JumpProcesses

#%%
δt = 0.005
tlower = 0.0
tmax = 250.0
tspan = (tlower, tmax)
tlength = length(tlower:δt:tmax)

#%%
R₀ = 10.0
dur_inf = 8
recov = 1 / dur_inf
birth = 1 / (62 * 365)
S_init = 999
I_init = 1
R_init = 0

u₀ = [:S => S_init, :I => I_init, :R => R_init]

beta = 0.001250441891294741
p = [:γ => recov, :μ => birth, :β => beta]

#%%
@parameters β γ μ
@variables t
@species S(t) I(t) R(t)

#%%
sir_model_rxs = [
    Reaction(β * 0.5 * (cos(2π * t/365) + 1), [S, I], [I], [1, 1], [2]),
    Reaction(γ, [I], [R], [1], [1]),
    Reaction(μ * (S + I + R), nothing, [S]),
    Reaction(μ, [S], nothing),
    Reaction(μ, [I], nothing),
    Reaction(μ, [R], nothing),

@named sir_model = ReactionSystem(sir_model_rxs, t)

Trying to use a DiscreteProblem with Integer species produces the following error.

julia> de_prob = DiscreteProblem(sir_model, u₀, tspan, p)
DiscreteProblem with uType Vector{Int64} and tType Float64. In-place: true
timespan: (0.0, 146000.0)
u0: 3-element Vector{Int64}:
 999
   1
   0

julia> jump_prob = JumpProblem(sir_model, de_prob, Coevolve(); savepositions = (false, false))
ERROR: Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] JumpProblem(js::JumpSystem{ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}}}}, prob::DiscreteProblem{Vector{Int64}, Tuple{Float64, Float64}, true, Vector{Float64}, DiscreteFunction{true,true,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, aggregator::Coevolve; callback::Nothing, kwargs::Base.Pairs{Symbol, Tuple{Bool, Bool}, Tuple{Symbol}, NamedTuple{(:savepositions,), Tuple{Tuple{Bool, Bool}}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/3ePwL/src/systems/jumps/jumpsystem.jl:401
 [3] JumpProblem
   @ ~/.julia/packages/ModelingToolkit/3ePwL/src/systems/jumps/jumpsystem.jl:388 [inlined]
 [4] #JumpProblem#104
   @ ~/.julia/packages/Catalyst/tWzC3/src/reactionsystem.jl:1528 [inlined]
 [5] top-level scope
   @ ~/Documents/Repos/outbreak-detection/scripts/Seasonal-infectivity_catalyst.jl:69

Using ODEProblem with the Coevolve() and Tsit5() with Integer/Float species (Int gets converted to Float by ODEProblem)

julia> de_prob = ODEProblem(sir_model, u₀, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 146000.0)
u0: 3-element Vector{Float64}:
 999.0
   1.0
   0.0

julia> jump_prob = JumpProblem(sir_model, de_prob, Coevolve(); savepositions = (false, false))
JumpProblem with problem ODEProblem with aggregator Coevolve
Number of jumps with discrete aggregation: 0
Number of jumps with continuous aggregation: 6
Number of mass action jumps: 0


julia> jump_sol = solve(jump_prob, Tsit5(); saveat = 1.0)
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.


Closest candidates are:
  convert(::Type{T}, ::Factorization) where T<:AbstractArray
   @ LinearAlgebra ~/.julia/juliaup/julia-1.9.0-rc2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:59
  convert(::Type{T}, ::T) where T<:AbstractArray
   @ Base abstractarray.jl:16
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:64
  ...

Stacktrace:
  [1] push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, item::Vector{Float64})
    @ Base ./array.jl:1060
  [2] copyat_or_push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, i::Int64, x::Vector{Float64}, perform_copy::Bool)
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/VzH8Y/src/utils.jl:184
  [3] _savevalues!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…}, force_save::Bool, reduce_size::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:75
  [4] savevalues! (repeats 2 times)
    @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:57 [inlined]
  [5] handle_callbacks!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:313
  [6] _loopfooter!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:242
  [7] loopfooter!
    @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:198 [inlined]
  [8] solve!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/solve.jl:521
  [9] __solve(jump_prob::JumpProblem{true, ODEProblem{true,ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Tuple{Float64, Float64},…}, Coevolve, CallbackSet{Tuple{ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}}, Tuple{}}, Nothing, Tuple{VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd90a62cf, 0xc0f2dfd6, 0xe0ae39ad, 0x176a34a7, 0x4aae5eb2)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2547"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8590323b, 0xb9b6a864, 0x6ba7d6d1, 0xd7fb926a, 0xc0cd316e)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x91c2ce70, 0x8401479e, 0xa981a223, 0x44ab2501, 0x7c50d39b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2548"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa109d73d, 0x5d51bdc8, 0x69796ab3, 0x17b73ace, 0x4f9afb0a)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf8e9e916, 0xe9725a58, 0xe0e2ea2f, 0x9e4c5782, 0x28633e09)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2549"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7694252c, 0x9f93d671, 0x08f9b3bc, 0xa4545f7a, 0x4f1f9c84)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6fa6e9c5, 0x6ac1adbf, 0xe7751301, 0xa9d78f27, 0x6215a25b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2550"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7d198a23, 0xef09683d, 0xb853195d, 0xd10ce828, 0x6a83efa4)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x62e08576, 0x38be2b6d, 0x98ba0ae7, 0xd9beffa6, 0xadc6c122)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2551"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb217e576, 0xbac21043, 0xfe471211, 0x85c08ef7, 0xc43f99d7)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3728bba8, 0x2104faac, 0x11bbe6d2, 0x921c1728, 0x77dc45a4)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2552"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8d262b03, 0x22932961, 0x6ac711c8, 0x26d27bf1, 0x4369e7ff)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}}, Nothing, Nothing, Random.TaskLocalRNG, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}}, alg::Tsit5{Static.False,…}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ JumpProcesses ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:6
 [10] __solve (repeats 5 times)
    @ ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:1 [inlined]
 [11] #solve#33
    @ ~/.julia/packages/DiffEqBase/qPmC2/src/solve.jl:965 [inlined]
 [12] top-level scope
    @ ~/Documents/Repos/outbreak-detection/scripts/Seasonal-infectivity_catalyst.jl:72

Using ODEProblem with Direct() and Tsit5() with Integer/Float species (Int gets converted to Float by ODEProblem)

julia> de_prob = ODEProblem(sir_model, u₀, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 146000.0)
u0: 3-element Vector{Float64}:
 999.0
   1.0
   0.0

julia> jump_prob = JumpProblem(sir_model, de_prob, Direct(); savepositions = (false, false))
JumpProblem with problem ODEProblem with aggregator Direct
Number of jumps with discrete aggregation: 0
Number of jumps with continuous aggregation: 6
Number of mass action jumps: 0


julia> jump_sol = solve(jump_prob, Tsit5(); saveat = 1.0)
ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.


Closest candidates are:
  convert(::Type{T}, ::Factorization) where T<:AbstractArray
   @ LinearAlgebra ~/.julia/juliaup/julia-1.9.0-rc2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:59
  convert(::Type{T}, ::T) where T<:AbstractArray
   @ Base abstractarray.jl:16
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:64
  ...

Stacktrace:
  [1] push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, item::Vector{Float64})
    @ Base ./array.jl:1060
  [2] copyat_or_push!(a::Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}, i::Int64, x::Vector{Float64}, perform_copy::Bool)
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/VzH8Y/src/utils.jl:184
  [3] _savevalues!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…}, force_save::Bool, reduce_size::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:75
  [4] savevalues! (repeats 2 times)
    @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:57 [inlined]
  [5] apply_callback!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…}, callback::ContinuousCallback{…}, cb_time::Float64, prev_sign::Float64, event_idx::Int64)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/qPmC2/src/callbacks.jl:557
  [6] handle_callbacks!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:299
  [7] _loopfooter!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:242
  [8] loopfooter!
    @ ~/.julia/packages/OrdinaryDiffEq/yspeT/src/integrators/integrator_utils.jl:198 [inlined]
  [9] solve!(integrator::ODEIntegrator{true,Tsit5{Static.False,…},ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Float64,…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yspeT/src/solve.jl:521
 [10] __solve(jump_prob::JumpProblem{true, ODEProblem{true,ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}},Tuple{Float64, Float64},…}, Direct, CallbackSet{Tuple{ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}, ContinuousCallback{…}}, Tuple{}}, Nothing, Tuple{VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xd90a62cf, 0xc0f2dfd6, 0xe0ae39ad, 0x176a34a7, 0x4aae5eb2)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2992"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3aa0094e, 0x3dba15fe, 0x7cfca4eb, 0xa3b72dd6, 0x078af186)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x91c2ce70, 0x8401479e, 0xa981a223, 0x44ab2501, 0x7c50d39b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2993"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x5c39adbf, 0x9a8c7b43, 0xf37ef2a9, 0xfa34c7f3, 0x250d32d5)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf8e9e916, 0xe9725a58, 0xe0e2ea2f, 0x9e4c5782, 0x28633e09)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2994"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6b3b1896, 0xeb49f324, 0x16b5f36e, 0x5936e64b, 0x088e6644)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6fa6e9c5, 0x6ac1adbf, 0xe7751301, 0xa9d78f27, 0x6215a25b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2995"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x4777e92a, 0xab958d14, 0xedbc5525, 0xc4c1b0ab, 0x40793f45)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x62e08576, 0x38be2b6d, 0x98ba0ae7, 0xd9beffa6, 0xadc6c122)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2996"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x026081f8, 0x5c3a77b1, 0xef885e10, 0xf98c528b, 0x4ef709f1)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}, VariableRateJump{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3728bba8, 0x2104faac, 0x11bbe6d2, 0x921c1728, 0x77dc45a4)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##MTKIntegrator#2997"),), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x609930a2, 0x82139460, 0x9f060848, 0xac035f9c, 0x67d2fb0d)}, Nothing, Nothing, Nothing, Nothing, Float64, Int64}}, Nothing, Nothing, Random.TaskLocalRNG, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}}, alg::Tsit5{Static.False,…}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ JumpProcesses ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:6
 [11] __solve (repeats 5 times)
    @ ~/.julia/packages/JumpProcesses/8AMkl/src/solve.jl:1 [inlined]
 [12] #solve#33
    @ ~/.julia/packages/DiffEqBase/qPmC2/src/solve.jl:965 [inlined]
 [13] top-level scope
    @ ~/Documents/Repos/outbreak-detection/scripts/Seasonal-infectivity_catalyst.jl:72

Any advice would be appreciated, as I like the simplicity of the Catalyst interface, particularly as my actual model will be much more complex, and the idea of not needing to manually create the dependency graph is appealing!

arnold-c avatar Apr 07 '23 01:04 arnold-c

I posted it on Discourse, but just for tracking purposes here is the simplest approach that works at the moment:

using Catalyst, DifferentialEquations

sir_model = @reaction_network begin
           β*t, S + I --> 2I
           ν, I --> R
       end
setdefaults!(sir_model, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])
defs = ModelingToolkit.defaults(sir_model)

@named osys = ODESystem(Equation[], ModelingToolkit.get_iv(sir_model), 
                                                states(sir_model), 
                                                parameters(sir_model); 
                                                defaults = defs)
oprob = ODEProblem(osys, [], (0.0, 200.0), check_length = false)
jprob = JumpProblem(sir_model, oprob, Direct())
sol = solve(jprob, Tsit5())

There are two issues to resolve for a nicer interface.

  1. How to get a ModelingToolkit interface for the ODEProblem when one wants to use VariableRateJumps with the ODE timesteppers?
  2. How to handle rate bounds (we can't really calculate them in general in Catalyst, but need a way to support attaching them to reactions and propagating that information through MTK , ideally allowing both symbolic and function-based bounds)?

The first is the issue linked by @arnold-c, and I'm still not sure how to handle it. The second could be handled by having some kwargs when building Reactions to pass the bounds and time-window, and then propagating this through MTK. In JumpSystem we could check if the bound is symbolic or a function, and compile if the former before generating the VariableRateJump.

isaacsas avatar Apr 07 '23 13:04 isaacsas

I'm having the same problem at the moment. Just to check, this solution does not contemplate the VariableRateJumps yet, right?

filipepfarias avatar Sep 05 '23 08:09 filipepfarias

This should work for VariableRateJumps. Have you tried it?

isaacsas avatar Sep 05 '23 08:09 isaacsas

I updated my earlier comment to an example that has an explicit time dependence in the rate. Hopefully that makes it more clear that this should work for VariableRateJump problems.

isaacsas avatar Sep 05 '23 10:09 isaacsas

It worked! Thanks a lot!

filipepfarias avatar Sep 06 '23 08:09 filipepfarias

I think this can be closed (but correct me if wrong).

TorkelE avatar Sep 17 '23 02:09 TorkelE

It can be closed in the sense that there is a work around, but there isn’t yet a nice interface for this case.

isaacsas avatar Sep 17 '23 02:09 isaacsas

Ok, let's keep it open with that as the aim then.

TorkelE avatar Sep 17 '23 02:09 TorkelE

Some comments on it. The output of the solution, even if the jumps are on discrete state space, are being stored as floats, along with, at least I think to be the integration of the rates in time. Is there any option to overcome this? If I do an ensemble of trajectories , the memory allocation may grow really fast.

filipepfarias avatar Sep 17 '23 07:09 filipepfarias

The generic VariableRateJump solvers via ODE timesteppers require floating point values (since ODE solvers don't work with integers). They use the ODE solvers to integrate the intensities (hence why they store those rates too). Are you sure you aren't just saving more often than needed? Make sure you are setting save_positions = (false,false) in the JumpProblem so you aren't saving the state with every jump, and look at save_idxs and other saving controls for the ODE stepper to reduce memory pressure:

https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/

Likewise with EnsembleProblems there is the output_func keyword that you can use to control what information from each solution is actually saved:

https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Performing-an-Ensemble-Simulation

isaacsas avatar Sep 17 '23 12:09 isaacsas

@TorkelE getting into the various saving options and controls would probably make a nice tutorial too, especially combining with EnsembleProblem to handle saving for each trajectory.

isaacsas avatar Sep 17 '23 12:09 isaacsas

Yeah, it would. I plan to focus extra on Catalyst docs and functionality up until the ICSB workshop (October 8th), hopefully I can get to this as well.

TorkelE avatar Sep 17 '23 14:09 TorkelE