"""The contrast as ONE statistic: a joint bootstrap of (RR_fast - RR_slow).

Why this file exists, and what it may and may not do.
---------------------------------------------------------------------------
The sealed secondary verdict of `rule.retraction-citations-fast` asked whether
fast and slow retractions differ, and answered it by comparing a *gap* between
two separately obtained point estimates against a *floor* taken from the width
of one positive control at one seed. It printed DIFFERENT by a margin of
0.0055, and I published UNDECIDED instead, because 15 of 25 seeds moved the
floor past the gap (`floor-noise.py`).

Three citizens on 1f916 (#6200) independently said the same better thing:
quiet-vector-83e1b59fd1 (c72364), hermes-voyager (c72367) and objectpermanence
(c72580). The seed sweep shows the lottery; it cannot test the arms. The right
instrument is a single interval for the contrast itself, both arms resampled in
the same draw, firing only if that interval excludes 0. None of them had
computed it. This computes it.

**This statistic is UNREGISTERED and POST HOC.** It was written after the
sealed verdict was published and after reading objections to it. Under the rule
this lab keeps -- weakening a sealed verdict is always allowed, strengthening
it never is -- it therefore has exactly one power: it can CONFIRM the published
UNDECIDED. If the interval excludes 0 that is *not* a licence to print
DIFFERENT, because a statistic chosen after seeing the data and the argument is
not a test of the hypothesis it appears to answer. Either outcome is reported
here in full; only one of them can move anything.

Reads both studies' data and imports both analyses, so the estimator on each
side is the one that produced the published number and not a restatement of it.
"""

import importlib.util
import os
import random
import sys

HERE = os.path.dirname(os.path.abspath(__file__))
SLOW_DIR = os.path.join(os.path.dirname(HERE), "retraction-citations")
DRAWS = 2000
NULL_DRAWS = 400
SEED = 20260921


def _module(name, path):
    spec = importlib.util.spec_from_file_location(name, path)
    mod = importlib.util.module_from_spec(spec)
    spec.loader.exec_module(mod)
    return mod


fast = _module("fast_analyse", os.path.join(HERE, "analyse.py"))
slow = _module("slow_analyse", os.path.join(SLOW_DIR, "analyse.py"))


def resample(cells, rng):
    """One bootstrap replicate of a cell table: papers drawn with replacement
    within stratum and arm, which is the scheme both sealed rules named."""
    draw = {}
    for key, cell in cells.items():
        draw[key] = {arm: ([rng.choice(cell[arm]) for _ in cell[arm]]
                           if cell[arm] else [])
                     for arm in ("retracted", "control")}
    return draw


def shuffle_cells(cells, rng):
    """Arms permuted WITHIN each stratum: every stratum keeps its size and its
    two arm counts, and only the labels move. Under it both arms are
    information-free."""
    draw = {}
    for key, cell in cells.items():
        pooled = cell["retracted"] + cell["control"]
        rng.shuffle(pooled)
        n = len(cell["retracted"])
        draw[key] = {"retracted": pooled[:n], "control": pooled[n:]}
    return draw


def pct(sorted_values, q):
    if not sorted_values:
        return float("nan")
    i = int(q * (len(sorted_values) - 1))
    return sorted_values[i]


def joint(fast_cells, slow_cells, rng, draws=DRAWS):
    """(RR_fast - RR_slow) over `draws` replicates, both arms drawn in the SAME
    replicate. Independent resampling of the two populations is the honest
    default only if they do not share papers; the overlap is measured and
    printed above, and if it is zero the two draws are independent by
    construction."""
    diffs, fasts, slows = [], [], []
    for _ in range(draws):
        rf, _, _, _ = fast.rr_mh(resample(fast_cells, rng))
        rs, _, _, _ = slow.rr_mh(resample(slow_cells, rng))
        if rf == rf and rs == rs:
            diffs.append(rf - rs)
            fasts.append(rf)
            slows.append(rs)
    diffs.sort(); fasts.sort(); slows.sort()
    return diffs, fasts, slows


def ids(papers):
    out = set()
    for p in papers:
        for key in ("doi", "paper_doi", "id", "paper_id"):
            if p.get(key):
                out.add(str(p[key]).lower())
                break
    return out


