"""
Three referee-requested reanalyses, all from frozen local data. No network access.

1. CITATION REMATCH on baseline citation visibility.
   The original design matched on publication year, retraction year and publisher, then
   discovered residual imbalance in citations accumulated through the retraction year
   (SMD -0.152) and handled it after the fact with calipers. The referee asked for the
   matching to be rebuilt with baseline citations built in. Because baseline citations are
   only observable after retrieval, this uses the frozen candidate pool (t39) rather than
   new API calls. Reports: exact-baseline matching, conditional-logistic odds ratio on
   discordant pairs, match-tier subgroups, and stability across repeated random assignments.

2. TERMINOLOGY AUDIT, broadened.
   Adds spacing, plural, hyphenation and version variants the referee listed. Emits the
   exact regular expressions so they can be published.

3. BOOTSTRAP STABILITY at 10,000 replicates with Monte Carlo error on the endpoints.
"""

# _REPRO_PATHS_ : resolve paths relative to the repository root, overridable via
# the REPRO_ROOT environment variable. No absolute paths are baked in.
import os as _os
_HERE = _os.path.dirname(_os.path.abspath(__file__))
_ROOT = _os.environ.get("REPRO_ROOT") or _os.path.dirname(_HERE)
def _p(*parts):
    return _os.path.join(_ROOT, *parts)
import csv
import os
import re
import sys
import math
import collections
from datetime import datetime

import numpy as np

csv.field_size_limit(10**9)

# CANONICAL date formats, identical in every script in this package. The AM/PM variant is
# required by exactly one record (ID 18930, a 1756 retraction); omitting it silently drops
# that record and yields 66,012 rather than 66,013 analysed retractions.
_DATE_FORMATS = ("%m/%d/%Y %I:%M:%S %p", "%m/%d/%Y %H:%M:%S",
                 "%m/%d/%Y %H:%M", "%m/%d/%Y", "%Y-%m-%d", "%d/%m/%Y")
SEED = 20260731
N_BOOT = 10000
N_ASSIGN = 20                      # repeated random control assignments
RAW = _p("data", "retraction_watch_gitlab.csv")
T = _p("runs", "ai-retractions-v3", "tables")
REASON_PAT = ("computer-aided content", "computer-generated content")


def rd(name):
    with open(os.path.join(T, name), encoding="utf-8-sig", newline="") as f:
        return list(csv.DictReader(f))


def wr(name, rows, header=None):
    if not rows:
        return
    keys = header or list(rows[0].keys())
    with open(os.path.join(T, name), "w", newline="", encoding="utf-8") as f:
        w = csv.DictWriter(f, fieldnames=keys, quoting=csv.QUOTE_ALL)
        w.writeheader()
        w.writerows(rows)
    print(f"  wrote {name}")


def wilson(k, n):
    if n == 0:
        return (0.0, 0.0)
    z = 1.959963984540054
    p = k / n
    d = 1 + z * z / n
    c = (p + z * z / (2 * n)) / d
    h = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d
    return (100 * max(0, c - h), 100 * min(1, c + h))


def mcnemar(b, c):
    """b = machine-only, c = control-only. Returns (chi2, p) with continuity correction."""
    if b + c == 0:
        return 0.0, 1.0
    chi2 = (abs(b - c) - 1) ** 2 / (b + c)
    # survival of chi2 with 1 df
    p = math.erfc(math.sqrt(max(chi2, 0) / 2))
    return chi2, p


def paired_diff_ci(b, c, n):
    """Agresti-Min style Wald interval for paired proportion difference."""
    d = (b - c) / n
    var = (b + c - (b - c) ** 2 / n) / (n * n)
    se = math.sqrt(max(var, 0))
    return 100 * d, 100 * (d - 1.959964 * se), 100 * (d + 1.959964 * se)


