add test 1.8 and minor fix

Author: xiaoming

.github/workflows/CI.yml

@@ -20,6 +20,7 @@ jobs:
         version:
           - '1.6'
           - '1.7'
+          - '1.8'
         os:
           - ubuntu-latest
         arch:

src/delayproblem.jl

@@ -71,44 +71,45 @@ delaysets = DelayJumpSet(delay_trigger,delay_complete,delay_interrupt)
 
 """
 mutable struct DelayJumpSet{T1<:Union{Function,Vector},T2<:Union{Function,Vector},T3<:Function}
-    """reactions in the Markovian part that trigger the change of the state of the delay channels or/and the state of the reactants upon initiation."""
-    delay_trigger::Dict{Int,T1}
-    """reactions in the Markovian part that change the state of the delay channels or/and the state of the reactants in the middle of on-going delay reactions."""
-    delay_complete::Dict{Int,T2}
-    """reactions that are initiated by delay trigger reactions and change the state of the delay channels or/and the state of the reactants upon completion."""
-    delay_interrupt::Dict{Int,T3}
-    """collection of indices of reactions that can interrupt the delay reactions. of `delay_trigger`."""
-    delay_trigger_set::Vector{Int}
-    """collection of indices of `delay_interrupt`."""
-    delay_interrupt_set::Vector{Int}    
+  """reactions in the Markovian part that trigger the change of the state of the delay channels or/and the state of the reactants upon initiation."""
+  delay_trigger::Dict{Int,T1}
+  """reactions in the Markovian part that change the state of the delay channels or/and the state of the reactants in the middle of on-going delay reactions."""
+  delay_complete::Dict{Int,T2}
+  """reactions that are initiated by delay trigger reactions and change the state of the delay channels or/and the state of the reactants upon completion."""
+  delay_interrupt::Dict{Int,T3}
+  """collection of indices of reactions that can interrupt the delay reactions. of `delay_trigger`."""
+  delay_trigger_set::Vector{Int}
+  """collection of indices of `delay_interrupt`."""
+  delay_interrupt_set::Vector{Int}
 end
-function DelayJumpSet(delay_trigger::Dict,delay_complete::Dict,delay_interrupt::Dict) 
-    delay_trigger_, delay_complete_, delay_interrupt_ = convert_delayset.([delay_trigger, delay_complete, delay_interrupt])
-    DelayJumpSet(delay_trigger_, delay_complete_, delay_interrupt_, collect(keys(delay_trigger_)), collect(keys(delay_interrupt_)))
+function DelayJumpSet(delay_trigger::Dict, delay_complete::Dict, delay_interrupt::Dict)
+  delay_trigger_, delay_complete_, delay_interrupt_ = convert_delayset.([delay_trigger, delay_complete, delay_interrupt])
+  DelayJumpSet(delay_trigger_, delay_complete_, delay_interrupt_, collect(keys(delay_trigger_)), collect(keys(delay_interrupt_)))
 end
 
-convert_delayset(delay_set::Dict) = isempty(delay_set) ? convert(Dict{Int,Function},delay_set) : delay_set
+convert_delayset(delay_set::Dict) = isempty(delay_set) ? convert(Dict{Int,Function}, delay_set) : delay_set
 
 #BEGIN DelayJump
-mutable struct DelayJumpProblem{iip,P,A,C,J<:Union{Nothing,AbstractJumpAggregator},J2,J3,J4,J5,deType} <: DiffEqBase.AbstractJumpProblem{P,J}
-    prob::P
-    aggregator::A
-    discrete_jump_aggregation::J
-    jump_callback::C
-    variable_jumps::J2
-    regular_jump::J3
-    massaction_jump::J4
-    delayjumpsets::J5
-    de_chan0::deType
-    save_delay_channel::Bool
+mutable struct DelayJumpProblem{iip,P,A,C,J<:Union{Nothing,AbstractJumpAggregator},J2,J3,J4,J5,deType,K} <: DiffEqBase.AbstractJumpProblem{P,J}
+  prob::P
+  aggregator::A
+  discrete_jump_aggregation::J
+  jump_callback::C
+  variable_jumps::J2
+  regular_jump::J3
+  massaction_jump::J4
+  delayjumpsets::J5
+  de_chan0::deType
+  save_delay_channel::Bool
+  kwargs::K
 end
 
-function DelayJumpProblem(p::P,a::A,dj::J,jc::C,vj::J2,rj::J3,mj::J4,djs::J5,de_chan0::deType,save_delay_channel::Bool) where {P,A,J,C,J2,J3,J4,J5,deType}
-    if !(typeof(a)<:AbstractDelayAggregatorAlgorithm)
-      error("To solve DelayJumpProblem, one has to use one of the delay aggregators.")
-    end
-    iip = isinplace_jump(p,rj)
-    DelayJumpProblem{iip,P,A,C,J,J2,J3,J4,J5,deType}(p,a,dj,jc,vj,rj,mj,djs,de_chan0,save_delay_channel)
+function DelayJumpProblem(p::P, a::A, dj::J, jc::C, vj::J2, rj::J3, mj::J4, djs::J5, de_chan0::deType, save_delay_channel::Bool, kwargs::K) where {P,A,J,C,J2,J3,J4,J5,deType,K}
+  if !(typeof(a) <: AbstractDelayAggregatorAlgorithm)
+    error("To solve DelayJumpProblem, one has to use one of the delay aggregators.")
+  end
+  iip = isinplace_jump(p, rj)
+  DelayJumpProblem{iip,P,A,C,J,J2,J3,J4,J5,deType,K}(p, a, dj, jc, vj, rj, mj, djs, de_chan0, save_delay_channel, kwargs)
 end
 
 
@@ -136,27 +137,27 @@ end
     The initial condition of the delay channel.
 """
 function DelayJumpProblem(prob, aggregator::AbstractDelayAggregatorAlgorithm, jumps::JumpSet, delayjumpsets::DelayJumpSet, de_chan0;
-                     save_positions = typeof(prob) <: DiffEqBase.AbstractDiscreteProblem ? (false,true) : (true,true),
-                     rng = Xorshifts.Xoroshiro128Star(rand(UInt64)), scale_rates = false, useiszero = true,  save_delay_channel = false, kwargs...)
+  save_positions=typeof(prob) <: DiffEqBase.AbstractDiscreteProblem ? (false, true) : (true, true),
+  rng=Xorshifts.Xoroshiro128Star(rand(UInt64)), scale_rates=false, useiszero=true, save_delay_channel=false, callback=nothing, kwargs...)
 
   # initialize the MassActionJump rate constants with the user parameters
