#!/usr/bin/env Rscript
# longevity_rci.R — Component 2 · RCI management analysis
#
# Part of the Seed Longevity Lab's downloadable lab (site v0.17.0,
# longevity.seedsofhawaii.org/scripts). Independent (embeds Component 1). Its accession-baseline stage is the one thing Component 3 depends on.
#
# Standalone by design: the shared cores below are embedded verbatim in every
# script of the suite — no source() of siblings, no packages beyond base R
# (stats::glm is the probit engine).
# R 4.0+. Input: a raw completed-test CSV (one row per completed germination
# test) OR the suite's normalized_assays.csv interchange file.
#
# Same scientific contract as the Python scripts and the website engine, and
# cross-validated against both by the site's test suite. Each component's
# page card at /scripts spells out exactly where behavior differs from the
# archival workflow behind the published catalog.
#
# Counts rule, everywhere: actual number sown/tested + number germinated are
# required. Nothing in this suite ever assumes a sample size.

# ---------------------------------------------------------------------------
# SHARED CORE (R) · statistics helpers — identical block in every script.
# The probit fits themselves use base stats::glm(binomial(link="probit")),
# i.e. the same IRLS the reference workflow used. Separation / extreme fits
# are FLAGGED (warning captured), never silently replaced with a different
# estimator — see the per-component notes on the site's /scripts page.
# ---------------------------------------------------------------------------
DAYS_PER_YEAR <- 365.25
DAYS_PER_MONTH <- 365.25 / 12

wilson_lo <- function(k, n, z = 1.959963984540054) {
  if (is.na(k) || is.na(n) || n <= 0) return(NA_real_)
  p <- k / n
  den <- 1 + z * z / n
  center <- (p + z * z / (2 * n)) / den
  half <- (z / den) * sqrt(p * (1 - p) / n + z * z / (4 * n * n))
  max(0, center - half)
}

glm_probit <- function(g, n, X) {
  # X: data.frame of predictors. Returns list(coef, vcov, deviance, converged,
  # extreme, fitted) or NULL. Extreme-probability warning is CAPTURED as a flag.
  dat <- cbind(data.frame(.g = g, .ng = n - g), X)
  extreme <- FALSE
  fml <- as.formula(paste("cbind(.g, .ng) ~", paste(names(X), collapse = " + ")))
  m <- tryCatch(
    withCallingHandlers(
      glm(fml, family = binomial(link = "probit"), data = dat,
          control = glm.control(epsilon = 1e-10, maxit = 200)),
      warning = function(w) {
        if (grepl("fitted probabilities numerically 0 or 1", conditionMessage(w), fixed = TRUE))
          extreme <<- TRUE
        invokeRestart("muffleWarning")
      }),
    error = function(e) NULL)
  if (is.null(m)) return(NULL)
  list(coef = unname(coef(m)), vcov = vcov(m), deviance = deviance(m),
       converged = isTRUE(m$converged), extreme = extreme,
       separated = !isTRUE(m$converged) || any(abs(coef(m)) > 40, na.rm = TRUE),
       fitted = fitted(m), model = m)
}

crossing_t <- function(intercept, slope, target_prop) {
  if (!is.finite(intercept) || !is.finite(slope) || slope >= 0) return(NA_real_)
  if (!is.finite(target_prop) || target_prop <= 0 || target_prop >= 1) return(NA_real_)
  (qnorm(target_prop) - intercept) / slope
}

crossing_display <- function(raw) {
  if (is.na(raw) || !is.finite(raw)) return("not interpretable")
  if (raw <= 0) return("not interpretable (crossing ≤0 y)")
  sprintf("%.2f y", raw)
}

fmt_rci_year <- function(y) {
  if (abs(y - round(y)) < 5e-3) as.character(round(y)) else sprintf("%.1f", y)
}

# ---------------------------------------------------------------------------
# SHARED CORE (R) · upload/normalization layer — identical block in every
# script; the same contract as the site and the Python lab:
#   one row = one completed germination test · counts REQUIRED (no assumed
#   sample size, ever) · exact time (days/months/years/dates, no rounding) ·
#   is_initial explicit (legacy storage="Initial" recognized) · protocol_key
#   from substrate + treatment with trivial normalization only · drying
#   optional · duplicates flagged, never merged.
# ---------------------------------------------------------------------------
ROLES <- list(
  family = c("family"),
  taxon = c("taxonname", "taxon", "species", "scientificname", "speciesname"),
  accession = c("acquisitionnum", "accession", "accessionnum", "accessionnumber", "accno", "lot", "lotid", "collectionid"),
  assay_id = c("assayid", "testid", "germinationtestid", "assay", "testnumber", "testno"),
  source_institution = c("sourceinstitution", "institution", "seedbank", "source"),
  is_initial = c("isinitial", "initial", "initialtest"),
  storage = c("storedrycode", "storagecondition", "storagecode", "storecode", "storage", "storagetemp", "storagetemperature", "storagetempc"),
  drying = c("dryingtreatment", "drying", "drytreatment", "dryingprotocol"),
  substrate = c("substrate", "medium"),
  treatment = c("germinationtreatment", "treatment", "pretreatment"),
  tested = c("numbersown", "numsown", "sown", "seedsown", "seedssown", "tested", "numbertested", "ntested", "nsown", "n"),
  germinated = c("numbergermin", "numbergerminated", "germinated", "numgerminated", "germin", "germ", "ngerminated"),
  percent = c("%germin", "percentgermination", "percentgerminated", "germinationpercent", "germinationpct", "pctgermin", "percent", "pct", "viability", "viabilitypct"),
  age_days = c("storageagedays", "storagedays", "days", "day", "ageindays"),
  age_months = c("storageagemonths", "storagemonths", "months", "month", "ageinmonths"),
  age_years = c("yearsdifference", "storageageyears", "storageyears", "years", "year", "yr", "yeardiff", "yearsdiff", "ageinyears"),
  age_ambiguous = c("storageage", "age", "time", "interval", "storagetime"),
  collection_date = c("collectiondate", "storagestartdate", "collecteddate", "datecollected", "collectdate"),
  test_date = c("teststartdate", "germinationteststartdate", "testdate", "germinationdate", "startdate", "sowdate", "datesown"),
  notes = c("notes", "comments", "note")
)
TRUTHY <- c("true", "t", "yes", "y", "1", "initial")
FALSY <- c("false", "f", "no", "n", "0", "")

