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

## ----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[, `:=`(
  a_gompertz = stats::runif(n = K, min = 0.0045, max = 0.0055),
  b_gompertz = stats::runif(n = K, min = 0.085, max = 0.095)
  )]

## ----function-lambda-gompertz-------------------------------------------------
l_gompertz <- function(t, a = pop$a_gompertz, b = pop$b_gompertz, ...) {
  a * b * exp(b * t)
}

L_gompertz <- function(t, a = pop$a_gompertz, b = pop$b_gompertz, ...) {
  a * exp(b * t) - a
}

Li_gompertz <- function(z, a = pop$a_gompertz, b = pop$b_gompertz, ...) {
  log(z/a +1) / b
}
 

## ----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_gompertz,
  breaks = breaks_mat,
  is_monotone = TRUE
)

pop[
  ,
  t_thinning := nhppp::vdraw_intensity(
    lambda = l_gompertz,
    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 = 10
  )
]
tictoc::toc(log = TRUE) # timer end

## ----method-2-----------------------------------------------------------------
tictoc::tic("Method 2 (vectorized, inversion)")
pop[
  ,
  t_inversion := nhppp::vdraw_cumulative_intensity(
    Lambda = L_gompertz,
    Lambda_inv = Li_gompertz,
    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)

## ----all_events---------------------------------------------------------------
Z <- nhppp::vdraw_cumulative_intensity(
    Lambda = L_gompertz,
    Lambda_inv = Li_gompertz,
    t_min = pop$T0,
    t_max = pop$T1,
    atmost1 = FALSE
  )

dim(Z) # a lot of events!

## ----histogram----------------------------------------------------------------
if(nchar(system.file(package = 'ggplot2'))>0) {
  pop2 <- setDT(list(id = 1:1000)) 
  pop2[, `:=`(
    T0 = 20, 
    T1 = 100, 
    a = 0.005, 
    b = 0.090
  )]
  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, b = pop2$b, ...) a * b * exp(b * t)
  L  <- function(t, a = pop2$a, b = pop2$b, ...) a * exp(b * t) - a
  Li <- function(z, a = pop2$a, b = pop2$b, ...) log(z/a +1) / b
  
  ## 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], b = pop2$b[1]) / (
        L(100, a = pop2$a[1], b = pop2$b[1]) - 
        L(20, a = pop2$a[1], b = pop2$b[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 Gompertz times", x = "Time", y = "Empirical and theoretical density")
}