-  if using_params(jumps.massaction_jump) 
+  if using_params(jumps.massaction_jump)
     rates = jumps.massaction_jump.param_mapper(prob.p)
-    maj = MassActionJump(rates, jumps.massaction_jump.reactant_stoch, jumps.massaction_jump.net_stoch, 
-                         jumps.massaction_jump.param_mapper; scale_rates=scale_rates, useiszero=useiszero, 
-                         nocopy=true)
+    maj = MassActionJump(rates, jumps.massaction_jump.reactant_stoch, jumps.massaction_jump.net_stoch,
+      jumps.massaction_jump.param_mapper; scale_rates=scale_rates, useiszero=useiszero,
+      nocopy=true)
   else
     maj = jumps.massaction_jump
   end
 
 
   ## Constant Rate Handling
-  t,end_time,u = prob.tspan[1],prob.tspan[2],prob.u0
+  t, end_time, u = prob.tspan[1], prob.tspan[2], prob.u0
   if (typeof(jumps.constant_jumps) <: Tuple{}) && (maj === nothing)  # check if there are no jumps
     disc = nothing
     constant_jump_callback = CallbackSet()
   else
-    disc = aggregate(aggregator,u,prob.p,t,end_time,jumps.constant_jumps,maj,save_positions,rng; kwargs...)
+    disc = aggregate(aggregator, u, prob.p, t, end_time, jumps.constant_jumps, maj, save_positions, rng; kwargs...)
     constant_jump_callback = DiscreteCallback(disc)
   end
 
