Renewal and delays

The discrete-time renewal recursion, the convolutions that carry infections through each delay, and the conversions between growth rate, doubling time and reproduction number. These are plain functions with no sampling in them, so they can be called directly and are covered by unit tests.

Index

• BVDOutbreakSize.ForecastHorizon
• BVDOutbreakSize.convolve_delay
• BVDOutbreakSize.convolve_pmf
• BVDOutbreakSize.deviation_knots
• BVDOutbreakSize.discretise_censored
• BVDOutbreakSize.doubling_time
• BVDOutbreakSize.euler_lotka_r
• BVDOutbreakSize.future_knot_days
• BVDOutbreakSize.horizon_days
• BVDOutbreakSize.implied_national_Rt
• BVDOutbreakSize.implied_national_Rt_at
• BVDOutbreakSize.importation_from_kernel
• BVDOutbreakSize.interpolate_knots
• BVDOutbreakSize.knot_days
• BVDOutbreakSize.lognormal_meansd
• BVDOutbreakSize.patch_infections
• BVDOutbreakSize.r_to_R0
• BVDOutbreakSize.relative_multiplier
• BVDOutbreakSize.renewal_infections
• BVDOutbreakSize.safe_rate
• BVDOutbreakSize.seed_at_renewal_start
• BVDOutbreakSize.seed_infections
• BVDOutbreakSize.seeding_age
• BVDOutbreakSize.sigmoid_ramp

Reference

BVDOutbreakSize.ForecastHorizon Type

struct ForecastHorizon

ForecastHorizon(days)

The number of days a model runs past its cut-off to forecast. Passed as a composer's forecast keyword (see with_horizon). The latent grid then runs to day n + days, the walks draw fresh innovations for the future knots, and each stream draws its future counts as missing observations for predict to generate. The default forecast = nothing is the fitted model. It is a type rather than an integer so the fitted model compiles with no forecast code in it.

---

Fields

• days::Int64

source

BVDOutbreakSize.convolve_delay Method

convolve_delay(
    x::AbstractVector,
    delay::AbstractVector
) -> Any

Convolve a daily trajectory x (infections or onsets) with a delay PMF delay (indexed from lag 0), returning the expected daily counts of the delayed event on the same daily grid: entry t sums x[t−d] · delay[d+1] over lags d that stay in range. Maps infections to onsets, onsets to deaths, onsets to reports and onsets to detected exports. Type-stable and AD-transparent.

Each lag adds one scaled, shifted copy of x with axpy!, a single BLAS call on float arrays. BLAS skips a lag whose weight is exactly zero, so on float arrays an Inf or NaN in x does not reach the days that lag feeds; other element types add 0 · x there and propagate it.

source

BVDOutbreakSize.convolve_pmf Method

convolve_pmf(a::AbstractVector, b::AbstractVector) -> Any

Discrete convolution of two delay PMFs a and b (both indexed from delay 0 at element 1), giving the PMF of the summed delay a ⊕ b. The result has length length(a) + length(b) - 1. Its mass equals sum(a) * sum(b), so normalised inputs give a normalised output. Used to build the infection→detection delay (incubation ⊕ onset-to-detection) and the infection→death delay (incubation ⊕ onset-to-death) for the exports streams from their component PMFs. Type-stable and AD-transparent.

source

BVDOutbreakSize.deviation_knots Method

deviation_knots(
    z_level::AbstractVector,
    z_drift::AbstractVector,
    σ_level::Real,
    σ_δ::AbstractVector,
    φ::Real,
    groups::AbstractVector{<:UnitRange},
    factors::AbstractVector,
    drift_factors::AbstractVector,
    walking::AbstractVector{Bool},
    walk_index::AbstractVector{<:Integer},
    n_walking::Integer,
    n_knots::Integer
) -> Any

Group-centred AR(1) deviation knots over the health zones of bvd_zone, the construction of the province deviations in patch_rt_model applied within each group. groups are the index ranges the deviations are centred within, one per patch.