norm_header <- function(h) gsub("[^a-z0-9%]", "", tolower(h))

detect_roles <- function(headers) {
  nh <- norm_header(headers)
  mapping <- list()
  for (i in seq_along(nh)) {
    for (role in names(ROLES)) {
      if (is.null(mapping[[role]]) && nh[i] %in% ROLES[[role]]) {
        mapping[[role]] <- i
        break
      }
    }
  }
  mapping
}

norm_protocol_text <- function(x) {
  s <- tolower(gsub("\\s+", " ", trimws(ifelse(is.na(x), "", x))))
  ifelse(s == "", "(missing)", s)
}
protocol_key <- function(substrate, treatment) {
  s <- norm_protocol_text(substrate); t <- norm_protocol_text(treatment)
  ifelse(s == "(missing)" & t == "(missing)", "unspecified", paste(s, "||", t))
}

normalize_storage <- function(raw) {
  s <- trimws(ifelse(is.na(raw), "", raw))
  l <- gsub("[−–—]", "-", gsub("[°º]", "", gsub("\\s+", "", tolower(s))))
  if (grepl("^init", s, ignore.case = TRUE)) return(list(code = "Initial", original = s, legacy = TRUE))
  if (grepl("^-18c?$", l) || l %in% c("minus18", "minus18c")) return(list(code = "-18C", original = s, legacy = FALSE))
  if (grepl("^\\+?5c?$", l)) return(list(code = "5C", original = s, legacy = FALSE))
  if (grepl("^amb", l)) return(list(code = "Ambient", original = s, legacy = FALSE))
  list(code = s, original = s, legacy = FALSE)
}

parse_date_loose <- function(v) {
  s <- trimws(ifelse(is.na(v), "", v))
  if (s == "") return(NA)
  if (grepl("^[0-9]+(\\.[0-9]+)?$", s)) {
    num <- as.numeric(s)  # Excel serial (days since 1899-12-30)
    if (num > 59 && num < 200000) return(as.Date(num, origin = "1899-12-30"))
    return(NA)
  }
  if (grepl("^\\d{4}-\\d{1,2}-\\d{1,2}", s)) {
    d <- suppressWarnings(as.Date(substr(s, 1, 10), format = "%Y-%m-%d"))
    return(d)
  }
  m <- regmatches(s, regexec("^(\\d{1,2})[/-](\\d{1,2})[/-](\\d{2,4})$", s))[[1]]
  if (length(m) == 4) {
    yr <- as.integer(m[4]); if (yr < 100) yr <- 2000 + yr
    d <- suppressWarnings(as.Date(sprintf("%04d-%02d-%02d", yr, as.integer(m[2]), as.integer(m[3]))))
    return(d)
  }
  suppressWarnings(as.Date(s))
}

numf <- function(v) suppressWarnings(as.numeric(gsub("%", "", ifelse(is.na(v), "", v))))

NORMALIZED_SIGNATURE <- c("row_id", "is_initial", "protocol_key", "storage_age_years", "tested", "germinated")

issue <- function(level, msg) message(sprintf("[%s] %s", toupper(level), msg))