@@ -167,21 +168,23 @@ function DelayJumpProblem(prob, aggregator::AbstractDelayAggregatorAlgorithm, ju
     new_prob = prob
     variable_jump_callback = CallbackSet()
   else
-    new_prob = extend_problem(prob,jumps)
-    variable_jump_callback = build_variable_callback(CallbackSet(),0,jumps.variable_jumps...)
+    new_prob = extend_problem(prob, jumps)
+    variable_jump_callback = build_variable_callback(CallbackSet(), 0, jumps.variable_jumps...)
   end
-  callbacks = CallbackSet(constant_jump_callback,variable_jump_callback)
+  callbacks = CallbackSet(constant_jump_callback, variable_jump_callback)
+
+  solkwargs = make_kwarg(; callback)
 
   DelayJumpProblem{iip,typeof(new_prob),typeof(aggregator),typeof(callbacks),
-              typeof(disc),typeof(jumps.variable_jumps),
-              typeof(jumps.regular_jump),typeof(maj),typeof(delayjumpsets),typeof(de_chan0)}(
-                        new_prob,aggregator,disc,
-                        callbacks,
-                        jumps.variable_jumps,
-                        jumps.regular_jump, maj, delayjumpsets, de_chan0, save_delay_channel)
+    typeof(disc),typeof(jumps.variable_jumps),
+    typeof(jumps.regular_jump),typeof(maj),typeof(delayjumpsets),typeof(de_chan0),typeof(solkwargs)}(
+    new_prob, aggregator, disc,
+    callbacks,
+    jumps.variable_jumps,
+    jumps.regular_jump, maj, delayjumpsets, de_chan0, save_delay_channel, solkwargs)
 end
 
-
+make_kwarg(; kwargs...) = kwargs
 """
     function DelayJumpProblem(js::JumpSystem, prob, aggregator, delayjumpset, de_chan0; kwargs...)
 # Fields
@@ -203,34 +206,37 @@ end
 
     The initial condition of the delay channel.
 """
