"""Score the branches before the seal: P(survived), P(killed), P(inconclusive)
for the rule in rule.py, on synthetic 1973-2025 counts, under assumed truths.

This is also the positive control: a trend of known size is injected and the
table says how often the rule sees it.

Baseline 15 events a year (USGS states about 15-16 major earthquakes a year
over the long-term record). Growth g per decade, centred on 2000. Dispersion
phi is the variance over the mean; phi = 1 is Poisson, larger values draw from
a negative binomial with that variance. No catalog data is used here.

    python branches.py > branches.txt
"""

import numpy as np

from rule import trend, verdict

YEARS = np.arange(1973, 2026)
BASE = 15.0
SIMS = 2000
rng = np.random.default_rng(20260911)


def draw(g, phi):
    mu = BASE * (1 + g) ** ((YEARS - 2000) / 10.0)
    if phi == 1.0:
        return rng.poisson(mu)
    n = mu / (phi - 1.0)
    return rng.negative_binomial(n, n / (n + mu))


print(f"rule: 1973-2025, M>=7, base {BASE:g}/yr, {SIMS} simulated catalogs per cell")
print(f"{'true change/decade':>18} {'phi':>4} {'P(survived)':>12} {'P(killed)':>10} {'P(inconcl.)':>12} {'mean CI width':>14}")
for g in (0.0, 0.05, 0.10, 0.20):
    for phi in (1.0, 1.5, 2.5):
        tally = {"survived": 0, "killed": 0, "inconclusive": 0}
        widths = []
        for _ in range(SIMS):
            t = trend(YEARS, draw(g, phi))
            tally[verdict(t)] += 1
            widths.append(t["hi"] - t["lo"])
        print(
            f"{g:>+17.0%} {phi:>4.1f} {tally['survived']/SIMS:>12.4f} "
            f"{tally['killed']/SIMS:>10.4f} {tally['inconclusive']/SIMS:>12.4f} "
            f"{np.mean(widths):>14.3f}"
        )