Catalyst.jl VariableRateJump

Specific Domains Modelling & Simulations
, ,
1

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!

1 Like
2

I’ll put an example together and get back to you later today or tomorrow. It isn’t completely automated so I need to check how the ODEProblem to JumpProblem coupling should be done such that indexing with symbols still works.

Longer-term this is definitely a good question and something we need to make easier with MTK and Catalyst. If you wouldn’t mind opening an issue on Catalyst about that to keep it visible there I’d appreciate it.

1 Like
3

Thanks very much @isaacsas. I have opened a GitHub issue on Catalyst.jl here

4

This seems to work

using Catalyst, DifferentialEquations

sir_model = @reaction_network begin
           β, 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)
u0 = map(s -> defs[s], states(sir_model))
p = map(s -> defs[s], parameters(sir_model))
ofun = ODEFunction((du,u,p,t) -> nothing; sys = sir_model, 
                                 syms = Symbol.(states(sir_model)),
                                 indepsym = Symbol(ModelingToolkit.get_iv(sir_model)),
                                 paramsyms = Symbol.(parameters(sir_model)))
oprob = ODEProblem(ofun, u0, (0.0, 200.0), p)
jprob = JumpProblem(sir_model, oprob, Direct())
sol = solve(jprob, Tsit5())

A few notes:

  1. There is no point using coevolve when going via Catalyst as Catalyst can’t auto-generate the rate bounds/windows (this seems like it would be quite hard to do in general). I’d like to add the ability for users to provide such bounds to Catalyst.Reactions, and then propagate them through to JumpProcesses eventually, but right now they can’t be used via Catalyst. So if you want the better performance of coevolve you need to manually create your jumps or see 3. below.
  2. This means one has to simulate VariableRateJump systems via the ODE integrator interface when coming via Catalyst.
  3. All that said, if you look at JumpSystem in ModelingToolkit, and are willing to dig into the internals a bit, you could always build the jumps you need via Catalyst/ModelingToolkit, along with the dependency graphs, using the ModelingToolkit functionality in JumpSystem, and then manually add in rate bounds for those that are VariableRateJumps. That would allow you to get the symbolically generated rate and affect functions, MassActionJumps, and dependency graphs, which you could then supplement with manually coded rate bounds for the VariableRateJumps.
5

Thanks!

6

Here is another approach, that is a bit simpler than my first suggestion. It involves making an empty ODESystem with all the species and parameters from Catalyst:

using Catalyst, DifferentialEquations

sir_model = @reaction_network begin
           β, 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())