-function DelayJumpProblem(js::JumpSystem, prob, aggregator, delayjumpset, de_chan0; scale_rates = false, save_delay_channel = false, kwargs...)
-    statetoid = Dict(value(state) => i for (i,state) in enumerate(states(js)))
-    eqs       = equations(js)
-    invttype  = prob.tspan[1] === nothing ? Float64 : typeof(1 / prob.tspan[2])
-
-    # handling parameter substition and empty param vecs
-    p = (prob.p isa DiffEqBase.NullParameters || prob.p === nothing) ? Num[] : prob.p
-
-    majpmapper = ModelingToolkit.JumpSysMajParamMapper(js, p; jseqs=eqs, rateconsttype=invttype)
-    majs = isempty(eqs.x[1]) ? nothing : assemble_maj(eqs.x[1], statetoid, majpmapper)
-    crjs = ConstantRateJump[assemble_crj(js, j, statetoid) for j in eqs.x[2]]
-    vrjs = VariableRateJump[assemble_vrj(js, j, statetoid) for j in eqs.x[3]]
-    ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")
-    jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs)
-
-    if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator)
-        jdeps = asgraph(js)
-        vdeps = variable_dependencies(js)
-        vtoj = jdeps.badjlist
-        jtov = vdeps.badjlist
-        jtoj = needs_depgraph(aggregator) ? eqeq_dependencies(jdeps, vdeps).fadjlist : nothing
-        dep_graph_delay = dep_gr_delay(delayjumpset, vtoj, length(eqs))
-    else
-        vtoj = nothing; jtov = nothing; jtoj = nothing; dep_graph_delay = nothing;
-    end
+function DelayJumpProblem(js::JumpSystem, prob, aggregator, delayjumpset, de_chan0; scale_rates=false, save_delay_channel=false, kwargs...)
+  statetoid = Dict(value(state) => i for (i, state) in enumerate(states(js)))
+  eqs = equations(js)
+  invttype = prob.tspan[1] === nothing ? Float64 : typeof(1 / prob.tspan[2])
+
+  # handling parameter substition and empty param vecs
+  p = (prob.p isa DiffEqBase.NullParameters || prob.p === nothing) ? Num[] : prob.p
+
+  majpmapper = ModelingToolkit.JumpSysMajParamMapper(js, p; jseqs=eqs, rateconsttype=invttype)
+  majs = isempty(eqs.x[1]) ? nothing : assemble_maj(eqs.x[1], statetoid, majpmapper)
+  crjs = ConstantRateJump[assemble_crj(js, j, statetoid) for j in eqs.x[2]]
+  vrjs = VariableRateJump[assemble_vrj(js, j, statetoid) for j in eqs.x[3]]
+  ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")
+  jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs)
+
+  if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator)
+    jdeps = asgraph(js)
+    vdeps = variable_dependencies(js)
+    vtoj = jdeps.badjlist
+    jtov = vdeps.badjlist
+    jtoj = needs_depgraph(aggregator) ? eqeq_dependencies(jdeps, vdeps).fadjlist : nothing
+    dep_graph_delay = dep_gr_delay(delayjumpset, vtoj, length(eqs))
+  else
+    vtoj = nothing
+    jtov = nothing
+    jtoj = nothing
+    dep_graph_delay = nothing
+  end
 
 
-    DelayJumpProblem(prob, aggregator, jset, delayjumpset, de_chan0; save_delay_channel=save_delay_channel, dep_graph=jtoj, vartojumps_map=vtoj, jumptovars_map=jtov, scale_rates=scale_rates, dep_graph_delay = dep_graph_delay, nocopy=true, kwargs...)
+  DelayJumpProblem(prob, aggregator, jset, delayjumpset, de_chan0; save_delay_channel=save_delay_channel, dep_graph=jtoj, vartojumps_map=vtoj, jumptovars_map=jtov, scale_rates=scale_rates, dep_graph_delay=dep_graph_delay, nocopy=true, kwargs...)
 end
 
 
@@ -242,47 +248,47 @@ function DiffEqBase.remake(thing::DelayJumpProblem; kwargs...)
   errmesg = """
   DelayJumpProblems can currently only be remade with new u0, de_chan0, p, tspan, delayjumpsets fields, prob fields. 
   """
-  !issubset(keys(kwargs),((:u0,:de_chan0,:p,:tspan,:prob)...,propertynames(thing.delayjumpsets)...)) && error(errmesg)
+  !issubset(keys(kwargs), ((:u0, :de_chan0, :p, :tspan, :prob)..., propertynames(thing.delayjumpsets)...)) && error(errmesg)
 
   if :prob ∉ keys(kwargs)
     dprob = DiffEqBase.remake(thing.prob; kwargs...)
     # if the parameters were changed we must remake the MassActionJump too
     if (:p ∈ keys(kwargs)) && JumpProcesses.using_params(thing.massaction_jump)
-        JumpProcesses.update_parameters!(thing.massaction_jump, dprob.p; kwargs...)
-    end      
+      JumpProcesses.update_parameters!(thing.massaction_jump, dprob.p; kwargs...)
+    end
   else
-    any(k -> k in keys(kwargs), (:u0,:p,:tspan)) && error("If remaking a DelayJumpProblem you can not pass both prob and any of u0, p, or tspan.")
+    any(k -> k in keys(kwargs), (:u0, :p, :tspan)) && error("If remaking a DelayJumpProblem you can not pass both prob and any of u0, p, or tspan.")
     dprob = kwargs[:prob]
 
     # we can't know if p was changed, so we must remake the MassActionJump
     if JumpProcesses.using_params(thing.massaction_jump)
       JumpProcesses.update_parameters!(thing.massaction_jump, dprob.p; kwargs...)
-    end 
+    end
   end
