#!/usr/bin/env Rscript
# longevity_formal.R — Component 3 · formal protocol-specific survival models
#
# Part of the Seed Longevity Lab's downloadable lab (site v0.17.0,
# longevity.seedsofhawaii.org/scripts). Embeds Components 1–2; the RCI-P70 bridge threshold is its genuine contingency on the RCI baseline stage (recomputed internally, or chained via --thresholds).
#
# 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)
}

# ---------------------------------------------------------------------------
# SHARED CORE (R) · formal protocol-specific survival — identical block where
# used. Trajectory = taxon × accession × storage × protocol_key (+ drying
# only when supplied AND variable); matched same-protocol Initial (≤1 y)
# pooled once at t = 0, never borrowed across protocols; late Initials
# excluded AND reported; same-time rows pooled by summed counts (identical
# binomial likelihood); plateau + early-rise QC; bridge tiers; crossings
# kept raw AND displayed safely (≤0 y is never "negative longevity").
# ---------------------------------------------------------------------------
INITIAL_MAX_YEAR <- 1

group_formal_trajectories <- function(assays) {
  counted <- assays[!is.na(assays$tested) & !is.na(assays$germinated), ]
  dry_var <- FALSE
  dd <- counted[counted$drying != "", ]
  if (nrow(dd)) {
    kk <- paste(dd$taxon, dd$accession, dd$storage_code, dd$protocol_key, sep = "\u0001")
    dry_var <- any(vapply(split(tolower(trimws(dd$drying)), kk), function(x) length(unique(x)) > 1, FALSE))
  }
  dryk <- function(d) if (dry_var) tolower(trimws(ifelse(d == "", "(none)", d))) else ""
  post <- counted[!counted$is_initial, ]
  init <- counted[counted$is_initial & counted$years <= INITIAL_MAX_YEAR, ]
  late <- counted[counted$is_initial & counted$years > INITIAL_MAX_YEAR, ]
  for (i in seq_len(nrow(late)))
    issue("note", sprintf("Initial at %.1f y (%s %s, %s) excluded from t=0 anchoring (older than 1 y).",
                          late$years[i], late$taxon[i], late$accession[i], late$row_id[i]))
  if (dry_var) issue("note", "drying treatment varies — added as a grouping dimension.")

  ikey <- paste(init$taxon, init$accession, init$protocol_key, dryk(init$drying), sep = "\u0001")
  init_pool <- if (nrow(init)) lapply(split(seq_len(nrow(init)), ikey), function(idx)
    list(sown = sum(init$tested[idx]), germ = sum(init$germinated[idx]),
         pct = 100 * sum(init$germinated[idx]) / sum(init$tested[idx]))) else list()

  tkey <- paste(post$taxon, post$accession, post$storage_code, post$protocol_key, dryk(post$drying), sep = "\u0001")
  trajs <- lapply(split(seq_len(nrow(post)), tkey), function(idx) {
    r1 <- post[idx[1], ]
    ini <- init_pool[[paste(r1$taxon, r1$accession, r1$protocol_key, dryk(r1$drying), sep = "\u0001")]]
    pts <- aggregate(cbind(n = post$tested[idx], g = post$germinated[idx]),
                     by = list(t = post$years[idx]), FUN = sum)
    pts$includes_initial <- FALSE
    if (!is.null(ini)) pts <- rbind(data.frame(t = 0, n = ini$sown, g = ini$germ, includes_initial = TRUE), pts)
    pts <- pts[order(pts$t), ]
    early <- post$pct[idx][post$years[idx] > 0 & post$years[idx] <= 5]
    rise <- if (!is.null(ini) && length(early)) max(early) - ini$pct else NA_real_
    list(taxon = r1$taxon, accession = r1$accession, family = r1$family,
         storage_code = r1$storage_code, protocol_key = r1$protocol_key,
         substrate = r1$substrate, treatment = r1$treatment,
         protocol_complete = r1$protocol_complete,
         drying = if (r1$drying == "") NA_character_ else r1$drying,
         matched_initial = ini, points = pts,
         post_rows = post[idx, ],
         n_post_timepoints = length(unique(post$years[idx])),
         n_seeds = sum(pts$n),
         early_rise_pp = rise,
         early_rise_ge5 = !is.na(rise) && rise >= 5,
         early_rise_ge10 = !is.na(rise) && rise >= 10,
         early_rise_class = if (is.null(ini)) "not assessable: no same-protocol Initial"
           else if (is.na(rise)) "not assessable: no post-storage test within 0–5 y"
           else if (rise >= 10) "substantial early rise (≥10 pp)"
           else if (rise >= 5) "early rise (≥5 to <10 pp)"
           else "no ≥5 pp early rise detected")
  })
  ord <- order(vapply(trajs, function(t) paste(t$taxon, t$accession, t$storage_code, t$protocol_key), ""))
  list(trajectories = trajs[ord], drying_keyed = dry_var, n_late_initials = nrow(late))
}

