"""Analyses added in the correction, after the results were known and in answer to review.

The registered analysis (code/analyze.py) counts an outcome in an arm as the sum, over the
outcome's terms, of participants affected, so a participant listed under two of its terms is
counted twice. The posted tables don't say which participants share terms, but they bound how
many distinct participants have the outcome in each adverse-event group:

- at least the largest of its terms' counts, and at least the group's participants with any
  serious event less the sum of its other serious terms' counts;
- at most the sum of its terms' counts, and at most the group's participants with any serious
  event.

Groups are disjoint, so an arm's bounds are the sums of its groups'. This script repeats the
registered estimates at every count the bounds allow, and again at the narrower bounds the
trials' public reports leave (data/resolved.json), and adds the other analyses deviations.json
lists. It imports code/analyze.py, so the comparisons, selections, codings and estimators are
the registered ones, and it writes results/R6.json and results/ambiguous.csv.

Usage: python3 code/correct.py <bundle root>
"""

import csv
import itertools
import json
import math
import os
import re

import numpy as np
from numpy.polynomial.hermite import hermgauss
from scipy import optimize, special, stats

import analyze as A

LIMIT = 1_000_000  # configurations enumerated for one analysis at most
CHI2_95 = float(stats.chi2.ppf(0.95, 1))


def load(path):
    with open(os.path.join(A.ROOT, path), encoding="utf-8") as f:
        return json.load(f)


# ---------------------------------------------------------------------------
# Bounds on distinct participants
# ---------------------------------------------------------------------------


def group_bounds(rec, gid, outcome, only=None):
    """Bounds on the participants in one adverse-event group with an outcome (or with one of its
    preferred terms, `only`): (least, most, the outcome's terms the group lists with their
    counts). A serious outcome's bounds also use the group's total of participants with any
    serious event, where the tables are consistent with it; a non-serious table lists only the
    terms above its threshold, so its total bounds nothing."""
    ae = rec["resultsSection"]["adverseEventsModule"]
    group = next(g for g in ae["eventGroups"] if g["id"] == gid)
    serious = A.OUTCOMES[outcome]["table"] == "serious"
    rows = ae.get("seriousEvents" if serious else "otherEvents", [])
    terms = A.outcome_terms(outcome)
    mine, others = [], []
    for e in rows:
        term = A.norm(e.get("term"))
        n = sum(s.get("numAffected") or 0 for s in e.get("stats", []) if s.get("groupId") == gid)
        if n:
            ours = term in terms and (only is None or A.CODING[term]["canonical"] == only)
            (mine if ours else others).append((e.get("term"), n))
    least, most = max((n for _, n in mine), default=0), sum(n for _, n in mine)
    total = group.get("seriousNumAffected")
    if serious and total is not None and least <= total <= most + sum(n for _, n in others):
        least, most = max(least, total - sum(n for _, n in others)), min(most, total)
    at_risk = A.group_at_risk(group, rows, "seriousNumAtRisk" if serious else "otherNumAtRisk")
    return least, min(most, at_risk), mine


def comparison_groups(s):
    nct, k = s["id"].split("#")
    comp = A.SCREENING[nct]["comparisons"][int(k) - 1]
    return comp["antidepressant_groups"], comp["placebo_groups"]


RECORDS = {}


def record(nct):
    if nct not in RECORDS:
        RECORDS[nct] = A.record(nct)
    return RECORDS[nct]


AMBIGUOUS = {}  # (comparison, group, outcome) -> what results/ambiguous.csv lists

# What the trials' public reports settle about groups the tables leave open (data/resolved.json):
# (comparison, group, outcome) -> (least, most, how).
RESOLVED = {(g["comparison"], g["group"], outcome): (bounds[0], bounds[1], g["basis"])
            for g in load("data/resolved.json")["groups"] for outcome, bounds in g["outcomes"].items()}