parse_completed_tests <- function(path, time_unit = NULL, strict = TRUE) {
  # Raw completed-test CSV/TSV -> canonical assay data.frame. Blocking issues
  # go to stderr; percent-only input stops analytical scripts (exit 2).
  first <- readLines(path, n = 1, warn = FALSE)
  sep <- if (grepl("\t", first)) "\t" else if (lengths(regmatches(first, gregexpr(";", first))) >
    lengths(regmatches(first, gregexpr(",", first)))) ";" else ","
  raw <- read.csv(path, sep = sep, colClasses = "character", check.names = FALSE,
                  fileEncoding = "UTF-8-BOM", strip.white = TRUE)
  headers <- names(raw)
  mp <- detect_roles(headers)
  err <- FALSE
  for (role in c("taxon", "accession")) {
    if (is.null(mp[[role]])) { issue("error", sprintf("Missing a %s column. Headers: %s", role, paste(headers, collapse = ", "))); err <- TRUE }
  }
  if (is.null(mp$storage) && is.null(mp$is_initial)) {
    issue("error", "Missing a storage-condition column (post-storage rows need one; Initial rows may use is_initial)."); err <- TRUE
  }
  has_counts <- !is.null(mp$tested) && !is.null(mp$germinated)
  has_pct <- !is.null(mp$percent)
  if (!has_counts && !has_pct) { issue("error", "No germination data found — need tested (number sown) + germinated counts."); err <- TRUE }
  percent_only <- (!has_counts) && has_pct
  if (percent_only) {
    issue("block", paste("Counts required for formal modelling: this file has germination percentages but no number",
      "sown/tested and number germinated. A binomial longevity model requires those counts because 80% based on",
      "8/10 seeds carries different statistical information from 80/100. Add the original counts — no sample",
      "size will be assumed."))
    if (strict) quit(status = 2)
  }
  has_days <- !is.null(mp$age_days); has_months <- !is.null(mp$age_months)
  has_years <- !is.null(mp$age_years); has_amb <- !is.null(mp$age_ambiguous)
  has_dates <- !is.null(mp$collection_date) && !is.null(mp$test_date)
  amb_unit <- NULL
  if (!has_days && !has_months && !has_years && has_amb) {
    if (!is.null(time_unit) && time_unit %in% c("days", "months", "years")) {
      amb_unit <- time_unit
      issue("info", sprintf("Ambiguous time column read as %s (your --time-unit selection).", time_unit))
    } else {
      issue("choice", sprintf("The time column \"%s\" doesn't say its unit. Pass --time-unit days|months|years — this script won't guess.", headers[mp$age_ambiguous]))
      quit(status = 1)
    }
  }
  if (!has_days && !has_months && !has_years && !has_amb && !has_dates) {
    issue("error", "No storage-age column found — need days / months / years, or CollectionDate + TestStartDate."); err <- TRUE
  }
  if (err) quit(status = 1)

  g <- function(row, role) if (is.null(mp[[role]])) "" else trimws(raw[row, mp[[role]]])
  out <- list(); rej <- 0; legacy_n <- 0; dated_n <- 0
  for (i in seq_len(nrow(raw))) {
    taxon <- g(i, "taxon"); accession <- g(i, "accession")
    if (taxon == "" || accession == "") { rej <- rej + 1; next }
    st <- normalize_storage(g(i, "storage"))
    is_init <- st$legacy
    if (!is.null(mp$is_initial)) {
      v <- tolower(g(i, "is_initial"))
      if (v %in% TRUTHY) is_init <- TRUE
    }
    if (st$legacy) legacy_n <- legacy_n + 1
    tested <- NA_real_; germ <- NA_real_; pct <- NA_real_
    if (has_counts) {
      tested <- numf(g(i, "tested")); germ <- numf(g(i, "germinated"))
      if (is.na(tested) || is.na(germ) || tested <= 0 || germ < 0 || germ > tested) { rej <- rej + 1; next }
      pct <- 100 * germ / tested
    } else if (has_pct) {
      pct <- min(100, max(0, numf(g(i, "percent"))))
      if (is.na(pct)) { rej <- rej + 1; next }
    }
    age_orig <- NA_real_; age_unit <- NA_character_; days <- NA_real_
    if (has_days && g(i, "age_days") != "") { age_orig <- numf(g(i, "age_days")); age_unit <- "days"; days <- age_orig }
    else if (has_years && g(i, "age_years") != "") { age_orig <- numf(g(i, "age_years")); age_unit <- "years"; days <- age_orig * DAYS_PER_YEAR }
    else if (has_months && g(i, "age_months") != "") { age_orig <- numf(g(i, "age_months")); age_unit <- "months"; days <- age_orig * DAYS_PER_MONTH }
    else if (!is.null(amb_unit) && g(i, "age_ambiguous") != "") {
      age_orig <- numf(g(i, "age_ambiguous")); age_unit <- amb_unit
      days <- age_orig * switch(amb_unit, days = 1, months = DAYS_PER_MONTH, years = DAYS_PER_YEAR)
    }
    cd <- parse_date_loose(g(i, "collection_date")); td <- parse_date_loose(g(i, "test_date"))
    if (is.na(days) && !is.na(cd) && !is.na(td)) {
      dd <- as.numeric(td - cd)  # exact calendar days, never rounded
      if (dd >= 0) { age_orig <- dd; age_unit <- "dates"; days <- dd; dated_n <- dated_n + 1 }
    }
    if (is.na(days) || days < 0) { rej <- rej + 1; next }
    substrate <- g(i, "substrate"); treatment <- g(i, "treatment")
    out[[length(out) + 1]] <- data.frame(
      row_id = sprintf("row-%d", i), assay_id = g(i, "assay_id"),
      family = g(i, "family"), taxon = taxon, accession = accession,
      is_initial = is_init,
      storage_code = if (is_init) "Initial" else st$code,
      storage_original = st$original,
      drying = g(i, "drying"), substrate = substrate, treatment = treatment,
      protocol_key = protocol_key(substrate, treatment),
      protocol_complete = norm_protocol_text(substrate) != "(missing)" && norm_protocol_text(treatment) != "(missing)",
      tested = tested, germinated = germ, pct = pct,
      age_original = age_orig, age_unit_original = age_unit,
      days = days, years = days / DAYS_PER_YEAR,
      collection_date = if (is.na(cd)) "" else format(cd), test_date = if (is.na(td)) "" else format(td),
      stringsAsFactors = FALSE)
  }
  assays <- if (length(out)) do.call(rbind, out) else
    stop("no valid assays parsed")
  # duplicate / same-date multiple tests: flag, never merge
  dkey <- with(assays, paste(taxon, accession, storage_code, protocol_key,
    ifelse(test_date != "", test_date, paste0("age:", days)), tested, germinated, sep = "\u0001"))
  for (k in unique(dkey[duplicated(dkey)])) {
    rows <- assays[dkey == k, ]
    ids <- rows$assay_id
    if (all(ids != "") && !anyDuplicated(ids)) next  # distinct assay IDs = distinct tests
    issue("warn", sprintf("Possible duplicate or multiple same-date tests (%s %s): rows %s preserved — confirm distinct tests vs duplicate records.",
      rows$taxon[1], rows$accession[1], paste(rows$row_id, collapse = ", ")))
  }
  if (legacy_n) issue("info", sprintf("%d row(s) used legacy storage=\"Initial\" — mapped to is_initial=TRUE.", legacy_n))
  if (dated_n) issue("info", sprintf("%d row(s): storage age derived from dates as exact elapsed days (no rounding).", dated_n))
  if (rej) issue("warn", sprintf("%d row(s) rejected (missing identity, invalid counts, or invalid age).", rej))
  assays
}