plateau_diag <- function(pts) {
  pct <- 100 * pts$g / pts$n
  lead <- 0; for (v in pct) { if (!is.na(v) && v >= 99) lead <- lead + 1 else break }
  trail <- 0; for (v in rev(pct)) { if (!is.na(v) && v <= 0) trail <- trail + 1 else break }
  list(candidate = lead >= 2 || trail >= 2,
       note = if (lead >= 2 && trail >= 2) "leading ~100% plateau and trailing 0% plateau"
         else if (lead >= 2) "leading ~100% plateau"
         else if (trail >= 2) "trailing 0% plateau" else "none detected")
}

fit_formal_survival <- function(traj, threshold_pct = NA_real_) {
  pts <- traj$points
  n_time <- length(unique(pts$t))
  max_t <- max(pts$t)
  pl <- plateau_diag(pts)
  base <- list(taxon = traj$taxon, accession = traj$accession, family = traj$family,
    storage_code = traj$storage_code, protocol_key = traj$protocol_key,
    substrate = traj$substrate, treatment = traj$treatment, drying = traj$drying,
    points = pts, fit_ok = FALSE, fit_method = NA_character_, fit_warning = NA_character_,
    n_timepoints = n_time, n_post_timepoints = traj$n_post_timepoints,
    n_seeds = traj$n_seeds, max_time = max_t,
    matched_initial_available = !is.null(traj$matched_initial),
    matched_initial_pct = if (is.null(traj$matched_initial)) NA_real_ else traj$matched_initial$pct,
    Ki_probit = NA_real_, Ki_percent = NA_real_, slope = NA_real_, slope_se = NA_real_, slope_p = NA_real_,
    sigma_years = NA_real_, p50_years = NA_real_, p70_absolute_years = NA_real_,
    threshold_pct = threshold_pct,
    fitted_rci_p70_years_raw = NA_real_, fitted_rci_p70_display = "not interpretable",
    crossing_status = "model unavailable", plateau_candidate = pl$candidate, plateau_note = pl$note,
    early_rise_pp = traj$early_rise_pp, early_rise_ge5 = traj$early_rise_ge5,
    early_rise_ge10 = traj$early_rise_ge10, early_rise_class = traj$early_rise_class,
    bridge_category = "model unavailable", bridge_evidence_tier = "do not interpret",
    display_class = "NO FORMAL FIT", qc_display_status = NA_character_,
    protocol_complete = traj$protocol_complete)
  qc <- function(b) {
    if (!b$protocol_complete) return("review: incomplete protocol metadata")
    if (b$early_rise_ge10) return("review: substantial early rise; monotonic survival assumption questionable")
    if (b$early_rise_ge5) return("review: early rise; interpret monotonic survival fit cautiously")
    if (!b$matched_initial_available && b$bridge_evidence_tier == "primary model evidence")
      return("caution: within-range crossing but no same-protocol Initial")
    if (!b$matched_initial_available && b$bridge_evidence_tier == "extrapolated model evidence")
      return("caution: extrapolated crossing and no same-protocol Initial")
    b$bridge_evidence_tier
  }
  if (nrow(pts) < 2 || n_time < 2) { base$qc_display_status <- qc(base); return(base) }
  m <- glm_probit(pts$g, pts$n, data.frame(t = pts$t))
  if (is.null(m)) { base$qc_display_status <- qc(base); return(base) }
  b0 <- m$coef[1]; b1 <- m$coef[2]
  ki_pct <- pnorm(b0) * 100
  declining <- is.finite(b1) && b1 < 0
  se <- sqrt(m$vcov[2, 2]); pz <- if (se > 0) 2 * pnorm(-abs(b1 / se)) else NA_real_
  ug <- if (declining && !is.na(threshold_pct) && threshold_pct > 0 && threshold_pct < 100)
    crossing_t(b0, b1, threshold_pct / 100) else NA_real_
  cat_ <- if (!is.finite(b1) || b1 >= 0) "non-declining slope"
    else if (is.na(ug) || !is.finite(ug)) "crossing unavailable"
    else if (ug < 0) "fitted viability below threshold at model time 0"
    else if (ug <= max_t) "crossing within observed time range"
    else "crossing extrapolated beyond observed time range"
  tier <- if (!is.finite(b1) || b1 >= 0 || n_time < 2 || is.na(ug) || !is.finite(ug) || ug < 0) "do not interpret"
    else if (n_time == 2) "exploratory only"
    else if (ug <= max_t) "primary model evidence"
    else "extrapolated model evidence"
  base$fit_ok <- TRUE; base$fit_method <- "glm-probit"
  base$fit_warning <- if (m$separated) "possible separation: fit did not stabilize at finite estimates"
    else if (m$extreme) "fitted probabilities numerically 0 or 1 occurred" else NA_character_
  base$Ki_probit <- b0; base$Ki_percent <- ki_pct; base$slope <- b1
  base$slope_se <- se; base$slope_p <- pz
  base$sigma_years <- if (declining) -1 / b1 else NA_real_
  base$p50_years <- if (declining && ki_pct >= 50) crossing_t(b0, b1, 0.5) else NA_real_
  base$p70_absolute_years <- if (declining && ki_pct >= 70) crossing_t(b0, b1, 0.7) else NA_real_
  base$fitted_rci_p70_years_raw <- ug
  base$fitted_rci_p70_display <- if (is.na(threshold_pct)) "no RCI-P70 threshold (accession baseline unavailable)" else crossing_display(ug)
  base$crossing_status <- if (!declining) "non-declining fitted slope; longevity crossing not interpreted"
    else if (is.na(threshold_pct)) "RCI-P70 threshold unavailable"
    else if (is.na(ug)) "threshold crossing unavailable"
    else if (ug < 0) "fitted viability already below RCI P70 at model time 0"
    else if (ug <= max_t) "crossing within observed storage-time range"
    else "crossing extrapolated beyond observed storage-time range"
  base$bridge_category <- cat_; base$bridge_evidence_tier <- tier
  base$display_class <- if (tier == "do not interpret") "DO NOT INTERPRET"
    else if (tier == "exploratory only") "EXPLORATORY" else "FORMAL SURVIVAL"
  base$extreme <- m$extreme
  base$qc_display_status <- qc(base)
  base
}

