#!/usr/bin/env python3
"""
14-model.py: the cumulative-case model for truth.logos.so (research area 14, integration).

Bayesian confirmation in log-odds form over six rival hypotheses (N, D, O, J, I, C), with
  - per-line likelihood-ratio RANGES (low / central / high) read from an input file,
  - a likelihood PROFILE per line: which hypotheses the evidence favours (full weight), which get
    part of the effect ("partial", an exponent 0..1 on the LR) and which get none (baseline),
  - a MODEL-ERROR cap: with probability eps the analysis behind a line is wrong and the evidence is
    as expected under the hypothesis it seemed to count against; this bounds any one line's LR at
    about 1/eps (and below at about eps), so no single "heroic number" can decide the case,
  - a DEPENDENCE discount: lines in the same cluster are not treated as independent; the strongest
    line in a cluster counts in full and every other line in it counts with weight w (0 = only the
    strongest counts, 1 = full independence). Clusters are assumed conditionally independent of
    each other given each hypothesis (stated, and testable by the reader),
  - counter-evidence as LRs < 1 on the hypotheses it counts against,
  - hierarchical priors (P(N), then shares within "a Mind exists"), in presets skeptic / neutral /
    generous, plus the reader's own,
and it prints: posteriors across presets and LR scenarios, a ledger of each line's weight in bans
(log10 units), break-even priors ("what would you need to believe for naturalism to stay ahead?"),
robustness checks (drop the strongest line / cluster / tier, double the counter-evidence), a
dependence x model-error sweep, and a Monte Carlo over the LR ranges.

Usage:
  python3 14-model.py [inputs.json|inputs.yaml] [--mc 20000] [--seed 1] [--json out.json]
Default input: 14-inputs.PLACEHOLDER.json next to this script. Standard library only (PyYAML optional).
"""
import json
import math
import os
import random
import sys

HERE = os.path.dirname(os.path.abspath(__file__))
SCEN = ("low", "central", "high")


# ----------------------------------------------------------------------------------------------
# input
# ----------------------------------------------------------------------------------------------
def load(path):
    with open(path) as f:
        if path.endswith((".yaml", ".yml")):
            import yaml  # optional dependency
            return yaml.safe_load(f)
        return json.load(f)


def flat_prior(preset, hyps):
    """Hierarchical preset -> flat prior over hypotheses."""
    pN = preset["N"]
    share = preset["within_mind"]
    tot = sum(share.values())
    pri = {"N": pN}
    for h in hyps:
        if h != "N":
            pri[h] = (1 - pN) * share.get(h, 0) / tot
    return pri


# ----------------------------------------------------------------------------------------------
# core arithmetic
# ----------------------------------------------------------------------------------------------
def capped_lr(lr, eps):
    """Model-error mixture. With prob. eps the line's analysis is wrong and E is as probable under the
    disfavoured hypothesis as under the favoured one. LR>=1 -> LR/((1-eps)+eps*LR) (< 1/eps);
    LR<1 -> (1-eps)*LR + eps (> eps)."""
    if eps <= 0:
        return lr
    if lr >= 1:
        return lr / ((1 - eps) + eps * lr)
    return (1 - eps) * lr + eps


def item_logs(item, hyps, lr, eps):
    """log10 likelihood multiplier for each hypothesis (relative to the baseline hypotheses)."""
    L = math.log10(capped_lr(lr, eps))
    out = {}
    for h in hyps:
        if h in item["favours"]:
            w = 1.0
        else:
            w = item.get("partial", {}).get(h, 0.0)
        out[h] = w * L
    return out


def combine(items, hyps, lrs, eps, w_dep):
    """Sum log-likelihoods with the within-cluster dependence discount.
    lrs: dict item_id -> LR value. Returns dict h -> total log10 multiplier, and per-item effective logs."""
    clusters = {}
    for it in items:
        clusters.setdefault(it["cluster"], []).append(it)
    total = {h: 0.0 for h in hyps}
    eff = {it["id"]: {h: 0.0 for h in hyps} for it in items}
    for cl, members in clusters.items():
        # One ranking per cluster, used for every hypothesis (so the ledger stays consistent):
        # pro lines (LR >= 1) and counter lines (LR < 1) are ranked separately by |log10 LR|
        # (ties: input order); the first of each counts in full, the rest with weight w_dep.
        order = {it["id"]: n for n, it in enumerate(members)}
        factor = {}
        for sign in (1, -1):
            grp = [it for it in members if (lrs[it["id"]] >= 1) == (sign == 1)]
            grp.sort(key=lambda it: (-abs(math.log10(lrs[it["id"]])), order[it["id"]]))
            for rank, it in enumerate(grp):
                factor[it["id"]] = 1.0 if rank == 0 else w_dep
        for it in members:
            raw = item_logs(it, hyps, lrs[it["id"]], eps)
            for h in hyps:
                eff[it["id"]][h] = factor[it["id"]] * raw[h]
                total[h] += eff[it["id"]][h]
    return total, eff


