#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""longevity_normalize.py — Component 1 · the shared upload/normalization layer

Part of the Seed Longevity Lab's downloadable Python lab (site v0.17.0,
longevity.seedsofhawaii.org/scripts). No upstream dependency — every other component builds on this one's output (or embeds it).

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

# ============================ command line ============================
def main():
    import argparse
    ap = argparse.ArgumentParser(description=(
        "Component 1 — normalize a completed-test CSV into the Longevity Lab "
        "interchange format (normalized_assays.csv). One row per completed "
        "germination test; counts required; time kept exact; Initial made "
        "explicit; protocol keys built; duplicates flagged, never merged. "
        "Every other script in this suite accepts either your raw CSV or "
        "this file."))
    ap.add_argument("input", help="raw completed-test CSV/TSV")
    ap.add_argument("-o", "--out", default="normalized_assays.csv")
    ap.add_argument("--time-unit", choices=["days", "months", "years"],
                    help="unit for an ambiguous time column (never guessed)")
    ap.add_argument("--allow-percent-only", action="store_true",
                    help="write the normalized file for INSPECTION even when "
                         "counts are missing (analytical scripts will still refuse it)")
    args = ap.parse_args()
    with open(args.input, newline="", encoding="utf-8-sig") as fh:
        parsed = parse_completed_tests(fh.read(), time_unit=args.time_unit)
    report_issues(parsed, strict=not args.allow_percent_only)
    rows = assays_to_normalized_rows(parsed["assays"])
    write_csv(args.out, NORMALIZED_COLUMNS, rows)
    for g in parsed["duplicate_groups"]:
        print(f"[REVIEW] possible duplicate/multiple same-date tests: {g['taxon']} {g['accession']} "
              f"{g['storage']} @ {g['date_key']} rows {', '.join(g['row_ids'])}", file=sys.stderr)
    print(f"{args.out}: {len(rows)} assays normalized "
          f"({len(parsed['rejected'])} rejected, {len(parsed['duplicate_groups'])} duplicate group(s) flagged)")

if __name__ == "__main__":
    main()