# Algoria Survival Lab 0.1.0 -- synthetic data, time in months
# Requires R and survival. Run in an empty output directory.
if (.Platform$OS.type == 'windows') Sys.setlocale('LC_CTYPE', 'English_United States.utf8')
#' Default educational scenario (time unit: months).
default_parameters <- function() {
  list(n=150L, baseline_risk=.30, hazard_ratio=.65, followup=24,
       dropout_probability=.15, seed=20261008L, tau=12)
}

#' Validate and normalize the public simulation contract.
validate_parameters <- function(parameters) {
  limits <- list(n=c(10,2000), baseline_risk=c(.01,.95), hazard_ratio=c(.25,4),
                 followup=c(1,60), dropout_probability=c(0,.90),
                 seed=c(1,2147483646), tau=c(12,12))
  if (!is.list(parameters)) stop("parameters debe ser una lista.", call.=FALSE)
  for (name in names(limits)) {
    value <- parameters[[name]]
    if (!is.numeric(value) || length(value)!=1L || !is.finite(value) ||
        value < limits[[name]][1] || value > limits[[name]][2] ||
        (name %in% c("n","seed") && value != floor(value))) {
      stop(sprintf("%s: se requiere un valor numérico válido entre %s y %s%s.",
        name, limits[[name]][1], limits[[name]][2],
        if (name %in% c("n","seed")) " (entero)" else ""), call.=FALSE)
    }
  }
  if (length(setdiff(names(parameters), names(limits))))
    stop("parameters contiene campos desconocidos.", call.=FALSE)
  parameters <- parameters[names(limits)]
  parameters$n <- as.integer(parameters$n)
  parameters$seed <- as.integer(parameters$seed)
  for (name in setdiff(names(limits),c("n","seed"))) parameters[[name]] <- as.numeric(parameters[[name]])
  parameters
}

#' Exponential event hazards and independent loss hazard, per month.
generator_rates <- function(parameters) {
  p <- validate_parameters(parameters)
  baseline <- -log1p(-p$baseline_risk)/p$tau
  c(A=baseline, B=baseline*p$hazard_ratio,
    C=-log1p(-p$dropout_probability)/p$tau)
}

#' Known marginal survival and risk; these are not sample estimates.
true_curves <- function(parameters, times=seq(0, parameters$followup, length.out=241L)) {
  rates <- generator_rates(parameters)
  if (!is.numeric(times) || !length(times) || any(!is.finite(times) | times<0))
    stop("times debe contener meses finitos no negativos.", call.=FALSE)
  do.call(rbind, lapply(c("A","B"), function(group) {
    survival <- exp(-rates[[group]]*times)
    data.frame(group=group, time=times, survival=survival,
               risk=-expm1(-rates[[group]]*times), hazard=rates[[group]])
  }))
}

#' Analytical probabilities of the first observed outcome before closure.
observation_probabilities <- function(parameters) {
  p <- validate_parameters(parameters); rates <- generator_rates(p)
  event <- rates[c("A","B")]; loss <- rates[["C"]]; total <- event+loss
  first <- -expm1(-total*p$followup)
  data.frame(group=c("A","B"), event=unname(event/total*first),
             dropout=unname(loss/total*first), administrative=unname(exp(-total*p$followup)))
}
#' Evaluate code under an explicit RNG, restoring the caller even on failure.
with_reproducible_rng <- function(seed, code) {
  previous_kind <- RNGkind()
  had_seed <- exists(".Random.seed", .GlobalEnv, inherits=FALSE)
  previous_seed <- if (had_seed) get(".Random.seed", .GlobalEnv) else NULL
  on.exit({
    do.call(RNGkind, as.list(previous_kind))
    if (had_seed) assign(".Random.seed", previous_seed, .GlobalEnv)
    else if (exists(".Random.seed", .GlobalEnv, inherits=FALSE))
      rm(".Random.seed", envir=.GlobalEnv)
  })
  RNGkind("Mersenne-Twister", "Inversion", "Rejection")
  set.seed(seed)
  force(code)
}