def posterior(prior, total):
    logs = {h: math.log10(prior[h]) + total[h] for h in prior if prior[h] > 0}
    m = max(logs.values())
    z = sum(10 ** (v - m) for v in logs.values())
    post = {h: 10 ** (logs[h] - m) / z for h in logs}
    for h in prior:
        post.setdefault(h, 0.0)
    return post


def run(data, preset_name, scen="central", w_dep=None, eps=None, drop=(), counter_mult=1.0,
        lr_override=None, prior_override=None, coarse=False):
    """scen: 'low' = every line at the end of its range least favourable to the biblical hypotheses
    (pro and counter lines at 'low', rival-tradition lines at 'high'); 'high' = the reverse.
    coarse=True: the pro lines are grouped one cluster per tier (worst-case dependence grouping for the
    case); counter and rival lines keep their declared clusters, so this is a deliberate stress test."""
    hyps = list(data["hypotheses"].keys())
    preset = data["prior_presets"][preset_name]
    prior = prior_override or flat_prior(preset, hyps)
    w = preset["dependence_weight"] if w_dep is None else w_dep
    e = preset["model_error"] if eps is None else eps
    items = [it for it in data["items"] if it["id"] not in drop and it["cluster"] not in drop
             and it.get("tier") not in drop]
    lrs = {}
    for it in items:
        sc = scen
        if it.get("kind") == "rival" and scen in ("low", "high"):
            sc = "high" if scen == "low" else "low"
        lr = (lr_override or {}).get(it["id"], it["lr"][sc])
        if it.get("kind") == "counter" and counter_mult != 1.0:
            lr = lr ** counter_mult
        lrs[it["id"]] = lr
    if coarse:
        items = [dict(it, cluster="tier-" + str(it.get("tier"))) if it.get("kind") == "pro" else it
                 for it in items]
    total, eff = combine(items, hyps, lrs, e, w)
    return posterior(prior, total), total, eff, prior


def group_p(post, members):
    return sum(post[h] for h in members)


def logodds(p):
    p = min(max(p, 1e-300), 1 - 1e-16)
    return math.log10(p / (1 - p))


def fmt_p(p):
    if p >= 0.9995:
        return "%.4f" % p if p < 0.99995 else ">0.9999"
    if p < 1e-4:
        return "%.1e" % p
    return "%.3f" % p


def break_even_prior(data, preset_name, scen, target, members):
    """Smallest prior on the `members` group (other hypotheses keeping their relative prior weights)
    at which the group's posterior reaches `target`. Found by bisection on log prior."""
    hyps = list(data["hypotheses"].keys())
    base = flat_prior(data["prior_presets"][preset_name], hyps)
    inside = sum(base[h] for h in members)
    outside = 1 - inside

    def post_at(q):
        pri = {}
        for h in hyps:
            pri[h] = (base[h] / inside * q) if h in members else (base[h] / outside * (1 - q))
        post, *_ = run(data, preset_name, scen, prior_override=pri)
        return group_p(post, members)

    lo, hi = -300.0, math.log10(0.999999)
    if post_at(10 ** lo) >= target:
        return 10 ** lo
    if post_at(10 ** hi) < target:
        return None
    for _ in range(200):
        mid = (lo + hi) / 2
        if post_at(10 ** mid) >= target:
            hi = mid
        else:
            lo = mid
    return 10 ** hi


def one_in(p):
    if p is None:
        return "not reachable"
    if p >= 0.5:
        return "%.3f" % p
    return "1 in %s" % (f"{1/p:,.0f}" if 1 / p < 1e7 else "%.1e" % (1 / p))