def bounded_rows(strata, outcome, only=None, resolved=False):
    """For each comparison an analysis covers: (stratum, (least, most) on antidepressants, at
    risk, (least, most) on placebo, at risk). With `resolved`, a group the trials' reports settle
    takes their bounds, which must lie within the tables'."""
    rows = []
    for s, a, n1, c, n0 in A.tables(strata, outcome, only):
        drug, placebo = comparison_groups(s)
        b = []
        for arm, groups in (("antidepressant", drug), ("placebo", placebo)):
            parts = []
            for g in groups:
                lo, hi, terms = group_bounds(record(s["nct"]), g, outcome, only)
                if lo != hi and only is None:
                    AMBIGUOUS[(s["id"], g, outcome)] = (arm, lo, hi, terms)
                key = (s["id"], g, outcome)
                if resolved and only is None and key in RESOLVED:
                    rlo, rhi, _ = RESOLVED[key]
                    if not lo <= rlo <= rhi <= hi:
                        raise ValueError(f"{key}: the reports' {rlo} to {rhi} is outside the tables' {lo} to {hi}")
                    lo, hi = rlo, rhi
                parts.append((lo, hi))
            b.append((sum(p[0] for p in parts), sum(p[1] for p in parts)))
        rows.append((s, b[0], n1, b[1], n0))
    return rows


def write_ambiguous():
    """results/ambiguous.csv: every adverse-event group whose distinct participants with an
    outcome the tables leave open, with the terms it lists and its bounds."""
    with open(os.path.join(A.OUT, "ambiguous.csv"), "w", newline="", encoding="utf-8") as f:
        w = csv.writer(f)
        w.writerow(["comparison", "group", "arm", "group_title", "outcome", "terms", "serious_any", "least", "most",
                    "reports_least", "reports_most", "reports_settle_by"])
        for (sid, gid, outcome), (arm, lo, hi, terms) in sorted(AMBIGUOUS.items()):
            rec = record(sid.split("#")[0])
            group = next(g for g in rec["resultsSection"]["adverseEventsModule"]["eventGroups"] if g["id"] == gid)
            rlo, rhi, how = RESOLVED.get((sid, gid, outcome), (lo, hi, ""))
            w.writerow([sid, gid, arm, group.get("title"), outcome, "; ".join(f"{t} {n}" for t, n in terms),
                        group.get("seriousNumAffected"), lo, hi, rlo, rhi, how])


# ---------------------------------------------------------------------------
# The registered estimator over every count the bounds allow
# ---------------------------------------------------------------------------


def mh_parts(a, n1, c, n0):
    """One comparison's terms of the Mantel-Haenszel risk ratio and the Greenland-Robins variance,
    as code/analyze.py's mantel_haenszel adds them."""
    if a + c == 0:
        return 0.0, 0.0, 0.0
    N = n1 + n0
    return a * n0 / N, c * n1 / N, (n1 * n0 * (a + c) - a * c * N) / N ** 2


def mh_result(R, S, V, drug, placebo, informative):
    out = {"drug_affected": drug, "placebo_affected": placebo, "comparisons_with_events": informative}
    if R == 0 or S == 0:
        out.update({"rr": None, "lo": None, "hi": None, "p": None})
        return out
    log_rr = math.log(R / S)
    se = math.sqrt(V / (R * S))
    z = log_rr / se
    out.update({"rr": math.exp(log_rr), "lo": math.exp(log_rr - A.Z95 * se), "hi": math.exp(log_rr + A.Z95 * se),
                "log_rr": log_rr, "se": se, "p": 2 * stats.norm.sf(abs(z))})
    return out


def at_counts(rows, pick):
    """The registered estimate with each comparison's counts chosen by pick((least, most)) for
    each arm, and the counts it used."""
    R = S = V = 0.0
    drug = placebo = informative = 0
    counts = []
    for s, bd, n1, bp, n0 in rows:
        a, c = pick(bd, "drug"), pick(bp, "placebo")
        r, q, v = mh_parts(a, n1, c, n0)
        R += r; S += q; V += v
        drug += a; placebo += c
        informative += a + c > 0
        counts.append((s, a, n1, c, n0))
    return mh_result(R, S, V, drug, placebo, informative), counts


