#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""longevity_rci.py — Component 2 · RCI management analysis

Part of the Seed Longevity Lab's downloadable Python lab (site v0.17.0,
longevity.seedsofhawaii.org/scripts). Independent (embeds Component 1). Its accession-baseline stage is the one thing Component 3 depends on.

Standalone by design: the shared cores below are embedded verbatim in every
script of the suite, 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}

# ============================ command line ============================
RCI_RESULT_COLUMNS = [
    "family", "taxon", "recommended_storage", "storage_behavior",
    "recommendation_reason", "overall_observed_rci_years_raw",
    "overall_observed_rci_display", "overall_driver_accession",
    "overall_driver_storage", "rci_at_recommended_years_raw",
    "rci_at_recommended_display", "rci_at_recommended_driver_accession",
    "rebound_failure_flags", "accession_baseline_percent",
    "rci_p70_threshold_percent", "fresh_seeds_baseline", "recent_germ_at_rci",
    "supporting_accessions", "potentially_supporting_accessions",
    "supporting_published", "rci_eligible_accessions", "raw_accessions",
    "override_to_minus18", "R5", "R18", "freeze_sensitive", "freeze_reason",
]
INTERVAL_COLUMNS = [
    "family", "taxon", "accession", "storage", "baseline_year",
    "baseline_viab", "threshold", "n_years", "max_year_meeting_p70",
    "glm_support", "wilson_support", "glm_lo", "wilson_lo",
]

def main():
    import argparse
    ap = argparse.ArgumentParser(description=(
        "Component 2 — RCI management analysis (raw CSV or "
        "normalized_assays.csv in). Empirical early-window baseline, RCI P70, "
        "observed intervals with GLM+Wilson support, species roll-up with "
        "Overall observed RCI SEPARATE from observed RCI at recommended "
        "storage, freeze flags, censoring-honest notation. Public upload rule "
        "(no archival two-point exception). Self-contained: embeds the "
        "normalize component."))
    ap.add_argument("input")
    ap.add_argument("-o", "--outdir", default=".")
    ap.add_argument("--time-unit", choices=["days", "months", "years"])
    args = ap.parse_args()
    assays, parsed, was_norm = load_input(args.input, time_unit=args.time_unit)
    report_issues(parsed)
    if not any(a["tested"] is not None for a in assays):
        print("[BLOCK] counts required — this input has no sown/germinated counts.", file=sys.stderr)
        sys.exit(2)
    res = run_rci_pipeline(assays)
    import os
    os.makedirs(args.outdir, exist_ok=True)
    p1 = os.path.join(args.outdir, "rci_results.csv")
    write_csv(p1, RCI_RESULT_COLUMNS, res["publication"])
    ivals = [{**{c: it.get(c) for c in INTERVAL_COLUMNS}} for it in res["intervals"]]
    p2 = os.path.join(args.outdir, "rci_accession_intervals.csv")
    write_csv(p2, INTERVAL_COLUMNS, ivals)
    print(f"{p1}: {len(res['publication'])} species")
    print(f"{p2}: {len(res['intervals'])} accession × storage series "
          "(baseline/threshold columns feed longevity_formal.py --thresholds)")

if __name__ == "__main__":
    main()