# ----------------------------------------------------------------------------------------------
# report
# ----------------------------------------------------------------------------------------------
def main(argv):
    path = os.path.join(HERE, "14-inputs.PLACEHOLDER.json")
    mc_n, seed, out_json = 20000, 1, None
    args = list(argv)
    while args:
        a = args.pop(0)
        if a == "--mc":
            mc_n = int(args.pop(0))
        elif a == "--seed":
            seed = int(args.pop(0))
        elif a == "--json":
            out_json = args.pop(0)
        else:
            path = a
    data = load(path)
    hyps = list(data["hypotheses"].keys())
    G = data["groups"]
    presets = list(data["prior_presets"].keys())
    results = {}

    if any(it.get("placeholder") for it in data["items"]) or "_PLACEHOLDER" in data:
        print("=" * 100)
        print("PLACEHOLDER INPUTS: the LR numbers are dummies, not estimates. Read the structure, not the numbers.")
        print("=" * 100)
    print("input:", os.path.relpath(path, HERE))
    print()

    print("## 1. Lines of evidence (likelihood profile)")
    print("%-34s %-6s %-18s %-14s %-22s %s" % ("id", "tier", "cluster", "LR low/c/high", "favours", "partial"))
    for it in data["items"]:
        r = it["lr"]
        print("%-34s %-6s %-18s %-14s %-22s %s" % (
            it["id"], it["tier"], it["cluster"], "%g/%g/%g" % (r["low"], r["central"], r["high"]),
            ",".join(it["favours"]), ",".join("%s:%g" % kv for kv in it.get("partial", {}).items()) or "-"))
    print()

    print("## 2. Posteriors by prior preset and LR scenario (low = every line at the end of its range least favourable to the biblical hypotheses, high = most favourable; preset's own w and eps)")
    hdr = "%-9s %-8s %5s %6s " % ("preset", "LRs", "w", "eps") + " ".join("%8s" % h for h in hyps) + "   P(Mind)  P(C)"
    print(hdr)
    for pn in presets:
        pr = data["prior_presets"][pn]
        pri = flat_prior(pr, hyps)
        print("%-9s %-8s %5s %6s " % (pn, "PRIOR", "", "") + " ".join("%8s" % fmt_p(pri[h]) for h in hyps)
              + "   %7s %6s" % (fmt_p(group_p(pri, G["Mind"])), fmt_p(pri["C"])))
        for sc in SCEN:
            post, *_ = run(data, pn, sc)
            results.setdefault(pn, {})[sc] = post
            print("%-9s %-8s %5g %6g " % (pn, sc, pr["dependence_weight"], pr["model_error"])
                  + " ".join("%8s" % fmt_p(post[h]) for h in hyps)
                  + "   %7s %6s" % (fmt_p(group_p(post, G["Mind"])), fmt_p(post["C"])))
    print()

    print("## 3. Ledger (neutral preset, central LRs): each line's effective weight in bans (log10), after caps and discounts")
    post, total, eff, pri = run(data, "neutral", "central")
    pairs = [("C", "N"), ("C", "J"), ("C", "I"), ("C", "D"), ("J", "N")]
    print("%-34s " % "line" + " ".join("%8s" % ("%s:%s" % p) for p in pairs))
    for it in data["items"]:
        e = eff[it["id"]]
        print("%-34s " % it["id"] + " ".join("%+8.3f" % (e[a] - e[b]) for a, b in pairs))
    print("%-34s " % "TOTAL (log10 Bayes factor)" + " ".join("%+8.3f" % (total[a] - total[b]) for a, b in pairs))
    print("(1 ban = a factor of 10; 0.1 ban = 1 deciban, a factor of about 1.26. Positive = favours the first hypothesis.)")
    print()

    print("## 4. Break-even priors: what you would need to believe beforehand")
    print("Smallest prior at which the evidence brings the group to 50% (other hypotheses keep their preset proportions).")
    print("%-9s %-8s %-34s %-34s %-28s" % ("preset", "LRs", "P(Mind) needed for Mind>50%", "P(C) needed for C>50%", "P(C) needed for C>95%"))
    for pn in presets:
        for sc in SCEN:
            a = break_even_prior(data, pn, sc, 0.5, G["Mind"])
            b = break_even_prior(data, pn, sc, 0.5, ["C"])
            c = break_even_prior(data, pn, sc, 0.95, ["C"])
            print("%-9s %-8s %-34s %-34s %-28s" % (pn, sc, one_in(a), one_in(b), one_in(c)))
    print("Read: naturalism stays ahead only if your prior on 'a Mind' is below the first column.")
    print()

    print("## 5. Robustness (neutral and skeptic presets, central LRs)")
    _, _, eff_n, _ = run(data, "neutral", "central")
    strongest = max(data["items"], key=lambda it: eff_n[it["id"]]["C"] - eff_n[it["id"]]["N"]
                    if it.get("kind") != "counter" else -1e9)
    clus_w = {}
    for it in data["items"]:
        if it.get("kind") != "counter":
            clus_w[it["cluster"]] = clus_w.get(it["cluster"], 0) + eff_n[it["id"]]["C"] - eff_n[it["id"]]["N"]
    strongest_cl = max(clus_w, key=clus_w.get)
    checks = [
        ("baseline", {}),
        ("drop strongest line: " + strongest["id"], {"drop": (strongest["id"],)}),
        ("drop strongest cluster: " + strongest_cl, {"drop": (strongest_cl,)}),
        ("drop all of Tier I", {"drop": ("I",)}),
        ("drop all of Tier III", {"drop": ("III",)}),
        ("counter-evidence doubled (LR^2)", {"counter_mult": 2.0}),
        ("full independence (w=1)", {"w_dep": 1.0}),
        ("only strongest per cluster (w=0)", {"w_dep": 0.0}),
        ("pro lines one cluster per tier, preset w", {"coarse": True}),
        ("worst case: pro lines per tier and w=0", {"coarse": True, "w_dep": 0.0}),
    ]
    print("%-48s %10s %8s %10s %8s" % ("check", "neu P(Mind)", "neu P(C)", "skp P(Mind)", "skp P(C)"))
    for name, kw in checks:
        pn_ = run(data, "neutral", "central", **kw)[0]
        ps_ = run(data, "skeptic", "central", **kw)[0]
        print("%-48s %10s %8s %10s %8s" % (name, fmt_p(group_p(pn_, G["Mind"])), fmt_p(pn_["C"]),
                                           fmt_p(group_p(ps_, G["Mind"])), fmt_p(ps_["C"])))
    print()

    print("## 6. Dependence weight w x model error eps (neutral prior, central LRs): P(C) [P(Mind)]")
    ws = (0.0, 0.25, 0.5, 0.75, 1.0)
    es = (0.1, 0.01, 0.001)
    print("%-8s" % "eps\\w" + "".join("%18s" % w for w in ws))
    for e in es:
        row = "%-8g" % e
        for w in ws:
            p = run(data, "neutral", "central", w_dep=w, eps=e)[0]
            row += "%18s" % ("%s [%s]" % (fmt_p(p["C"]), fmt_p(group_p(p, G["Mind"]))))
        print(row)
    print()

    if mc_n > 0:
        print("## 7. Monte Carlo over the LR ranges (%d draws; log10 LR ~ triangular(low, central, high))" % mc_n)
        rng = random.Random(seed)
        print("%-9s %-24s %-24s" % ("preset", "P(Mind) 5% / 50% / 95%", "P(C) 5% / 50% / 95%"))
        for pn in presets:
            pm, pc = [], []
            for _ in range(mc_n):
                ov = {}
                for it in data["items"]:
                    lo, c, hi = (math.log10(it["lr"][s]) for s in SCEN)
                    ov[it["id"]] = 10 ** (rng.triangular(lo, hi, c) if hi > lo else c)
                p = run(data, pn, "central", lr_override=ov)[0]
                pm.append(group_p(p, G["Mind"]))
                pc.append(p["C"])
            pm.sort()
            pc.sort()
            q = lambda xs, f: xs[min(len(xs) - 1, int(f * len(xs)))]
            print("%-9s %-24s %-24s" % (pn, " / ".join(fmt_p(q(pm, f)) for f in (.05, .5, .95)),
                                         " / ".join(fmt_p(q(pc, f)) for f in (.05, .5, .95))))
        print("(The spread here reflects only the stated LR ranges; it does not cover errors in the structure.)")
        print()

    if out_json:
        with open(out_json, "w") as f:
            json.dump(results, f, indent=1)


if __name__ == "__main__":
    main(sys.argv[1:])