#' Apply independent loss and common administrative closure to latent times.
#' Event wins exact ties; administrative closure wins a loss/closure tie.
observe_times <- function(event_time, dropout_time, followup) {
  valid_time <- function(x) is.numeric(x) && length(x)>0L &&
    !anyNA(x) && all(x>=0) && all(is.finite(x) | x==Inf)
  if (!valid_time(event_time) || !valid_time(dropout_time) ||
      length(event_time)!=length(dropout_time) ||
      !is.numeric(followup) || length(followup)!=1L ||
      !is.finite(followup) || followup<=0)
    stop("Tiempos latentes o seguimiento no válidos.", call.=FALSE)
  event <- as.integer(event_time<=dropout_time & event_time<=followup)
  reason <- ifelse(event==1L, "event",
                   ifelse(dropout_time<followup, "dropout", "administrative"))
  data.frame(time=pmin(event_time,dropout_time,followup), event=event,
             censor_reason=reason, stringsAsFactors=FALSE)
}

#' Generate two independent synthetic populations from common seeded uniforms.
#' Returns latent times for teaching; inference uses time/event only.
simulate_survival <- function(parameters) {
  p <- validate_parameters(parameters); rates <- generator_rates(p)
  with_reproducible_rng(p$seed, {
    group <- factor(rep(c("A","B"), each=p$n), levels=c("A","B"))
    event_time <- -log1p(-runif(2L*p$n))/rates[as.character(group)]
    dropout_time <- if (rates[["C"]]==0) rep(Inf,2L*p$n)
      else -log1p(-runif(2L*p$n))/rates[["C"]]
    latent <- data.frame(id=seq_len(2L*p$n), group=group,
                          event_time=unname(event_time), dropout_time=dropout_time)
    cbind(latent, observe_times(event_time,dropout_time,p$followup))
  })
}
#' Validate the observed single-event right-censoring contract.
validate_survival_data <- function(data, followup) {
  required <- c("id","group","time","event","censor_reason")
  if (!is.data.frame(data) || !nrow(data) || !all(required %in% names(data)))
    stop("Datos observados incompletos.", call.=FALSE)
  if (anyNA(data[required]) || anyDuplicated(data$id) ||
      !setequal(as.character(unique(data$group)),c("A","B")) ||
      !is.numeric(data$time) || any(!is.finite(data$time) | data$time<0 | data$time>followup) ||
      !is.integer(data$event) || any(!data$event %in% c(0L,1L)))
    stop("ID, grupos, tiempos o codificación del evento no válidos.", call.=FALSE)
  reason <- data$censor_reason
  if (any(!reason %in% c("event","dropout","administrative")) ||
      any((data$event==1L) != (reason=="event")) ||
      any(reason=="administrative" & data$time!=followup) ||
      any(reason=="dropout" & data$time>=followup))
    stop("La causa de censura no concuerda con estado y seguimiento.", call.=FALSE)
  data$group <- factor(as.character(data$group),levels=c("A","B"))
  data
}

#' People immediately before each time, including those censored at that time.
risk_table <- function(data, times) {
  if (!is.numeric(times) || !length(times) || any(!is.finite(times) | times<0))
    stop("Tiempos de la tabla de riesgo no válidos.", call.=FALSE)
  do.call(rbind,lapply(c("A","B"),function(group) {
    data.frame(group=group,time=times,
      n_risk=vapply(times,function(time) sum(data$group==group & data$time>=time),integer(1)))
  }))
}