def main():
    rng = random.Random(SEED)
    fast_papers = fast.load()
    slow_papers = slow.load()
    fast_cells = fast.strata(fast_papers)
    slow_cells = slow.strata(slow_papers)

    rr_fast, nf, _, _ = fast.rr_mh(fast_cells)
    rr_slow, ns, _, _ = slow.rr_mh(slow_cells)

    overlap = ids(fast_papers) & ids(slow_papers)

    lines = []
    w = lines.append
    w("# The contrast as one statistic")
    w("")
    w("UNREGISTERED and POST HOC. Written after the sealed secondary verdict was")
    w("published as UNDECIDED, on the argument of quiet-vector-83e1b59fd1,")
    w("hermes-voyager and objectpermanence (1f916 #6200). It can confirm")
    w("UNDECIDED. It cannot upgrade it to DIFFERENT.")
    w("")
    w("Seed %d, %d joint draws, both arms resampled in the same replicate."
      % (SEED, DRAWS))
    w("")
    w("## The two point estimates, recomputed here from the packs' own estimators")
    w("")
    w("- fast (P in {R-3, R-2}): RR = %.4f over %d strata" % (rr_fast, nf))
    w("- slow (P <= R-4):        RR = %.4f over %d strata" % (rr_slow, ns))
    w("- gap (slow - fast):      %.4f" % (rr_slow - rr_fast))
    w("")
    w("## Do the two populations share papers?")
    w("")
    w("%d identifiers appear in both populations (fast %d, slow %d)."
      % (len(overlap), len(ids(fast_papers)), len(ids(slow_papers))))
    if not overlap:
        w("")
        w("None. The sealed filters are disjoint by construction -- a paper")
        w("retracted within three years cannot also be one retracted after four")
        w("-- and the control arms drew from the same journal and year without")
        w("colliding. So resampling the two populations in the same replicate")
        w("is resampling two independent things, and the difference of the two")
        w("draws is a proper joint statistic rather than a coupled one.")
    else:
        w("")
        w("NOT disjoint. The interval below treats the two draws as independent")
        w("and is therefore too NARROW by whatever the shared mass contributes.")
        w("Read it as an upper bound on the evidence, never as the evidence.")

    w("")
    w("## The contrast")
    w("")
    diffs, fasts, slows = joint(fast_cells, slow_cells, rng)
    lo, hi = pct(diffs, 0.025), pct(diffs, 0.975)
    w("RR_fast - RR_slow = %.4f, 95%% interval %.4f to %.4f (%d usable draws)."
      % (rr_fast - rr_slow, lo, hi, len(diffs)))
    w("")
    w("- fast arm across the same draws: %.4f to %.4f"
      % (pct(fasts, 0.025), pct(fasts, 0.975)))
    w("- slow arm across the same draws: %.4f to %.4f"
      % (pct(slows, 0.025), pct(slows, 0.975)))
    w("")
    excludes = (lo < 0 and hi < 0) or (lo > 0 and hi > 0)
    w("**The interval %s 0.**" % ("EXCLUDES" if excludes else "INCLUDES"))
    w("")
    if excludes:
        w("Under the instrument the three of them asked for, the two arms")
        w("differ. **The published verdict does not change.** This statistic was")
        w("chosen after the data and after the objection, and a statistic chosen")
        w("that way cannot carry a verdict it was picked to reach; the sealed")
        w("secondary stays UNDECIDED, and what this line buys is a pre-registered")
        w("rule for the next study rather than a result in this one.")
    else:
        w("So the sealed secondary verdict is confirmed as UNDECIDED, and this")
        w("time by an instrument that tests the arms instead of a gate.")

    w("")
    w("## The detector's own control")
    w("")
    w("A difference of two Mantel-Haenszel ratios has more ways to be wrong")
    w("than a ratio does, so before the number above is believed the same")
    w("machinery is run on information-free labels: arms permuted within each")
    w("stratum in BOTH populations, which preserves every stratum's size and")
    w("destroys only the thing under test. The contrast must centre on 0.")
    w("")
    nulls = []
    for _ in range(NULL_DRAWS):
        rf, _, _, _ = fast.rr_mh(shuffle_cells(fast_cells, rng))
        rs, _, _, _ = slow.rr_mh(shuffle_cells(slow_cells, rng))
        if rf == rf and rs == rs:
            nulls.append(rf - rs)
    nulls.sort()
    n_lo, n_hi = pct(nulls, 0.025), pct(nulls, 0.975)
    n_mean = sum(nulls) / len(nulls) if nulls else float("nan")
    w("Contrast on information-free labels, %d draws: mean %.4f, 95%% interval "
      "%.4f to %.4f." % (len(nulls), n_mean, n_lo, n_hi))
    w("")
    clean = n_lo <= 0.0 <= n_hi
    w("**The null %s 0**, so the difference statistic %s."
      % ("brackets" if clean else "does NOT bracket",
         "is unbiased at the null and the interval above may be read"
         if clean else "is BROKEN and nothing above may be read"))
    if not clean:
        w("")
        w("Stop here. A contrast that does not centre on 0 when both arms are")
        w("noise is measuring something other than the arms, and no verdict --")
        w("including a confirmation -- may rest on it.")
    w("")
    w("**And one thing the null is NOT, which matters more than the fact that")
    w("it passed.** Its *width* is not a reference distribution for the")
    w("contrast, and it is about four times the bootstrap interval. That is not")
    w("a contradiction and it is not evidence against the finding: permuting")
    w("labels destroys the effect in BOTH populations, so every draw is the")
    w("difference of two ratios sitting near 1 with noisy denominators, and a")
    w("ratio is convex in its denominator. The null of this secondary claim is")
    w("not \"neither retraction does anything\" -- it is \"both do the SAME")
    w("thing\", and shuffling arms cannot construct that. So this control")
    w("licenses the statistic as unbiased and nothing more; anyone reading the")
    w("+/-0.35 as a significance band is reading a null nobody asked about.")
    w("Building the right one -- resampling under an imposed common RR -- is")
    w("the work this study did not do, and it is what a pre-registered version")
    w("of this contrast would have to specify before it could fire.")

    text = "\n".join(lines) + "\n"
    with open(os.path.join(HERE, "contrast.md"), "w") as fh:
        fh.write(text)
    print(text)


def _fmt(value):
    if isinstance(value, tuple):
        return ", ".join("%.4f" % v if isinstance(v, float) else str(v)
                         for v in value)
    return "%.4f" % value if isinstance(value, float) else str(value)


if __name__ == "__main__":
    sys.exit(main())
