#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""longevity_formal.py — Component 3 · formal protocol-specific survival models

Part of the Seed Longevity Lab's downloadable Python 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, so this file has NO imports from its siblings and no
third-party dependencies.
Python 3.9+. Input: a raw completed-test CSV (one row per completed
germination test) OR the suite's normalized_assays.csv interchange file.

Scientific contract: the Longevity Lab developer handoff
(specification/02_ANALYSIS_CONTRACT). Numerically cross-validated against
the site's JS engine, against R-equivalent references, and against the
matching R script (this component ships in both languages) 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.
"""

import sys

# ---------------------------------------------------------------------------
# SHARED CORE · statistics (identical block in every Longevity Lab script)
# Pure standard library. The IRLS binomial-probit GLM here is the same
# algorithm as the site engine (src/lib/formal.js) and is cross-validated
# against it, against R-equivalent statsmodels references, and against the
# viabilitymetrics R package lineage. No sample size is ever assumed.
# ---------------------------------------------------------------------------
import math

def normcdf(z):
    return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))

def dnorm(z):
    return math.exp(-0.5 * z * z) / math.sqrt(2.0 * math.pi)

def probit(p):
    """Inverse normal CDF (Acklam) polished with one Newton step."""
    if p <= 0.0 or p >= 1.0:
        raise ValueError("probit needs 0<p<1")
    a = [-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
         1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00]
    b = [-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
         6.680131188771972e+01, -1.328068155288572e+01]
    c = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
         -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00]
    d = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
         3.754408661907416e+00]
    plow, phigh = 0.02425, 1 - 0.02425
    if p < plow:
        q = math.sqrt(-2 * math.log(p))
        x = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1)
    elif p <= phigh:
        q = p - 0.5; r = q * q
        x = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1)
    else:
        q = math.sqrt(-2 * math.log(1 - p))
        x = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1)
    e = normcdf(x) - p
    x -= e * math.sqrt(2.0 * math.pi) * math.exp(0.5 * x * x)
    return x

def wilson(k, n, z=1.959963984540054):
    if n <= 0:
        return (float("nan"), float("nan"))
    p = k / n
    den = 1 + z * z / n
    center = (p + z * z / (2 * n)) / den
    half = (z / den) * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n))
    return (max(0.0, center - half), min(1.0, center + half))

def _solve(A, b):
    n = len(b)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for col in range(n):
        piv = max(range(col, n), key=lambda r: abs(M[r][col]))
        if abs(M[piv][col]) < 1e-12:
            return None
        M[col], M[piv] = M[piv], M[col]
        for r in range(n):
            if r == col:
                continue
            f = M[r][col] / M[col][col]
            for k2 in range(col, n + 1):
                M[r][k2] -= f * M[col][k2]
    return [M[i][n] / M[i][i] for i in range(n)]

def _invert(A):
    n = len(A)
    cols = []
    for j in range(n):
        e = [0.0] * n
        e[j] = 1.0
        x = _solve([row[:] for row in A], e)
        if x is None:
            return None
        cols.append(x)
    return [[cols[j][i] for j in range(n)] for i in range(n)]

def _binom_deviance(rows, mu):
    dev = 0.0
    for i, (x, n, g) in enumerate(rows):
        m = min(1 - 1e-12, max(1e-12, mu[i]))
        if g > 0:
            dev += g * math.log(g / (n * m))
        if n - g > 0:
            dev += (n - g) * math.log((n - g) / (n * (1 - m)))
    return 2 * dev

def glm_probit(rows):
    """Binomial GLM, probit link, general design matrix (IRLS).
    rows: list of (x_vector, n, g). Returns dict or None.
    Separation/extreme fits are FLAGGED, never bias-reduced (the site's
    'flagged, not faked' rule; see docs/DIFFERENCES_FROM_R_REFERENCE.md)."""
    usable = [(list(x), float(n), float(g)) for x, n, g in rows
              if n > 0 and 0 <= g <= n and all(math.isfinite(v) for v in x)]
    if not usable:
        return None
    p = len(usable[0][0])
    totn = sum(r[1] for r in usable)
    totg = sum(r[2] for r in usable)
    beta = [0.0] * p
    beta[0] = probit(min(0.999, max(0.001, totg / max(1.0, totn))))
    dev = float("inf")
    converged = False
    for _ in range(200):
        eta = [sum(xi * beta[j] for j, xi in enumerate(x)) for x, n, g in usable]
        mu = [min(1 - 1e-12, max(1e-12, normcdf(e))) for e in eta]
        XtWX = [[0.0] * p for _ in range(p)]
        XtWz = [0.0] * p
        for i, (x, n, g) in enumerate(usable):
            phi = dnorm(eta[i])
            v = mu[i] * (1 - mu[i])
            w = n * phi * phi / max(v, 1e-12)
            z = eta[i] + (g / n - mu[i]) / max(phi, 1e-12)
            for a2 in range(p):
                XtWz[a2] += w * x[a2] * z
                for b2 in range(p):
                    XtWX[a2][b2] += w * x[a2] * x[b2]
        nxt = _solve([row[:] for row in XtWX], XtWz)
        if nxt is None or any(not math.isfinite(v) for v in nxt):
            break
        cand = nxt
        for _h in range(21):
            ceta = [sum(xi * cand[j] for j, xi in enumerate(x)) for x, n, g in usable]
            cmu = [min(1 - 1e-12, max(1e-12, normcdf(e))) for e in ceta]
            cdev = _binom_deviance(usable, cmu)
            if not (math.isfinite(dev) and cdev > dev + 1e-8):
                break
            cand = [(v + beta[j]) / 2 for j, v in enumerate(cand)]
        delta = max(abs(v - beta[j]) for j, v in enumerate(cand))
        beta = cand
        rel = abs(cdev - dev) / (abs(cdev) + 0.1)
        dev = cdev
        if delta < 1e-10 or rel < 1e-13:
            converged = True
            break
    eta = [sum(xi * beta[j] for j, xi in enumerate(x)) for x, n, g in usable]
    mu = [min(1 - 1e-12, max(1e-12, normcdf(e))) for e in eta]
    extreme = any(abs(e) > 6.3 for e in eta)
    pinned = all(m > 0.999 or m < 0.001 for m in mu)
    separated = (not converged) or pinned or any(abs(b3) > 40 for b3 in beta)
    XtWX = [[0.0] * p for _ in range(p)]
    for i, (x, n, g) in enumerate(usable):
        phi = dnorm(eta[i])
        v = mu[i] * (1 - mu[i])
        w = n * phi * phi / max(v, 1e-12)
        for a2 in range(p):
            for b2 in range(p):
                XtWX[a2][b2] += w * x[a2] * x[b2]
    cov = _invert(XtWX)
    loglik = sum(g * math.log(mu[i]) + (n - g) * math.log(1 - mu[i])
                 for i, (x, n, g) in enumerate(usable))
    return {"coef": beta, "cov": cov, "deviance": dev,
            "df_residual": len(usable) - p, "converged": converged,
            "extreme_prob": extreme, "separated": separated,
            "fitted": mu, "loglik": loglik}

def crossing(intercept, slope, target_prop):
    if not (math.isfinite(intercept) and math.isfinite(slope)) or slope >= 0:
        return None
    if not (0 < target_prop < 1):
        return None
    return (probit(target_prop) - intercept) / slope

def crossing_display(raw):
    if raw is None or not math.isfinite(raw):
        return "not interpretable"
    if raw <= 0:
        return "not interpretable (crossing ≤0 y)"
    return f"{raw:.2f} y"

def chisq1_sf(x):
    if not math.isfinite(x) or x < 0:
        return float("nan")
    return 2 * (1 - normcdf(math.sqrt(x)))

# ---------------------------------------------------------------------------
# SHARED CORE · upload/normalization layer (identical block in every script)
# Port of the site's src/lib/upload.js — the handoff contract:
#   · one row = one completed germination test (same age may repeat)
#   · actual tested + germinated counts required; percent recomputed from
#     counts; NOTHING here assumes a sample size for percentage-only data
#   · time in days / months / years / dates; original value + unit retained
#     alongside exact normalized days and years; no rounding, ever
#   · Initial is a baseline state (is_initial preferred; legacy
#     storage="Initial" recognized), not a storage treatment
#   · protocol_key = substrate + germination treatment, trivial case and
#     whitespace normalization ONLY — never synonym inference
#   · drying treatment optional; internal row IDs; ambiguous same-date
#     duplicates flagged for review, never merged or deleted
# ---------------------------------------------------------------------------
import csv, io, re, sys
from datetime import date, timedelta

DAYS_PER_YEAR = 365.25
DAYS_PER_MONTH = 365.25 / 12.0

ROLES = {
    "family": ["family"],
    "taxon": ["taxonname", "taxon", "species", "scientificname", "speciesname"],
    "accession": ["acquisitionnum", "accession", "accessionnum", "accessionnumber", "accno", "lot", "lotid", "collectionid"],
    "assay_id": ["assayid", "testid", "germinationtestid", "assay", "testnumber", "testno"],
    "source_institution": ["sourceinstitution", "institution", "seedbank", "source"],
    "is_initial": ["isinitial", "initial", "initialtest"],
    "storage": ["storedrycode", "storagecondition", "storagecode", "storecode", "storage", "storagetemp", "storagetemperature", "storagetempc"],
    "drying": ["dryingtreatment", "drying", "drytreatment", "dryingprotocol"],
    "substrate": ["substrate", "medium"],
    "treatment": ["germinationtreatment", "treatment", "pretreatment"],
    "tested": ["numbersown", "numsown", "sown", "seedsown", "seedssown", "tested", "numbertested", "ntested", "nsown", "n"],
    "germinated": ["numbergermin", "numbergerminated", "germinated", "numgerminated", "germin", "germ", "ngerminated"],
    "percent": ["%germin", "percentgermination", "percentgerminated", "germinationpercent", "germinationpct", "pctgermin", "percent", "pct", "viability", "viabilitypct"],
    "age_days": ["storageagedays", "storagedays", "days", "day", "ageindays"],
    "age_months": ["storageagemonths", "storagemonths", "months", "month", "ageinmonths"],
    "age_years": ["yearsdifference", "storageageyears", "storageyears", "years", "year", "yr", "yeardiff", "yearsdiff", "ageinyears"],
    "age_ambiguous": ["storageage", "age", "time", "interval", "storagetime"],
    "collection_date": ["collectiondate", "storagestartdate", "collecteddate", "datecollected", "collectdate"],
    "test_date": ["teststartdate", "germinationteststartdate", "testdate", "germinationdate", "startdate", "sowdate", "datesown"],
    "notes": ["notes", "comments", "note"],
}

TRUTHY = {"true", "t", "yes", "y", "1", "initial"}
FALSY = {"false", "f", "no", "n", "0", ""}

def _norm_header(h):
    return re.sub(r"[^a-z0-9%]", "", str(h or "").lower())

def detect_roles(headers):
    mapping = {}
    for i, h in enumerate(headers):
        n = _norm_header(h)
        for role, names in ROLES.items():
            if role not in mapping and n in names:
                mapping[role] = i
                break
    return mapping

def norm_protocol_text(x):
    s = re.sub(r"\s+", " ", str(x or "").strip()).lower()
    return s if s else "(missing)"

def protocol_key(substrate, treatment):
    s, t = norm_protocol_text(substrate), norm_protocol_text(treatment)
    if s == "(missing)" and t == "(missing)":
        return "unspecified"
    return f"{s} || {t}"

def normalize_storage(raw):
    s = str(raw or "").strip()
    l = re.sub(r"[°º]", "", re.sub(r"\s+", "", s.lower()))
    l = re.sub(r"[−–—]", "-", l)
    if re.match(r"^init", s, re.I):
        return {"code": "Initial", "original": s, "legacy_initial": True}
    if re.match(r"^-18c?$", l) or l in ("minus18", "minus18c"):
        return {"code": "-18C", "original": s, "legacy_initial": False}
    if re.match(r"^\+?5c?$", l):
        return {"code": "5C", "original": s, "legacy_initial": False}
    if l.startswith("amb"):
        return {"code": "Ambient", "original": s, "legacy_initial": False}
    return {"code": s, "original": s, "legacy_initial": False}

def _parse_date_loose(v):
    s = str(v or "").strip()
    if not s:
        return None
    if re.match(r"^\d+(\.\d+)?$", s):
        num = float(s)  # Excel serial (days since 1899-12-30)
        if 59 < num < 200000:
            return date(1899, 12, 30) + timedelta(days=num)
        return None
    m = re.match(r"^(\d{4})-(\d{1,2})-(\d{1,2})", s)
    if m:
        try:
            return date(int(m.group(1)), int(m.group(2)), int(m.group(3)))
        except ValueError:
            return None
    m = re.match(r"^(\d{1,2})[/\-](\d{1,2})[/\-](\d{2,4})$", s)
    if m:
        yr = int(m.group(3))
        yr = 2000 + yr if yr < 100 else yr
        try:
            return date(yr, int(m.group(1)), int(m.group(2)))
        except ValueError:
            return None
    return None

def _f(v):
    try:
        x = float(str(v).replace("%", "").strip())
        return x
    except (ValueError, TypeError):
        return None

NORMALIZED_SIGNATURE = {"row_id", "is_initial", "protocol_key", "storage_age_years", "tested", "germinated"}

def parse_completed_tests(text, time_unit=None):
    """Raw completed-test CSV/TSV → canonical assay dicts.
    Returns dict with assays, issues, mapping, headers, needs_time_unit,
    percent_only, duplicate_groups, rejected."""
    issues, assays, rejected = [], [], []
    lines = [l for l in text.lstrip("").strip().splitlines() if l.strip()]
    out = {"assays": assays, "issues": issues, "mapping": {}, "headers": [],
           "needs_time_unit": False, "percent_only": False,
           "duplicate_groups": [], "rejected": rejected}
    if not lines:
        issues.append(("error", "No data."))
        return out
    delim = "\t" if "\t" in lines[0] else (";" if lines[0].count(";") > lines[0].count(",") else ",")
    rows = list(csv.reader(io.StringIO("\n".join(lines)), delimiter=delim))
    headers = [h.strip() for h in rows[0]]
    out["headers"] = headers
    mapping = detect_roles(headers)
    out["mapping"] = mapping

    for role in ("taxon", "accession"):
        if role not in mapping:
            issues.append(("error", f"Missing a {role} column. Headers seen: {', '.join(headers)}"))
    if "storage" not in mapping and "is_initial" not in mapping:
        issues.append(("error", "Missing a storage-condition column (post-storage rows need one; Initial rows may use is_initial)."))

    has_counts = "tested" in mapping and "germinated" in mapping
    has_pct = "percent" in mapping
    if not has_counts and not has_pct:
        issues.append(("error", 'No germination data found — need "tested" (number sown) + "germinated" counts.'))
    out["percent_only"] = (not has_counts) and has_pct
    if out["percent_only"]:
        issues.append(("block",
            "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 and upload again — the site will not assume a sample size."))

    has_days = "age_days" in mapping
    has_months = "age_months" in mapping
    has_years = "age_years" in mapping
    has_amb = "age_ambiguous" in mapping
    has_dates = "collection_date" in mapping and "test_date" in mapping
    amb_unit = None
    if not (has_days or has_months or has_years) and has_amb:
        if time_unit in ("days", "months", "years"):
            amb_unit = time_unit
            issues.append(("info", f'Ambiguous time column read as {time_unit} (your --time-unit selection).'))
        else:
            out["needs_time_unit"] = True
            issues.append(("choice",
                f'The time column "{headers[mapping["age_ambiguous"]]}" doesn\'t say its unit. '
                "Pass --time-unit days|months|years — this script won't guess."))
    if not (has_days or has_months or has_years or has_amb or has_dates):
        issues.append(("error", "No storage-age column found — need days / months / years, or CollectionDate + TestStartDate."))

    if any(l == "error" for l, _ in issues) or out["needs_time_unit"]:
        return out

    n_legacy_initial = n_date_derived = n_discrepant = n_pct_mismatch = 0
    for li, cells in enumerate(rows[1:], start=1):
        if not any(str(c).strip() for c in cells):
            continue
        row_id = f"row-{li}"
        warnings = []
        def get(role):
            i = mapping.get(role)
            return (cells[i].strip() if i is not None and i < len(cells) else "")
        taxon, accession = get("taxon"), get("accession")
        if not taxon or not accession:
            rejected.append({"row_id": row_id, "line": li + 1, "reason": "missing taxon or accession"})
            continue
        st = normalize_storage(get("storage"))
        is_initial = st["legacy_initial"]
        if "is_initial" in mapping:
            v = get("is_initial").lower()
            if v in TRUTHY:
                is_initial = True
            elif v not in FALSY:
                warnings.append(f'unrecognized is_initial value "{get("is_initial")}" — treated as FALSE')
        if st["legacy_initial"]:
            n_legacy_initial += 1

        tested = germinated = pct = None
        pct_supplied = _f(get("percent")) if has_pct else None
        if has_counts:
            tested, germinated = _f(get("tested")), _f(get("germinated"))
            if tested is None or germinated is None:
                rejected.append({"row_id": row_id, "line": li + 1, "reason": "missing counts"}); continue
            if tested <= 0:
                rejected.append({"row_id": row_id, "line": li + 1, "reason": f"number sown must be > 0 (got {tested:g})"}); continue
            if germinated < 0:
                rejected.append({"row_id": row_id, "line": li + 1, "reason": "germinated < 0"}); continue
            if germinated > tested:
                rejected.append({"row_id": row_id, "line": li + 1, "reason": f"germinated ({germinated:g}) > sown ({tested:g})"}); continue
            pct = 100.0 * germinated / tested
            if pct_supplied is not None and abs(pct_supplied - pct) > 2:
                warnings.append(f"supplied percent {pct_supplied:g}% ≠ recomputed {pct:.1f}% — counts are authoritative")
                n_pct_mismatch += 1
        elif pct_supplied is not None:
            pct = min(100.0, max(0.0, pct_supplied))  # inspection only — no invented n
        else:
            rejected.append({"row_id": row_id, "line": li + 1, "reason": "no germination value"}); continue

        age_original = age_unit = None
        days = None
        if has_days and get("age_days") != "":
            age_original, age_unit = _f(get("age_days")), "days"; days = age_original
        elif has_years and get("age_years") != "":
            age_original, age_unit = _f(get("age_years")), "years"
            days = None if age_original is None else age_original * DAYS_PER_YEAR
        elif has_months and get("age_months") != "":
            age_original, age_unit = _f(get("age_months")), "months"
            days = None if age_original is None else age_original * DAYS_PER_MONTH
        elif amb_unit and get("age_ambiguous") != "":
            age_original = _f(get("age_ambiguous")); age_unit = amb_unit
            scale = 1.0 if amb_unit == "days" else (DAYS_PER_MONTH if amb_unit == "months" else DAYS_PER_YEAR)
            days = None if age_original is None else age_original * scale
        date_derived_years = None
        cd, td = _parse_date_loose(get("collection_date")), _parse_date_loose(get("test_date"))
        if cd and td:
            dd = (td - cd).days  # exact calendar-day difference, never rounded
            if dd >= 0:
                date_derived_years = dd / DAYS_PER_YEAR
                if days is None:
                    age_original, age_unit, days = dd, "dates", float(dd)
                    n_date_derived += 1
                elif abs(days / DAYS_PER_YEAR - date_derived_years) > 0.15:
                    warnings.append(f"explicit age {days / DAYS_PER_YEAR:.2f} y disagrees with date-derived {date_derived_years:.2f} y — explicit age used")
                    n_discrepant += 1
        if days is None or days < 0:
            rejected.append({"row_id": row_id, "line": li + 1, "reason": "missing or invalid storage age"}); continue

        substrate, treatment = get("substrate"), get("treatment")
        s_norm, t_norm = norm_protocol_text(substrate), norm_protocol_text(treatment)
        assays.append({
            "row_id": row_id,
            "assay_id": get("assay_id") or None,
            "family": get("family"), "taxon": taxon, "accession": accession,
            "is_initial": is_initial,
            "storage_code": "Initial" if is_initial else st["code"],
            "storage_original": st["original"],
            "drying": get("drying") or None,
            "substrate": substrate or None, "treatment": treatment or None,
            "substrate_norm": s_norm, "treatment_norm": t_norm,
            "protocol_key": protocol_key(substrate, treatment),
            "protocol_complete": s_norm != "(missing)" and t_norm != "(missing)",
            "tested": tested, "germinated": germinated,
            "pct": pct, "pct_supplied": pct_supplied,
            "age_original": age_original, "age_unit_original": age_unit,
            "days": days, "years": days / DAYS_PER_YEAR,
            "date_derived_years": date_derived_years,
            "collection_date": cd.isoformat() if cd else None,
            "test_date": td.isoformat() if td else None,
            "source_institution": get("source_institution") or None,
            "notes": get("notes") or None,
            "warnings": warnings,
        })

    # ambiguous duplicate / multiple same-date tests — flag, never delete
    dup = {}
    for a in assays:
        date_key = a["test_date"] or f"age:{a['days']}"
        k = (a["taxon"], a["accession"], a["storage_code"], a["protocol_key"],
             date_key, a["tested"], a["germinated"])
        dup.setdefault(k, []).append(a)
    for k, group in dup.items():
        if len(group) < 2:
            continue
        aids = [a["assay_id"] for a in group]
        if all(x is not None for x in aids) and len(set(aids)) == len(aids):
            continue  # distinct institutional IDs establish distinct tests
        out["duplicate_groups"].append({
            "taxon": k[0], "accession": k[1], "storage": k[2], "protocol": k[3],
            "date_key": k[4], "row_ids": [a["row_id"] for a in group]})
    if out["duplicate_groups"]:
        issues.append(("warn",
            f"Possible duplicate or multiple same-date tests: {len(out['duplicate_groups'])} group(s) of "
            "identical rows preserved — confirm whether they are distinct completed tests or duplicate records."))
    if n_legacy_initial:
        issues.append(("info", f'{n_legacy_initial} row(s) used legacy storage="Initial" — mapped to is_initial=TRUE.'))
    if n_date_derived:
        issues.append(("info", f"{n_date_derived} row(s): storage age derived from dates as exact elapsed days (no rounding)."))
    if n_discrepant:
        issues.append(("warn", f"{n_discrepant} row(s): explicit age and date-derived age disagree by >0.15 y (explicit age used)."))
    if n_pct_mismatch:
        issues.append(("warn", f"{n_pct_mismatch} row(s): supplied percent disagrees with counts by >2 points — counts are authoritative."))
    if rejected:
        reasons = "; ".join(sorted({r["reason"] for r in rejected}))
        issues.append(("warn", f"{len(rejected)} row(s) rejected: {reasons}."))
    return out

NORMALIZED_COLUMNS = [
    "row_id", "assay_id", "family", "taxon", "accession", "is_initial",
    "storage_condition", "storage_condition_original", "drying_treatment",
    "substrate", "germination_treatment", "protocol_key",
    "tested", "germinated", "germination_percent_recomputed",
    "germination_percent_supplied", "storage_age_original",
    "storage_age_unit_original", "storage_age_days", "storage_age_years",
    "date_derived_years", "collection_date", "test_date",
    "source_institution", "notes", "warnings",
]

def assays_to_normalized_rows(assays):
    rows = []
    for a in assays:
        rows.append({
            "row_id": a["row_id"], "assay_id": a["assay_id"] or "",
            "family": a["family"], "taxon": a["taxon"], "accession": a["accession"],
            "is_initial": "TRUE" if a["is_initial"] else "FALSE",
            "storage_condition": "" if a["is_initial"] else a["storage_code"],
            "storage_condition_original": a["storage_original"],
            "drying_treatment": a["drying"] or "",
            "substrate": a["substrate"] or "", "germination_treatment": a["treatment"] or "",
            "protocol_key": a["protocol_key"],
            "tested": "" if a["tested"] is None else f"{a['tested']:g}",
            "germinated": "" if a["germinated"] is None else f"{a['germinated']:g}",
            "germination_percent_recomputed": "" if a["pct"] is None else f"{a['pct']:.4f}",
            "germination_percent_supplied": "" if a["pct_supplied"] is None else f"{a['pct_supplied']:g}",
            "storage_age_original": "" if a["age_original"] is None else f"{a['age_original']:g}",
            "storage_age_unit_original": a["age_unit_original"] or "",
            "storage_age_days": f"{a['days']:.4f}",
            "storage_age_years": f"{a['years']:.6f}",
            "date_derived_years": "" if a["date_derived_years"] is None else f"{a['date_derived_years']:.6f}",
            "collection_date": a["collection_date"] or "", "test_date": a["test_date"] or "",
            "source_institution": a["source_institution"] or "", "notes": a["notes"] or "",
            "warnings": " | ".join(a["warnings"]),
        })
    return rows

def normalized_rows_to_assays(dict_rows):
    assays = []
    for r in dict_rows:
        tested = _f(r.get("tested"))
        germinated = _f(r.get("germinated"))
        years = _f(r.get("storage_age_years"))
        days = _f(r.get("storage_age_days"))
        if years is None and days is not None:
            years = days / DAYS_PER_YEAR
        if years is None:
            continue
        is_initial = str(r.get("is_initial", "")).strip().lower() in TRUTHY
        substrate = r.get("substrate") or None
        treatment = r.get("germination_treatment") or None
        s_norm, t_norm = norm_protocol_text(substrate), norm_protocol_text(treatment)
        pct = _f(r.get("germination_percent_recomputed"))
        if pct is None and tested and germinated is not None and tested > 0:
            pct = 100.0 * germinated / tested
        assays.append({
            "row_id": r.get("row_id") or f"row-{len(assays)+1}",
            "assay_id": r.get("assay_id") or None,
            "family": r.get("family") or "", "taxon": r.get("taxon") or "",
            "accession": r.get("accession") or "",
            "is_initial": is_initial,
            "storage_code": "Initial" if is_initial else (r.get("storage_condition") or r.get("storage_condition_original") or ""),
            "storage_original": r.get("storage_condition_original") or "",
            "drying": r.get("drying_treatment") or None,
            "substrate": substrate, "treatment": treatment,
            "substrate_norm": s_norm, "treatment_norm": t_norm,
            "protocol_key": r.get("protocol_key") or protocol_key(substrate, treatment),
            "protocol_complete": s_norm != "(missing)" and t_norm != "(missing)",
            "tested": tested, "germinated": germinated,
            "pct": pct, "pct_supplied": _f(r.get("germination_percent_supplied")),
            "age_original": _f(r.get("storage_age_original")),
            "age_unit_original": r.get("storage_age_unit_original") or None,
            "days": days if days is not None else years * DAYS_PER_YEAR,
            "years": years, "date_derived_years": _f(r.get("date_derived_years")),
            "collection_date": r.get("collection_date") or None,
            "test_date": r.get("test_date") or None,
            "source_institution": r.get("source_institution") or None,
            "notes": r.get("notes") or None, "warnings": [],
        })
    return assays

def load_input(path, time_unit=None):
    """Read EITHER a raw completed-test CSV or a normalized_assays.csv
    (this suite's interchange format — also what the website exports).
    Returns (assays, parse_result_or_None, was_normalized)."""
    with open(path, newline="", encoding="utf-8-sig") as fh:
        text = fh.read()
    first = text.lstrip("").splitlines()[0] if text.strip() else ""
    header_set = {h.strip() for h in first.split(",")}
    if NORMALIZED_SIGNATURE.issubset(header_set):
        with io.StringIO(text) as fh2:
            dict_rows = list(csv.DictReader(fh2))
        return normalized_rows_to_assays(dict_rows), None, True
    parsed = parse_completed_tests(text, time_unit=time_unit)
    return parsed["assays"], parsed, False

def write_csv(path, columns, dict_rows):
    with open(path, "w", newline="", encoding="utf-8") as fh:
        w = csv.DictWriter(fh, fieldnames=columns, extrasaction="ignore")
        w.writeheader()
        for r in dict_rows:
            w.writerow(r)

def report_issues(parsed, strict=True):
    """Print parser issues to stderr. Exit 2 on blocking/percent-only when
    strict (analytical scripts); exit 1 on hard errors."""
    if parsed is None:
        return
    for level, msg in parsed["issues"]:
        print(f"[{level.upper()}] {msg}", file=sys.stderr)
    if any(l == "error" for l, _ in parsed["issues"]) or parsed["needs_time_unit"]:
        sys.exit(1)
    if strict and parsed["percent_only"]:
        sys.exit(2)

# ---------------------------------------------------------------------------
# SHARED CORE · RCI management analysis (identical block in RCI-dependent scripts)
# Port of the site's src/lib/rci.js under the PUBLIC upload rule:
#   · empirical early-window baseline per accession = highest observed
#     germination from the earliest assay age through +5 years — even when
#     only two distinct assay ages exist (the archival two-point exception is
#     NOT ported; handoff fixture 11 guards this)
#   · RCI P70 = 0.70 × baseline; max-at-age replicate rule (RCI branch ONLY)
#   · per accession × storage latest passing age; Initial excluded from
#     storage evidence; baseline-adds-a-new-year inclusion rule
#   · GLM + Wilson support at the latest passing age
#   · species roll-up: ≥2 distinct test years, max across accessions;
#     Overall observed RCI and observed RCI at recommended storage are
#     SEPARATE outputs — a recommendation override never rewrites history
#   · freeze-sensitivity: ratio ≤ 0.5 and CI-asymmetry triggers, −18C-only
#     NA→FALSE override, −18C publication override
#   · censoring-honest a–b / a+ / <b notation, fractional precision kept
# ---------------------------------------------------------------------------
EPS = 1e-8
BASELINE_WINDOW_YEARS = 5.0
THRESHOLD_FRACTION = 0.70
FREEZE_RATIO = 0.5
MARGIN_YEARS = 5.0
TEMP_LABEL = {"-18C": "−18 °C", "5C": "5 °C", "Ambient": "Ambient"}

def _rci_rows(assays):
    rows = []
    for a in assays:
        rows.append({
            "family": a["family"] or "", "taxon": a["taxon"], "accession": a["accession"],
            "storage": "Initial" if a["is_initial"] else a["storage_code"],
            "sown": a["tested"], "germ": a["germinated"], "pct": a["pct"],
            "years": a["years"], "protocol_key": a["protocol_key"],
        })
    return [r for r in rows if r["pct"] is not None]

def collapse_by_year(rows, by_storage=True, strategy="max"):
    m = {}
    for r in rows:
        if r["sown"] is None or r["sown"] <= 0 or r["germ"] is None or r["germ"] < 0:
            continue
        k = (r["taxon"], r["accession"], r["storage"] if by_storage else None, r["years"])
        cur = m.get(k)
        if strategy == "pooled":
            if cur is None:
                m[k] = dict(r)
            else:
                cur["sown"] += r["sown"]; cur["germ"] += r["germ"]
        else:
            p = r["germ"] / r["sown"]
            if cur is None or p > cur["germ"] / cur["sown"] or (
                    p == cur["germ"] / cur["sown"] and r["sown"] > cur["sown"]):
                m[k] = dict(r)
    out = list(m.values())
    for r in out:
        r["pct"] = 100.0 * r["germ"] / r["sown"]
    return out

def early_baselines(rows):
    """PUBLIC rule: highest observed germination in [earliest, earliest+5y]."""
    by_acc = {}
    for r in rows:
        k = (r["family"], r["taxon"], r["accession"])
        yr = by_acc.setdefault(k, {})
        if r["years"] not in yr or r["pct"] > yr[r["years"]]:
            yr[r["years"]] = r["pct"]
    out = {}
    for k, yr in by_acc.items():
        years = sorted(yr)
        if not years:
            continue
        upper = years[0] + BASELINE_WINDOW_YEARS
        pick_y, pick_v = None, None
        for y in years:
            if y > upper:
                break
            v = yr[y]
            if pick_v is None or v > pick_v or (v == pick_v and y < pick_y):
                pick_y, pick_v = y, v
        out[k] = {"baseline_year": pick_y, "baseline_viab": pick_v,
                  "threshold": THRESHOLD_FRACTION * pick_v}
    return out

def _glm_year(points):
    """Probit GLM of germination ~ years on collapsed points; returns
    (fit, predict_ci) or (None, None). Mirrors the site's support-CI GLM:
    ≥3 rows, ≥2 distinct years, separation flagged not faked."""
    pts = [(p["year"], p["sown"], p["germ"]) for p in points
           if p["sown"] and p["sown"] > 0 and p["germ"] is not None and 0 <= p["germ"] <= p["sown"]]
    if len(pts) < 3 or len({p[0] for p in pts}) < 2:
        return None, None
    m = glm_probit([([1.0, y], n, g) for y, n, g in pts])
    if m is None or m["separated"]:
        return m, None
    b0, b1 = m["coef"]
    V = m["cov"]
    def predict_ci(t, z=1.959963984540054):
        if V is None:
            return None
        var_eta = V[0][0] + 2 * t * V[0][1] + t * t * V[1][1]
        if var_eta < 0:
            return None
        se = math.sqrt(var_eta)
        eta = b0 + b1 * t
        return {"fit": normcdf(eta), "lo": normcdf(eta - z * se), "hi": normcdf(eta + z * se)}
    return m, predict_ci

def per_storage_intervals(rows, baselines):
    collapsed = collapse_by_year([r for r in rows if r["storage"] != "Initial"], by_storage=True)
    by_series_c, by_series_raw = {}, {}
    for r in collapsed:
        by_series_c.setdefault((r["family"], r["taxon"], r["accession"], r["storage"]), []).append(r)
    for r in rows:
        if r["storage"] == "Initial":
            continue
        by_series_raw.setdefault((r["family"], r["taxon"], r["accession"], r["storage"]), []).append(r)
    out = []
    for k, series in by_series_raw.items():
        family, taxon, accession, storage = k
        bl = baselines.get((family, taxon, accession))
        if not bl or not (bl["baseline_viab"] > 0):
            continue
        years_set = {r["years"] for r in series}
        adds_new = bl["baseline_year"] not in years_set
        if len(years_set) + (1 if adds_new else 0) < 2:
            continue
        latest = None
        for r in series:
            if r["pct"] >= bl["threshold"] and (latest is None or r["years"] > latest):
                latest = r["years"]
        sub = sorted(({"year": r["years"], "sown": r["sown"], "germ": r["germ"], "pct": r["pct"]}
                      for r in by_series_c.get(k, [])), key=lambda p: p["year"])
        glm_support = wilson_support = None
        glm_lo = wilson_lo = None
        if latest is not None:
            _, predict_ci = _glm_year(sub)
            if predict_ci:
                ci = predict_ci(latest)
                if ci:
                    glm_lo = ci["lo"] * 100
                    glm_support = glm_lo >= bl["threshold"]
            at = [p for p in sub if p["year"] == latest]
            if not at and sub:
                at = [min(sub, key=lambda p: abs(p["year"] - latest))]
            if len(at) == 1 and at[0]["sown"] and at[0]["sown"] > 0:
                lo, _hi = wilson(at[0]["germ"], at[0]["sown"])
                wilson_lo = lo * 100
                wilson_support = wilson_lo >= bl["threshold"]
        out.append({"family": family, "taxon": taxon, "accession": accession,
                    "storage": storage, "baseline_year": bl["baseline_year"],
                    "baseline_viab": bl["baseline_viab"], "threshold": bl["threshold"],
                    "n_years": len(years_set), "max_year_meeting_p70": latest,
                    "glm_support": glm_support, "wilson_support": wilson_support,
                    "glm_lo": glm_lo, "wilson_lo": wilson_lo, "series": sub})
    return out

def species_rollup(rows, baselines):
    by_acc = {}
    for r in rows:
        by_acc.setdefault((r["family"], r["taxon"], r["accession"]), []).append(r)
    acc_max = []
    for k, series in by_acc.items():
        family, taxon, accession = k
        bl = baselines.get(k)
        if not bl or not (bl["baseline_viab"] > 0):
            continue
        if len({r["years"] for r in series}) < 2:
            continue
        latest = None
        for r in series:
            if r["pct"] >= bl["threshold"] and (latest is None or r["years"] > latest):
                latest = r["years"]
        acc_max.append({"family": family, "taxon": taxon, "accession": accession,
                        "baseline_year": bl["baseline_year"], "baseline_viab": bl["baseline_viab"],
                        "threshold": bl["threshold"], "acc_interval": latest})
    by_sp = {}
    for a in acc_max:
        by_sp.setdefault((a["family"], a["taxon"]), []).append(a)
    species = []
    for k, accs in by_sp.items():
        family, taxon = k
        with_int = [a for a in accs if a["acc_interval"] is not None]
        if not with_int:
            species.append({"family": family, "taxon": taxon, "interval": None, "accs": accs})
            continue
        interval = max(a["acc_interval"] for a in with_int)
        chosen = sorted(with_int, key=lambda a: (-a["acc_interval"], -a["baseline_viab"], str(a["accession"])))[0]
        species.append({
            "family": family, "taxon": taxon, "interval": interval,
            "supporting": sum(1 for a in with_int if a["acc_interval"] + EPS >= interval),
            "total_qualifying": len({a["accession"] for a in accs}),
            "chosen_accession": chosen["accession"], "chosen_threshold": chosen["threshold"],
            "chosen_baseline_year": chosen["baseline_year"], "chosen_baseline_viab": chosen["baseline_viab"],
            "accs": accs})
    return acc_max, species

def freeze_flags(intervals):
    by_sp = {}
    for it in intervals:
        by_sp.setdefault((it["family"], it["taxon"]), []).append(it)
    out = {}
    for k, rows in by_sp.items():
        def smax(vals):
            vals = [v for v in vals if v is not None]
            return max(vals) if vals else None
        R5 = smax([r["max_year_meeting_p70"] for r in rows if r["storage"] == "5C"])
        R18 = smax([r["max_year_meeting_p70"] for r in rows if r["storage"] == "-18C"])
        ci5 = any(r["storage"] == "5C" and (r["glm_support"] is True or r["wilson_support"] is True) for r in rows)
        ci18 = any(r["storage"] == "-18C" and (r["glm_support"] is True or r["wilson_support"] is True) for r in rows)
        have5, have18 = R5 is not None, R18 is not None
        ratio = have5 and have18 and R5 > 0 and R18 / R5 <= FREEZE_RATIO
        ci_trig = ci5 and not ci18 and have5 and have18 and R5 > R18
        if not have5 and not have18:
            flag, reason = None, "no 5C or -18C data"
        elif not have5:
            flag, reason = None, "no 5C data"
        elif not have18:
            flag, reason = None, "no -18C data"
        elif ratio:
            flag, reason = True, f"-18C interval ({fmt_rci_year(R18)} y) <= 50% of 5C ({fmt_rci_year(R5)} y)"
        elif ci_trig:
            flag, reason = True, "CI supports 5C at its latest; 5C interval > -18C"
        else:
            flag, reason = False, "no evidence of freeze sensitivity"
        out[k] = {"R5": R5, "R18": R18, "ci5": ci5, "ci18": ci18, "flag": flag, "reason": reason}
    return out

def fmt_rci_year(y):
    if abs(y - round(y)) < 5e-3:
        return str(int(round(y)))
    return f"{y:.1f}"

def rci_notation(series, threshold):
    s = sorted(series, key=lambda p: p["year"])
    hits = [p["year"] for p in s if p["pct"] >= threshold]
    if not hits:
        fails = [p for p in s if p["pct"] < threshold]
        return f"<{fmt_rci_year(fails[0]['year'])}" if fails else "<1"
    a = max(hits)
    after = [p["year"] for p in s if p["year"] > a and p["pct"] < threshold]
    if not after:
        return f"{fmt_rci_year(a)}+"
    return f"{fmt_rci_year(a)}–{fmt_rci_year(min(after))}"

def run_rci_pipeline(assays):
    rows = _rci_rows(assays)
    baselines = early_baselines(rows)
    intervals = per_storage_intervals(rows, baselines)
    acc_max, species = species_rollup(rows, baselines)
    freeze = freeze_flags(intervals)

    totals_raw = {}
    for r in rows:
        totals_raw.setdefault((r["family"], r["taxon"]), set()).add(r["accession"])

    by_storage_species = {}
    for it in intervals:
        v = it["max_year_meeting_p70"]
        if v is None or not v > 0:
            continue
        k = (it["family"], it["taxon"], it["storage"])
        if k not in by_storage_species or v > by_storage_species[k]:
            by_storage_species[k] = v

    collapsed_all = collapse_by_year(rows, by_storage=True)
    series_by_as = {}
    for r in collapsed_all:
        series_by_as.setdefault((r["taxon"], r["accession"], r["storage"]), []).append(
            {"year": r["years"], "sown": r["sown"], "germ": r["germ"]})

    potential_by_sp = {}
    for sp in species:
        if sp["interval"] is None:
            continue
        pot = set()
        for a in sp["accs"]:
            if a["acc_interval"] is not None and a["acc_interval"] + EPS >= sp["interval"]:
                continue
            if a["acc_interval"] is not None and a["acc_interval"] + MARGIN_YEARS >= sp["interval"]:
                pot.add(a["accession"]); continue
            pred_best = lo_best = None
            for st in ("-18C", "5C", "Ambient"):
                ser = series_by_as.get((sp["taxon"], a["accession"], st))
                if not ser:
                    continue
                pts = sorted(({"year": p["year"], "sown": p["sown"], "germ": p["germ"],
                               "pct": 100.0 * p["germ"] / p["sown"]} for p in ser), key=lambda p: p["year"])
                m, predict_ci = _glm_year(pts)
                if m is None or m["separated"] or predict_ci is None:
                    continue
                b0, b1 = m["coef"]
                pr = normcdf(b0 + b1 * sp["interval"]) * 100
                if pred_best is None or pr > pred_best:
                    pred_best = pr
                ci = predict_ci(sp["interval"])
                if ci and (lo_best is None or ci["lo"] * 100 > lo_best):
                    lo_best = ci["lo"] * 100
            if (pred_best is not None and pred_best >= a["threshold"]) or \
               (lo_best is not None and lo_best >= a["threshold"]):
                pot.add(a["accession"])
        potential_by_sp[(sp["family"], sp["taxon"])] = len(pot)

    publication = []
    for sp in species:
        if sp["interval"] is None:
            continue
        k = (sp["family"], sp["taxon"])
        fz = dict(freeze.get(k) or {"flag": None, "reason": "no per-storage intervals", "R5": None, "R18": None})
        order = {"-18C": 0, "5C": 1, "Ambient": 2}
        anchor = anchor_years = None
        for st in ("-18C", "5C", "Ambient"):
            v = by_storage_species.get((sp["family"], sp["taxon"], st))
            if v is None:
                continue
            if anchor_years is None or v > anchor_years or (v == anchor_years and order[st] < order[anchor]):
                anchor, anchor_years = st, v
        if anchor is None:
            cand = [r for r in rows if r["taxon"] == sp["taxon"] and r["accession"] == sp["chosen_accession"]
                    and r["storage"] != "Initial" and r["years"] == sp["interval"]]
            if cand:
                anchor = sorted(cand, key=lambda r: -r["pct"])[0]["storage"]
        if anchor is None:
            fb = [r for r in rows if r["taxon"] == sp["taxon"] and r["accession"] == sp["chosen_accession"]
                  and r["storage"] != "Initial" and r["years"] <= sp["interval"] and r["pct"] >= sp["chosen_threshold"]]
            if fb:
                anchor = sorted(fb, key=lambda r: (-r["years"], -r["pct"], str(r["storage"])))[0]["storage"]
        # R-reference NA→FALSE override for −18C-only species
        if fz["flag"] is None and anchor == "-18C":
            fz["flag"] = False
            fz["reason"] = f"only -18C data; species interval = {fmt_rci_year(sp['interval'])} y at -18C"

        override = fz["flag"] is False and fz["reason"] == "no evidence of freeze sensitivity"
        rec_storage = "-18C" if override else anchor
        behavior = "Orthodox" if rec_storage == "-18C" else ("Freeze-sensitive" if rec_storage == "5C" else "Unknown")

        anchor_series = [{"year": r["years"], "pct": r["pct"]} for r in collapse_by_year(
            [r for r in rows if r["taxon"] == sp["taxon"] and r["accession"] == sp["chosen_accession"]
             and r["storage"] == anchor], by_storage=False)]
        notation = rci_notation(anchor_series, sp["chosen_threshold"]) if anchor_series else None

        rec_interval = by_storage_species.get((sp["family"], sp["taxon"], rec_storage)) if rec_storage else None
        rec_notation = rec_driver = None
        if rec_storage and rec_interval is not None:
            rec_its = sorted((it for it in intervals if it["family"] == sp["family"] and it["taxon"] == sp["taxon"]
                              and it["storage"] == rec_storage and it["max_year_meeting_p70"] is not None),
                             key=lambda it: (-it["max_year_meeting_p70"], -it["baseline_viab"]))
            if rec_its:
                rec_driver = rec_its[0]["accession"]
                rec_notation = rci_notation(rec_its[0]["series"], rec_its[0]["threshold"])

        rebound = sorted({f"{it['accession']} @ {TEMP_LABEL.get(it['storage'], it['storage'])}"
                          for it in intervals
                          if it["family"] == sp["family"] and it["taxon"] == sp["taxon"]
                          and it["max_year_meeting_p70"] is not None
                          and any(p["year"] < it["max_year_meeting_p70"] and p["pct"] < it["threshold"]
                                  for p in it["series"])})

        recent = None
        exact = [p for p in anchor_series if abs(p["year"] - sp["interval"]) < 1e-9]
        if len(exact) == 1:
            recent = exact[0]["pct"]
        else:
            ok2 = sorted((p for p in anchor_series if p["year"] <= sp["interval"] and p["pct"] >= sp["chosen_threshold"]),
                         key=lambda p: (-p["year"], -p["pct"]))
            if ok2:
                recent = ok2[0]["pct"]

        if override:
            rec_reason = f"no evidence of freeze sensitivity (R5={fz.get('R5')} y, R18={fz.get('R18')} y) → −18 °C publication rule"
        elif fz["flag"] is True:
            rec_reason = f"freeze-sensitivity evidence: {fz['reason']}"
        elif anchor:
            rec_reason = f"longest demonstrated interval under {TEMP_LABEL.get(anchor, anchor)}"
        else:
            rec_reason = "attribution fallback from the driving accession"

        publication.append({
            "family": sp["family"], "taxon": sp["taxon"],
            "recommended_storage": TEMP_LABEL.get(rec_storage, rec_storage) if rec_storage else None,
            "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": TEMP_LABEL.get(anchor, anchor) if anchor else None,
            "rci_at_recommended_years_raw": rec_interval,
            "rci_at_recommended_display": (rec_notation or (f"{fmt_rci_year(rec_interval)}+" if rec_interval is not None else None))
                if rec_interval is not None else
                (f"not yet demonstrated at {TEMP_LABEL.get(rec_storage, rec_storage)}" if rec_storage else None),
            "rci_at_recommended_driver_accession": rec_driver,
            "rebound_failure_flags": "; ".join(rebound) if rebound else "",
            "accession_baseline_percent": round(sp["chosen_baseline_viab"], 1),
            "rci_p70_threshold_percent": round(sp["chosen_threshold"], 2),
            "fresh_seeds_baseline": ("Yes" if sp["chosen_baseline_year"] == 0 else "No")
                if sp["chosen_baseline_year"] is not None else None,
            "recent_germ_at_rci": round(recent, 1) if recent is not None else None,
            "supporting_accessions": sp["supporting"],
            "potentially_supporting_accessions": potential_by_sp.get(k, 0),
            "supporting_published": sp["supporting"] + potential_by_sp.get(k, 0),
            "rci_eligible_accessions": sp["total_qualifying"],
            "raw_accessions": len(totals_raw.get(k, set())),
            "override_to_minus18": override,
            "R5": fz.get("R5"), "R18": fz.get("R18"),
            "freeze_sensitive": fz["flag"], "freeze_reason": fz["reason"],
        })
    publication.sort(key=lambda p: (p["family"], p["taxon"]))
    return {"rows": rows, "baselines": baselines, "intervals": intervals,
            "acc_max": acc_max, "species": species, "freeze": freeze,
            "publication": publication}

# ---------------------------------------------------------------------------
# SHARED CORE · formal protocol-specific survival (identical block in
# formal-dependent scripts). Port of the site's src/lib/formal.js, itself
# built to the handoff's 13_protocol_specific_survival_v1_6.R +
# 10_survival_plateau_qc semantics:
#   · trajectory = taxon × accession × storage × protocol_key (+ drying only
#     when supplied and actually variable) — never pooled across protocols
#   · matched same-protocol Initial (≤ 1 y) pooled once at t = 0, never
#     borrowed across protocols; older "Initials" excluded AND reported
#   · same-time rows pool by summed counts (identical binomial likelihood)
#   · plateau + early-rise QC; bridge category/tier; raw crossing kept,
#     ≤0 y displayed as "not interpretable (crossing ≤0 y)"
# ---------------------------------------------------------------------------
INITIAL_MAX_YEAR = 1.0

def group_formal_trajectories(assays):
    counted = [a for a in assays if a["tested"] is not None and a["germinated"] is not None]
    probe = {}
    for a in counted:
        if not a["drying"]:
            continue
        k = (a["taxon"], a["accession"], a["storage_code"], a["protocol_key"])
        probe.setdefault(k, set()).add(a["drying"].strip().lower())
    drying_keyed = any(len(s) > 1 for s in probe.values())
    dry = (lambda a: (a["drying"] or "(none)").strip().lower()) if drying_keyed else (lambda a: "")

    post = [a for a in counted if not a["is_initial"]]
    initial = [a for a in counted if a["is_initial"] and a["years"] <= INITIAL_MAX_YEAR]
    excluded_late = [{"row_id": a["row_id"], "taxon": a["taxon"], "accession": a["accession"], "years": a["years"]}
                     for a in counted if a["is_initial"] and a["years"] > INITIAL_MAX_YEAR]

    init_by = {}
    for a in initial:
        k = (a["taxon"], a["accession"], a["protocol_key"], dry(a))
        e = init_by.setdefault(k, {"n_assays": 0, "sown": 0.0, "germ": 0.0})
        e["n_assays"] += 1; e["sown"] += a["tested"]; e["germ"] += a["germinated"]

    traj_map = {}
    for a in post:
        k = (a["taxon"], a["accession"], a["storage_code"], a["protocol_key"], dry(a))
        t = traj_map.setdefault(k, {
            "key": k, "taxon": a["taxon"], "accession": a["accession"], "family": a["family"],
            "storage_code": a["storage_code"], "protocol_key": a["protocol_key"],
            "substrate": a["substrate"], "treatment": a["treatment"],
            "protocol_complete": a["protocol_complete"], "drying": a["drying"] or None,
            "post_rows": []})
        t["post_rows"].append(a)

    trajectories = []
    for t in traj_map.values():
        ik = (t["taxon"], t["accession"], t["protocol_key"],
              (t["drying"] or "(none)").strip().lower() if drying_keyed else "")
        init = init_by.get(ik)
        t["matched_initial"] = None if not init else {
            "n_assays": init["n_assays"], "sown": init["sown"], "germ": init["germ"],
            "pct": 100.0 * init["germ"] / init["sown"]}
        by_t = {}
        def add(tt, n, g, is_init):
            e = by_t.setdefault(tt, {"t": tt, "n": 0.0, "g": 0.0, "n_assays": 0, "includes_initial": False})
            e["n"] += n; e["g"] += g; e["n_assays"] += 1
            e["includes_initial"] = e["includes_initial"] or is_init
        if init:
            add(0.0, init["sown"], init["germ"], True)
        for a in t["post_rows"]:
            add(a["years"], a["tested"], a["germinated"], False)
        t["points"] = sorted(by_t.values(), key=lambda p: p["t"])
        post_years = sorted({a["years"] for a in t["post_rows"]})
        t["n_post_timepoints"] = len(post_years)
        t["max_post_years"] = post_years[-1] if post_years else None
        t["n_seeds"] = sum(p["n"] for p in t["points"])
        early = [a["pct"] for a in t["post_rows"] if 0 < a["years"] <= 5]
        max_early = max(early) if early else None
        if t["matched_initial"] and max_early is not None:
            rise = max_early - t["matched_initial"]["pct"]
            t["early_rise_pp"] = rise
            t["early_rise_ge5"], t["early_rise_ge10"] = rise >= 5, rise >= 10
            t["early_rise_class"] = ("substantial early rise (≥10 pp)" if rise >= 10
                                     else "early rise (≥5 to <10 pp)" if rise >= 5
                                     else "no ≥5 pp early rise detected")
        else:
            t["early_rise_pp"] = None
            t["early_rise_ge5"] = t["early_rise_ge10"] = False
            t["early_rise_class"] = ("not assessable: no same-protocol Initial" if not t["matched_initial"]
                                     else "not assessable: no post-storage test within 0–5 y")
        trajectories.append(t)
    trajectories.sort(key=lambda t: (t["taxon"], str(t["accession"]), t["storage_code"], t["protocol_key"]))
    return {"trajectories": trajectories, "drying_keyed": drying_keyed,
            "excluded_late_initials": excluded_late}

def plateau_diagnostics(points, near100=99.0):
    s = sorted(points, key=lambda p: p["t"])
    pct = [100.0 * p["g"] / p["n"] for p in s]
    lead = 0
    for v in pct:
        if v >= near100: lead += 1
        else: break
    trail = 0
    for v in reversed(pct):
        if v <= 0: trail += 1
        else: break
    note = ("leading ~100% plateau and trailing 0% plateau" if lead >= 2 and trail >= 2
            else "leading ~100% plateau" if lead >= 2
            else "trailing 0% plateau" if trail >= 2 else "none detected")
    return {"leading_near100": lead, "trailing_zero": trail,
            "plateau_candidate": lead >= 2 or trail >= 2, "plateau_note": note}

def fit_formal_survival(traj, threshold_pct=None):
    pts = traj["points"]
    n_time = len({p["t"] for p in pts})
    max_t = max((p["t"] for p in pts), default=None)
    plateau = plateau_diagnostics(pts)
    series_tier = ("insufficient (<2 distinct times)" if n_time < 2
                   else "exploratory (2 distinct times)" if n_time == 2
                   else "interpretable candidate (≥3 distinct times)")
    base = {
        "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"],
        "protocol_complete": traj["protocol_complete"], "drying": traj["drying"],
        "points": pts, "fit_ok": False, "fit_method": None, "fit_reason": None,
        "fit_warning": None, "extreme_prob_warning": False, "converged": None,
        "series_tier": series_tier, "n_model_rows": len(pts), "n_timepoints": n_time,
        "n_post_timepoints": traj["n_post_timepoints"], "n_seeds": traj["n_seeds"],
        "max_time": max_t,
        "matched_initial_available": bool(traj["matched_initial"]),
        "matched_initial_pct": traj["matched_initial"]["pct"] if traj["matched_initial"] else None,
        "Ki_probit": None, "Ki_percent": None, "slope": None, "slope_se": None, "slope_p": None,
        "sigma_years": None, "p50_years": None, "p70_absolute_years": None,
        "threshold_pct": threshold_pct,
        "fitted_rci_p70_years_raw": None, "fitted_rci_p70_display": "not interpretable",
        "crossing_status": "model unavailable", "sparse_series": n_time < 3, **plateau,
        "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": None,
    }
    def qc_display(res):
        if not res["protocol_complete"]:
            return "review: incomplete protocol metadata"
        if res["early_rise_ge10"]:
            return "review: substantial early rise; monotonic survival assumption questionable"
        if res["early_rise_ge5"]:
            return "review: early rise; interpret monotonic survival fit cautiously"
        if not res["matched_initial_available"] and res["bridge_evidence_tier"] == "primary model evidence":
            return "caution: within-range crossing but no same-protocol Initial"
        if not res["matched_initial_available"] and res["bridge_evidence_tier"] == "extrapolated model evidence":
            return "caution: extrapolated crossing and no same-protocol Initial"
        return res["bridge_evidence_tier"]

    if len(pts) < 2 or n_time < 2:
        base["fit_reason"] = "fewer than 2 distinct storage times"
        base["qc_display_status"] = qc_display(base)
        return base
    m = glm_probit([([1.0, p["t"]], p["n"], p["g"]) for p in pts])
    if m is None:
        base["fit_reason"] = "probit fit failed"
        base["qc_display_status"] = qc_display(base)
        return base
    b0, b1 = m["coef"]
    ki_pct = normcdf(b0) * 100
    declining = math.isfinite(b1) and b1 < 0
    slope_se = math.sqrt(max(0.0, m["cov"][1][1])) if m["cov"] else None
    slope_p = 2 * (1 - normcdf(abs(b1 / slope_se))) if slope_se and slope_se > 0 else None
    sigma = -1 / b1 if declining else None
    p50 = crossing(b0, b1, 0.5) if declining and ki_pct >= 50 else None
    p70a = crossing(b0, b1, 0.7) if declining and ki_pct >= 70 else None
    ug = (crossing(b0, b1, threshold_pct / 100)
          if declining and threshold_pct is not None and 0 < threshold_pct < 100 else None)
    status = ("non-declining fitted slope; longevity crossing not interpreted" if not declining
              else "RCI-P70 threshold unavailable" if threshold_pct is None
              else "threshold crossing unavailable" if ug is None or not math.isfinite(ug)
              else "fitted viability already below RCI P70 at model time 0" if ug < 0
              else "crossing within observed storage-time range" if ug <= max_t
              else "crossing extrapolated beyond observed storage-time range")
    cat = ("non-declining slope" if not math.isfinite(b1) or b1 >= 0
           else "insufficient time points" if n_time < 2
           else "crossing unavailable" if ug is None or not math.isfinite(ug)
           else "fitted viability below threshold at model time 0" if ug < 0
           else "crossing within observed time range" if ug <= max_t
           else "crossing extrapolated beyond observed time range")
    tier = ("do not interpret" if not math.isfinite(b1) or b1 >= 0 or n_time < 2
            or ug is None or not math.isfinite(ug) or ug < 0
            else "exploratory only" if n_time == 2
            else "primary model evidence" if ug <= max_t
            else "extrapolated model evidence")
    display = ("DO NOT INTERPRET" if tier == "do not interpret"
               else "EXPLORATORY" if tier == "exploratory only"
               else "FORMAL SURVIVAL")
    base.update({
        "fit_ok": True, "fit_method": "irls-probit",
        "fit_reason": "ok" if declining else "non-declining slope",
        "fit_warning": ("possible separation: fit did not stabilize at finite estimates" if m["separated"]
                        else "fitted probabilities numerically 0 or 1 occurred" if m["extreme_prob"] else None),
        "extreme_prob_warning": m["extreme_prob"], "converged": m["converged"],
        "Ki_probit": b0, "Ki_percent": ki_pct, "slope": b1,
        "slope_se": slope_se, "slope_p": slope_p,
        "sigma_years": sigma, "p50_years": p50, "p70_absolute_years": p70a,
        "fitted_rci_p70_years_raw": ug,
        "fitted_rci_p70_display": ("no RCI-P70 threshold (accession baseline unavailable)"
                                   if threshold_pct is None else crossing_display(ug)),
        "crossing_status": status, "bridge_category": cat,
        "bridge_evidence_tier": tier, "display_class": display,
    })
    base["qc_display_status"] = qc_display(base)
    return base

def thresholds_from_baselines(baselines):
    """(family, taxon, accession) → RCI-P70 threshold %, from the RCI branch's
    accession early baselines — the ONE genuine contingency of the formal
    workflow on the RCI stage."""
    return {k: v["threshold"] for k, v in baselines.items()}

# ============================ command line ============================
FORMAL_COLUMNS = [
    "taxon", "accession", "storage_code", "drying", "protocol_key",
    "substrate", "treatment", "display_class", "fit_method", "fit_warning",
    "qc_display_status", "bridge_evidence_tier", "bridge_category",
    "n_model_rows", "n_timepoints", "n_post_timepoints", "n_seeds",
    "matched_initial_available", "Ki_probit", "Ki_percent", "slope",
    "slope_se", "slope_p", "sigma_years", "p50_years", "p70_absolute_years",
    "threshold_pct", "fitted_rci_p70_years_raw", "fitted_rci_p70_display",
    "crossing_status", "early_rise_class", "plateau_note",
]

def _load_thresholds_csv(path):
    import csv as _csv
    thr = {}
    with open(path, newline="", encoding="utf-8-sig") as fh:
        for r in _csv.DictReader(fh):
            try:
                thr[(r.get("family") or "", r["taxon"], r["accession"])] = float(r["threshold"])
            except (KeyError, ValueError, TypeError):
                continue
    return thr

def main():
    import argparse, os
    ap = argparse.ArgumentParser(description=(
        "Component 3 — formal protocol-specific survival models (raw CSV or "
        "normalized_assays.csv in). Trajectory = taxon × accession × storage "
        "× germination protocol; matched same-protocol Initial anchors t=0; "
        "actual counts only; Formal/Exploratory/Do-not-interpret/No-fit "
        "classes with early-rise, plateau and separation QC; fitted RCI-P70 "
        "crossing kept raw AND displayed safely. CONTINGENCY: the RCI-P70 "
        "bridge threshold comes from the RCI component's accession baselines "
        "— computed internally by default, or supplied via --thresholds "
        "rci_accession_intervals.csv (from longevity_rci.py)."))
    ap.add_argument("input")
    ap.add_argument("-o", "--outdir", default=".")
    ap.add_argument("--time-unit", choices=["days", "months", "years"])
    ap.add_argument("--thresholds", help="rci_accession_intervals.csv from longevity_rci.py "
                                         "(otherwise baselines are recomputed from this input)")
    args = ap.parse_args()
    assays, parsed, was_norm = load_input(args.input, time_unit=args.time_unit)
    report_issues(parsed)
    counted = [a for a in assays if a["tested"] is not None]
    if not counted:
        print("[BLOCK] counts required for formal/binomial modelling — no sample size is ever assumed.", file=sys.stderr)
        sys.exit(2)
    if args.thresholds:
        thr = _load_thresholds_csv(args.thresholds)
    else:
        thr = thresholds_from_baselines(early_baselines(_rci_rows(assays)))
    g = group_formal_trajectories(assays)
    for x in g["excluded_late_initials"]:
        print(f"[NOTE] Initial at {x['years']:.1f} y ({x['taxon']} {x['accession']}, {x['row_id']}) "
              "excluded from t=0 anchoring (older than 1 y).", file=sys.stderr)
    if g["drying_keyed"]:
        print("[NOTE] drying treatment varies — added as a grouping dimension.", file=sys.stderr)
    results = [fit_formal_survival(t, thr.get((t["family"] or "", t["taxon"], t["accession"])))
               for t in g["trajectories"]]
    os.makedirs(args.outdir, exist_ok=True)
    p1 = os.path.join(args.outdir, "formal_survival_results.csv")
    write_csv(p1, FORMAL_COLUMNS, results)
    n_formal = sum(1 for r in results if r["display_class"] == "FORMAL SURVIVAL")
    print(f"{p1}: {len(results)} trajectories ({n_formal} FORMAL SURVIVAL)")

if __name__ == "__main__":
    main()