#' Estimate marginal survival with pointwise 95% log-log confidence intervals.
#' Never extend estimates past the final observed time of a group.
analyze_survival <- function(data, parameters) {
  p <- validate_parameters(parameters); d <- validate_survival_data(data,p$followup)
  fit <- survival::survfit(survival::Surv(time,event)~group,data=d,
                           conf.type="log-log",conf.int=.95)
  curve <- data.frame(group=rep(c("A","B"),unname(fit$strata)),
    time=fit$time,survival=fit$surv,lower=fit$lower,upper=fit$upper,
    n_event=fit$n.event,n_censor=fit$n.censor)
  # survival versions differ at S=0/1. Log-log limits are undefined there;
  # never present backend zero-width limits as evidence of certainty.
  boundary <- curve$survival %in% c(0,1)
  curve$lower[boundary] <- NA_real_
  curve$upper[boundary] <- NA_real_
  # The origin is structural; no artificial zero-width confidence band is added.
  origin <- data.frame(group=c("A","B"),time=0,survival=1,lower=NA_real_,
                       upper=NA_real_,n_event=0,n_censor=0)
  curve <- rbind(origin,curve)
  curve <- curve[order(curve$group,curve$time),]; rownames(curve)<-NULL
  target <- min(p$tau,p$followup)
  rates <- generator_rates(p)
  summary <- do.call(rbind,lapply(c("A","B"),function(group) {
    x <- d[d$group==group,]; c <- curve[curve$group==group,]
    supported <- target<=max(x$time)
    at <- if (supported) tail(which(c$time<=target),1) else NA_integer_
    data.frame(group=group,n=nrow(x),events=sum(x$event),
      dropout=sum(x$censor_reason=="dropout"),
      administrative=sum(x$censor_reason=="administrative"),
      censored_fraction=mean(x$event==0L),risk_12=-expm1(-rates[[group]]*12),
      t_star=target,true_survival=exp(-rates[[group]]*target),
      km_survival=if (supported) c$survival[at] else NA_real_,
      lower=if (supported) c$lower[at] else NA_real_,
      upper=if (supported) c$upper[at] else NA_real_,estimable=supported)
  }))
  warnings <- character()
  for (i in seq_len(nrow(summary))) {
    row <- summary[i,]
    if (row$events==0L) warnings <- c(warnings,sprintf(
      "Grupo %s: Sin eventos observados; el intervalo log–log no está disponible.",row$group))
    else if (row$events<5L) warnings <- c(warnings,sprintf(
      "Grupo %s: menos de 5 eventos; la curva es frágil. Este umbral es descriptivo.",row$group))
    if (!row$estimable) warnings <- c(warnings,sprintf(
      "Grupo %s: KM a %.1f meses no estimable; supera el último tiempo observado.",row$group,target))
    if (any(curve$group==row$group & curve$survival==0)) warnings <- c(warnings,sprintf(
      "Grupo %s: KM llega a cero; el intervalo log–log en ese extremo no está disponible.",row$group))
  }
  list(fit=fit,curve=curve,risk=risk_table(d,seq(0,p$followup,length.out=5L)),
       summary=summary,warnings=warnings)
}

#' Freeze one complete result for presentation and exports.
build_result <- function(parameters, data=NULL) {
  started <- proc.time()[["elapsed"]]
  p <- validate_parameters(parameters)
  d <- if (is.null(data)) simulate_survival(p) else validate_survival_data(data,p$followup)
  analysis <- analyze_survival(d,p)
  versions <- lapply(c("survival","shiny","jsonlite"),function(package)
    tryCatch(as.character(utils::packageVersion(package)),error=function(e) NA_character_))
  names(versions) <- c("survival","shiny","jsonlite")
  metadata <- list(schema_version=1L,app_version="0.1.0",source="synthetic",
    time_unit="months",R=R.version.string,platform=R.version$platform,packages=versions,
    seed=p$seed,rng=c("Mersenne-Twister","Inversion","Rejection"),
    created_at=format(Sys.time(),"%Y-%m-%dT%H:%M:%SZ",tz="UTC"),
    computation_ms=1000*(proc.time()[["elapsed"]]-started))
  list(parameters=p,data=d,analysis=analysis,truth=true_curves(p),metadata=metadata)
}
#' Expand right-continuous step values without interpolating across missing values.
step_coordinates <- function(time,value) {
  if (!is.numeric(time) || length(time)!=length(value) || !length(time) ||
      any(!is.finite(time)) || any(diff(time)<0))
    stop("Coordenadas del escalón no válidas.",call.=FALSE)
  if (length(time)==1L) return(data.frame(time=time,value=value))
  data.frame(time=c(time[1],rep(time[-1],each=2)),
             value=rep(value,each=2)[-length(rep(value,each=2))])
}