load_input <- function(path, time_unit = NULL, strict = TRUE) {
  # Raw CSV or the suite's normalized_assays.csv interchange file (the site
  # exports the identical format) — autodetected by signature columns.
  first <- readLines(path, n = 1, warn = FALSE)
  cols <- trimws(strsplit(gsub("^", "", first), ",")[[1]])
  if (all(NORMALIZED_SIGNATURE %in% cols)) {
    d <- read.csv(path, colClasses = "character", check.names = FALSE, fileEncoding = "UTF-8-BOM")
    years <- numf(d$storage_age_years)
    keep <- !is.na(years)
    d <- d[keep, ]; years <- years[keep]
    tested <- numf(d$tested); germ <- numf(d$germinated)
    pct <- numf(d$germination_percent_recomputed)
    pct[is.na(pct) & !is.na(tested) & tested > 0] <-
      (100 * germ / tested)[is.na(pct) & !is.na(tested) & tested > 0]
    return(data.frame(
      row_id = d$row_id, assay_id = d$assay_id, family = d$family,
      taxon = d$taxon, accession = d$accession,
      is_initial = tolower(d$is_initial) %in% TRUTHY,
      storage_code = ifelse(tolower(d$is_initial) %in% TRUTHY, "Initial",
        ifelse(d$storage_condition != "", d$storage_condition, d$storage_condition_original)),
      storage_original = d$storage_condition_original,
      drying = d$drying_treatment, substrate = d$substrate, treatment = d$germination_treatment,
      protocol_key = ifelse(d$protocol_key != "", d$protocol_key, protocol_key(d$substrate, d$germination_treatment)),
      protocol_complete = norm_protocol_text(d$substrate) != "(missing)" & norm_protocol_text(d$germination_treatment) != "(missing)",
      tested = tested, germinated = germ, pct = pct,
      age_original = numf(d$storage_age_original), age_unit_original = d$storage_age_unit_original,
      days = ifelse(is.na(numf(d$storage_age_days)), years * DAYS_PER_YEAR, numf(d$storage_age_days)),
      years = years, collection_date = d$collection_date, test_date = d$test_date,
      stringsAsFactors = FALSE))
  }
  parse_completed_tests(path, time_unit = time_unit, strict = strict)
}

require_counts <- function(assays) {
  if (!any(!is.na(assays$tested))) {
    issue("block", "Counts required — this input has no sown/germinated counts; no sample size is ever assumed.")
    quit(status = 2)
  }
}

# ---------------------------------------------------------------------------
# SHARED CORE (R) · RCI management analysis — identical block where used.
# The PUBLIC upload rule (the site's specification), matching the Python lab
# and the website engine exactly:
#   baseline = highest observed germination in [earliest age, +5 y] (no
#   archival two-point exception) · RCI P70 = 0.70 × baseline · max-at-age
#   replicate rule (RCI branch only) · per-storage latest passing age with
#   GLM + Wilson support · Overall observed RCI SEPARATE from RCI at
#   recommended storage · freeze flags incl. the only-−18C FALSE override ·
#   censoring-honest a–b / a+ / <b notation, fractional precision kept.
# ---------------------------------------------------------------------------
BASELINE_WINDOW_YEARS <- 5; THRESHOLD_FRACTION <- 0.7
FREEZE_RATIO <- 0.5; MARGIN_YEARS <- 5; RCI_EPS <- 1e-8
TEMP_LABEL <- c("-18C" = "−18 °C", "5C" = "5 °C", "Ambient" = "Ambient")
tlab <- function(x) ifelse(x %in% names(TEMP_LABEL), TEMP_LABEL[x], x)