The first knot is a level and the later knots revert toward zero at retention φ,

so every group sums to zero at every knot and no unit is privileged. factors[i] is a lower-triangular correlation factor over the units of group i, and drift_factors[i] the same over its walking units; nothing in either place leaves those draws independent. Scaling comes before centring, so a per-unit σ_δ still leaves the group sum at zero.

Only a walking unit carries an innovation and a level-only unit decays along the mean path φ^{k-1} δ_u(1). walk_index[u] is the unit's position among the n_walking walking units, zero for a level-only unit, and z_drift holds the innovations knot by knot at (k - 2) n_walking + walk_index[u].

Returns an (n_units × n_knots) matrix.

source

BVDOutbreakSize.discretise_censored Method

discretise_censored(dist, nmax::Integer) -> Any

Daily probability mass function for the continuous delay dist over lags 0, 1, …, nmax, discretised by double interval censoring (uniform primary event over a one-day window, then unit-interval censoring of the secondary event, truncated at nmax), as CensoredDistributions.double_interval_censored defines it. The truncation holds the CDF at one from nmax on, so the lag nmax entry is zero. For a LogNormal or Gamma delay the CDF differentiates cleanly under Mooncake, so this is AD-safe. Extreme warmup proposals that drive the total to a non-finite or zero value fall back to a uniform PMF, so the downstream convolution stays finite (the proposal is still rejected through its low log-likelihood). Returns a vector whose element type follows the delay parameters.

source

BVDOutbreakSize.doubling_time Method

doubling_time(r) -> Any

Doubling time log(2) / r implied by an exponential growth rate r. Returns a non-finite value as r crosses zero, matching the limit of an unbounded doubling time at zero growth.

source

BVDOutbreakSize.euler_lotka_r Method

euler_lotka_r(R, g::AbstractVector) -> Any

Exponential growth rate r implied by a reproduction number R and a generation-interval PMF g (indexed from lag 1): the root of the Euler–Lotka identity R · Σ_s g_s e^{−r s} = 1. Starts from the small-r approximation r ≈ (R − 1) / (R · ḡ) with ḡ the mean generation time, or below R = 1 from log(R) / ḡ, then takes Newton steps on log(R Σ_s g_s e^{−r s}) = 0 until a step is below 1e-12, at most 20. That function is convex and falling in r, and both starts sit below the root, so Newton climbs to it without overshooting. The second start is closer as R falls, and the log form stays close to linear there, so a near-zero R from a depleted pool gives a finite rate. R = 0 has no finite root and gives -Inf. Its Mooncake rule differentiates the root through the implicit function theorem. Mirrors the R_to_r seeding helper in EpiAware.jl and the implied-growth initialisation in the EpiNow2 Stan model.

source

BVDOutbreakSize.future_knot_days Method

future_knot_days(n::Integer, horizon::Integer; week) -> Any

Knot days after the cut-off n for a forecast of horizon days: one every week days and one on the last day. knot_days pins a knot to the cut-off, so these continue the fitted knots without moving any of them. They are also the future vintages the weekly forecast quantities are binned to. Empty when horizon is zero.

source

BVDOutbreakSize.horizon_days Method

horizon_days(_::Nothing) -> Int64

Days a model runs past its cut-off: 0 for the fitted model.

source

BVDOutbreakSize.implied_national_Rt Method

implied_national_Rt(
    infections_total::AbstractVector,
    g::AbstractVector
) -> Any

Derive the implied national reproduction number from a summed infection trajectory by inverting the renewal equation:

Rt_national(t) = I_total(t) / sum_s I_total(t-s) * g_s

This reconstructs what a single-patch model would estimate as the national Rt from the aggregated infection count. The first day is set to zero (no prior infections to divide by). Days where the force of infection is zero (no prior infections) also return zero. AD-transparent under Mooncake (only arithmetic and @inbounds loops).

A model that needs the reproduction number on one day only should call implied_national_Rt_at, which this is the trajectory form of.

source

BVDOutbreakSize.implied_national_Rt_at Method

