#!/usr/bin/env Rscript
# longevity_normalize.R — Component 1 · the shared upload/normalization layer
#
# Part of the Seed Longevity Lab's downloadable lab (site v0.17.0,
# longevity.seedsofhawaii.org/scripts). No upstream dependency — every other component builds on this one's output (or embeds it).
#
# 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)
  }
}

# ============================ command line ============================
main <- function() {
  args <- commandArgs(trailingOnly = TRUE)
  if (!length(args)) { cat("usage: Rscript longevity_normalize.R INPUT.csv [-o normalized_assays.csv] [--time-unit days|months|years]\n"); quit(status = 1) }
  input <- args[1]
  out <- if ("-o" %in% args) args[which(args == "-o") + 1] else "normalized_assays.csv"
  tu <- if ("--time-unit" %in% args) args[which(args == "--time-unit") + 1] else NULL
  a <- parse_completed_tests(input, time_unit = tu, strict = !("--allow-percent-only" %in% args))
  d <- data.frame(
    row_id = a$row_id, assay_id = a$assay_id, family = a$family, taxon = a$taxon,
    accession = a$accession, is_initial = ifelse(a$is_initial, "TRUE", "FALSE"),
    storage_condition = ifelse(a$is_initial, "", a$storage_code),
    storage_condition_original = a$storage_original,
    drying_treatment = a$drying, substrate = a$substrate, germination_treatment = a$treatment,
    protocol_key = a$protocol_key, tested = a$tested, germinated = a$germinated,
    germination_percent_recomputed = round(a$pct, 4),
    germination_percent_supplied = "", storage_age_original = a$age_original,
    storage_age_unit_original = a$age_unit_original,
    storage_age_days = round(a$days, 4), storage_age_years = round(a$years, 6),
    date_derived_years = "", collection_date = a$collection_date, test_date = a$test_date,
    source_institution = "", notes = "", warnings = "", check.names = FALSE)
  write.csv(d, out, row.names = FALSE, na = "")
  cat(sprintf("%s: %d assays normalized\n", out, nrow(d)))
}
main()