def enumerate_counts(rows, against=()):
    """The registered estimate at every count the bounds allow: each ambiguous arm takes every
    count from its least to its most, independently. Returns the range of the estimate, of its
    interval's ends and of its p-value, how many configurations give p below alpha, and for each
    earlier estimate named in `against`, how many configurations' intervals contain it."""
    fixed = [(s, bd[0], n1, bp[0], n0) for s, bd, n1, bp, n0 in rows if bd[0] == bd[1] and bp[0] == bp[1]]
    open_rows = [(s, bd, n1, bp, n0) for s, bd, n1, bp, n0 in rows if not (bd[0] == bd[1] and bp[0] == bp[1])]
    R0 = S0 = V0 = 0.0
    d0 = p0 = i0 = 0
    for s, a, n1, c, n0 in fixed:
        r, q, v = mh_parts(a, n1, c, n0)
        R0 += r; S0 += q; V0 += v; d0 += a; p0 += c; i0 += a + c > 0
    choices = [[(a, c) for a in range(bd[0], bd[1] + 1) for c in range(bp[0], bp[1] + 1)] for _, bd, _, bp, _ in open_rows]
    size = math.prod(len(ch) for ch in choices)
    if size > LIMIT:
        raise ValueError(f"{size} configurations is more than {LIMIT}")
    parts = [[mh_parts(a, n1, c, n0) + (a, c) for a, c in ch] for ch, (_, _, n1, _, n0) in zip(choices, open_rows)]
    found = {"configurations": size, "significant": 0, "rr_min": math.inf, "rr_max": -math.inf,
             "lo_min": math.inf, "lo_max": -math.inf, "hi_min": math.inf, "hi_max": -math.inf,
             "p_min": math.inf, "p_max": -math.inf, "contains": {name: 0 for name in against}}
    values = {name: EARLIER["estimates"][name]["value"] for name in against}
    for combo in itertools.product(*parts):
        R, S, V = R0, S0, V0
        d, p, inf = d0, p0, i0
        for r, q, v, a, c in combo:
            R += r; S += q; V += v; d += a; p += c; inf += a + c > 0
        m = mh_result(R, S, V, d, p, inf)
        if m["rr"] is None:
            continue
        found["rr_min"] = min(found["rr_min"], m["rr"]); found["rr_max"] = max(found["rr_max"], m["rr"])
        found["lo_min"] = min(found["lo_min"], m["lo"]); found["lo_max"] = max(found["lo_max"], m["lo"])
        found["hi_min"] = min(found["hi_min"], m["hi"]); found["hi_max"] = max(found["hi_max"], m["hi"])
        found["p_min"] = min(found["p_min"], m["p"]); found["p_max"] = max(found["p_max"], m["p"])
        found["significant"] += m["p"] < A.P["alpha"]
        for name, value in values.items():
            found["contains"][name] += m["lo"] <= value <= m["hi"]
    return found


def least(b, arm):
    return b[0]


def most(b, arm):
    return b[1]


def lowest(b, arm):
    """The count that makes the risk ratio lowest: fewest on antidepressants, most on placebo."""
    return b[0] if arm == "drug" else b[1]


def highest(b, arm):
    return b[1] if arm == "drug" else b[0]


def ambiguous_arms(rows):
    out = []
    for s, bd, n1, bp, n0 in rows:
        for arm, b in (("antidepressant", bd), ("placebo", bp)):
            if b[0] != b[1]:
                out.append({"comparison": s["id"], "arm": arm, "least": b[0], "most": b[1]})
    return out