rci_rows <- function(assays) {
  d <- assays[!is.na(assays$pct), ]
  data.frame(family = d$family, taxon = d$taxon, accession = d$accession,
             storage = ifelse(d$is_initial, "Initial", d$storage_code),
             sown = d$tested, germ = d$germinated, pct = d$pct, years = d$years,
             protocol_key = d$protocol_key, stringsAsFactors = FALSE)
}

collapse_by_year <- function(rows, by_storage = TRUE) {
  # max-at-age: highest germination proportion; ties -> larger n (RCI branch only)
  d <- rows[!is.na(rows$sown) & rows$sown > 0 & !is.na(rows$germ) & rows$germ >= 0, ]
  if (!nrow(d)) return(d)
  key <- paste(d$taxon, d$accession, if (by_storage) d$storage else "", d$years, sep = "\u0001")
  keep <- unlist(lapply(split(seq_len(nrow(d)), key), function(idx) {
    p <- d$germ[idx] / d$sown[idx]
    idx[order(-p, -d$sown[idx])][1]
  }), use.names = FALSE)
  out <- d[sort(keep), ]
  out$pct <- 100 * out$germ / out$sown
  out
}

early_baselines <- function(rows) {
  out <- list()
  for (k in unique(paste(rows$family, rows$taxon, rows$accession, sep = "\u0001"))) {
    p <- strsplit(k, "\u0001")[[1]]
    d <- rows[rows$family == p[1] & rows$taxon == p[2] & rows$accession == p[3], ]
    agg <- aggregate(pct ~ years, d, max)
    agg <- agg[order(agg$years), ]
    upper <- agg$years[1] + BASELINE_WINDOW_YEARS
    cand <- agg[agg$years <= upper, ]
    pick <- cand[order(-cand$pct, cand$years), ][1, ]
    out[[k]] <- list(baseline_year = pick$years, baseline_viab = pick$pct,
                     threshold = THRESHOLD_FRACTION * pick$pct)
  }
  out
}

glm_year_ci <- function(sub) {
  # site support-GLM: ≥3 rows, ≥2 distinct years; separation flagged not faked
  if (nrow(sub) < 3 || length(unique(sub$year)) < 2) return(NULL)
  m <- glm_probit(sub$germ, sub$sown, data.frame(year = sub$year))
  if (is.null(m) || m$separated) return(NULL)
  m
}
predict_lo <- function(m, t, z = 1.959963984540054) {
  V <- m$vcov
  var_eta <- V[1, 1] + 2 * t * V[1, 2] + t * t * V[2, 2]
  if (is.na(var_eta) || var_eta < 0) return(NA_real_)
  eta <- m$coef[1] + m$coef[2] * t
  pnorm(eta - z * sqrt(var_eta))
}

per_storage_intervals <- function(rows, baselines) {
  post <- rows[rows$storage != "Initial", ]
  coll <- collapse_by_year(post, by_storage = TRUE)
  out <- list()
  for (k in unique(paste(post$family, post$taxon, post$accession, post$storage, sep = "\u0001"))) {
    p <- strsplit(k, "\u0001")[[1]]
    bl <- baselines[[paste(p[1], p[2], p[3], sep = "\u0001")]]
    if (is.null(bl) || !(bl$baseline_viab > 0)) next
    ser <- post[post$family == p[1] & post$taxon == p[2] & post$accession == p[3] & post$storage == p[4], ]
    yrs <- unique(ser$years)
    if (length(yrs) + as.integer(!(bl$baseline_year %in% yrs)) < 2) next
    pass <- ser$years[ser$pct >= bl$threshold]
    latest <- if (length(pass)) max(pass) else NA_real_
    sub <- coll[coll$family == p[1] & coll$taxon == p[2] & coll$accession == p[3] & coll$storage == p[4], ]
    sub <- data.frame(year = sub$years, sown = sub$sown, germ = sub$germ, pct = sub$pct)
    sub <- sub[order(sub$year), ]
    glm_support <- NA; wilson_support <- NA; glm_lo_v <- NA_real_; wilson_lo_v <- NA_real_
    if (!is.na(latest)) {
      m <- glm_year_ci(sub)
      if (!is.null(m)) {
        lo <- predict_lo(m, latest)
        if (!is.na(lo)) { glm_lo_v <- lo * 100; glm_support <- glm_lo_v >= bl$threshold }
      }
      at <- sub[sub$year == latest, ]
      if (!nrow(at) && nrow(sub)) at <- sub[which.min(abs(sub$year - latest)), ]
      if (nrow(at) == 1 && !is.na(at$sown) && at$sown > 0) {
        wilson_lo_v <- wilson_lo(at$germ, at$sown) * 100
        wilson_support <- wilson_lo_v >= bl$threshold
      }
    }
    out[[length(out) + 1]] <- list(family = p[1], taxon = p[2], accession = p[3], storage = p[4],
      baseline_year = bl$baseline_year, baseline_viab = bl$baseline_viab, threshold = bl$threshold,
      n_years = length(yrs), max_year_meeting_p70 = latest,
      glm_support = glm_support, wilson_support = wilson_support,
      glm_lo = glm_lo_v, wilson_lo = wilson_lo_v, series = sub)
  }
  out
}

