L -> rateinterval

Author: isaacsas

docs/src/tutorials/discrete_stochastic_example.md

@@ -510,28 +510,29 @@ define an infection reaction as a bounded `VariableRateJump`, requiring us to
 again provide `rate` and `affect` functions, but also give functions that
 calculate an upper-bound on the rate (`urate(u,p,t)`), an optional lower-bound
 on the rate (`lrate(u,p,t)`), and a time window over which the bounds are valid
-as long as any states these three rates depend on are unchanged (`L(u,p,t)`).
+as long as any states these three rates depend on are unchanged
+(`rateinterval(u,p,t)`).
 The lower- and upper-bounds of the rate should be valid from the time they are
-computed `t` until `t + L(u, p, t)`:
+computed `t` until `t + rateinterval(u, p, t)`:
 
 ```@example tut2
 H = zeros(Float64, 10)
 rate3(u, p, t) = p[1]*u[1]*u[2] + p[3]*u[1]*sum(exp(-p[4]*(t - _t)) for _t in H)
 lrate = rate1              # β*S*I
 urate = rate3
-L(u, p, t) = 1 / (2*urate(u, p, t))
+rateinterval(u, p, t) = 1 / (2*urate(u, p, t))
 function affect3!(integrator)
   integrator.u[1] -= 1     # S -> S - 1
   integrator.u[2] += 1     # I -> I + 1
   push!(H, integrator.t)
   nothing
 end
-jump3 = VariableRateJump(rate3, affect3!; lrate=lrate, urate=urate, L=L)
+jump3 = VariableRateJump(rate3, affect3!; lrate, urate, rateinterval)
 ```
 Note that here we set the lower bound rate to be the normal SIR infection rate,
 and set the upper bound rate equal to the new rate of infection (`rate3`). As
 long as `u[1]` and `u[2]` are unchanged by another jump, for any `s` in `[t,t +
-L(u,p,t)]` we have that `lrate(u,p,t) <= rate3(u,p,s) <= urate(u,p,t)`.
+rateinterval(u,p,t)]` we have that `lrate(u,p,t) <= rate3(u,p,s) <= urate(u,p,t)`.
 
 Next, we redefine the recovery jump's `affect!` such that a random infection is
 removed from `H` for every recovery.

src/aggregators/coevolve.jl

@@ -18,14 +18,13 @@ mutable struct CoevolveJumpAggregation{T, S, F1, F2, RNG, GR, PQ} <:
     pq::PQ                            # priority queue of next time
     lrates::F1                        # vector of rate lower bound functions
     urates::F1                        # vector of rate upper bound functions
-    Ls::F1                            # vector of interval length functions
+    rateintervals::F1                 # vector of interval length functions
 end
 
 function CoevolveJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::Nothing,
                                  maj::S, rs::F1, affs!::F2, sps::Tuple{Bool, Bool},
-                                 rng::RNG; u::U,
-                                 dep_graph = nothing,
-                                 lrates, urates, Ls) where {T, S, F1, F2, RNG, U}
+                                 rng::RNG; u::U, dep_graph = nothing, lrates, urates,
+                                 rateintervals) where {T, S, F1, F2, RNG, U}
     if dep_graph === nothing
         if (get_num_majumps(maj) == 0) || !isempty(rs)
             error("To use Coevolve a dependency graph between jumps must be supplied.")
@@ -47,30 +46,22 @@ function CoevolveJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::Not
 
     pq = MutableBinaryMinHeap{T}()
     CoevolveJumpAggregation{T, S, F1, F2, RNG, typeof(dg),
-                            typeof(pq)
-                            }(nj, nj, njt,
-                              et,
-                              crs, sr, maj,
-                              rs,
-                              affs!, sps,
-                              rng,
-                              dg, pq,
-                              lrates, urates, Ls)
+                            typeof(pq)}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng,
+                                        dg, pq, lrates, urates, rateintervals)
 end
 
 # creating the JumpAggregation structure (tuple-based variable jumps)
 function aggregate(aggregator::Coevolve, u, p, t, end_time, constant_jumps,
-                   ma_jumps, save_positions, rng;
-                   dep_graph = nothing, variable_jumps = nothing,
-                   kwargs...)
+                   ma_jumps, save_positions, rng; dep_graph = nothing,
+                   variable_jumps = nothing, kwargs...)
     AffectWrapper = FunctionWrappers.FunctionWrapper{Nothing, Tuple{Any}}
     RateWrapper = FunctionWrappers.FunctionWrapper{typeof(t),
                                                    Tuple{typeof(u), typeof(p), typeof(t)}}
     affects! = Vector{AffectWrapper}()
     rates = Vector{RateWrapper}()
     lrates = Vector{RateWrapper}()
     urates = Vector{RateWrapper}()
-    Ls = Vector{RateWrapper}()
+    rateintervals = Vector{RateWrapper}()
 
     if (constant_jumps !== nothing) && !isempty(constant_jumps)
         append!(affects!,
@@ -86,7 +77,7 @@ function aggregate(aggregator::Coevolve, u, p, t, end_time, constant_jumps,
         append!(rates, [RateWrapper(j.rate) for j in variable_jumps])
         append!(lrates, [RateWrapper(j.lrate) for j in variable_jumps])
         append!(urates, [RateWrapper(j.urate) for j in variable_jumps])
-        append!(Ls, [RateWrapper(j.L) for j in variable_jumps])
+        append!(rateintervals, [RateWrapper(j.rateinterval) for j in variable_jumps])
     end
 
     num_jumps = get_num_majumps(ma_jumps) + length(urates)
@@ -96,9 +87,7 @@ function aggregate(aggregator::Coevolve, u, p, t, end_time, constant_jumps,
     next_jump_time = typemax(t)
     CoevolveJumpAggregation(next_jump, next_jump_time, end_time, cur_rates, sum_rate,
                             ma_jumps, rates, affects!, save_positions, rng;
-                            u = u,
-                            dep_graph = dep_graph,
-                            lrates = lrates, urates = urates, Ls = Ls)
+                            u, dep_graph, lrates, urates, rateintervals)
 end
 
 # set up a new simulation and calculate the first jump / jump time
@@ -144,8 +133,8 @@ end
     @inbounds return p.urates[uidx](u, params, t)
 end
 
-@inline function get_L(p::CoevolveJumpAggregation, lidx, u, params, t)
-    @inbounds return p.Ls[lidx](u, params, t)
+@inline function get_rateinterval(p::CoevolveJumpAggregation, lidx, u, params, t)
+    @inbounds return p.rateintervals[lidx](u, params, t)
 end
 
 @inline function get_lrate(p::CoevolveJumpAggregation, lidx, u, params, t)
@@ -173,9 +162,9 @@ function next_time(p::CoevolveJumpAggregation{T}, u, params, t, i, tstop::T) whe
     _t = t + s
     if lidx > 0
         while t < tstop
-            L = get_L(p, lidx, u, params, t)
-            if s > L
-                t = t + L
+            rateinterval = get_rateinterval(p, lidx, u, params, t)
+            if s > rateinterval
+                t = t + rateinterval
                 urate = get_urate(p, uidx, u, params, t)
                 s = urate == zero(t) ? typemax(t) : randexp(rng) / urate
                 _t = t + s

test/hawkes_test.jl

@@ -36,7 +36,7 @@ function hawkes_jump(i::Int, g, h)
     rate = hawkes_rate(i, g, h)
     lrate(u, p, t) = p[1]
     urate = rate
-    function L(u, p, t)
+    function rateinterval(u, p, t)
         _lrate = lrate(u, p, t)
         _urate = urate(u, p, t)
         return _urate == _lrate ? typemax(t) : 1 / (2 * _urate)
@@ -45,7 +45,7 @@ function hawkes_jump(i::Int, g, h)
         push!(h[i], integrator.t)
         integrator.u[i] += 1
     end
-    return VariableRateJump(rate, affect!; lrate = lrate, urate = urate, L = L)
+    return VariableRateJump(rate, affect!; lrate, urate, rateinterval)
 end
 
 function hawkes_jump(u, g, h)
@@ -57,8 +57,7 @@ function hawkes_problem(p, agg::Coevolve; u = [0.0], tspan = (0.0, 50.0),
                         g = [[1]], h = [[]])
     dprob = DiscreteProblem(u, tspan, p)
     jumps = hawkes_jump(u, g, h)
-    jprob = JumpProblem(dprob, agg, jumps...;
-                        dep_graph = g, save_positions = save_positions, rng = rng)
+    jprob = JumpProblem(dprob, agg, jumps...; dep_graph = g, save_positions, rng)
     return jprob
 end
 
@@ -72,7 +71,7 @@ function hawkes_problem(p, agg; u = [0.0], tspan = (0.0, 50.0),
                         g = [[1]], h = [[]])
     oprob = ODEProblem(f!, u, tspan, p)
     jumps = hawkes_jump(u, g, h)
-    jprob = JumpProblem(oprob, agg, jumps...; save_positions = save_positions, rng = rng)
+    jprob = JumpProblem(oprob, agg, jumps...; save_positions, rng)
     return jprob
 end
 
@@ -102,7 +101,7 @@ algs = (Direct(), Coevolve())
 Nsims = 250
 
 for alg in algs
-    jump_prob = hawkes_problem(p, alg; u = u0, tspan = tspan, g = g, h = h)
+    jump_prob = hawkes_problem(p, alg; u = u0, tspan, g, h)
     if typeof(alg) <: Coevolve
         stepper = SSAStepper()
     else