implied_national_Rt_at(
    infections_total::AbstractVector,
    g::AbstractVector,
    t::Integer
) -> Any

implied_national_Rt_at(infections_total, g, t)

The implied national reproduction number on a single day t, the one entry implied_national_Rt would put at index t. Day one and any day whose force of infection is zero give zero, as they do there.

The model reports the aggregate reproduction number at the cut-off alone, and building the whole trajectory to read its last entry would put n divisions and n force sums on the gradient tape for one number.

source

BVDOutbreakSize.importation_from_kernel Method

importation_from_kernel(
    K::AbstractMatrix,
    I_prev::AbstractVector,
    epsilon::Real
) -> Any

importation_from_kernel(K, I_prev, epsilon)

Per-patch importation into each of n_patches patches on a single day, given the n_patches x n_patches importation kernel K, the previous day's infections per patch I_prev (length n_patches), and the importation intensity epsilon.

K[p, q] is the per-capita daily travel rate from patch q to patch p (the first index is the destination). Diagonal entries should be zero (no self-importation). Each entry is unitless (a rate per day per traveller in the source patch).

Returns a length-n_patches vector of imported infections expected on the current day. AD-transparent under Mooncake.

source

BVDOutbreakSize.interpolate_knots Method

interpolate_knots(
    knot_vals::AbstractVector,
    days::AbstractVector{<:Integer},
    n::Integer
) -> Any

Linearly interpolate the knot values knot_vals, placed on the day indices days, onto the full daily grid 1:n, returning the length-n series. Piecewise-linear between bracketing knots, so the series bends only at the knots and is otherwise straight. Applied on the log-R_t scale, this gives weekly random-walk knots with within-week linear interpolation. Type-stable and AD-transparent (the output element type follows knot_vals).

Outside the knot span the series is held flat at the nearest knot value rather than extrapolated (the interpolation fraction is clamped to [0, 1]). This lets the reproduction-number walk start at a day > 1 and hold R_t flat at the established R0 over every earlier day, rather than running the first segment's slope backwards off the start of the grid.

source

BVDOutbreakSize.knot_days Method

knot_days(n::Integer; week, start) -> Any

