# Simulated randomized trials (plan section 8, D18).
#
# Realistic RCT: hospital-based (OBL = 0) physicians in the top third of skill and the
# bottom half of style; 30% of their eligible patients consent; n = 240.
# Ideal pragmatic RCT: all physicians, everyone consents, large n.
# Assignment A: PVI at a scheduled month 0..2 (uniform); 10% refuse and follow usual care.
# Assignment B: no protocol PVI in the window; crossover happens through the DGP's own
# PVI process (usual care), and usual care continues after the window.

RCT_DESIGNS <- list(
  realistic = list(n = 240L, consent = 0.30, select = TRUE),
  ideal = list(n = NA, consent = 1.00, select = FALSE)
)

#' Physicians a realistic trial would recruit.
trial_physicians <- function(phys) {
  phys[obl == 0 & skill >= quantile(skill, 2 / 3) & style <= median(style), phys_id]
}

#' Draw a trial population from `pop` (patients at diagnosis) and its random assignment.
#' @return list(pop = population restricted to enrolled patients, arm, refuse, sched).
enrol_trial <- function(pop, design = c("realistic", "ideal"), n = NULL, seed = 1L) {
  design <- match.arg(design)
  spec <- RCT_DESIGNS[[design]]
  set.seed(seed)
  pat <- pop$pat
  if (spec$select) pat <- pat[phys_id %in% trial_physicians(pop$phys)]
  pat <- pat[runif(.N) < spec$consent]
  n <- if (!is.null(n)) n else if (is.na(spec$n)) nrow(pat) else spec$n
  if (nrow(pat) < n) stop("not enough consenting patients: ", nrow(pat), " < ", n)
  pat <- pat[sort(sample.int(.N, n))]
  pat[, id := seq_len(.N)]
  arm <- sample(rep(c("A", "B"), length.out = n))
  list(pop = list(phys = pop$phys, pat = pat), arm = arm,
       refuse = runif(n) < 0.10, sched = sample(0:2, n, replace = TRUE))
}

#' The trial protocol as a patient-specific policy.
trial_policy <- function(trial, assign = NULL) {
  arm <- if (is.null(assign)) trial$arm else rep(assign, length(trial$arm))
  force <- arm == "A" & !trial$refuse
  policy("custom", decide = function(t, natural, idx) {
    go <- natural
    f <- force[idx]
    go[f] <- t >= trial$sched[idx][f]       # PVI at the scheduled month, or the first month eligible after it
    go
  })
}

#' Run one simulated trial and analyse it: ITT (by assignment, unadjusted) and naive
#' as-treated (PVI in the window vs not, ignoring assignment).
run_rct <- function(par, design = "realistic", rep = 1L, outcomes = c(male = "t_male", disease = "t_disease"),
                    source_n = 60000L) {
  seed <- replicate_seed(par$world, paste0(par$scenario, "|rct|", design), rep)
  pop <- make_population(par, n = source_n, seed = seed)
  trial <- enrol_trial(pop, design, n = if (design == "ideal") par$size$n_patients else NULL, seed = seed + 1L)
  ev <- simulate_course(trial$pop, par, trial_policy(trial), seed = seed + 2L, record = FALSE)$events
  W <- par$size$window
  ev[, `:=`(phys_id = trial$pop$pat$phys_id, assigned = trial$arm,
            treated_win = ifelse(!is.na(t_pvi) & t_pvi < W, "early", "none"))]
  rbindlist(lapply(names(outcomes), function(o) {
    p <- ev[, .(id, phys_id, assigned, treated_win, t_out = get(outcomes[[o]]), t_death)]
    itt <- rd_from_aj(aj_cif(expand_followup(copy(p)[, `:=`(arm = assigned, w = 1)], 0L)), "A", "B")
    at <- rd_from_aj(aj_cif(expand_followup(copy(p)[, `:=`(arm = treated_win, w = 1)], 0L)), "early", "none")
    rbind(cbind(analysis = "itt", itt), cbind(analysis = "as_treated", at))[
      , `:=`(outcome = o, design = paste0("rct_", design), n = nrow(ev),
             crossover_B = mean(ev$treated_win[ev$assigned == "B"] == "early"),
             nonadherence_A = mean(ev$treated_win[ev$assigned == "A"] == "none"))]
  }))
}

#' ITT truth for a trial design: everyone in a large trial-eligible population assigned
#' A vs everyone assigned B, same protocol (refusal, scheduling, natural crossover).
#' Consent is random, so it does not change the trial population and is skipped here.
rct_truth <- function(par, design = "realistic", n = 1e6, seed = 20261029L, horizon = 60L) {
  pop <- truth_population(par, n = n, seed = seed)
  pat <- if (design == "realistic") pop$pat[phys_id %in% trial_physicians(pop$phys)] else pop$pat
  set.seed(seed + 1L)
  k <- nrow(pat)
  trial <- list(pop = list(phys = pop$phys, pat = copy(pat)[, id := seq_len(.N)]),
                arm = rep("A", k), refuse = runif(k) < 0.10, sched = sample(0:2, k, replace = TRUE))
  runs <- lapply(c("A", "B"), function(a)
    simulate_course(trial$pop, par, trial_policy(trial, assign = a), seed = seed + 2L, record = FALSE)$events)
  rbindlist(lapply(names(OUTCOMES), function(o) {
    d <- as.numeric(by_h(runs[[1]][[OUTCOMES[[o]]]], horizon)) - as.numeric(by_h(runs[[2]][[OUTCOMES[[o]]]], horizon))
    data.table(world = par$world, scenario = par$scenario, outcome = o,
               estimand = paste0("rct_", design, "_itt"), value = mean(d), mcse = sd(d) / sqrt(length(d)),
               n = length(d))
  }))
}