# Truths by simulating the interventions in a large population (plan section 5.5).
#
# The truth population keeps the replicate's panel structure (patients per physician)
# but has proportionally more physicians, so the truth averages over the physician
# distribution rather than over one draw of 1,500 physicians. Every policy is run on the
# same population with the same course seed (common random numbers), so risk
# differences are paired and their Monte Carlo SE is sd(paired difference) / sqrt(n).
#
# Landmark, naive and RCT estimands are not interventional in this sense; they are added
# with the estimators (oracle large-sample limits).

OUTCOMES <- c(limb = "t_limb", male = "t_male", disease = "t_disease", clti = "t_clti_code", amp = "t_amp",
              bypass = "t_bypass", rr = "t_rr", death = "t_death", nco1 = "t_nco1", nco2 = "t_nco2")

truth_population <- function(par, n, seed) {
  scale <- n / par$size$n_patients
  par$size$n_physicians <- max(1L, round(par$size$n_physicians * scale))
  make_population(par, n = n, seed = seed)
}

# event indicator by month h (event in months 0..h-1)
by_h <- function(t, h) !is.na(t) & t < h

#' Truths for one world x scenario.
#'
#' @return list(truths = long table: outcome, estimand, value, mcse, n;
#'   curves = cumulative incidence by month under natural, A, B1, B2 for MALE and disease).
compute_truths <- function(world, scenario, n = 1e6, seed = 20261027L, horizon = 60L,
                           window = NULL, verbose = TRUE) {
  par <- load_params(world, scenario)
  W <- if (is.null(window)) par$size$window else window
  par$size$window <- W
  t0 <- Sys.time()
  pop <- truth_population(par, n, seed)
  cseed <- seed + 1L
  run <- function(p) simulate_course(pop, par, p, seed = cseed, record = FALSE)$events

  ev <- list(natural = run(policy("natural")), A = run(policy("A")),
             B1 = run(policy("B1")), B2 = run(policy("B2")))
  rows <- list()
  add <- function(outcome, estimand, d, nn = length(d)) {
    rows[[length(rows) + 1L]] <<- data.table(world = world, scenario = scenario, outcome = outcome,
                                             estimand = estimand, value = mean(d),
                                             mcse = sd(d) / sqrt(nn), n = nn)
  }
  for (o in names(OUTCOMES)) {
    col <- OUTCOMES[[o]]
    y <- lapply(ev, function(e) as.numeric(by_h(e[[col]], horizon)))
    for (p in names(y)) add(o, paste0("risk_", p), y[[p]])
    add(o, "rd_A_B1", y$A - y$B1)
    add(o, "rd_A_B2", y$A - y$B2)
  }
  # total procedures per patient (secondary outcome, D7)
  np <- lapply(ev, function(e) e$n_proc)
  for (p in names(np)) add("n_proc", paste0("mean_", p), np[[p]])
  add("n_proc", "diff_A_B1", np$A - np$B1)

  # sequential trials, ITT (D14): among patients eligible at month m under the natural
  # course, PVI forced vs forbidden at m, natural afterwards; risk by month `horizon`
  # since diagnosis; pooled over m = 0..W-1 weighted by the number eligible.
  nat <- ev$natural
  seq_rows <- list()
  for (m in 0:(W - 1)) {
    elig <- (is.na(nat$t_death) | nat$t_death > m) & (is.na(nat$t_pvi) | nat$t_pvi >= m) &
      (is.na(nat$t_male) | nat$t_male > m)
    f <- run(policy("seq_force", m)); b <- run(policy("seq_forbid", m))
    for (o in names(OUTCOMES)) {
      col <- OUTCOMES[[o]]
      d <- as.numeric(by_h(f[[col]][elig], horizon) & f[[col]][elig] > m) -
        as.numeric(by_h(b[[col]][elig], horizon) & b[[col]][elig] > m)
      seq_rows[[length(seq_rows) + 1L]] <- data.table(outcome = o, m = m, n_elig = sum(elig),
                                                      rd = mean(d), var = var(d))
    }
  }
  seq_dt <- rbindlist(seq_rows)
  pooled <- seq_dt[, .(value = sum(rd * n_elig) / sum(n_elig),
                       mcse = sqrt(sum(n_elig^2 * var / n_elig)) / sum(n_elig),
                       n = sum(n_elig)), by = outcome]
  rows[[length(rows) + 1L]] <- pooled[, .(world, scenario, outcome, estimand = "seq_itt", value, mcse, n)]
  for (mm in 0:(W - 1)) {
    x <- seq_dt[m == mm]
    rows[[length(rows) + 1L]] <- x[, .(world, scenario, outcome, estimand = paste0("seq_itt_m", mm),
                                       value = rd, mcse = sqrt(var / n_elig), n = n_elig)]
  }

  curves <- rbindlist(lapply(names(ev), function(p) rbindlist(lapply(c("limb", "male", "disease"), function(o)
    data.table(world = world, scenario = scenario, policy = p, outcome = o, month = seq_len(horizon),
               cuminc = vapply(seq_len(horizon), function(h) mean(by_h(ev[[p]][[OUTCOMES[[o]]]], h)), 0))))))
  if (verbose) message(sprintf("%s/%s: %.1f min", world, scenario,
                               as.numeric(Sys.time() - t0, units = "mins")))
  list(truths = rbindlist(rows), curves = curves)
}