Day indices of the weekly reproduction-number knots over an n-day grid. The first knot sits on day start (default 1) and the last knot on day n, with regular knots every week days, so a knot is pinned to start and to the end of the grid. With start > 1 the reproduction number is held flat (at the first knot's value) for all days before start (interpolate_knots clamps below the first knot), so the random walk only varies R_t from start onward. This fixes R_t over the pre-establishment seeding window before the genetic TMRCA bound. Returns a sorted vector of unique day indices.

source

BVDOutbreakSize.lognormal_meansd Method

lognormal_meansd(mean, sd) -> Distributions.LogNormal

LogNormal with the given mean and standard deviation sd, by moment matching var = mean^2 (exp(σ^2) − 1). The inputs are passed through safe_rate first so a NaN-prone warmup proposal cannot push σ = sqrt(log1p(·)) into NaN territory and trip the LogNormal domain check. Used by every delay submodel so a delay is parameterised by its mean and SD rather than the log-scale parameters.

source

BVDOutbreakSize.patch_infections Method

patch_infections(
    Rt_matrix::AbstractMatrix,
    g::AbstractVector,
    seeds_matrix::AbstractMatrix,
    importation_kernel::AbstractMatrix,
    epsilon::Union{Real, AbstractMatrix},
    N::AbstractVector
) -> NamedTuple{(:infections, :importation), <:Tuple{Any, Any}}

patch_infections(Rt_matrix, g, seeds_matrix, importation_kernel, epsilon, N)

Multi-patch (meta-population) renewal with between-patch importation and weak susceptible depletion. Each patch p generates infections from its own renewal force on a shared daily grid, and importation relocates a share of each day's generated infections through a kernel K:

y_{p,t} is taken as a rate on the patch's own pool of N[p] as in renewal_infections: I_{p,t} = S_{p,t−1}(1 − e^{−y_{p,t} / N_p}) and S_{p,t} = S_{p,t−1} e^{−y_{p,t} / N_p}.

Arguments

• Rt_matrix: n_patches x n_days matrix whose [p, t] entry is the reproduction number in patch p on day t. Each row is one patch's daily R_t trajectory.

• g: shared generation-interval PMF (indexed from lag 1, so g[1] is the probability of a one-day generation interval). Same for all patches.

• seeds_matrix: n_patches x L matrix whose [p, :] row is the pre-computed seed infection trajectory for patch p (see seed_infections). The seed fills days 1 ... L and the renewal recursion begins on day L+1.

• importation_kernel: n_patches x n_patches matrix K where K[p, q] is the share of patch q's transmission that lands in patch p rather than at home. The diagonal should be zero, and each column's off-diagonal sum times epsilon must be at most one, so a patch cannot export more transmission than it generates. Both hold for province_importation_kernel at any epsilon in [0, 1].

• epsilon: importation intensity of each origin, a scalar for every origin and day or an n_patches x n_days matrix. Importation is a transfer on the day it happens, so the origin patch is debited exactly what the destination patches are credited and coupling is never a source of infections. It is not conserved across days. The destination grows at its own reproduction number, so relocating infections from a fast patch to a slow one lowers the national total and the reverse raises it. With one shared reproduction number the transfer cancels exactly. Depletion acts on the destination after the transfer, so where a pool binds the destination takes in less than the origin sent.

• N: population of each patch, the pool it depletes. The seed is drawn from it first.

Returns

(; infections, importation). infections is a matrix of shape (n_patches, n_days) where row p is the daily infection trajectory for patch p. The first L days are copied from seeds_matrix and the remaining days are the renewal recursion with importation. importation is the matching matrix of infections each patch received from the others, the arrivals term alone before depletion rather than the net of arrivals and departures. The element type is promoted from all input types. AD-transparent under Mooncake.

source

BVDOutbreakSize.r_to_R0 Method

r_to_R0(r, g::AbstractVector) -> Any

Reproduction number R implied by an exponential growth rate r and a generation-interval PMF g (indexed from lag 1), the forward Euler–Lotka relation R = 1 / Σ_s g_s e^{−r s}. The inverse of euler_lotka_r, so a prior can be placed on the growth rate and the reproduction number derived from it under the model's generation interval. e^{−r s} is a running product of one exp, so the gradient tapes one exp, not one per lag.

source

BVDOutbreakSize.relative_multiplier Method

relative_multiplier(
    z::AbstractVector,
    σ::Real,
    groups::AbstractVector{<:UnitRange}
) -> Any

Partially pooled relative multiplier over the units of each group, with log contrasts that sum to zero within the group,

Q_g is the sum-to-zero basis of group g (sum_to_zero_basis), so the log multipliers have the distribution of n_g independent N(0, σ²) draws centred within the group, from the n_g - 1 directions a composition over the group can see. The level stays with the group and as σ shrinks the composition weights its units by the modelled quantity alone. This is the construction of province_composition_model, applied with one group per patch over its health zones in bvd_zone.

z holds each group's draws in turn, relative_multiplier_dims in all. A group of one unit takes no draw and its multiplier is one.

source

BVDOutbreakSize.renewal_infections Method

renewal_infections(
    Rt::AbstractVector,
    g::AbstractVector,
    seed::AbstractVector,
    N::Real
) -> Any

Daily latent infections from the renewal equation with weak susceptible depletion. Each day's renewal force R_t Σ_{s ≥ 1} I_{t−s} g_s, with generation-interval PMF g (indexed from lag 1) and per-day reproduction numbers Rt (length n), is taken as a rate on a pool of N:

While the pool is large against the outbreak, I_t ≈ (S_{t−1} / N) R_t f_t, the renewal scaled by the susceptible fraction. I_t never exceeds the pool left, so an overflowing force takes the pool rather than giving Inf. N is a population size, not an estimate of who can be reached. The models pass census counts (PROVINCE_POPULATIONS), and at that scale it is a light bound. Rt is therefore the reproduction number in a fully susceptible population.

A pre-computed seed of length L < n fills the first L days (see seed_infections) and is drawn from the pool before the recursion runs for days L+1 … n. Returns the length-n infection trajectory. The output element type is promoted from Rt, g, seed and N.

Multi-patch analogue

See patch_infections for the meta-population extension with between-patch importation.

source

BVDOutbreakSize.safe_rate Method

safe_rate(x) -> Any

NaN / Inf-safe positive rate. The renewal recursion can transiently overflow on extreme NUTS warmup proposals, giving a non-finite expected count. A plain max(x, eps) would propagate the NaN (max(NaN, eps) = NaN) and trip the Poisson / NegativeBinomial domain check.

source

BVDOutbreakSize.seed_at_renewal_start Method

seed_at_renewal_start(C_T_prior) -> Any

Renewal-start seed magnitude for the two-phase renewal: the daily infection incidence the analytic cryptic phase reaches on the renewal-start day.

It is an incidence, not a cumulative count. The cryptic total over the renewal_start grid days is larger by more than an order of magnitude, so reading the seed as cumulative would shift m by several generations.

The cryptic phase runs from the origin to the renewal start (≈ the genetic TMRCA day, off the renewal grid), spanning T = m · G days for m generations at mean generation interval G. A single daily infection at the origin grows over it at the cryptic rate r to

which the renewal grows forward under R_t, so the realised cut-off size stays data-driven while the prior fixes only the renewal-start scale.

The magnitude is referenced to the origin, not the cut-off. A cut-off-referenced size would put r into the seed and the renewal growth in opposing directions, cancelling for a fixed realised size and leaving a flat ridge along which R0 could slide. Referenced to the origin they compound instead, so r does enter the seed.

C_T_prior is returned unchanged. The helper names the quantity at the seeding call site.

source

BVDOutbreakSize.seed_infections Method

seed_infections(I0, r, len::Integer) -> Any

Seed the first len days of the infection trajectory as exponential growth I_t = I0 · e^{r (t − len)} at the implied growth rate r (see euler_lotka_r), so the seeding window is pinned at I0 on its last day and tails off backwards. This is the initialisation used by EpiNow2 and EpiAware.jl. Placing the whole seed on a single day instead injects a transient the renewal recursion has to relax away from. Returns a length-len vector whose element type follows I0 and r.

source

BVDOutbreakSize.seeding_age Method

seeding_age(cumulative::AbstractVector, n::Integer) -> Any

Outbreak age in days: the elapsed time from the model-implied seeding day to the cut-off (day n), where the seeding day is the smooth crossing at which cumulative infections first reach one. The crossing is linearly interpolated between the two days that bracket a cumulative of one, so it is a continuous function of the trajectory. Before the trajectory reaches one it returns n (the full grid). Used for the seeding-date plots and the genetic-TMRCA bound.

source

BVDOutbreakSize.sigmoid_ramp Method

sigmoid_ramp(
    n::Integer,
    day::Missing;
    ramp
) -> Vector{Float64}

Smooth intervention ramp over an n-day grid: the logistic curve 1 / (1 + e^{−(t − day) / ramp}) for each day t, rising from ≈0 well before day to ≈1 well after, with ramp setting the transition width in days. Multiplied by a sampled effect size and added to log-R_t, this gives an intervention (e.g. the first WHO situation report) a gradual ramped effect on transmission rather than an instantaneous step. Returns a length-n Float64 vector. day = missing gives an all-zero ramp (no intervention). Type-stable and AD-transparent in the effect size it multiplies.

Split into two dispatches on day's concrete type rather than one method with a runtime ismissing(day) branch. Mooncake cannot build a reverse rule for the branching form (TypeError: non-boolean (Missing) used in boolean context). Dispatch resolves day's type at the call site, so each method body is differentiated on its own concretely-typed slot.

source