def bounded(strata, outcome, only=None, against=(), resolved=False):
    """An analysis at every count its bounds allow."""
    rows = bounded_rows(strata, outcome, only, resolved)
    out = {
        "drug_least": sum(r[1][0] for r in rows), "drug_most": sum(r[1][1] for r in rows),
        "placebo_least": sum(r[3][0] for r in rows), "placebo_most": sum(r[3][1] for r in rows),
        "drug_at_risk": sum(r[2] for r in rows), "placebo_at_risk": sum(r[4] for r in rows),
        "comparisons": len(rows), "trials": len({r[0]["nct"] for r in rows}), "ambiguous": ambiguous_arms(rows),
    }
    for name, pick in (("lowest", lowest), ("highest", highest), ("shared", least), ("separate", most)):
        out[name] = at_counts(rows, pick)[0]
    out["every"] = enumerate_counts(rows, against)
    out["testable"] = out["drug_least"] + out["placebo_least"] >= A.P["min_events"]
    return out


# ---------------------------------------------------------------------------
# The earlier analyses' trials and estimates
# ---------------------------------------------------------------------------

EARLIER = load("data/earlier.json")
REANALYSED = {t["nct"] for t in EARLIER["datasets"]["khan_2018"]["trials"]}


def in_reanalyses(s):
    """Whether the re-analyses' data hold the comparison's trial, by the FDA reviews' study lists."""
    return s["nct"] in REANALYSED


def possibly_in_fda_2006(s):
    """Whether the comparison could be among the FDA's 2006 adult trials, by drug and date."""
    d = EARLIER["datasets"]["stone_2009"]
    return bool(s["completion"]) and s["completion"] <= d["completed_by"] and bool(set(s["molecules"]) & set(d["molecules"]))


AGAINST = {
    "P1": ("stone_2009_behaviour_25_64", "stone_2009_behaviour_25_and_over", "hengartner_2019_attempts", "kaminski_2020_attempts"),
    "P2": ("stone_2009_suicidality_25_64", "kaminski_2020_attempts"),
    "A19": ("hammad_2006_pediatric_suicidality",),
}


def independence(strata, analysis):
    """An adult co-primary analysis split by whether the earlier analyses' data hold its trials."""
    sel = A.select(strata, analysis["select"])
    held = [s for s in sel if in_reanalyses(s)]
    maybe = [s for s in sel if possibly_in_fda_2006(s) and not in_reanalyses(s)]
    rest = [s for s in sel if not in_reanalyses(s) and not possibly_in_fda_2006(s)]
    against = AGAINST.get(analysis["id"], ())
    part = lambda ss: {"participants": sum(s["drug"].at_risk + s["placebo"].at_risk for s in ss),
                       "ids": sorted(s["id"] for s in ss), **bounded(ss, analysis["outcome"], against=against),
                       "resolved": bounded(ss, analysis["outcome"], against=against, resolved=True)}
    return {"in_reanalyses": part(held), "possibly_fda_2006": part(maybe), "neither": part(rest)}


# ---------------------------------------------------------------------------
# Random effects without a continuity correction
# ---------------------------------------------------------------------------

NODES, WEIGHTS = hermgauss(64)
Z_NODES = math.sqrt(2) * NODES
LOG_W = np.log(WEIGHTS / math.sqrt(math.pi))


def conditional_loglik(mu, tau, a, m, offset):
    """Log-likelihood of the random-effects Poisson model conditional on each comparison's events:
    a_i given m_i is binomial with log odds log(n1_i / n0_i) + theta_i, theta_i ~ N(mu, tau^2),
    integrated over theta_i by 64-point Gauss-Hermite quadrature. Constants are left out."""
    eta = offset[:, None] + mu + tau * Z_NODES[None, :]
    ll = a[:, None] * eta - m[:, None] * np.logaddexp(0.0, eta)
    return float(np.sum(special.logsumexp(ll + LOG_W[None, :], axis=1)))


def tau_at(mu, a, m, offset, tau_max):
    """The tau in [0, tau_max] that maximizes the likelihood at mu, and that likelihood."""
    res = optimize.minimize_scalar(lambda t: -conditional_loglik(mu, t, a, m, offset), bounds=(0.0, tau_max),
                                   method="bounded", options={"xatol": 1e-9})
    at_zero = conditional_loglik(mu, 0.0, a, m, offset)
    return (0.0, at_zero) if at_zero >= -float(res.fun) else (float(res.x), -float(res.fun))


