## ----include = FALSE----------------------------------------------------------
set.seed(20250820)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
K <- 10^6

## ----setup--------------------------------------------------------------------
library(data.table)
library(nhppp)

## ----population---------------------------------------------------------------
pop <- setDT(
  list(
    id = 1:K,
    T0 = rep(20, K),
    T1 = stats::runif(n = K, min = 60, max = 100)
  )
)
setindex(pop, id)
pop

## -----------------------------------------------------------------------------
pop[, `:=`(
  alpha_weibull = stats::runif(n = K, min = 2.28, max = 3.28),
  sigma_weibull = stats::runif(n = K, min = 64, max = 84)
  )]

## ----function-lambda-weibull--------------------------------------------------
l_weibull <- function(t, alpha = pop$alpha_weibull, sigma = pop$sigma_weibull, ...) {
  alpha/sigma * (t/sigma)^(alpha-1)
}

L_weibull <- function(t, alpha = pop$alpha_weibull, sigma = pop$sigma_weibull, ...) {
  (t/sigma)^alpha
}

Li_weibull<- function(z, alpha = pop$alpha_weibull, sigma = pop$sigma_weibull, ...) {
  sigma * z^(1/alpha)
} 

## ----method-1-----------------------------------------------------------------
tictoc::tic("Method 1 (vectorized, thinning)")
M <- 5
break_points <- seq.int(from = 20, to = 100, length.out = M + 1)
breaks_mat <- matrix(rep(break_points, each = K), nrow = K)

lmaj_mat <- nhppp::get_step_majorizer(
  fun = l_weibull,
  breaks = breaks_mat,
  is_monotone = TRUE
)

pop[
  ,
  t_thinning := nhppp::vdraw_intensity(
    lambda = l_weibull,
    lambda_maj_matrix = lmaj_mat,
    rate_matrix_t_min = 20,
    rate_matrix_t_max = 100,
    t_min = pop$T0,
    t_max = pop$T1,
    atmost1 = TRUE,
    atmostB = 5
  )
]
tictoc::toc(log = TRUE) # timer end

## ----method-2-----------------------------------------------------------------
tictoc::tic("Method 2 (vectorized, inversion)")
pop[
  ,
  t_inversion := nhppp::vdraw_cumulative_intensity(
    Lambda = L_weibull,
    Lambda_inv = Li_weibull,
    t_min = pop$T0,
    t_max = pop$T1,
    atmost1 = TRUE
  )
]
tictoc::toc(log = TRUE) # timer end

## ----qq-plot, fig.alt="QQ plot comparing simulated times with the two methods. The QQ plot indicates excellent agreement."----
qqplot(pop$t_thinning, pop$t_inversion)

## ----histogram----------------------------------------------------------------
if(nchar(system.file(package = 'ggplot2'))>0) {
  pop2 <- setDT(list(id = 1:10000)) 
  pop2[, `:=`(
    T0 = 20, 
    T1 = 100, 
    a = 2.78, 
    s = 74
  )]
  setindex(pop2, id)
  pop2
  
  ## re-define functions to use `pop2` by default 
  ##    clunckiness of `nhppp` -- help us improve it!
  l  <- function(t, a = pop2$a, s = pop2$s, ...) a/s * (t/s)^(a-1)
  L  <- function(t, a = pop2$a, s = pop2$s, ...) (t/s)^a
  Li <- function(z, a = pop2$a, s = pop2$s, ...)  s * z^(1/a)
  
  ## get all times (full trajectories)
  Z <- nhppp::vdraw_cumulative_intensity(
      Lambda = L,
      Lambda_inv = Li,
      t_min = pop2$T0,
      t_max = pop2$T1,
      atmost1 = FALSE
    )
  
  ## ignore non-realized times 
  Z <- as.vector(Z)
  Z <- Z[!is.na(Z)]
  
  x <- seq(from = 20, to = 100, length.out = 100)
  ## normalize the hazard rate to have area 1 between 20 and 100
  y <- l(x, a = pop2$a[1], s = pop2$s[1]) / (
        L(100, a = pop2$a[1], s = pop2$s[1]) - 
        L(20, a = pop2$a[1], s = pop2$s[1])
       )
  
  
  library(ggplot2)
  ggplot(xlim = c(0, 110)) +
    geom_histogram(
      aes(x = Z, y = after_stat(density)), 
      bins = 100, 
      fill = "gray", 
      color = "white") +
    geom_line(
      aes(x = x, y = y), 
      color = "red"
    ) +
    labs(title = "Histogram of all simulated Weibull times", x = "Time", y = "Empirical and theoretical density")
}