# =====================================================================
# 1. CITATION REMATCH
# =====================================================================
def citation_rematch():
    print("\n=== 1. CITATION REMATCH ON BASELINE CITATIONS ===")
    cites = rd("t19_citations_full.csv")
    cand = rd("t39_control_candidates_openalex.csv")
    topics = {r["Doi"].lower(): r.get("Field", "") for r in rd("t42_openalex_topics.csv")}

    # Publisher for every machine-attributed record, taken from the Retraction Watch
    # snapshot itself. An earlier implementation read this from the initial matched-pair
    # table, which contains only records that were successfully paired in the weaker first
    # design; that made the repaired design partly conditional on success in the design it
    # was meant to replace. Reading from the snapshot removes that dependence.
    pub_of = {}
    with open(RAW, encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            if (row.get("RetractionNature") or "").strip() != "Retraction":
                continue
            if not any(p in (row.get("Reason") or "").lower() for p in REASON_PAT):
                continue
            doi = (row.get("OriginalPaperDOI") or "").strip().lower()
            if doi and doi != "unavailable":
                pub_of[doi] = (row.get("Publisher") or "").strip()
    print(f"publisher map built from snapshot for {len(pub_of)} machine-attributed DOIs "
          f"(previously {len(rd('t18_matched_pairs.csv'))} from the matched-pair table)")

    def num(x, d=0):
        try:
            return int(float(x))
        except (TypeError, ValueError):
            return d

    machine = []
    seen_m = set()
    dup_m = 0
    for r in cites:
        if r.get("Cohort") != "machine":
            continue
        doi = r["Doi"].lower()
        # Retain only the first record per DOI, matching the deduplication applied in the
        # initial design. The frozen extract contains repeated rows where a work has more
        # than one OpenAlex entry.
        if doi in seen_m:
            dup_m += 1
            continue
        seen_m.add(doi)
        machine.append({
            "doi": doi, "py": num(r["PubYear"]), "ry": num(r["RetractionYear"]),
            "exp": num(r["ExposureYears"]), "base": num(r["UpToRetraction"]),
            "after": num(r["AfterRetraction"]), "y1": num(r["CitesY1"]),
            "pub": pub_of.get(doi, ""), "field": topics.get(doi, ""),
        })

    controls = []
    seen_c = set()
    dup_c = 0
    for r in cand:
        doi = r["Doi"].lower()
        if doi in seen_c:
            dup_c += 1
            continue
        seen_c.add(doi)
        controls.append({
            "doi": doi, "py": num(r["PubYear"]), "ry": num(r["RetractionYear"]),
            "exp": num(r["ExposureYears"]), "base": num(r["UpToRetraction"]),
            "after": num(r["AfterRetraction"]), "y1": num(r["CitesY1"]),
            "pub": (r.get("Publisher") or "").strip(),
            "field": topics.get(doi, r.get("Field", "") or ""),
        })
    mdoi = {m["doi"] for m in machine}
    controls = [c for c in controls if c["doi"] not in mdoi]
    print(f"machine={len(machine)} (dropped {dup_m} duplicate DOI rows)  "
          f"control candidates={len(controls)} (dropped {dup_c} duplicate DOI rows)")

    # Diagnostic: how many machine records now carry a publisher, and how the records that
    # remain without one differ. This is reported so the eligible population for the
    # publisher-exact design can be judged independently of the initial matching.
    with_pub = [m for m in machine if m["pub"]]
    without = [m for m in machine if not m["pub"]]
    print(f"  machine records with a publisher value: {len(with_pub)}; without: {len(without)}")
    if without:
        def summarise(rs, label):
            if not rs:
                return
            b = np.array([r["base"] for r in rs], float)
            print(f"    {label:14s} n={len(rs):5d}  median pub year "
                  f"{int(np.median([r['py'] for r in rs]))}  median ret year "
                  f"{int(np.median([r['ry'] for r in rs]))}  mean baseline cites {b.mean():.2f}")
        summarise(with_pub, "with publisher")
        summarise(without, "without")

    def run_assignment(seed, require_publisher):
        """Exact match on (py, ry, baseline citations); publisher exact if required."""
        rng = np.random.default_rng(seed)
        pool = collections.defaultdict(list)
        for c in controls:
            key = (c["py"], c["ry"], c["base"], c["pub"]) if require_publisher \
                else (c["py"], c["ry"], c["base"])
            pool[key].append(c)
        for k in pool:
            rng.shuffle(pool[k])
        out = []
        for m in machine:
            key = (m["py"], m["ry"], m["base"], m["pub"]) if require_publisher \
                else (m["py"], m["ry"], m["base"])
            p = pool.get(key)
            if p:
                out.append((m, p.pop()))
        return out

    results = []
    for label, reqpub in (("baseline + year, publisher exact", True),
                          ("baseline + year, publisher free", False)):
        mp = run_assignment(SEED, reqpub)
        if not mp:
            continue
        elig = [(m, c) for m, c in mp if m["exp"] >= 1 and c["exp"] >= 1]
        b = sum(1 for m, c in elig if m["after"] > 0 and c["after"] == 0)
        cc = sum(1 for m, c in elig if m["after"] == 0 and c["after"] > 0)
        both = sum(1 for m, c in elig if m["after"] > 0 and c["after"] > 0)
        nei = len(elig) - b - cc - both
        n = len(elig)
        pm = 100 * (b + both) / n
        pc = 100 * (cc + both) / n
        diff, lo, hi = paired_diff_ci(b, cc, n)
        chi2, p = mcnemar(b, cc)
        # conditional logistic OR for 1:1 matched binary data = b/c on discordant pairs
        orr = (b / cc) if cc else float("inf")
        se_lnor = math.sqrt(1 / b + 1 / cc) if b and cc else float("nan")
        or_lo = orr * math.exp(-1.959964 * se_lnor) if b and cc else float("nan")
        or_hi = orr * math.exp(1.959964 * se_lnor) if b and cc else float("nan")
        # baseline SMD check
        mb = np.array([m["base"] for m, _ in elig], float)
        cb = np.array([c["base"] for _, c in elig], float)
        sd = math.sqrt((mb.var(ddof=1) + cb.var(ddof=1)) / 2)
        smd = (mb.mean() - cb.mean()) / sd if sd > 0 else 0.0
        print(f"\n  [{label}]")
        print(f"    pairs={n}  machine={pm:.2f}%  control={pc:.2f}%")
        print(f"    paired diff={diff:+.2f} pp (95% CI {lo:.2f} to {hi:.2f})  "
              f"McNemar p={p:.3g}")
        print(f"    conditional OR={orr:.3f} (95% CI {or_lo:.3f} to {or_hi:.3f})  "
              f"discordant b={b} c={cc}")
        print(f"    baseline SMD={smd:+.4f}  (machine mean {mb.mean():.2f}, "
              f"control mean {cb.mean():.2f})")
        results.append({
            "MatchingSpecification": label, "Pairs": n,
            "MachinePct": round(pm, 2), "ControlPct": round(pc, 2),
            "PairedDiffPP": round(diff, 2), "CiLo": round(lo, 2), "CiHi": round(hi, 2),
            "McNemarP": f"{p:.3g}",
            "ConditionalOR": round(orr, 3) if math.isfinite(orr) else "",
            "OrCiLo": round(or_lo, 3) if math.isfinite(or_lo) else "",
            "OrCiHi": round(or_hi, 3) if math.isfinite(or_hi) else "",
            "DiscordantMachineOnly": b, "DiscordantControlOnly": cc,
            "BaselineSMD": round(smd, 4),
        })

    # stability across repeated random control assignments
    print(f"\n  [stability across {N_ASSIGN} random control assignments]")
    diffs = []
    npairs = set()
    for i in range(N_ASSIGN):
        mp = run_assignment(SEED + 1000 + i, True)
        elig = [(m, c) for m, c in mp if m["exp"] >= 1 and c["exp"] >= 1]
        if not elig:
            continue
        npairs.add(len(elig))
        b = sum(1 for m, c in elig if m["after"] > 0 and c["after"] == 0)
        cc = sum(1 for m, c in elig if m["after"] == 0 and c["after"] > 0)
        d, _, _ = paired_diff_ci(b, cc, len(elig))
        diffs.append(d)
    diffs = np.array(diffs)
    # The eligible pair count is fixed by the match keys; only which eligible control is
    # drawn varies, so report the exact count rather than an average.
    pair_note = (str(npairs.pop()) if len(npairs) == 1
                 else f"{min(npairs)}-{max(npairs)}")
    print(f"    eligible pairs: {pair_note} in every assignment")
    print(f"    paired difference: mean={diffs.mean():+.2f} pp  sd={diffs.std(ddof=1):.2f}  "
          f"range {diffs.min():+.2f} to {diffs.max():+.2f}")
    results.append({
        "MatchingSpecification": f"stability, {N_ASSIGN} random assignments "
                                 f"(publisher exact)",
        "Pairs": pair_note, "MachinePct": "", "ControlPct": "",
        "PairedDiffPP": round(float(diffs.mean()), 2),
        "CiLo": round(float(diffs.min()), 2), "CiHi": round(float(diffs.max()), 2),
        "McNemarP": f"sd={diffs.std(ddof=1):.2f}", "ConditionalOR": "",
        "OrCiLo": "", "OrCiHi": "", "DiscordantMachineOnly": "",
        "DiscordantControlOnly": "", "BaselineSMD": "",
    })
    wr("t47_citation_rematched.csv", results)


# =====================================================================
# 2. TERMINOLOGY AUDIT
# =====================================================================
PATTERNS = collections.OrderedDict([
    ("ChatGPT (incl. spaced/hyphenated)", r"\bchat\s*[-]?\s*gpt\b"),
    ("GPT with version number",           r"\bgpt[\s\-]?[2-9](?:\.\d)?\b"),
    ("GPT bare token",                    r"\bgpt\b"),
    ("LLM / LLMs",                        r"\bllms?\b"),
    ("large language model(s)",           r"\blarge[\s\-]language[\s\-]models?\b"),
    ("language model(s)",                 r"\blanguage[\s\-]models?\b"),
    ("generative AI / GenAI",             r"\b(?:generative\s+a\.?i\.?|gen[\s\-]?ai)\b"),
    ("AI-generated / AI generated",       r"\bai[\s\-]generated\b"),
    ("machine-generated / machine written", r"\bmachine[\s\-](?:generated|written)\b"),
    ("automatically generated",           r"\bautomatic(?:ally)?[\s\-]generated\b"),
    ("automated (manuscript )generation", r"\bautomated\s+(?:\w+\s+)?generation\b"),
    ("synthetic text",                    r"\bsynthetic\s+text\b"),
    ("neural text generator",             r"\bneural\s+text\s+generat\w+\b"),
    ("hallucinated citations/references", r"\bhallucinat\w+\s+(?:citation|reference)s?\b"),
    ("generated references",              r"\bgenerated\s+references?\b"),
    ("named assistants",                  r"\b(?:bard|gemini|copilot|claude|llama|mistral"
                                          r"|deepseek|qwen|grok)\b"),
    ("SCIgen / Mathgen",                  r"\b(?:sci\s*gen|math\s*gen)\b"),
    ("tortured phrase(s)",                r"\btortured\s+phrase\w*\b"),
    ("paraphrasing tool / spinner",       r"\b(?:paraphras\w+|spinbot|article\s+spinner"
                                          r"|text\s+spinn\w+)\b"),
    ("artificial intelligence (generic)", r"\bartificial\s+intelligence\b"),
])


def terminology():
    print("\n=== 2. BROADENED TERMINOLOGY AUDIT ===")
    comp = {k: re.compile(v, re.I) for k, v in PATTERNS.items()}
    specific = [k for k in PATTERNS if k != "artificial intelligence (generic)"]

    counts = collections.Counter()
    machine_counts = collections.Counter()
    any_hit = set()
    any_specific = set()
    detail = []

    with open(RAW, encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            if (row.get("RetractionNature") or "").strip() != "Retraction":
                continue
            rid = row.get("Record ID") or ""
            hay = ((row.get("Title") or "") + " || " + (row.get("Notes") or ""))
            is_m = any(p in (row.get("Reason") or "").lower() for p in REASON_PAT)
            hits = [k for k, rx in comp.items() if rx.search(hay)]
            if not hits:
                continue
            for k in hits:
                counts[k] += 1
                if is_m:
                    machine_counts[k] += 1
            any_hit.add(rid)
            sp = [k for k in hits if k in specific]
            if sp:
                any_specific.add(rid)
                detail.append({
                    "RecordId": rid, "MachineAttributed": int(is_m),
                    "TermsMatched": "; ".join(sp),
                    "Title": (row.get("Title") or "")[:300],
                    "NotesSnippet": (row.get("Notes") or "")[:300],
                })

    print(f"records matching ANY term: {len(any_hit)}")
    print(f"records matching a SPECIFIC generative term: {len(any_specific)}")
    print("\n  per-pattern counts (all retractions / machine-attributed):")
    out = []
    for k in PATTERNS:
        print(f"    {counts[k]:6d} / {machine_counts[k]:5d}   {k}")
        out.append({"Term": k, "Regex": PATTERNS[k],
                    "RecordsAllRetractions": counts[k],
                    "RecordsMachineAttributed": machine_counts[k]})
    out.append({"Term": "ANY term above", "Regex": "(union)",
                "RecordsAllRetractions": len(any_hit), "RecordsMachineAttributed": ""})
    out.append({"Term": "ANY specific generative term", "Regex": "(union, excl. generic AI)",
                "RecordsAllRetractions": len(any_specific), "RecordsMachineAttributed": ""})
    wr("t48_terminology_broadened.csv", out)
    wr("t49_terminology_specific_matches.csv", detail)


# =====================================================================
# 3. BOOTSTRAP STABILITY
# =====================================================================
def bootstrap():
    print(f"\n=== 3. CLUSTER BOOTSTRAP AT {N_BOOT} REPLICATES ===")

    def parse_date(s):
        s = (s or "").strip()
        for fmt in _DATE_FORMATS:
            try:
                return datetime.strptime(s, fmt)
            except ValueError:
                continue
        return None

    recs = []          # (cluster_key, is_machine, lag, mill)
    with open(RAW, encoding="utf-8-sig", newline="") as f:
        for row in csv.DictReader(f):
            if (row.get("RetractionNature") or "").strip() != "Retraction":
                continue
            od = parse_date(row.get("OriginalPaperDate"))
            rdt = parse_date(row.get("RetractionDate"))
            if not od or not rdt:
                continue
            lag = (rdt - od).days
            if lag < 0:
                continue
            reason = (row.get("Reason") or "").lower()
            # A retraction DOI of blank, 'unavailable' or 'Unavailable' carries no notice
            # identity, so each such record is its own singleton cluster. Treating the
            # literal string as a shared key would pool thousands of unrelated records
            # into one cluster and badly distort the cluster bootstrap.
            raw = (row.get("RetractionDOI") or "").strip()
            if raw == "" or raw.lower() == "unavailable":
                key = f"__singleton_{row.get('Record ID')}"
            else:
                key = raw.lower()
            recs.append((key, any(p in reason for p in REASON_PAT), lag,
                         "paper mill" in reason))

    keys = sorted({r[0] for r in recs})
    kidx = {k: i for i, k in enumerate(keys)}
    K = len(keys)
    print(f"clusters (notices) = {K}")

    # Cluster bootstrap == weighting each cluster by its multinomial draw count.
    # Pre-sort each arm by lag so the weighted median is a searchsorted lookup.
    def arm(is_m):
        sub = [(r[2], kidx[r[0]], r[3]) for r in recs if r[1] == is_m]
        sub.sort(key=lambda t: t[0])
        return (np.array([t[0] for t in sub], float),
                np.array([t[1] for t in sub], np.int64),
                np.array([t[2] for t in sub], float))

    lag_m, clu_m, mill_m = arm(True)
    lag_c, clu_c, mill_c = arm(False)
    print(f"machine={lag_m.size}  comparison={lag_c.size}")

    def wstats(lag, clu, mill, cnt):
        w = cnt[clu].astype(float)
        tot = w.sum()
        if tot <= 0:
            return None
        cw = np.cumsum(w)
        med = lag[min(int(np.searchsorted(cw, tot / 2.0)), lag.size - 1)]
        mad = float(np.dot(w, np.abs(lag - med)) / tot)
        return float(med), mad, float(np.dot(w, mill) / tot)

    def triple(cnt):
        a = wstats(lag_m, clu_m, mill_m, cnt)
        b = wstats(lag_c, clu_c, mill_c, cnt)
        if a is None or b is None:
            return (np.nan, np.nan, np.nan)
        return (a[0] - b[0], a[1] - b[1], 100.0 * (a[2] - b[2]))

    obs = triple(np.ones(K, np.int64))
    rng = np.random.default_rng(SEED)
    probs = np.full(K, 1.0 / K)
    boot = np.empty((N_BOOT, 3))
    for i in range(N_BOOT):
        boot[i] = triple(rng.multinomial(K, probs))

    names = ["Median lag difference (days)",
             "Mean-absolute-deviation difference (days)",
             "Paper-mill co-coding difference (percentage points)"]
    rows = []
    for j, nm in enumerate(names):
        col = boot[~np.isnan(boot[:, j]), j]
        lo, hi = np.percentile(col, [2.5, 97.5])
        # Monte Carlo SE of each endpoint, from 4 independent quarters of the replicates
        q = [np.percentile(col[k::4], [2.5, 97.5]) for k in range(4)]
        mc_lo = float(np.std([x[0] for x in q], ddof=1)) / 2.0
        mc_hi = float(np.std([x[1] for x in q], ddof=1)) / 2.0
        print(f"  {nm}")
        print(f"    point={obs[j]:.2f}  95% CI {lo:.2f} to {hi:.2f}  "
              f"(Monte Carlo SE of endpoints: {mc_lo:.3f}, {mc_hi:.3f})")
        rows.append({"Estimate": nm, "PointEstimate": round(float(obs[j]), 2),
                     "ClusterCiLo": round(float(lo), 2), "ClusterCiHi": round(float(hi), 2),
                     "McSeCiLo": round(mc_lo, 3), "McSeCiHi": round(mc_hi, 3),
                     "Replicates": N_BOOT})
    wr("t50_bootstrap_10k.csv", rows)


if __name__ == "__main__":
    which = sys.argv[1] if len(sys.argv) > 1 else "all"
    if which in ("all", "cite"):
        citation_rematch()
    if which in ("all", "term"):
        terminology()
    if which in ("all", "boot"):
        bootstrap()
    print("\ndone")