def poisson_normal(counts, tau_max=3.0):
    """Maximum-likelihood fit of that model over the comparisons with events: the summary risk
    ratio exp(mu), its profile-likelihood 95% interval and likelihood-ratio p, and tau; and the
    same model with tau fixed at 0 (a common risk ratio, no continuity correction). Both are
    found by Brent's method in one dimension at a time, so the fit is the same on any machine
    to the digits reported."""
    rows = [(a, n1, c, n0) for _, a, n1, c, n0 in counts if a + c > 0]
    if len(rows) < 2:
        return None
    a = np.array([r[0] for r in rows], dtype=float)
    m = np.array([r[0] + r[2] for r in rows], dtype=float)
    offset = np.array([math.log(r[1] / r[3]) for r in rows])
    start = math.log(max(a.sum(), 0.5) / max((m - a).sum(), 0.5)) - math.log(sum(r[1] for r in rows) / sum(r[3] for r in rows))
    models = {
        "random": lambda u: tau_at(u, a, m, offset, tau_max)[1],
        "common": lambda u: conditional_loglik(u, 0.0, a, m, offset),
    }
    out = {"comparisons": len(rows)}
    for name, ll_at in models.items():
        res = optimize.minimize_scalar(lambda u: -ll_at(u), bracket=(start - 1, start + 1), tol=1e-10)
        mu, ll = float(res.x), -float(res.fun)
        tau = tau_at(mu, a, m, offset, tau_max)[0] if name == "random" else 0.0
        gap = lambda u: 2 * (ll - ll_at(u)) - CHI2_95
        lo = optimize.brentq(gap, mu - 6, mu, xtol=1e-10)
        hi = optimize.brentq(gap, mu, mu + 6, xtol=1e-10)
        lr = max(0.0, 2 * (ll - ll_at(0.0)))
        out[name] = {"rr": math.exp(mu), "lo": math.exp(lo), "hi": math.exp(hi), "p": float(stats.chi2.sf(lr, 1)), "tau": round(tau, 3)}
    return out


def inverse_variance_corrected(counts):
    """The inverse-variance fixed-effect risk ratio with 0.5 added to every cell of a comparison
    with a zero cell, as the registered random-effects estimate adds it."""
    y, w = [], []
    for _, a, n1, c, n0 in counts:
        if a + c == 0:
            continue
        if a == 0 or c == 0 or a == n1 or c == n0:
            a, c, n1, n0 = a + 0.5, c + 0.5, n1 + 1, n0 + 1
        y.append(math.log((a / n1) / (c / n0)))
        w.append(1 / (1 / a - 1 / n1 + 1 / c - 1 / n0))
    mu = sum(wi * yi for wi, yi in zip(w, y)) / sum(w)
    se = 1 / math.sqrt(sum(w))
    return {"rr": math.exp(mu), "lo": math.exp(mu - A.Z95 * se), "hi": math.exp(mu + A.Z95 * se),
            "p": 2 * stats.norm.sf(abs(mu / se))}


def estimators(strata, outcome):
    """The registered estimates and those that need no continuity correction, at the registered
    counts and at the counts that give the lowest and highest risk ratio."""
    rows = bounded_rows(strata, outcome)
    registered = A.tables(strata, outcome)
    zero_cell = sum(1 for _, a, _, c, _ in registered if a + c > 0 and (a == 0 or c == 0))
    out = {"comparisons_with_events": sum(1 for _, a, _, c, _ in registered if a + c > 0), "with_a_zero_cell": zero_cell,
           "registered_random_effects": A.dersimonian_laird(registered),
           "inverse_variance_corrected": inverse_variance_corrected(registered)}
    for name, counts in (("registered", registered), ("lowest", at_counts(rows, lowest)[1]), ("highest", at_counts(rows, highest)[1])):
        out[name] = poisson_normal(counts)
    return out


