# Aggregate FDA counts and deterministic synthetic examples; no patient data.
response_summary <- function(x, n) {
  integer_scalar <- function(z) is.numeric(z) && length(z) == 1L &&
    is.finite(z) && z == floor(z)
  if (!integer_scalar(x) || !integer_scalar(n) || n <= 0 || x < 0 || x > n) {
    stop("x and n must be finite integer counts with 0 <= x <= n and n > 0.")
  }
  ci <- stats::binom.test(x, n, conf.level = 0.95)$conf.int
  data.frame(responses = x, total = n, estimate = x / n,
             lower95 = unname(ci[1]), upper95 = unname(ci[2]))
}

missing_scenarios <- function(response) {
  if (!is.numeric(response) || length(response) == 0L ||
      any(!is.na(response) & !(response %in% c(0, 1)))) {
    stop("response must be a nonempty numeric vector containing only 0, 1 or NA.")
  }
  known <- sum(response == 1, na.rm = TRUE)
  missing <- sum(is.na(response))
  completed <- 0:missing
  data.frame(missing_responders_assumed = completed, known_responders = known,
             missing = missing, total = length(response),
             estimate = (known + completed) / length(response))
}

run_exercise <- function(out = "response-results") {
  dir.create(out, recursive = TRUE, showWarnings = FALSE)
  published <- response_summary(48, 77)
  published$provenance <- "FDA CIS efficacy population, 2024; not papillary-only"
  write.csv(published, file.path(out, "fda-response-recomputed.csv"), row.names = FALSE)
  denominators <- data.frame(
    quantity = c("CR_any_time", "duration_ge12_among_responders",
                 "duration_ge24_among_responders", "documented_ge12_over77",
                 "documented_ge24_over77"),
    numerator = c(48, 28, 19, 28, 19), denominator = c(77, 48, 48, 77, 77))
  denominators$proportion <- denominators$numerator / denominators$denominator
  denominators$provenance <- "FDA 2024 aggregate counts; last two rows are descriptive arithmetic, not KM"
  write.csv(denominators, file.path(out, "fda-denominators.csv"), row.names = FALSE)
  synthetic <- data.frame(id = sprintf("SYN-%03d", 1:100),
                          response = c(rep(1, 60), rep(0, 30), rep(NA_real_, 10)))
  write.csv(synthetic, file.path(out, "synthetic-response.csv"), row.names = FALSE, na = "")
  sensitivity <- missing_scenarios(synthetic$response)
  write.csv(sensitivity, file.path(out, "synthetic-missingness.csv"), row.names = FALSE)
  precision <- do.call(rbind, lapply(c(20, 80, 200), function(n) response_summary(0.6 * n, n)))
  write.csv(precision, file.path(out, "synthetic-precision.csv"), row.names = FALSE)
  for (lang in c("es", "en")) {
    es <- lang == "es"
    grDevices::png(file.path(out, paste0("precision-", lang, ".png")),
                   width = 1500, height = 850, res = 150)
    graphics::par(mar = c(6.5, 6, 4, 2), family = "sans")
    graphics::plot(100 * precision$estimate, 3:1, xlim = c(0, 100), ylim = c(0.5, 3.5),
                   yaxt = "n", pch = 19, col = "#0b5563", cex = 1.5,
                   xlab = "", ylab = "",
                   main = if (es) "Misma proporci\u00f3n, distinta precisi\u00f3n" else "Same proportion, different precision")
    graphics::axis(2, at = 3:1, labels = paste0(precision$responses, "/", precision$total), las = 1)
    graphics::segments(100 * precision$lower95, 3:1, 100 * precision$upper95, 3:1,
                      col = "#0b5563", lwd = 3)
    graphics::abline(v = 60, lty = 3, col = "#687d87")
    graphics::mtext(if (es) "Respuesta e IC exacto del 95% (%)" else
                     "Response and exact 95% CI (%)", side = 1, line = 3)
    # Unicode escapes keep source() portable when Windows uses a non-UTF-8 locale.
    graphics::mtext(if (es) "Tres muestras sint\u00e9ticas independientes; no son resultados de Anktiva." else
                     "Three independent synthetic samples; these are not Anktiva results.", side = 1, line = 5, cex = 0.75)
    grDevices::dev.off()
  }
  capture.output(sessionInfo(), file = file.path(out, "sessionInfo.txt"))
  cat("FDA aggregate recomputation:\n")
  print(published)
  cat("Synthetic complete-case proportion:", 60 / 90, "\n")
  cat("Synthetic sample completion bounds:", range(sensitivity$estimate), "\n")
  invisible(list(published = published, denominators = denominators,
                 sensitivity = sensitivity, precision = precision))
}