## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 5
)

## ----setup--------------------------------------------------------------------
library(ringbp)
library(data.table)
library(tinyplot)

## -----------------------------------------------------------------------------
offspring <- offspring_opts(
  community = \(n) rnbinom(n = n, mu = 2.5, size = 0.16),
  isolated = \(n) rgeom(n = n, prob = 0.8), 
  asymptomatic = \(n) rpois(n = n, lambda = 1.5)
)

## -----------------------------------------------------------------------------
delays <- delay_opts(
  incubation_period = \(n) rlnorm(n, meanlog = 0.9, sdlog = 0.5),
  onset_to_isolation = \(n) rgamma(n, shape = 3, scale = 2)
)

## -----------------------------------------------------------------------------
event_probs <- event_prob_opts(
  asymptomatic = 0.1,
  presymptomatic_transmission = 0.5,
  symptomatic_traced = 0.2
)

## -----------------------------------------------------------------------------
interventions <- intervention_opts(quarantine = FALSE)

## -----------------------------------------------------------------------------
sim <- sim_opts(
  cap_max_days = 350,
  cap_cases = 4500
)

## -----------------------------------------------------------------------------
set.seed(1)

## -----------------------------------------------------------------------------
outbreak <- scenario_sim(
  n = 100,
  initial_cases = 1,
  offspring = offspring,
  delays = delays,
  event_probs = event_probs,
  interventions = interventions,
  sim = sim
)

## -----------------------------------------------------------------------------
tinyplot(
  cumulative ~ week | as.factor(sim), 
  data = outbreak, 
  type = "l", 
  lwd = 3,
  ylab = "Cumulative number of cases", 
  xlab = "Week",
  legend = FALSE,
  theme = "clean"
)

## -----------------------------------------------------------------------------
tinyplot(
  weekly_cases ~ week | as.factor(sim), 
  data = outbreak, 
  type = "l", 
  lwd = 3,
  ylab = "Weekly number of cases", 
  xlab = "Week",
  legend = FALSE,
  theme = "clean"
)

## -----------------------------------------------------------------------------
extinct_prob(outbreak)

## -----------------------------------------------------------------------------
scenarios <- data.table(
  expand.grid(
    initial_cases = c(5, 20),
    r0_community = c(1.5, 2.5),
    symptomatic_traced = c(0, 0.5, 1)
  )
)
scenarios[, scenario :=  1:.N]
scenarios <- scenarios[, list(data = list(.SD)), by = scenario]

## -----------------------------------------------------------------------------
n <- 10

scenarios[, sims := lapply(data, \(x, n) {
  scenario_sim(
    n = n,
    initial_cases = x$initial_cases,
    offspring = offspring_opts(
      community = \(n) rnbinom(n = n, mu = x$r0_community, size = 0.16),
      isolated = \(n) rnbinom(n = n, mu = 0, size = 1)
    ),
    delays = delay_opts(
      incubation_period = \(n) rweibull(n, shape = 1.65, scale = 4.28),
      onset_to_isolation = \(n) rweibull(n, shape = 1.65, scale = 2.31)
    ),
    event_probs = event_prob_opts(
       asymptomatic = 0.1,
       presymptomatic_transmission = 0.1,
       symptomatic_traced = x$symptomatic_traced
    ),
    interventions = intervention_opts(quarantine = FALSE),
    sim = sim_opts(cap_max_days = 365, cap_cases = 5000)
  )
  },
  n = n
)]

## -----------------------------------------------------------------------------

scenarios[, 
  pext := extinct_prob(sims[[1]], extinction_week = 12:16), 
  by = scenario
]

# probability of extinction for each scenario
scenarios$pext