# ---------------------------------------------------------------------------
# Multiplicity at the bounds
# ---------------------------------------------------------------------------


def adjusted(results, ids, method, key):
    """Holm or Benjamini-Hochberg over `ids`, from each analysis's p-value at `key`."""
    ps = [results[i]["every"][key] if results[i]["testable"] else None for i in ids]
    q = A.holm(ps) if method == "holm" else A.benjamini_hochberg(ps)
    return dict(zip(ids, q))


# ---------------------------------------------------------------------------
# Age at the bounds
# ---------------------------------------------------------------------------


def age_kind(rec):
    """Which posted statistic code/analyze.py's mean_age took the trial's age from: "mean",
    "median", or None. It follows mean_age step by step."""
    bl = rec.get("resultsSection", {}).get("baselineCharacteristicsModule", {})
    groups = {g["id"]: A.norm(g.get("title")) for g in bl.get("groups", [])}
    counts = {}
    for d in bl.get("denoms", []):
        for c in d.get("counts", []):
            try:
                counts[c["groupId"]] = float(c["value"])
            except (TypeError, ValueError):
                pass
    for kind in ("MEAN", "MEDIAN"):
        for m in bl.get("measures", []):
            if re.sub(r"[^a-z]", "", m.get("title", "").lower()) != "agecontinuous" or m.get("paramType") != kind:
                continue
            if A.norm(m.get("unitOfMeasure")) not in ("years", "year"):
                continue
            values = {}
            for cl in m.get("classes", []):
                for cat in cl.get("categories", []):
                    for x in cat.get("measurements", []):
                        try:
                            values[x["groupId"]] = float(x["value"])
                        except (TypeError, ValueError, KeyError):
                            pass
            if any(t == "total" and g in values for g, t in groups.items()):
                return kind.lower()
            if sum(counts.get(g, 0) for g in values if groups.get(g) != "total") > 0:
                return kind.lower()
    return None


def age_over_counts(strata, outcome):
    """The registered Peto meta-regression on trial age at every count the bounds allow: the
    range of its slope (as % per year) and of its p-value."""
    rows = [r for r in bounded_rows(strata, outcome) if r[0]["mean_age"] is not None]
    acc = [np.zeros(1) for _ in range(5)]  # sum V, sum V x, sum V x^2, sum (O-E), sum x (O-E)
    for s, bd, n1, bp, n0 in rows:
        x = s["mean_age"]
        opts = []
        for a in range(bd[0], bd[1] + 1):
            for c in range(bp[0], bp[1] + 1):
                t = A.peto_terms(a, n1, c, n0)
                opts.append((0.0, 0.0, 0.0, 0.0, 0.0) if not t else (t[1], t[1] * x, t[1] * x * x, t[0], x * t[0]))
        opts = np.array(opts)
        acc = [(acc[k][:, None] + opts[None, :, k]).ravel() for k in range(5)]
        if acc[0].size > LIMIT:
            raise ValueError("too many configurations")
    sv, svx, svx2, soe, sxoe = acc
    sxx = svx2 - svx ** 2 / sv
    slope = (sxoe - svx / sv * soe) / sxx
    p = 2 * stats.norm.sf(np.abs(slope * np.sqrt(sxx)))
    pct = 100 * (np.exp(slope) - 1)
    with_events = sum(1 for _, bd, _, bp, _ in rows if bd[0] + bp[0] > 0)
    return {"configurations": int(sv.size), "pct_min": float(pct.min()), "pct_max": float(pct.max()),
            "p_min": float(p.min()), "p_max": float(p.max()), "comparisons_with_events": with_events,
            "median_age_trials": sorted({s["nct"] for s, bd, _, bp, _ in rows if bd[1] + bp[1] > 0 and age_kind(record(s["nct"])) == "median"})}


# ---------------------------------------------------------------------------