-  if any(k -> k in keys(kwargs), propertynames(thing.delayjumpsets)) 
-      delayjumpsets = update_delayjumpsets(thing.delayjumpsets; kwargs...)
+  if any(k -> k in keys(kwargs), propertynames(thing.delayjumpsets))
+    delayjumpsets = update_delayjumpsets(thing.delayjumpsets; kwargs...)
   else
-      delayjumpsets = thing.delayjumpsets
+    delayjumpsets = thing.delayjumpsets
   end
   de_chan0 = :de_chan0 ∈ keys(kwargs) ? kwargs[:de_chan0] : thing.de_chan0
 
   DelayJumpProblem(dprob, thing.aggregator, thing.discrete_jump_aggregation, thing.jump_callback,
-     thing.variable_jumps, thing.regular_jump, thing.massaction_jump, delayjumpsets, de_chan0, thing.save_delay_channel)
+    thing.variable_jumps, thing.regular_jump, thing.massaction_jump, delayjumpsets, de_chan0, thing.save_delay_channel, thing.kwargs)
 end
 
 function update_delayjumpsets(delayjumpsets::DelayJumpSet; kwargs...)
-    @unpack delay_trigger, delay_complete, delay_interrupt = delayjumpsets
-    for (key, value) in kwargs
-        exp_key = toexpr(key)  
-        if exp_key == :delay_trigger
-            delay_trigger = value
-        elseif exp_key == :delay_complete
-            delay_complete = value
-        elseif exp_key == :delay_interrupt
-            delay_interrupt = value
-        end
+  @unpack delay_trigger, delay_complete, delay_interrupt = delayjumpsets
+  for (key, value) in kwargs
+    exp_key = toexpr(key)
+    if exp_key == :delay_trigger
+      delay_trigger = value
+    elseif exp_key == :delay_complete
+      delay_complete = value
+    elseif exp_key == :delay_interrupt
+      delay_interrupt = value
     end
-    DelayJumpSet(delay_trigger, delay_complete, delay_interrupt) 
+  end
+  DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
 end
 
 # function update_delayjumpsets(delayjumpsets::DelayJumpSet; kwargs...)
@@ -299,20 +305,20 @@ end
 #   return delayjumpsets_
 # end
 
-Base.summary(io::IO, prob::DelayJumpProblem) = string(DiffEqBase.parameterless_type(prob)," with problem ",DiffEqBase.parameterless_type(prob.prob)," and aggregator ",typeof(prob.aggregator))
+Base.summary(io::IO, prob::DelayJumpProblem) = string(DiffEqBase.parameterless_type(prob), " with problem ", DiffEqBase.parameterless_type(prob.prob), " and aggregator ", typeof(prob.aggregator))
 function Base.show(io::IO, mime::MIME"text/plain", A::DelayJumpProblem)
-  println(io,summary(A))
-  println(io,"Number of constant rate jumps: ",A.discrete_jump_aggregation === nothing ? 0 : num_constant_rate_jumps(A.discrete_jump_aggregation))
-  println(io,"Number of variable rate jumps: ",length(A.variable_jumps))
+  println(io, summary(A))
+  println(io, "Number of constant rate jumps: ", A.discrete_jump_aggregation === nothing ? 0 : num_constant_rate_jumps(A.discrete_jump_aggregation))
+  println(io, "Number of variable rate jumps: ", length(A.variable_jumps))
   if A.regular_jump !== nothing
-    println(io,"Have a regular jump")
+    println(io, "Have a regular jump")
   end
   if (A.massaction_jump !== nothing) && (get_num_majumps(A.massaction_jump) > 0)
-    println(io,"Have a mass action jump")
+    println(io, "Have a mass action jump")
   end
   if A.delayjumpsets !== nothing
-    println(io,"Number of delay trigger reactions: ",length(A.delayjumpsets.delay_trigger))
-    println(io,"Number of delay interrupt reactions: ",length(A.delayjumpsets.delay_interrupt))
+    println(io, "Number of delay trigger reactions: ", length(A.delayjumpsets.delay_trigger))
+    println(io, "Number of delay interrupt reactions: ", length(A.delayjumpsets.delay_interrupt))
   end
 end
 #END DelayJump