species_rollup <- function(rows, baselines) {
  acc_max <- list()
  for (k in unique(paste(rows$family, rows$taxon, rows$accession, sep = "\u0001"))) {
    p <- strsplit(k, "\u0001")[[1]]
    bl <- baselines[[k]]
    if (is.null(bl) || !(bl$baseline_viab > 0)) next
    d <- rows[rows$family == p[1] & rows$taxon == p[2] & rows$accession == p[3], ]
    if (length(unique(d$years)) < 2) next
    pass <- d$years[d$pct >= bl$threshold]
    acc_max[[length(acc_max) + 1]] <- list(family = p[1], taxon = p[2], accession = p[3],
      baseline_year = bl$baseline_year, baseline_viab = bl$baseline_viab,
      threshold = bl$threshold, acc_interval = if (length(pass)) max(pass) else NA_real_)
  }
  species <- list()
  spk <- unique(vapply(acc_max, function(a) paste(a$family, a$taxon, sep = "\u0001"), ""))
  for (k in spk) {
    accs <- Filter(function(a) paste(a$family, a$taxon, sep = "\u0001") == k, acc_max)
    withi <- Filter(function(a) !is.na(a$acc_interval), accs)
    p <- strsplit(k, "\u0001")[[1]]
    if (!length(withi)) { species[[length(species) + 1]] <- list(family = p[1], taxon = p[2], interval = NA_real_, accs = accs); next }
    interval <- max(vapply(withi, function(a) a$acc_interval, 0))
    ord <- order(-vapply(withi, function(a) a$acc_interval, 0),
                 -vapply(withi, function(a) a$baseline_viab, 0),
                 vapply(withi, function(a) as.character(a$accession), ""))
    chosen <- withi[[ord[1]]]
    species[[length(species) + 1]] <- list(family = p[1], taxon = p[2], interval = interval,
      supporting = sum(vapply(withi, function(a) a$acc_interval + RCI_EPS >= interval, FALSE)),
      total_qualifying = length(unique(vapply(accs, function(a) a$accession, ""))),
      chosen_accession = chosen$accession, chosen_threshold = chosen$threshold,
      chosen_baseline_year = chosen$baseline_year, chosen_baseline_viab = chosen$baseline_viab,
      accs = accs)
  }
  list(acc_max = acc_max, species = species)
}

rci_notation <- function(series, threshold) {
  s <- series[order(series$year), ]
  hits <- s$year[s$pct >= threshold]
  if (!length(hits)) {
    fails <- s$year[s$pct < threshold]
    return(if (length(fails)) paste0("<", fmt_rci_year(fails[1])) else "<1")
  }
  a <- max(hits)
  after <- s$year[s$year > a & s$pct < threshold]
  if (!length(after)) return(paste0(fmt_rci_year(a), "+"))
  paste0(fmt_rci_year(a), "–", fmt_rci_year(min(after)))
}