def main():
    strata = A.build_strata()
    analyses = {a["id"]: a for a in A.ANALYSES}
    serious = [a["id"] for a in A.ANALYSES if a["role"] in ("primary", "secondary") and A.OUTCOMES[a["outcome"]]["table"] == "serious"]
    primaries = [a["id"] for a in A.ANALYSES if a["role"] == "primary"]
    secondaries = [a["id"] for a in A.ANALYSES if a["role"] == "secondary" and a["id"] in serious]
    out = {}
    # The analyses at the tables' bounds, and at the bounds the trials' reports leave.
    for family, resolved in (("bounds", False), ("resolved", True)):
        out[family] = {}
        for i in serious:
            a = analyses[i]
            out[family][i] = {"name": a["name"], **bounded(A.select(strata, a["select"]), a["outcome"],
                                                          against=AGAINST.get(i, ()), resolved=resolved)}
        for key in ("p_max", "p_min"):
            for i, q in adjusted(out[family], primaries, "holm", key).items():
                out[family][i][f"holm_at_{key}"] = q
            for i, q in adjusted(out[family], secondaries, "bh", key).items():
                out[family][i][f"q_at_{key}"] = q
        for name in ("shared", "separate"):
            ps = [out[family][i][name]["p"] for i in primaries]
            for i, q in zip(primaries, A.holm(ps)):
                out[family][i][f"holm_{name}"] = q

    # Each co-primary analysis's most frequent preferred term alone, which a group lists once.
    out["single_term"] = {}
    for i in primaries:
        a = analyses[i]
        sel = A.select(strata, a["select"])
        term = A.dominant_term(sel, a["outcome"])
        out["single_term"][i] = {"term": term, **bounded(sel, a["outcome"], term)}

    out["estimators"] = {i: estimators(A.select(strata, analyses[i]["select"]), analyses[i]["outcome"])
                         for i in primaries + ["A6", "A19"]}
    out["age"] = {i: age_over_counts(A.select(strata, analyses[i]["select"]), analyses[i]["outcome"])
                  for i in ("G1", "G2")}
    # The nausea check: one adverse-event group lists two of its terms.
    out["checks"] = {"K1": bounded(A.select(strata, analyses["K1"]["select"]), analyses["K1"]["outcome"])}
    # Which trials the earlier analyses' data hold or could hold.
    out["independence"] = {i: independence(strata, analyses[i]) for i in primaries}
    pediatric = A.select(strata, analyses["A19"]["select"])
    cutoff = EARLIER["datasets"]["hammad_2006"]["completed_by"]
    out["independence"]["pediatric"] = {"comparisons": len(pediatric),
                                        "completed_after_fda_2004": sum(1 for s in pediatric if s["completion"] and s["completion"] > cutoff),
                                        "first_completion": min(s["completion"] for s in pediatric if s["completion"])}
    # How much of what the tables leave open the reports settle, for the serious outcomes.
    open_pairs = {k for k in AMBIGUOUS if A.OUTCOMES[k[2]]["table"] == "serious"}
    settled = {k for k in open_pairs if k in RESOLVED and RESOLVED[k][0] == RESOLVED[k][1]}
    narrowed = {k for k in open_pairs if k in RESOLVED and RESOLVED[k][0] != RESOLVED[k][1]}
    out["resolution"] = {"groups": len({k[:2] for k in open_pairs}), "groups_settled": len({k[:2] for k in open_pairs if all(
                             j not in open_pairs or j in settled for j in [(k[0], k[1], o) for o in A.OUTCOMES])}),
                         "groups_unresolved": len({k[:2] for k in open_pairs if not any(j in RESOLVED for j in open_pairs if j[:2] == k[:2])}),
                         "outcome_bounds": len(open_pairs), "settled": len(settled), "narrowed": len(narrowed)}
    out["resolution"]["groups_narrowed"] = (out["resolution"]["groups"] - out["resolution"]["groups_settled"]
                                            - out["resolution"]["groups_unresolved"])
    A.write("R6.json", out)
    write_ambiguous()


if __name__ == "__main__":
    main()