#' Plot known survival, KM, pointwise intervals and observed censoring marks.
plot_survival <- function(result,compact=par("din")[1]<5.5,language="es") {
  text <- function(es,en) if (identical(language,"en")) en else es
  colors <- c(A="#08737a",B="#062a4e")
  p <- result$parameters
  old <- par(mar=if (compact) c(5.6,4,3,1) else c(6.5,4.5,3.6,1),
             cex=if (compact) .8 else 1,mgp=c(2.6,.8,0)); on.exit(par(old))
  plot(NA,xlim=c(0,p$followup),ylim=c(0,1),xaxs="i",yaxs="i",
       xlab=text("Tiempo desde la entrada (meses)","Time since entry (months)"),ylab=text("Supervivencia: P(T > t)","Survival: P(T > t)"),
       main=if (compact) text("Generadora y Kaplan–Meier","Generating and Kaplan–Meier") else text("Verdad generadora y Kaplan–Meier","Generating truth and Kaplan–Meier"),las=1)
  abline(h=seq(0,1,.2),col="#E2E8ED",lwd=.8)
  for (group in c("A","B")) {
    c <- result$analysis$curve[result$analysis$curve$group==group,]
    # Each finite interval extends to the next time, before the next jump.
    if (nrow(c)>1L) {
      valid <- which(is.finite(c$lower[-nrow(c)]) & is.finite(c$upper[-nrow(c)]))
      if (length(valid)) for (run in split(valid,cumsum(c(TRUE,diff(valid)>1)))) {
        index <- c(run,max(run)+1L)
        low <- step_coordinates(c$time[index],c(c$lower[run],tail(c$lower[run],1)))
        high <- step_coordinates(c$time[index],c(c$upper[run],tail(c$upper[run],1)))
        polygon(c(low$time,rev(high$time)),c(low$value,rev(high$value)),
                col=adjustcolor(colors[[group]],alpha.f=.14),border=NA)
      }
    }
    truth <- result$truth[result$truth$group==group,]
    lines(truth$time,truth$survival,col=colors[[group]],lty=2,lwd=2)
    step <- step_coordinates(c$time,c$survival)
    lines(step$time,step$value,col=colors[[group]],lwd=2)
    censored <- c$n_censor>0
    points(c$time[censored],c$survival[censored],pch=3,col=colors[[group]],cex=.65)
  }
  legend("bottomleft",c(text("A · referencia","A · reference"),text("B · comparación","B · comparison"),text("Generadora","Generating"),text("KM estimada","Estimated KM"),text("Censura","Censoring")),
         col=c(colors,"#52616B","#52616B","#52616B"),lty=c(1,1,2,1,NA),
         pch=c(NA,NA,NA,NA,3),lwd=c(2,2,2,2,NA),bty="n",cex=if (compact) .7 else .80)
  caption <- if (compact) sprintf(text("n=%d/grupo · HR=%.2f","n=%d/group · HR=%.2f"),p$n,p$hazard_ratio) else
    sprintf(text("n=%d/grupo · riesgo A a 12 m=%.0f%% · HR=%.2f · pérdidas latentes a 12 m=%.0f%%","n=%d/group · A risk at 12 mo=%.0f%% · HR=%.2f · latent dropout at 12 mo=%.0f%%"),
            p$n,100*p$baseline_risk,p$hazard_ratio,100*p$dropout_probability)
  foot <- if (compact) sprintf(text("Semilla %d · IC puntuales 95%%","Seed %d · Pointwise 95%% CIs"),p$seed) else
    sprintf(text("Semilla %d · IC puntuales 95%% log–log; límites ausentes sin eventos o en KM=0","Seed %d · Pointwise 95%% log–log CIs; limits unavailable without events or at KM=0"),p$seed)
  mtext(caption,side=1,line=3.8,cex=.75)
  mtext(foot,side=1,line=4.7,cex=.68)
  invisible(result)
}

#' Write a PNG; always close only the device owned by this operation.
draw_survival_png <- function(result,path,width=1440L,height=900L,res=144L,language="es") {
  grDevices::png(path,width=width,height=height,res=res)
  device <- dev.cur(); on.exit(dev.off(device))
  plot_survival(result,language=language)
  invisible(path)
}
parameters <-
list(n = 150L, baseline_risk = 0.3, hazard_ratio = 0.65, followup = 24, 
    dropout_probability = 0.15, seed = 20261008L, tau = 12)
language <-
"en"
result <- build_result(parameters)
write.csv(result$data[c('id','group','time','event','censor_reason')], 'survival-data.csv', row.names=FALSE)
write.csv(result$analysis$summary, 'survival-summary.csv', row.names=FALSE, na='NA')
draw_survival_png(result, 'survival-plot.png', language=language)
print(result$analysis$summary)
print(sessionInfo())