run_rci_pipeline <- function(assays) {
  rows <- rci_rows(assays)
  baselines <- early_baselines(rows)
  intervals <- per_storage_intervals(rows, baselines)
  ru <- species_rollup(rows, baselines)
  species <- ru$species

  iv <- function(f) vapply(intervals, function(it) it[[f]], if (f %in% c("family", "taxon", "accession", "storage")) "" else 0)
  # freeze flags per species
  freeze <- list()
  for (k in unique(paste(iv("family"), iv("taxon"), sep = "\u0001"))) {
    p <- strsplit(k, "\u0001")[[1]]
    its <- Filter(function(it) it$family == p[1] && it$taxon == p[2], intervals)
    m18 <- vapply(Filter(function(it) it$storage == "-18C" && !is.na(it$max_year_meeting_p70), its), function(it) it$max_year_meeting_p70, 0)
    m5 <- vapply(Filter(function(it) it$storage == "5C" && !is.na(it$max_year_meeting_p70), its), function(it) it$max_year_meeting_p70, 0)
    R5 <- if (length(m5)) max(m5) else NA_real_; R18 <- if (length(m18)) max(m18) else NA_real_
    ci5 <- any(vapply(its, function(it) it$storage == "5C" && (isTRUE(it$glm_support) || isTRUE(it$wilson_support)), FALSE))
    ci18 <- any(vapply(its, function(it) it$storage == "-18C" && (isTRUE(it$glm_support) || isTRUE(it$wilson_support)), FALSE))
    ratio <- !is.na(R5) && !is.na(R18) && R5 > 0 && R18 / R5 <= FREEZE_RATIO
    citrig <- ci5 && !ci18 && !is.na(R5) && !is.na(R18) && R5 > R18
    fz <- if (is.na(R5) && is.na(R18)) list(flag = NA, reason = "no 5C or -18C data")
      else if (is.na(R5)) list(flag = NA, reason = "no 5C data")
      else if (is.na(R18)) list(flag = NA, reason = "no -18C data")
      else if (ratio) list(flag = TRUE, reason = sprintf("-18C interval (%s y) <= 50%% of 5C (%s y)", fmt_rci_year(R18), fmt_rci_year(R5)))
      else if (citrig) list(flag = TRUE, reason = "CI supports 5C at its latest; 5C interval > -18C")
      else list(flag = FALSE, reason = "no evidence of freeze sensitivity")
    fz$R5 <- R5; fz$R18 <- R18
    freeze[[k]] <- fz
  }

  by_storage_species <- list()
  for (it in intervals) {
    if (is.na(it$max_year_meeting_p70) || !(it$max_year_meeting_p70 > 0)) next
    k <- paste(it$family, it$taxon, it$storage, sep = "\u0001")
    if (is.null(by_storage_species[[k]]) || it$max_year_meeting_p70 > by_storage_species[[k]])
      by_storage_species[[k]] <- it$max_year_meeting_p70
  }

  coll_all <- collapse_by_year(rows, by_storage = TRUE)

  # potential supporters: within-5y margin + per-storage GLM point/CI probes
  # at the species interval (the union rule behind the published counts)
  potential_n <- list()
  for (sp in species) {
    if (is.na(sp$interval)) next
    pot <- character(0)
    for (a in sp$accs) {
      if (!is.na(a$acc_interval) && a$acc_interval + RCI_EPS >= sp$interval) next
      if (!is.na(a$acc_interval) && a$acc_interval + MARGIN_YEARS >= sp$interval) { pot <- c(pot, a$accession); next }
      hit <- FALSE
      for (st in c("-18C", "5C", "Ambient")) {
        ser <- coll_all[coll_all$taxon == sp$taxon & coll_all$accession == a$accession & coll_all$storage == st, ]
        if (!nrow(ser)) next
        sub <- data.frame(year = ser$years, sown = ser$sown, germ = ser$germ)
        sub <- sub[order(sub$year), ]
        m <- glm_year_ci(sub)
        if (is.null(m)) next
        pr <- pnorm(m$coef[1] + m$coef[2] * sp$interval) * 100
        lo <- predict_lo(m, sp$interval)
        if ((!is.na(pr) && pr >= a$threshold) || (!is.na(lo) && lo * 100 >= a$threshold)) hit <- TRUE
      }
      if (hit) pot <- c(pot, a$accession)
    }
    potential_n[[paste(sp$family, sp$taxon, sep = "\u0001")]] <- length(unique(pot))
  }
  raw_acc_n <- tapply(rows$accession, paste(rows$family, rows$taxon, sep = "\u0001"),
                      function(x) length(unique(x)))

  pub <- list()
  for (sp in species) {
    if (is.na(sp$interval)) next
    k <- paste(sp$family, sp$taxon, sep = "\u0001")
    fz <- freeze[[k]]; if (is.null(fz)) fz <- list(flag = NA, reason = "no per-storage intervals", R5 = NA, R18 = NA)
    anchor <- NULL; anchor_years <- NA_real_
    for (st in c("-18C", "5C", "Ambient")) {
      v <- by_storage_species[[paste(sp$family, sp$taxon, st, sep = "\u0001")]]
      if (is.null(v)) next
      if (is.na(anchor_years) || v > anchor_years) { anchor <- st; anchor_years <- v }
    }
    if (is.null(anchor)) {
      cand <- rows[rows$taxon == sp$taxon & rows$accession == sp$chosen_accession &
                   rows$storage != "Initial" & rows$years == sp$interval, ]
      if (nrow(cand)) anchor <- cand$storage[order(-cand$pct)][1]
    }
    if (is.null(anchor)) {
      fb <- rows[rows$taxon == sp$taxon & rows$accession == sp$chosen_accession &
                 rows$storage != "Initial" & rows$years <= sp$interval & rows$pct >= sp$chosen_threshold, ]
      if (nrow(fb)) anchor <- fb$storage[order(-fb$years, -fb$pct, fb$storage)][1]
    }
    # only-−18C NA→FALSE override (matches the reference workflow)
    if (is.na(fz$flag) && identical(anchor, "-18C")) {
      fz$flag <- FALSE
      fz$reason <- sprintf("only -18C data; species interval = %s y at -18C", fmt_rci_year(sp$interval))
    }
    override <- isTRUE(!fz$flag) && fz$reason == "no evidence of freeze sensitivity"
    rec <- if (override) "-18C" else anchor
    behavior <- if (identical(rec, "-18C")) "Orthodox" else if (identical(rec, "5C")) "Freeze-sensitive" else "Unknown"
    aser <- collapse_by_year(rows[rows$taxon == sp$taxon & rows$accession == sp$chosen_accession &
                                  rows$storage == (if (is.null(anchor)) "" else anchor), ], by_storage = FALSE)
    aser <- data.frame(year = aser$years, pct = aser$pct)
    notation <- if (nrow(aser)) rci_notation(aser, sp$chosen_threshold) else NA_character_
    rec_int <- if (!is.null(rec)) by_storage_species[[paste(sp$family, sp$taxon, rec, sep = "\u0001")]] else NULL
    rec_driver <- NA_character_
    rec_disp <- if (!is.null(rec_int)) {
      rits <- Filter(function(it) it$family == sp$family && it$taxon == sp$taxon &&
                       it$storage == rec && !is.na(it$max_year_meeting_p70), intervals)
      if (length(rits)) {
        rord <- order(-vapply(rits, function(it) it$max_year_meeting_p70, 0),
                      -vapply(rits, function(it) it$baseline_viab, 0))
        rec_driver <- rits[[rord[1]]]$accession
        rci_notation(rits[[rord[1]]]$series, rits[[rord[1]]]$threshold)
      } else paste0(fmt_rci_year(rec_int), "+")
    } else if (!is.null(rec)) sprintf("not yet demonstrated at %s", tlab(rec)) else NA_character_

    # rebound/failure transparency: a failure before the latest pass stays visible
    reb <- vapply(Filter(function(it) it$family == sp$family && it$taxon == sp$taxon &&
                           !is.na(it$max_year_meeting_p70) &&
                           any(it$series$year < it$max_year_meeting_p70 & it$series$pct < it$threshold),
                         intervals),
                  function(it) paste0(it$accession, " @ ", tlab(it$storage)), "")
    # most recent germ at the RCI (exact year, else latest earlier passing)
    recent <- NA_real_
    ex <- aser[abs(aser$year - sp$interval) < 1e-9, ]
    if (nrow(ex) == 1) recent <- ex$pct else {
      okr <- aser[aser$year <= sp$interval & aser$pct >= sp$chosen_threshold, ]
      if (nrow(okr)) recent <- okr$pct[order(-okr$year, -okr$pct)][1]
    }
    rec_reason <- if (override) sprintf("no evidence of freeze sensitivity (R5=%s y, R18=%s y) → −18 °C publication rule",
        ifelse(is.na(fz$R5), "–", fmt_rci_year(fz$R5)), ifelse(is.na(fz$R18), "–", fmt_rci_year(fz$R18)))
      else if (isTRUE(fz$flag)) sprintf("freeze-sensitivity evidence: %s", fz$reason)
      else if (!is.null(anchor)) sprintf("longest demonstrated interval under %s", tlab(anchor))
      else "attribution fallback from the driving accession"
    pot_n <- potential_n[[k]]; if (is.null(pot_n)) pot_n <- 0L
    pub[[length(pub) + 1]] <- data.frame(
      family = sp$family, taxon = sp$taxon,
      recommended_storage = if (is.null(rec)) NA else tlab(rec),
      storage_behavior = behavior,
      recommendation_reason = rec_reason,
      overall_observed_rci_years_raw = sp$interval,
      overall_observed_rci_display = notation,
      overall_driver_accession = sp$chosen_accession,
      overall_driver_storage = if (is.null(anchor)) NA else tlab(anchor),
      rci_at_recommended_years_raw = if (is.null(rec_int)) NA_real_ else rec_int,
      rci_at_recommended_display = rec_disp,
      rci_at_recommended_driver_accession = rec_driver,
      rebound_failure_flags = paste(sort(unique(reb)), collapse = "; "),
      accession_baseline_percent = round(sp$chosen_baseline_viab, 1),
      rci_p70_threshold_percent = round(sp$chosen_threshold, 2),
      fresh_seeds_baseline = if (sp$chosen_baseline_year == 0) "Yes" else "No",
      recent_germ_at_rci = if (is.na(recent)) NA_real_ else round(recent, 1),
      supporting_accessions = sp$supporting,
      potentially_supporting_accessions = pot_n,
      supporting_published = sp$supporting + pot_n,
      rci_eligible_accessions = sp$total_qualifying,
      raw_accessions = as.integer(raw_acc_n[[k]]),
      override_to_minus18 = override,
      R5 = if (is.na(fz$R5)) NA_real_ else fz$R5, R18 = if (is.na(fz$R18)) NA_real_ else fz$R18,
      freeze_sensitive = fz$flag, freeze_reason = fz$reason,
      stringsAsFactors = FALSE)
  }
  publication <- if (length(pub)) do.call(rbind, pub) else data.frame()
  if (nrow(publication)) publication <- publication[order(publication$family, publication$taxon), ]
  list(rows = rows, baselines = baselines, intervals = intervals,
       species = species, publication = publication)
}