thresholds_from_baselines <- function(baselines) {
  # (family|taxon|accession) -> RCI-P70 threshold % — the ONE genuine
  # contingency of the formal workflow on the RCI baseline stage
  setNames(vapply(baselines, function(b) b$threshold, 0), names(baselines))
}
threshold_for <- function(thr, family, taxon, accession) {
  v <- thr[[paste(family, taxon, accession, sep = "\u0001")]]
  if (is.null(v) || is.na(v)) NA_real_ else v
}

# ============================ command line ============================
main <- function() {
  args <- commandArgs(trailingOnly = TRUE)
  if (!length(args)) { cat("usage: Rscript longevity_formal.R INPUT.csv [-o OUTDIR] [--thresholds rci_accession_intervals.csv] [--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)
  thr <- if ("--thresholds" %in% args) {
    tf <- read.csv(args[which(args == "--thresholds") + 1], colClasses = "character")
    setNames(as.numeric(tf$threshold), paste(tf$family, tf$taxon, tf$accession, sep = "\u0001"))
  } else thresholds_from_baselines(early_baselines(rci_rows(a)))
  g <- group_formal_trajectories(a)
  res <- lapply(g$trajectories, function(t) fit_formal_survival(t, threshold_for(thr, t$family, t$taxon, t$accession)))
  d <- do.call(rbind, lapply(res, function(r) data.frame(
    taxon = r$taxon, accession = r$accession, storage_code = r$storage_code,
    drying = r$drying, protocol_key = r$protocol_key, substrate = r$substrate,
    treatment = r$treatment, display_class = r$display_class, fit_method = r$fit_method,
    fit_warning = r$fit_warning, qc_display_status = r$qc_display_status,
    bridge_evidence_tier = r$bridge_evidence_tier, bridge_category = r$bridge_category,
    n_timepoints = r$n_timepoints, n_post_timepoints = r$n_post_timepoints, n_seeds = r$n_seeds,
    matched_initial_available = r$matched_initial_available,
    Ki_probit = r$Ki_probit, Ki_percent = r$Ki_percent, slope = r$slope,
    slope_se = r$slope_se, slope_p = r$slope_p, sigma_years = r$sigma_years,
    p50_years = r$p50_years, p70_absolute_years = r$p70_absolute_years,
    threshold_pct = r$threshold_pct, fitted_rci_p70_years_raw = r$fitted_rci_p70_years_raw,
    fitted_rci_p70_display = r$fitted_rci_p70_display, crossing_status = r$crossing_status,
    early_rise_class = r$early_rise_class, plateau_note = r$plateau_note)))
  dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
  p1 <- file.path(outdir, "formal_survival_results.csv")
  write.csv(d, p1, row.names = FALSE, na = "")
  cat(sprintf("%s: %d trajectories (%d FORMAL SURVIVAL)\n", p1, nrow(d), sum(d$display_class == "FORMAL SURVIVAL")))
}
main()