# ============================ command line ============================
main <- function() {
  args <- commandArgs(trailingOnly = TRUE)
  if (!length(args)) { cat("usage: Rscript longevity_rci.R INPUT.csv [-o OUTDIR] [--time-unit U]\n"); quit(status = 1) }
  outdir <- if ("-o" %in% args) args[which(args == "-o") + 1] else "."
  tu <- if ("--time-unit" %in% args) args[which(args == "--time-unit") + 1] else NULL
  a <- load_input(args[1], time_unit = tu)
  require_counts(a)
  res <- run_rci_pipeline(a)
  dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
  p1 <- file.path(outdir, "rci_results.csv")
  write.csv(res$publication, p1, row.names = FALSE, na = "")
  iv <- do.call(rbind, lapply(res$intervals, function(it) data.frame(
    family = it$family, taxon = it$taxon, accession = it$accession, storage = it$storage,
    baseline_year = it$baseline_year, baseline_viab = it$baseline_viab, threshold = it$threshold,
    n_years = it$n_years, max_year_meeting_p70 = it$max_year_meeting_p70,
    glm_support = it$glm_support, wilson_support = it$wilson_support,
    glm_lo = it$glm_lo, wilson_lo = it$wilson_lo)))
  p2 <- file.path(outdir, "rci_accession_intervals.csv")
  write.csv(iv, p2, row.names = FALSE, na = "")
  cat(sprintf("%s: %d species\n%s: %d accession x storage series (feeds longevity_formal.R --thresholds)\n",
      p1, nrow(res$publication), p2, if (is.null(iv)) 0L else nrow(iv)))
}
main()