igneum/sim/finality_v2.py

#!/usr/bin/env python3
"""Igneum finality rule V2: checkpoint-level simulation with latency, partitions and eclipses.

Rule under test (CLAUDE.md, FINALITY RULE V2, review round 2, 3 October 2026):
  * Vote weight of a key = its blue blocks over a flat trailing 30-day window (DAA time,
    86,400 blocks a day). No damping. Dust threshold DUST blocks: a key below it is not a
    voter and is not counted in any denominator.
  * A checkpoint forms every 30 blocks (blue score 30i). Every voter signs every checkpoint
    it sees while online. Aggregators collect votes; a certificate (a lock) exists once the
    collected signatures reach QUORUM (2/3) of the denominator.
  * ACTIVE denominator: sum over eligible keys of weight x participation, where participation
    = min(1, certified votes of that key over the last PRESENCE (240) checkpoint indices / 240).
    A key whose first block is less than 240 checkpoints old counts as participation 1.
  * TOTAL denominator (the alternative): sum of weight over every eligible key, no factor.
  * Equivocation (one key signing two different checkpoint blocks at one index) zeroes the
    key's weight for the rest of the window once the evidence is seen.

Participation bookkeeping has three readings, selected with --pmode. The difference only shows
when a checkpoint index gets no certificate:
  cert    (the brief, literal) an uncertified index contributes 0 certified votes to EVERY key,
          so a stall decays everyone's participation and the denominator shrinks until the
          present signers reach 2/3 of it. Self-healing, and the partition hazard.
  seen    an uncertified index credits the keys whose votes the node observed on the wire.
          Silent keys decay, present keys do not. Not objective: each node counts its own view.
  frozen  the window is over the last 240 CERTIFIED indices, so a stall freezes participation.
          Fails safe and never recovers liveness within the window.

Model (all assumptions, repeated in results_v2.md):
  * Time step = one 30-s slot. The network mines Poisson(30) blocks per slot, split per key by
    hashrate share (perfect difficulty retarget). Every block is blue (GHOSTDAG abstracted).
  * Weight window = 720 hourly buckets of blocks per key (30 days), rolled every hour.
  * A checkpoint index forms on a side each time that side's blue score passes a multiple of
    30, so a partition side mining a fraction s of the hashrate forms checkpoints at s per slot.
    Checkpoint blocks on different sides are different blocks at the same index.
  * Regions: a one-way delay matrix, DELAY between regions, INTRA inside one, lognormal jitter
    per hop. A vote for checkpoint i issued by a key in region r reaches the aggregator in
    region a at t0 + d(src, r) + d(r, a) (+ a per-key extra delay for an eclipsed key).
    One aggregator per region on the side; the certificate forms at the first one to reach
    quorum; its signer set is every vote that arrived there by max(threshold time, t0 + GRACE).
  * Honest keys follow a two-state uptime chain: availability UPTIME (UPTIME_BIG for keys at or
    above 1% of hashrate, a pool with redundant infrastructure), mean outage OUTAGE slots.
  * A partition splits regions into sides, each with its own view (certificates, participation
    ring, blue score). At the heal the views merge: certificates are unioned, two certificates
    at one index with different blocks count as a CONFLICTING LOCK, participation rings are
    OR-merged by index, and any key that signed on two sides at a common index is stripped.
  * Weights are global (a side does not see the other side's blocks for the partition's
    duration; at most 150 minutes of a 30-day window, ignored).

Standard library plus numpy. Deterministic for a given --seed.

Usage:
  python3 finality_v2.py                       # every scenario, markdown on stdout
  python3 finality_v2.py --scenarios A,E --delay 5
  python3 finality_v2.py --quick               # shortened runs for development
"""

import argparse
import sys
import time

import numpy as np

SLOT_S = 30.0
BLOCKS_PER_CP = 30
SLOTS_PER_HOUR = 120
SLOTS_PER_DAY = 2_880
WINDOW_HOURS = 720
WINDOW_BLOCKS = 86_400 * 30
N_HONEST = 1_000
PARETO_SHAPE = 1.0
GEOGRAPHY = (0.45, 0.35, 0.20)  # honest hashrate per region, model assumption
TWO_THIRDS = 2.0 / 3.0


class P:
    """Parameters. Keyword overrides."""

    def __init__(self, **kw):
        self.delay = 2.0        # one-way inter-region delay, seconds
        self.intra = 0.1        # one-way intra-region delay, seconds
        self.jitter = 0.25      # lognormal sigma per hop
        self.grace = 15.0       # seconds after t0 during which late votes still enter the certificate
        self.uptime = 0.97      # availability of an honest key under 1% of hashrate
        self.uptime_big = 0.995  # availability of a key at or above 1% of hashrate (redundant pool infrastructure)
        self.big_share = 0.01
        self.outage = 20        # mean outage length, slots (10 minutes)
        self.denom = "active"   # "active" or "total"
        self.pmode = "cert"     # "cert", "seen", "frozen"
        self.presence = 240     # checkpoints in the participation window
        self.dust = 100         # blocks
        self.quorum = TWO_THIRDS
        self.floor = 0.0        # hybrid: active denominator never below floor x total weight (0 = off)
        self.daa = "none"       # "none": a partition side mines at its hashrate share; "full": each side retargets to 30 blocks/slot at once
        self.__dict__.update(kw)

    def label(self):
        if self.denom == "total":
            return "total"
        s = "active/" + self.pmode
        if self.floor > 0:
            s += "+floor%.2f" % self.floor
        if self.daa != "none":
            s += "+daa"
        return s


class View:
    """One side's view: participation ring, checkpoint counter, certificates."""

    def __init__(self, sid, regions, n, presence, next_idx):
        self.sid = sid
        self.regions = list(regions)
        self.ring = np.zeros((presence, n), dtype=np.int8)
        self.ring_idx = np.full(presence, -1, dtype=np.int64)
        self.pcount = np.zeros(n, dtype=np.int32)
        self.next_idx = int(next_idx)
        self.blue = 0.0
        self.certs = {}      # idx -> (sid, slot, latency_s)
        self.stalls = []     # (idx, slot)
        self.key_mask = None  # keys whose region is on this side
        self.mine_mask = None

    def clone_ring_from(self, other):
        self.ring[:] = other.ring
        self.ring_idx[:] = other.ring_idx
        self.pcount[:] = other.pcount


class Sim:
    def __init__(self, p, rng, n_regions=3):
        self.p = p
        self.rng = rng
        self.R = int(n_regions)
        self.D = np.full((self.R, self.R), p.delay)
        np.fill_diagonal(self.D, p.intra)
        self.N = 0
        self.hash = np.zeros(0)
        self.region = np.zeros(0, dtype=np.int64)
        self.online = np.zeros(0, dtype=bool)
        self.flaky = np.zeros(0, dtype=bool)
        self.signs = np.zeros(0, dtype=bool)
        self.equiv = np.zeros(0, dtype=bool)
        self.stripped = np.zeros(0, dtype=bool)
        self.extra_delay = np.zeros(0)
        self.first_idx = np.zeros(0, dtype=np.int64)
        self.buckets = np.zeros((WINDOW_HOURS, 0))
        self.weight = np.zeros(0)
        self.slot = 0
        self.views = []
        self.next_sid = 1
        self.records = []    # (slot, idx, sid, latency_s or -1, signed/denom, denom/total)
        self.conflicts = []  # (idx, slot the second certificate existed)
        self.q_on = 1.0 / p.outage
        self.q_off = np.zeros(0)

    # ------------------------------------------------------------- keys
    def add_keys(self, hashes, regions, flaky=True, equiv=False, signs=True):
        hashes = np.asarray(hashes, dtype=float)
        regions = np.asarray(regions, dtype=np.int64)
        n = hashes.size
        self.hash = np.concatenate([self.hash, hashes])
        self.region = np.concatenate([self.region, regions])
        self.online = np.concatenate([self.online, np.ones(n, dtype=bool)])
        self.flaky = np.concatenate([self.flaky, np.full(n, flaky, dtype=bool)])
        self.signs = np.concatenate([self.signs, np.full(n, signs, dtype=bool)])
        self.equiv = np.concatenate([self.equiv, np.full(n, equiv, dtype=bool)])
        self.stripped = np.concatenate([self.stripped, np.zeros(n, dtype=bool)])
        self.extra_delay = np.concatenate([self.extra_delay, np.zeros(n)])
        self.first_idx = np.concatenate([self.first_idx, np.full(n, -1, dtype=np.int64)])
        self.buckets = np.concatenate([self.buckets, np.zeros((WINDOW_HOURS, n))], axis=1)
        self.weight = np.concatenate([self.weight, np.zeros(n)])
        first = self.N
        self.N += n
        self.q_off = np.concatenate([self.q_off, np.zeros(n)])
        self.set_uptime()
        return np.arange(first, self.N)

    def set_uptime(self):
        """Per-key off-switch probability per slot from the size-dependent availability. Call after hashrate changes."""
        tot = self.hash.sum()
        share = self.hash / tot if tot > 0 else np.zeros(self.N)
        u = np.where(share >= self.p.big_share, self.p.uptime_big, self.p.uptime)
        self.q_off = self.q_on * (1.0 - u) / u

    def warm_start(self):
        """Fill the 30-day window as if the mining keys had been mining at their share for 30 days."""
        share = self.hash / self.hash.sum()
        lam = 3_600.0 * share
        for h in range(WINDOW_HOURS):
            self.buckets[h] = self.rng.poisson(lam)
        self.weight = self.buckets.sum(axis=0)
        self.first_idx[self.hash > 0] = -10 ** 9
        self.slot = 0

    def init_views(self, warm):
        v = View(0, range(self.R), self.N, self.p.presence, next_idx=(WINDOW_BLOCKS // BLOCKS_PER_CP if warm else 0))
        self._set_masks(v)
        if warm:
            elig = (self.weight >= self.p.dust) & self.signs
            v.ring[:] = elig.astype(np.int8)[None, :]
            v.ring_idx[:] = np.arange(v.next_idx - self.p.presence, v.next_idx)
            v.pcount[:] = v.ring.sum(axis=0)
        self.views = [v]
        self.next_sid = 1

    def _set_masks(self, v):
        v.key_mask = np.isin(self.region, v.regions)
        v.mine_mask = v.key_mask.copy()

    # ------------------------------------------------------------- partitions
    def split(self, groups):
        base = self.views[0]
        new = []
        for g in groups:
            v = View(self.next_sid, g, self.N, self.p.presence, next_idx=base.next_idx)
            self.next_sid += 1
            v.clone_ring_from(base)
            self._set_masks(v)
            new.append(v)
        self.split_idx = base.next_idx
        self.views = new

    def heal(self):
        views = self.views
        P_ = self.p.presence
        nxt = max(v.next_idx for v in views)
        m = View(self.next_sid, range(self.R), self.N, P_, next_idx=nxt)
        self.next_sid += 1
        self._set_masks(m)
        for i in range(max(0, nxt - P_), nxt):
            row = i % P_
            acc = np.zeros(self.N, dtype=np.int8)
            hit = False
            for v in views:
                if v.ring_idx[row] == i:
                    acc |= v.ring[row]
                    hit = True
            if hit:
                m.ring[row] = acc
                m.ring_idx[row] = i
        m.pcount[:] = m.ring.sum(axis=0)
        # certificates and conflicts
        for v in views:
            for idx, (sid, slot, lat) in v.certs.items():
                if idx in m.certs and m.certs[idx][0] != sid:
                    self.conflicts.append((idx, max(slot, m.certs[idx][1])))
                else:
                    m.certs.setdefault(idx, (sid, slot, lat))
            m.stalls.extend(v.stalls)
        # equivocation evidence: an equivocating key signed every side's checkpoint at every common index
        common = min(v.next_idx for v in views) > self.split_idx
        if common and self.equiv.any():
            self.stripped |= self.equiv
            self.strip_slot = self.slot
        self.views = [m]

    # ------------------------------------------------------------- stepping
    def run(self, n_slots, snap=None):
        """snap = (every_n_slots, fn(sim)) called before the step on matching slots."""
        for _ in range(int(n_slots)):
            if snap is not None and self.slot % snap[0] == 0:
                snap[1](self)
            self.step()

    def step(self):
        p = self.p
        slot = self.slot
        tot = self.hash.sum()
        if p.daa == "full" and len(self.views) > 1:
            blocks = np.zeros(self.N)
            for v in self.views:
                side_tot = self.hash[v.mine_mask].sum()
                if side_tot > 0:
                    blocks[v.mine_mask] = self.rng.poisson(BLOCKS_PER_CP * self.hash[v.mine_mask] / side_tot)
        else:
            blocks = self.rng.poisson(BLOCKS_PER_CP * self.hash / tot) if tot > 0 else np.zeros(self.N)
        hp = (slot // SLOTS_PER_HOUR) % WINDOW_HOURS
        if slot % SLOTS_PER_HOUR == 0:
            self.weight -= self.buckets[hp]
            self.buckets[hp] = 0.0
        self.buckets[hp] += blocks
        self.weight += blocks
        gidx = max(v.next_idx for v in self.views)
        newly = (blocks > 0) & (self.first_idx == -1)
        if newly.any():
            self.first_idx[newly] = gidx
        r = self.rng.random(self.N)
        flip = np.where(self.online, r < self.q_off, r < self.q_on) & self.flaky
        self.online ^= flip
        for v in self.views:
            v.blue += blocks[v.mine_mask].sum()
            while v.blue >= BLOCKS_PER_CP:
                v.blue -= BLOCKS_PER_CP
                self.checkpoint(v, blocks)
        self.slot += 1

    def participation(self, v, idx):
        part = np.minimum(1.0, v.pcount / float(self.p.presence))
        new = (self.first_idx > idx - self.p.presence) & (self.first_idx >= 0)
        part[new] = 1.0
        return part

    def checkpoint(self, v, blocks):
        p = self.p
        idx = v.next_idx
        v.next_idx += 1
        t0 = self.slot * SLOT_S
        # miner of the checkpoint block: a key on this side weighted by its blocks this slot
        bm = blocks * v.mine_mask
        cum = np.cumsum(bm)
        if cum[-1] > 0:
            j = int(np.searchsorted(cum, self.rng.random() * cum[-1], side="right"))
            src = int(self.region[min(j, self.N - 1)])
        else:
            src = v.regions[0]
        jit1 = np.exp(self.rng.normal(0.0, p.jitter, self.R))
        jit2 = np.exp(self.rng.normal(0.0, p.jitter, (self.R, self.R)))
        elig = (self.weight >= p.dust) & ~self.stripped
        voters = elig & self.online & self.signs & (v.key_mask | self.equiv)
        if p.denom == "active":
            denom = float((self.weight * self.participation(v, idx))[elig].sum())
        else:
            denom = float(self.weight[elig].sum())
        total = float(self.weight[elig].sum())
        if p.denom == "active" and p.floor > 0:
            denom = max(denom, p.floor * total)
        need = p.quorum * denom
        vi = np.flatnonzero(voters)
        signed_w = float(self.weight[vi].sum())
        ratio = signed_w / denom if denom > 0 else 0.0
        best = None
        if denom > 0 and vi.size > 0:
            w = self.weight[vi]
            reg = self.region[vi]
            hop1 = np.where(self.equiv[vi], p.intra, self.D[src, reg] * jit1[reg])
            issue = t0 + hop1 + self.extra_delay[vi]
            for a in v.regions:
                hop2 = np.where(self.equiv[vi], p.intra, self.D[reg, a] * jit2[reg, a])
                arr = issue + hop2
                order = np.argsort(arr, kind="stable")
                cw = np.cumsum(w[order])
                k = int(np.searchsorted(cw, need))
                if k < cw.size:
                    thr = float(arr[order[k]])
                    if best is None or thr < best[0]:
                        best = (thr, arr)
        if best is not None:
            thr, arr = best
            lat = thr - t0
            seal = max(thr, t0 + p.grace)
            mask = np.zeros(self.N, dtype=np.int8)
            mask[vi[arr <= seal]] = 1
            v.certs[idx] = (v.sid, self.slot, lat)
            self._push(v, idx, mask)
            self.records.append((self.slot, idx, v.sid, lat, ratio, denom / total if total > 0 else 0.0))
        else:
            v.stalls.append((idx, self.slot))
            if p.pmode == "cert":
                self._push(v, idx, np.zeros(self.N, dtype=np.int8))
            elif p.pmode == "seen":
                self._push(v, idx, voters.astype(np.int8))
            # frozen: no ring update
            self.records.append((self.slot, idx, v.sid, -1.0, ratio, denom / total if total > 0 else 0.0))

    def _push(self, v, idx, mask):
        row = idx % self.p.presence
        v.pcount -= v.ring[row]
        v.ring[row] = mask
        v.pcount += mask
        v.ring_idx[row] = idx

    # ------------------------------------------------------------- queries
    def recs(self):
        return np.array(self.records, dtype=float).reshape(-1, 6)

    def elig_mask(self):
        return (self.weight >= self.p.dust) & ~self.stripped

    def share(self, keys):
        e = self.elig_mask()
        tot = self.weight[e].sum()
        if tot <= 0:
            return 0.0
        m = np.zeros(self.N, dtype=bool)
        m[keys] = True
        return float(self.weight[m & e].sum() / tot)


# ---------------------------------------------------------------- helpers

def pareto_hashrates(rng, n, total=1.0):
    h = rng.pareto(PARETO_SHAPE, n) + 1.0
    return h * (total / h.sum())


def assign_regions(hashes, targets):
    """Largest key first, each to the region with the largest remaining hashrate deficit."""
    targets = np.asarray(targets, dtype=float)
    tot = hashes.sum()
    filled = np.zeros(targets.size)
    reg = np.zeros(hashes.size, dtype=np.int64)
    for i in np.argsort(-hashes):
        r = int(np.argmax(targets * tot - filled))
        reg[i] = r
        filled[r] += hashes[i]
    return reg


def pick_weight_subset(rng, weights, frac):
    """Random keys whose weight sums to about frac of total, never over."""
    target = frac * weights.sum()
    mask = np.zeros(weights.size, dtype=bool)
    acc = 0.0
    for i in rng.permutation(weights.size):
        if acc + weights[i] <= target:
            mask[i] = True
            acc += weights[i]
    return mask, acc / weights.sum()


def gini(x):
    x = np.sort(np.asarray(x, dtype=float))
    n = x.size
    if n == 0 or x.sum() == 0:
        return 0.0
    cum = np.cumsum(x)
    return (n + 1 - 2.0 * cum.sum() / cum[-1]) / n


def pct(x, d=1):
    return ("%%.%df%%%%" % d) % (100.0 * x)


def md_table(headers, rows):
    out = ["| " + " | ".join(str(h) for h in headers) + " |", "|" + "---|" * len(headers)]
    for r in rows:
        out.append("| " + " | ".join(str(c) for c in r) + " |")
    return "\n".join(out)


def lat_stats(recs, s_from=0, s_to=None):
    if s_to is None:
        s_to = np.inf
    sel = recs[(recs[:, 0] >= s_from) & (recs[:, 0] < s_to)]
    lat = sel[:, 3]
    ok = lat >= 0
    n = int(sel.shape[0])
    st = int((~ok).sum())
    cps = np.floor(lat[ok] / SLOT_S)
    d = dict(n=n, stalls=st, c0=int((cps == 0).sum()), c1=int((cps == 1).sum()), c2=int((cps >= 2).sum()))
    if ok.any():
        d.update(med=float(np.median(lat[ok])), p99=float(np.percentile(lat[ok], 99)), mx=float(lat[ok].max()),
                 ratio_min=float(sel[ok, 4].min()), ratio_med=float(np.median(sel[ok, 4])))
    else:
        d.update(med=float("nan"), p99=float("nan"), mx=float("nan"), ratio_min=float("nan"), ratio_med=float("nan"))
    return d


def lat_row(label, d):
    return [label, d["n"], d["c0"], d["c1"], d["c2"], d["stalls"], "%.1f" % d["med"], "%.1f" % d["p99"], "%.1f" % d["mx"],
            "%.3f" % d["ratio_min"]]


LAT_HEADERS = ["run", "checkpoints", "locked in slot 0", "slot 1", "slot 2+", "stalled", "median s", "p99 s", "max s",
               "min signed/denominator"]


def first_lock_after(recs, slot, sid=None):
    sel = recs[(recs[:, 0] >= slot) & (recs[:, 3] >= 0)]
    if sid is not None:
        sel = sel[sel[:, 2] == sid]
    if sel.shape[0] == 0:
        return None
    return int(sel[0, 0])


def stalls_between(recs, s_from, s_to, sid=None):
    sel = recs[(recs[:, 0] >= s_from) & (recs[:, 0] < s_to) & (recs[:, 3] < 0)]
    if sid is not None:
        sel = sel[sel[:, 2] == sid]
    return int(sel.shape[0])


def fmt_min(slots):
    if slots is None:
        return "never"
    m = slots * SLOT_S / 60.0
    if m < 120:
        return "%.0f min" % m
    if m < 48 * 60:
        return "%.1f h" % (m / 60.0)
    return "%.1f d" % (m / 1440.0)


def fmt_days(slots):
    return "never" if slots is None else "%.1f" % (slots / SLOTS_PER_DAY)


def build_honest(p, seed, n_regions=3, geography=GEOGRAPHY, big_pool=None):
    rng = np.random.default_rng(seed)
    sim = Sim(p, rng, n_regions=n_regions)
    if big_pool is None:
        h = pareto_hashrates(rng, N_HONEST)
        reg = assign_regions(h, geography)
        sim.add_keys(h, reg, flaky=True)
        pool = None
    else:
        h = pareto_hashrates(rng, N_HONEST - 1, total=1.0 - big_pool)
        reg = assign_regions(h, geography)
        sim.add_keys(h, reg, flaky=True)
        pool = int(sim.add_keys([big_pool], [3], flaky=False)[0])
    return sim, rng, pool


# ---------------------------------------------------------------- scenarios

def scenario_a(args):
    out = ["### A. Steady state from zero history, 1,000 honest keys, 3 regions %s, %d days" % (
        "/".join(pct(g, 0) for g in GEOGRAPHY), args.days_a), ""]
    days = args.days_a
    snaps = {}
    lat_rows = []
    prop_rows = []
    for denom in ("active", "total"):
        p = P(denom=denom, delay=args.delay)
        sim, rng, _ = build_honest(p, args.seed)
        sim.init_views(warm=False)
        series = []

        def snap(s, series=series):
            e = s.elig_mask()
            series.append((s.slot // SLOTS_PER_DAY, s.weight.sum() / WINDOW_BLOCKS, int((~e).sum())))

        sim.run(days * SLOTS_PER_DAY, snap=(SLOTS_PER_DAY, snap))
        snap(sim)
        recs = sim.recs()
        snaps[denom] = series
        lat_rows.append(lat_row("%s, days 0 to 1" % denom, lat_stats(recs, 0, SLOTS_PER_DAY)))
        lat_rows.append(lat_row("%s, days 1 to 30" % denom, lat_stats(recs, SLOTS_PER_DAY, 30 * SLOTS_PER_DAY)))
        lat_rows.append(lat_row("%s, days 30 to %d" % (denom, days), lat_stats(recs, 30 * SLOTS_PER_DAY, None)))
        w = sim.weight
        ws = w / w.sum()
        hs = sim.hash / sim.hash.sum()
        order = np.argsort(-hs)
        rel = ws / hs
        prop_rows.append([denom, "%.5f" % np.corrcoef(hs, ws)[0, 1], "%.3f / %.3f" % (gini(hs), gini(ws)),
                          pct(hs[order[0]]) + " / " + pct(ws[order[0]]),
                          pct(hs[order[:10]].sum()) + " / " + pct(ws[order[:10]].sum()),
                          pct(hs[order[500:]].sum()) + " / " + pct(ws[order[500:]].sum()),
                          "%.3f / %.3f" % (rel.min(), rel.max()), int((rel < 0.9).sum()), int((~sim.elig_mask()).sum())])
        # genesis stall run
        first = first_lock_after(recs, 0)
        snaps[denom + "_first"] = first
    s = snaps["active"]
    ramp = []
    for d in (1, 2, 5, 10, 20, 30, 31, 45, days):
        row = [x for x in s if x[0] == d]
        if row:
            ramp.append([d, pct(row[0][1]), row[0][2]])
    out.append("Weight ramp (active run; the total run mines the same blocks):")
    out.append("")
    out.append(md_table(["day", "total weight / full window", "keys under dust (100 blocks)"], ramp))
    out.append("")
    out.append("Weight against hashrate at day %d:" % days)
    out.append("")
    out.append(md_table(["denominator", "corr(hash, weight)", "Gini hash / weight", "top-1 hash / weight",
                         "top-10 hash / weight", "bottom-500 hash / weight", "min / max weight:hash", "keys under 0.9x",
                         "dust keys"], prop_rows))
    out.append("")
    out.append("Lock latency (seconds from checkpoint block to quorum; slot = 30 s), inter-region delay %.1f s:" % args.delay)
    out.append("")
    out.append(md_table(LAT_HEADERS, lat_rows))
    out.append("")
    out.append("First lock from genesis: active at slot %s, total at slot %s (the first checkpoints have no key above dust)." % (
        snaps["active_first"], snaps["total_first"]))
    out.append("")
    # delay comparison, short warm-started runs
    rows = []
    for delay in (0.5, 2.0, 5.0):
        for grace in (args.grace,):
            p = P(denom="active", delay=delay, grace=grace)
            sim, rng, _ = build_honest(p, args.seed)
            sim.warm_start()
            sim.init_views(warm=True)
            sim.run(args.days_delay * SLOTS_PER_DAY)
            recs = sim.recs()
            d = lat_stats(recs)
            # participation of the slowest region
            v = sim.views[0]
            part = sim.participation(v, v.next_idx)
            by_region = ["%.3f" % np.average(part[sim.region == r], weights=sim.weight[sim.region == r]) for r in range(3)]
            rows.append(lat_row("delay %.1f s, grace %.0f s, %d days warm" % (delay, grace, args.days_delay), d) + ["/".join(by_region)])
    out.append("Delay sweep, warm-started (steady-state weights), active denominator. Last column: weight-averaged participation per region.")
    out.append("")
    out.append(md_table(LAT_HEADERS + ["participation r0/r1/r2"], rows))
    return "\n".join(out)


def scenario_b(args):
    out = ["### B. Rental burst at day 60, one public key with a x honest hashrate, signs every checkpoint", ""]
    mults = (1, 2, 4, 9)
    series = {}
    cross = {}
    for a in mults:
        p = P(denom="active", delay=args.delay)
        sim, rng, _ = build_honest(p, args.seed, n_regions=4)
        att = int(sim.add_keys([0.0], [3], flaky=False)[0])
        sim.warm_start()
        sim.init_views(warm=True)
        sim.run(SLOTS_PER_HOUR)  # one hour of baseline
        t_event = sim.slot
        sim.hash[att] = float(a)
        s = []
        c13 = c23 = None
        for day in range(1, args.days_b + 1):
            for _ in range(SLOTS_PER_DAY // 24):
                sim.run(24)
                sh = sim.share([att])
                if c13 is None and sh >= 1.0 / 3.0:
                    c13 = sim.slot - t_event
                if c23 is None and sh >= TWO_THIRDS:
                    c23 = sim.slot - t_event
            s.append((day, sim.share([att])))
        recs = sim.recs()
        series[a] = s
        cross[a] = (c13, c23, lat_stats(recs, t_event, None))
    pick = [d for d in (1, 5, 10, 15, 20, 25, 30, 35) if d <= args.days_b]
    rows = []
    for d in pick:
        row = ["+%d" % d]
        for a in mults:
            sim_v = dict(series[a])[d]
            formula = min(d / 30.0, 1.0) * a / (1.0 + a)
            row.append("%s / %s" % (pct(sim_v), pct(formula)))
        rows.append(row)
    out.append("Attacker weight share, simulated / formula (t/30) x a/(1+a):")
    out.append("")
    out.append(md_table(["day after burst"] + ["a=%d (%s of hashrate)" % (a, pct(a / (1.0 + a), 0)) for a in mults], rows))
    out.append("")
    ev = [
        ["crosses 1/3 (honest alone can no longer lock), sim day"] + [fmt_days(cross[a][0]) for a in mults],
        ["crosses 1/3, formula 10(1+a)/a"] + ["%.1f" % (10.0 * (1 + a) / a) for a in mults],
        ["crosses 2/3 (locks alone), sim day"] + [fmt_days(cross[a][1]) for a in mults],
        ["crosses 2/3, formula 20(1+a)/a"] + ["%.1f" % (20.0 * (1 + a) / a) if 20.0 * (1 + a) / a <= 30 else "never (ceiling %s)" % pct(a / (1 + a), 0) for a in mults],
        ["max abs deviation sim vs formula, days 1 to 30, points"] + [
            "%.2f" % (100 * max(abs(v - min(d / 30.0, 1.0) * a / (1.0 + a)) for d, v in series[a] if d <= 30)) for a in mults],
        ["stalled checkpoints after the burst"] + [cross[a][2]["stalls"] for a in mults],
    ]
    out.append(md_table(["event"] + ["a=%d" % a for a in mults], ev))
    return "\n".join(out)


def scenario_c(args):
    out = ["### C. Silent set: a random set holding x of weight stops signing at day 60 (hour 3 of the run) and keeps mining", ""]
    fracs = (0.34, 0.40, 0.45, 0.50, 0.55)
    configs = [("active", "cert", args.hours_c, 0.0, 240), ("active", "seen", args.hours_c, 0.0, 240),
               ("active", "frozen", args.hours_c, 0.0, 240), ("active", "cert", args.hours_c_total, 0.0, 2880),
               ("active", "cert", args.hours_c_total, 0.8, 240), ("active", "cert", args.hours_c_total, 0.85, 240),
               ("total", "cert", args.hours_c_total, 0.0, 240)]
    labels = []
    for denom, pmode, hours, floor, presence in configs:
        lab = P(denom=denom, pmode=pmode, floor=floor, presence=presence).label()
        if presence != 240:
            lab += ", presence %d" % presence
        labels.append(lab)
    rows = []
    detail = []
    for frac in fracs:
        row = [pct(frac, 0)]
        for denom, pmode, hours, floor, presence in configs:
            p = P(denom=denom, pmode=pmode, delay=args.delay, floor=floor, presence=presence)
            sim, rng, _ = build_honest(p, args.seed)
            sim.warm_start()
            sim.init_views(warm=True)
            sim.run(3 * SLOTS_PER_HOUR)
            silent, got = pick_weight_subset(rng, sim.weight, frac)
            sim.signs[silent] = False
            t_event = sim.slot
            v0 = sim.views[0]
            e0 = sim.elig_mask()
            s_on = float(sim.weight[e0 & sim.online & ~silent].sum() / sim.weight[e0].sum())
            p0 = float((sim.weight * sim.participation(v0, v0.next_idx))[e0].sum() / sim.weight[e0].sum())
            predicted = max(0, int(round(p.presence * (p0 - 1.5 * s_on))))
            sim.run(int(hours * SLOTS_PER_HOUR))
            recs = sim.recs()
            fl = first_lock_after(recs, t_event)
            st = stalls_between(recs, t_event, sim.slot)
            v = sim.views[0]
            part = sim.participation(v, v.next_idx)
            sil_part = float(np.average(part[silent], weights=sim.weight[silent]))
            # stalls after the first lock (does it stay locked?)
            later = stalls_between(recs, fl, sim.slot) if fl is not None else st
            fl_txt = "never (%s)" % fmt_min(sim.slot - t_event) if fl is None else fmt_min(fl - t_event)
            row.append("%s, %d stalled" % (fl_txt, st))
            if denom == "active" and pmode == "cert" and floor == 0 and presence == 240:
                tail = recs[(recs[:, 0] >= sim.slot - SLOTS_PER_HOUR) & (recs[:, 3] >= 0)]
                detail.append([pct(frac, 0), pct(got), int(silent.sum()), pct(s_on), "%.3f" % p0, predicted, fl_txt, st, later, "%.3f" % sil_part,
                               "%.3f" % (np.median(tail[:, 4]) if tail.shape[0] else float("nan")),
                               "%.3f" % (np.median(tail[:, 5]) if tail.shape[0] else float("nan"))])
        rows.append(row)
    out.append("Time from the event to the first lock, and checkpoints stalled in the run (%d h for the first three columns, %d h for the rest):" % (
        args.hours_c, args.hours_c_total))
    out.append("")
    out.append(md_table(["silent weight"] + labels, rows))
    out.append("")
    out.append("Active/cert detail. Predicted stalls = presence x (p0 - 1.5 x online signing share), where p0 is the weight-averaged participation at the event:")
    out.append("")
    out.append(md_table(["silent weight", "picked", "keys", "online signing share at event", "p0", "predicted stalls", "first lock", "stalled", "stalled after first lock",
                         "silent participation at end", "signed/denominator at end (median, last hour)",
                         "active/total at end"], detail))
    return "\n".join(out)


def scenario_d(args):
    out = ["### D. Churn: a random set holding x of weight stops mining and signing at day 60 (hour 3 of the run)", ""]
    fracs = (0.35, 0.50)
    rows = []
    dconfigs = [P(denom="active", delay=args.delay), P(denom="active", presence=2880, delay=args.delay),
                P(denom="active", floor=0.8, delay=args.delay), P(denom="active", floor=0.85, delay=args.delay),
                P(denom="total", delay=args.delay)]
    for frac in fracs:
        for p in dconfigs:
            days = args.days_d35 if frac < 0.4 else args.days_d50
            denom = p.label() + (", presence 2880" if p.presence == 2880 else "")
            sim, rng, _ = build_honest(p, args.seed)
            sim.warm_start()
            sim.init_views(warm=True)
            sim.run(3 * SLOTS_PER_HOUR)
            gone, got = pick_weight_subset(rng, sim.weight, frac)
            sim.signs[gone] = False
            sim.online[gone] = False
            sim.flaky[gone] = False
            sim.hash[gone] = 0.0
            t_event = sim.slot
            live = ~gone
            sim.run(int(days * SLOTS_PER_DAY))
            recs = sim.recs()
            fl = first_lock_after(recs, t_event)
            st = stalls_between(recs, t_event, sim.slot)
            later = stalls_between(recs, fl, sim.slot) if fl is not None else 0
            live_share_at = None
            rows.append([pct(frac, 0), pct(got), int(gone.sum()), denom,
                         "never in %s" % fmt_min(sim.slot - t_event) if fl is None else fmt_min(fl - t_event),
                         st, later, pct(sim.share(np.flatnonzero(live)))])
    out.append(md_table(["churn weight", "picked", "keys", "denominator", "first lock after the event", "stalled checkpoints",
                         "stalled after first lock", "live share of total weight at end of run"], rows))
    out.append("")
    out.append("Analytic (perfect retarget): live share of total weight on day t = 1 - x(30 - t)/30, so the total "
               "denominator recovers at t = 30(1 - 1/(3x)): never for x <= 1/3, day 1.4 at 35%, day 10 at 50%.")
    return "\n".join(out)


def run_partition(seed, fracs, att_share, dur_min, p, pre_min=60, post_min=180):
    rng = np.random.default_rng(seed)
    nR = max(3, len(fracs) + 1)
    sim = Sim(p, rng, n_regions=nR)
    h = pareto_hashrates(rng, N_HONEST)
    reg = assign_regions(h, fracs)
    sim.add_keys(h, reg, flaky=True)
    groups = [[i] for i in range(len(fracs))]
    att = None
    if att_share > 0:
        att = int(sim.add_keys([att_share / (1.0 - att_share)], [len(fracs)], flaky=False, equiv=True)[0])
        groups[0].append(len(fracs))
    sim.warm_start()
    sim.init_views(warm=True)
    sim.run(pre_min * 2)
    t_split = sim.slot
    sim.split(groups)
    sides = [v.sid for v in sim.views]
    sim.run(dur_min * 2)
    t_heal = sim.slot
    sim.heal()
    sim.run(post_min * 2)
    recs = sim.recs()
    res = dict(conflicts=len(sim.conflicts), t_split=t_split, t_heal=t_heal)
    res["first_conflict_min"] = (min(c[1] for c in sim.conflicts) - t_split) / 2.0 if sim.conflicts else None
    res["side_first_lock"] = []
    res["side_locks"] = []
    res["side_stalls"] = []
    for sid in sides:
        fl = first_lock_after(recs[recs[:, 0] < t_heal], t_split, sid=sid)
        res["side_first_lock"].append(None if fl is None else (fl - t_split) / 2.0)
        sel = recs[(recs[:, 0] >= t_split) & (recs[:, 0] < t_heal) & (recs[:, 2] == sid)]
        res["side_locks"].append(int((sel[:, 3] >= 0).sum()))
        res["side_stalls"].append(int((sel[:, 3] < 0).sum()))
    res["post_stalls"] = stalls_between(recs, t_heal, sim.slot)
    fl = first_lock_after(recs, t_heal)
    res["post_first_lock_min"] = None if fl is None else (fl - t_heal) / 2.0
    res["att_share"] = sim.share([att]) if att is not None else 0.0
    return res


def fm(m):
    return "never" if m is None else "%.0f" % m


def scenario_e(args):
    out = ["### E. Partition: honest weight split across sides for a set time, then healed", ""]
    configs = [P(denom="active", pmode="cert", delay=args.delay), P(denom="active", pmode="seen", delay=args.delay),
               P(denom="total", delay=args.delay)]
    splits = [("50/50", [0.5, 0.5]), ("33/33/34", [0.33, 0.33, 0.34])]
    rows_conf = []
    rows_first = []
    for name, fr in splits:
        for att in (0.0, 0.34):
            for dur in (30, 90, 150):
                rc = [name, pct(att, 0), dur]
                rf = [name, pct(att, 0), dur]
                for p in configs:
                    r = run_partition(args.seed, fr, att, dur, p)
                    rc.append("%d%s" % (r["conflicts"], "" if r["conflicts"] == 0 else " (first at %s min)" % fm(r["first_conflict_min"])))
                    rf.append(" / ".join(fm(x) for x in r["side_first_lock"]) + " ; post-heal stalls %d" % r["post_stalls"])
                rows_conf.append(rc)
                rows_first.append(rf)
    labels = [p.label() for p in configs]
    out.append("Conflicting locks (two certificates at one index, different blocks). Attacker = equivocating key holding the stated share of total weight, honest weight split as stated. Pass needs 0 at 0% attacker for 30, 90 and 150 min.")
    out.append("")
    out.append(md_table(["honest split", "attacker", "partition min"] + labels, rows_conf))
    out.append("")
    out.append("Minutes after the split until each side's first lock (side order as in the split; attacker mines on the first side) and stalls in the 3 hours after the heal:")
    out.append("")
    out.append(md_table(["honest split", "attacker", "partition min"] + labels, rows_first))
    out.append("")
    # supplementary: more splits, 0% attacker, plus the DAA-retarget and floor variants
    configs2 = configs + [P(denom="active", pmode="cert", daa="full", delay=args.delay),
                          P(denom="active", pmode="seen", daa="full", delay=args.delay),
                          P(denom="active", pmode="cert", floor=0.8, daa="full", delay=args.delay),
                          P(denom="active", pmode="cert", floor=0.85, daa="full", delay=args.delay)]
    labels2 = [p.label() for p in configs2]
    rows = []
    for name, fr in (("50/50", [0.5, 0.5]), ("60/40", [0.6, 0.4]), ("67/33", [0.67, 0.33]), ("80/20", [0.8, 0.2]),
                     ("33/33/34", [0.33, 0.33, 0.34])):
        for dur in (150, 360):
            rc = [name, dur]
            for p in configs2:
                r = run_partition(args.seed, fr, 0.0, dur, p, post_min=120)
                rc.append("%d; %s" % (r["conflicts"], " / ".join(fm(x) for x in r["side_first_lock"])))
            rows.append(rc)
    out.append("Supplementary, 0% attacker, 150 and 360 min. Cell = conflicting locks; minutes to each side's first lock. "
               "'+daa' = each side retargets to 1 block/s at once (worst case for the presence clock); "
               "'+floor0.80' = active denominator never below 80% of total weight (a lock needs at least 53.3% of total), "
               "'+floor0.85' needs 56.7%.")
    out.append("")
    out.append(md_table(["honest split", "partition min"] + labels2, rows))
    out.append("")
    # presence window sweep, active/cert, 0% attacker
    rows = []
    for presence in (240, 720, 2880):
        for name, fr in (("50/50", [0.5, 0.5]), ("60/40", [0.6, 0.4]), ("33/33/34", [0.33, 0.33, 0.34])):
            dur = {240: 360, 720: 720, 2880: 1500}[presence]
            p = P(denom="active", pmode="cert", presence=presence, delay=args.delay)
            r = run_partition(args.seed, fr, 0.0, dur, p, post_min=120)
            s = min(fr)
            pred = presence * (1.0 - 1.5 * s) / s / 2.0 if s < TWO_THIRDS else None
            rows.append([presence, "%.1f h" % (presence / 120.0), name, dur, r["conflicts"],
                         " / ".join(fm(x) for x in r["side_first_lock"]), "%.0f" % pred if pred is not None else "no"])
    out.append("Presence window sweep, active/cert, 0% attacker. Prediction for the smallest side (share s): "
               "first lock after presence x (1 - 1.5 s) / s slots, in minutes = that / 2.")
    out.append("")
    out.append(md_table(["presence (checkpoints)", "presence (hours)", "honest split", "partition min", "conflicts",
                         "first lock per side, min", "predicted smallest-side lock, min"], rows))
    return "\n".join(out)


def run_eclipse(seed, dur_h, poisoned, p, pre_min=60, post_h=5):
    rng = np.random.default_rng(seed)
    sim = Sim(p, rng, n_regions=5)
    # shares of TOTAL weight: pool 20%, attacker 34% when poisoned, the rest honest
    others = 0.8 if not poisoned else 0.46
    h = pareto_hashrates(rng, N_HONEST - 1, total=others)
    reg = assign_regions(h, GEOGRAPHY)
    sim.add_keys(h, reg, flaky=True)
    K = int(sim.add_keys([0.2], [3], flaky=False)[0])
    att = None
    if poisoned:
        att = int(sim.add_keys([0.34], [4], flaky=False, equiv=True)[0])
    sim.warm_start()
    sim.init_views(warm=True)
    track = []

    def snap(s):
        v = [x for x in s.views if 0 in x.regions][0]
        track.append((s.slot, float(min(1.0, v.pcount[K] / p.presence))))

    sim.run(pre_min * 2, snap=(10, snap))
    t0 = sim.slot
    if poisoned:
        sim.split([[0, 1, 2], [3, 4]])
    else:
        sim.extra_delay[K] = dur_h * 3600.0
    sim.run(int(dur_h * SLOTS_PER_HOUR), snap=(10, snap))
    t_end = sim.slot
    if poisoned:
        sim.heal()
    else:
        sim.extra_delay[K] = 0.0
    sim.run(post_h * SLOTS_PER_HOUR, snap=(10, snap))
    snap(sim)
    recs = sim.recs()
    tr = np.array(track)
    during = tr[(tr[:, 0] >= t0) & (tr[:, 0] <= t_end)]
    after = tr[tr[:, 0] >= t_end]
    res = dict(min_part=float(tr[:, 1].min()), part_at_end=float(during[-1, 1]) if during.shape[0] else 1.0)
    below = tr[(tr[:, 0] >= t0) & (tr[:, 1] < 0.999)]
    res["drop_min"] = None if below.shape[0] == 0 else (below[0, 0] - t0) / 2.0
    rec = after[after[:, 1] >= 0.999]
    res["recover_min"] = None if rec.shape[0] == 0 else (rec[0, 0] - t_end) / 2.0
    res["conflicts"] = len(sim.conflicts)
    res["first_conflict_min"] = (min(c[1] for c in sim.conflicts) - t0) / 2.0 if sim.conflicts else None
    honest_sid = [v for v in [1]] if poisoned else [0]
    res["honest_stalls"] = stalls_between(recs, t0, t_end, sid=(1 if poisoned else 0))
    res["ecl_locks"] = int(((recs[:, 2] == 2) & (recs[:, 3] >= 0) & (recs[:, 0] >= t0) & (recs[:, 0] < t_end)).sum()) if poisoned else 0
    res["post_stalls"] = stalls_between(recs, t_end, sim.slot)
    return res


def scenario_f(args):
    out = ["### F. Eclipse of one pool holding 20% of weight", ""]
    configs = [P(denom="active", pmode="cert", delay=args.delay), P(denom="active", pmode="seen", delay=args.delay),
               P(denom="total", delay=args.delay)]
    rows = []
    for dur in (1, 2, 4):
        for p in configs:
            r = run_eclipse(args.seed, dur, False, p)
            rows.append([dur, p.label(), "%.3f" % r["min_part"], fm(r["drop_min"]), fm(r["recover_min"]), r["conflicts"],
                         r["honest_stalls"], r["post_stalls"]])
    out.append("F1. Delayed view: the pool receives every block and vote %s late, votes for the right blocks, late. Participation as the rest of the network computes it." % "D hours")
    out.append("")
    out.append(md_table(["eclipse h", "denominator", "pool participation, minimum", "first drop below 1, min after start",
                         "back to 1, min after end", "conflicting locks", "stalls during", "stalls after"], rows))
    out.append("")
    rows = []
    configs2 = configs + [P(denom="active", pmode="cert", floor=0.8, delay=args.delay),
                          P(denom="active", pmode="cert", floor=0.85, delay=args.delay)]
    for dur in (1, 2, 4):
        for p in configs2:
            r = run_eclipse(args.seed, dur, True, p)
            rows.append([dur, p.label(), "%.3f" % r["min_part"], fm(r["recover_min"]), r["conflicts"], fm(r["first_conflict_min"]),
                         r["ecl_locks"], r["honest_stalls"], r["post_stalls"]])
    out.append("F2. Poisoned view: an attacker holding 34% of weight feeds the pool a private fork for the eclipse, signs both forks, and is stripped at the heal. The pool votes for the attacker's checkpoints. Honest side = the other 46%.")
    out.append("")
    out.append(md_table(["eclipse h", "denominator", "pool participation, minimum (honest view)", "back to 1, min after end",
                         "conflicting locks", "first conflict, min after start", "locks on the eclipsed side",
                         "honest-side stalls during", "stalls after heal"], rows))
    return "\n".join(out)


def scenario_g(args):
    out = ["### G. Honest doubling overnight at day 60: 1,000 new keys, fresh Pareto draw, same total hashrate as the old 1,000", ""]
    rows_share = []
    rows_lat = []
    ev = []
    for denom in ("active", "total"):
        p = P(denom=denom, delay=args.delay)
        sim, rng, _ = build_honest(p, args.seed)
        h_new = pareto_hashrates(rng, N_HONEST, total=1.0)
        new = sim.add_keys(np.zeros(N_HONEST), assign_regions(h_new, GEOGRAPHY), flaky=True)
        old = np.arange(N_HONEST)
        sim.warm_start()
        sim.init_views(warm=True)
        sim.run(SLOTS_PER_HOUR)
        t_event = sim.slot
        sim.hash[new] = h_new
        sim.set_uptime()
        series = []
        last_23 = None
        for day in range(1, args.days_g + 1):
            sim.run(SLOTS_PER_DAY)
            so = sim.share(old)
            sn = sim.share(new)
            dust_new = int((sim.weight[new] < p.dust).sum())
            series.append((day, so, sn, dust_new))
            if so >= TWO_THIRDS:
                last_23 = day
        recs = sim.recs()
        if denom == "active":
            for d in (1, 5, 10, 15, 20, 21, 25):
                r = [x for x in series if x[0] == d]
                if r:
                    rows_share.append(["+%d" % d, pct(r[0][1]), pct(r[0][2]), pct(min(d / 60.0, 0.5)), r[0][3]])
        rows_lat.append(lat_row("%s, day before" % denom, lat_stats(recs, 0, t_event)))
        rows_lat.append(lat_row("%s, days 1 to 5 after" % denom, lat_stats(recs, t_event, t_event + 5 * SLOTS_PER_DAY)))
        rows_lat.append(lat_row("%s, days 5 to %d after" % (denom, args.days_g), lat_stats(recs, t_event + 5 * SLOTS_PER_DAY, None)))
        r45 = next((d for d, so, sn, _ in series if sn >= 0.45), None)
        r49 = next((d for d, so, sn, _ in series if sn >= 0.49), None)
        ev.append([denom, last_23, "not in run" if r45 is None else r45, "not in run" if r49 is None else r49,
                   stalls_between(recs, t_event, sim.slot)])
    out.append("Weight shares (active run):")
    out.append("")
    out.append(md_table(["day after doubling", "old cohort", "new cohort", "new cohort formula t/60", "new keys under dust"], rows_share))
    out.append("")
    out.append(md_table(["denominator", "last day old cohort holds 2/3 (locks alone)", "new cohort reaches 45%", "new cohort reaches 49%",
                         "stalled checkpoints"], ev))
    out.append("")
    out.append(md_table(LAT_HEADERS, rows_lat))
    return "\n".join(out)


# ---------------------------------------------------------------- additions, 3 October 2026, for section 3.11 Guarantees
# Scenarios H to K run the rule exactly as Q3 specifies it (active denominator, cert reading, floor 0.85, so a lock needs
# 2/3 of active and 17/30 of total) over several seeds. Nothing above this line changed.

NEW_SEEDS = (7, 11, 13, 17, 19)
P_FLOOR = 17.0 / 30.0


def rule_p(delay, **kw):
    """Q3 as specified: active/cert with the 0.85 floor."""
    base = dict(denom="active", pmode="cert", floor=0.85, delay=delay)
    base.update(kw)
    return P(**base)


def seeds_of(args):
    s = getattr(args, "seeds", "") or ""
    out = tuple(int(x) for x in s.split(",") if x.strip())
    return out or NEW_SEEDS


def span(vals, fmt="%d"):
    """'a' when every seed agrees, else 'a to b'."""
    vals = list(vals)
    lo, hi = min(vals), max(vals)
    return (fmt % lo) if lo == hi else (fmt % lo) + " to " + (fmt % hi)


def span_min(vals):
    """Minutes over seeds, 'never' when no seed produced the event, 'never in k of n' when some did."""
    vals = list(vals)
    got = [v for v in vals if v is not None]
    if not got:
        return "never"
    s = span(got, "%.0f")
    if len(got) < len(vals):
        s += " (never in %d of %d)" % (len(vals) - len(got), len(vals))
    return s


def lock_gaps_min(recs, s_from, s_to):
    """Longest run of consecutive slots without a lock inside [s_from, s_to), in minutes; the 'finality unavailable' interval."""
    locks = recs[(recs[:, 0] >= s_from) & (recs[:, 0] < s_to) & (recs[:, 3] >= 0), 0]
    edges = np.concatenate([[s_from], np.unique(locks), [s_to]])
    return float(np.max(np.diff(edges)) * SLOT_S / 60.0)


def run_partition2(seed, fracs, att_share, dur_min, p, pre_min=60, post_min=180):
    """run_partition with the heal check of section 3.11 item 4: every certificate any side held before the heal
    is in the merged view afterwards (a lock is never reversed), plus per-side lock counts."""
    rng = np.random.default_rng(seed)
    nR = max(3, len(fracs) + 1)
    sim = Sim(p, rng, n_regions=nR)
    h = pareto_hashrates(rng, N_HONEST)
    reg = assign_regions(h, fracs)
    sim.add_keys(h, reg, flaky=True)
    groups = [[i] for i in range(len(fracs))]
    att = None
    if att_share > 0:
        att = int(sim.add_keys([att_share / (1.0 - att_share)], [len(fracs)], flaky=False, equiv=True)[0])
        groups[0].append(len(fracs))
    sim.warm_start()
    sim.init_views(warm=True)
    sim.run(pre_min * 2)
    t_split = sim.slot
    sim.split(groups)
    sides = [v.sid for v in sim.views]
    sim.run(dur_min * 2)
    t_heal = sim.slot
    before = {}
    for v in sim.views:
        for idx, c in v.certs.items():
            before.setdefault(idx, set()).add(c[0])
    sim.heal()
    merged = sim.views[0].certs
    kept = all(idx in merged for idx in before)
    sim.run(post_min * 2)
    recs = sim.recs()
    res = dict(conflicts=len(sim.conflicts), kept=kept, pre_locks=len(before),
               att_share=sim.share([att]) if att is not None else 0.0)
    res["first_conflict_min"] = (min(c[1] for c in sim.conflicts) - t_split) / 2.0 if sim.conflicts else None
    res["side_first_lock"] = []
    res["side_locks"] = []
    for sid in sides:
        fl = first_lock_after(recs[recs[:, 0] < t_heal], t_split, sid=sid)
        res["side_first_lock"].append(None if fl is None else (fl - t_split) / 2.0)
        sel = recs[(recs[:, 0] >= t_split) & (recs[:, 0] < t_heal) & (recs[:, 2] == sid)]
        res["side_locks"].append(int((sel[:, 3] >= 0).sum()))
    fl = first_lock_after(recs, t_heal)
    res["post_first_lock_min"] = None if fl is None else (fl - t_heal) / 2.0
    res["post_stalls"] = stalls_between(recs, t_heal, sim.slot)
    return res


def scenario_h(args):
    """Partition with an equivocator under the rule as specified: the 4/30 bound of section 3.11."""
    seeds = seeds_of(args)
    q = getattr(args, "quick", False)
    durs = (60,) if q else (150, 360)
    atts = (0.0, 0.10, 0.13, 0.14, 0.20, 0.34)
    out = ["### H. Partition of a 50/50 honest network with an equivocating attacker, rule as specified (active/cert + floor 0.85), "
           "each side retargets at once (+daa, the median-time clock of Q1), seeds %s" % ",".join(str(s) for s in seeds), ""]
    out.append("Attacker = one key holding the stated share of TOTAL weight, mining on the first side, voting on both. Each side holds "
               "(1 - a)/2 + a of total. Prediction (section 3.11 item 2): a side holding s of total locks alone once s >= 17/30 and its "
               "view's active weight has decayed to 1.5 s, which under the cert reading is P x (1 - 1.5 s) slots after the split "
               "(P = 240, so 2 h x (1 - 1.5 s)); two sides over 17/30 need a >= 4/30 = 13.3%.")
    out.append("")
    rows = []
    for a in atts:
        s = (1.0 - a) / 2.0 + a
        pred = "no (side holds %s < 56.7%%)" % pct(s) if s < P_FLOOR else "%.0f min" % (120.0 * max(0.0, 1.0 - 1.5 * s))
        for dur in durs:
            rs = [run_partition2(sd, [0.5, 0.5], a, dur, rule_p(args.delay, daa="full")) for sd in seeds]
            rows.append([pct(a, 0), pct(s), dur, pred, span(r["conflicts"] for r in rs), span_min(r["first_conflict_min"] for r in rs),
                         " / ".join(span_min([r["side_first_lock"][i] for r in rs]) for i in range(2)),
                         "yes" if all(r["kept"] for r in rs) else "NO", span(r["post_stalls"] for r in rs)])
    out.append(md_table(["attacker (of total)", "each side holds", "partition min", "predicted first conflict", "conflicting locks",
                         "first conflict, min", "first lock per side, min", "every pre-heal lock kept at the heal", "stalls in 3 h after heal"], rows))
    return "\n".join(out)


def scenario_i(args):
    """The 40/40/20 split, three readings."""
    seeds = seeds_of(args)
    q = getattr(args, "quick", False)
    durs = (60,) if q else (150, 360)
    out = ["### I. The 40/40/20 split, rule as specified (active/cert + floor 0.85), +daa, seeds %s" % ",".join(str(s) for s in seeds), ""]
    cases = [("honest 40/40/20, no attacker", [0.4, 0.4, 0.2], 0.0),
             ("honest 40/40/20 plus a 10% equivocator (sides 36+10 / 36+10 / 18+10)", [0.4, 0.4, 0.2], 0.10),
             ("honest 40/40/20 plus a 20% equivocator (sides 32+20 / 32+20 / 16+20)", [0.4, 0.4, 0.2], 0.20),
             ("honest 40/40 plus a 20% equivocator reaching both (sides 40+20 / 40+20)", [0.5, 0.5], 0.20)]
    rows = []
    for name, fr, a in cases:
        for dur in durs:
            rs = [run_partition2(sd, fr, a, dur, rule_p(args.delay, daa="full")) for sd in seeds]
            n = len(fr)
            rows.append([name, dur, span(r["conflicts"] for r in rs), span_min(r["first_conflict_min"] for r in rs),
                         " / ".join(span([r["side_locks"][i] for r in rs]) for i in range(n)),
                         "yes" if all(r["kept"] for r in rs) else "NO",
                         span_min(r["post_first_lock_min"] for r in rs), span(r["post_stalls"] for r in rs)])
    out.append(md_table(["case", "partition min", "conflicting locks", "first conflict, min", "locks per side during",
                         "every pre-heal lock kept", "first lock after heal, min", "stalls in 3 h after heal"], rows))
    return "\n".join(out)


def run_silent_resume(seed, frac, hours, p, pre_h=3, post_h=3):
    sim, rng, _ = build_honest(p, seed)
    sim.warm_start()
    sim.init_views(warm=True)
    sim.run(pre_h * SLOTS_PER_HOUR)
    silent, got = pick_weight_subset(rng, sim.weight, frac)
    sim.signs[silent] = False
    t0 = sim.slot
    sim.run(int(hours * SLOTS_PER_HOUR))
    t1 = sim.slot
    sim.signs[silent] = True
    sim.run(post_h * SLOTS_PER_HOUR)
    recs = sim.recs()
    fl = first_lock_after(recs, t0)
    fr = first_lock_after(recs, t1)
    return dict(got=got, stalls=stalls_between(recs, t0, t1), first_lock=None if fl is None or fl >= t1 else (fl - t0) / 2.0,
                gap=lock_gaps_min(recs, t0, t1), resume=None if fr is None else (fr - t1) / 2.0,
                post_stalls=stalls_between(recs, t1, sim.slot), conflicts=len(sim.conflicts),
                locked_share=float(((recs[:, 0] >= t0) & (recs[:, 0] < t1) & (recs[:, 3] >= 0)).sum()) / max(1, int(((recs[:, 0] >= t0) & (recs[:, 0] < t1)).sum())))


def scenario_j(args):
    """Signing stops while mining continues, for 1, 6 and 24 hours, then resumes."""
    seeds = seeds_of(args)
    q = getattr(args, "quick", False)
    hours = (1,) if q else (1, 6, 24)
    fracs = (0.34, 0.40, 0.45)
    out = ["### J. Signing stops while mining continues for 1, 6 and 24 hours, then resumes; rule as specified, seeds %s" % ",".join(str(s) for s in seeds), ""]
    out.append("A random set holding x of weight stops signing at hour 3 and resumes after the stated time. Its weight never ages out because it keeps mining. "
               "'Longest gap' is the longest interval without a lock while they are silent: the time the node reports finality unavailable.")
    out.append("")
    rows = []
    for frac in fracs:
        for h in hours:
            rs = [run_silent_resume(sd, frac, h, rule_p(args.delay)) for sd in seeds]
            rows.append([pct(frac, 0), h, span_min(r["first_lock"] for r in rs), span(r["stalls"] for r in rs),
                         span(int(round(100 * r["locked_share"])) for r in rs) + "%",
                         span((r["gap"] for r in rs), "%.0f"), span_min(r["resume"] for r in rs), span(r["post_stalls"] for r in rs),
                         span(r["conflicts"] for r in rs)])
    out.append(md_table(["silent weight", "silent hours", "first lock after the stop, min", "stalled checkpoints while silent",
                         "checkpoints locked while silent", "longest gap without a lock, min", "first lock after resume, min",
                         "stalls in 3 h after resume", "conflicting locks"], rows))
    return "\n".join(out)


def run_acquired(seed, bought, att_hash, signs, days, p):
    """An attacker buys keys holding `bought` of window weight and starts mining at `att_hash` of network hashrate.
    The sellers keep their rigs and mine on under fresh keys."""
    sim, rng, _ = build_honest(p, seed)
    # the bought set is picked by hashrate, which the warm start turns into weight at the same share
    mask, got = pick_weight_subset(rng, sim.hash, bought)
    bidx = np.flatnonzero(mask)
    seller_hash = sim.hash[bidx].copy()
    fresh = sim.add_keys(np.zeros(bidx.size), sim.region[bidx].copy(), flaky=True)
    att = int(sim.add_keys([0.0], [0], flaky=False, signs=signs)[0])
    sim.warm_start()
    sim.init_views(warm=True)
    sim.run(SLOTS_PER_HOUR)
    honest = sim.hash.sum()
    sim.hash[fresh] = seller_hash
    sim.hash[bidx] = 0.0
    sim.flaky[bidx] = False
    sim.online[bidx] = True
    sim.signs[bidx] = signs
    sim.hash[att] = honest * att_hash / (1.0 - att_hash)
    sim.set_uptime()
    got = float(sim.share(bidx))
    owned = np.concatenate([bidx, [att]])
    t0 = sim.slot
    series = []
    above13 = None
    last13 = None
    for day in range(1, days + 1):
        s0 = sim.slot
        sim.run(SLOTS_PER_DAY)
        sh = sim.share(owned)
        recs = sim.recs()
        series.append((day, sh, stalls_between(recs, s0, sim.slot)))
        if sh >= 1.0 / 3.0:
            last13 = day
            if above13 is None:
                above13 = day
    recs = sim.recs()
    return dict(got=got, series=series, above13=above13, last13=last13, stalls=stalls_between(recs, t0, sim.slot),
                max_share=max(s for _, s, _ in series), end_share=series[-1][1], conflicts=len(sim.conflicts))


def scenario_k(args):
    """Acquired old keys against fresh hashrate."""
    seeds = seeds_of(args)
    q = getattr(args, "quick", False)
    days = 3 if q else 30
    att_hash = 0.30
    out = ["### K. Acquired keys: an attacker buys keys holding 20%% or 40%% of window weight and mines at 30%% of network hashrate from day 0; "
           "rule as specified, %d days, seeds %s" % (days, ",".join(str(s) for s in seeds)), ""]
    out.append("The sellers keep their rigs and mine on under fresh keys (so honest hashrate is unchanged and the fresh keys start under dust). "
               "Formula (section 3.11 item 5): share(t) = b (1 - t/30) + 0.30 t/30, the bought blocks age out of the window as the attacker's own blocks enter it. "
               "'Fresh hash only' is b = 0: share(t) = 0.30 t/30.")
    out.append("")
    configs = [(0.0, True), (0.20, True), (0.40, True), (0.20, False), (0.40, False)]
    results = {}
    for b, signs in configs:
        results[(b, signs)] = [run_acquired(sd, b, att_hash, signs, days, rule_p(args.delay)) for sd in seeds]
    pick = [d for d in (1, 5, 10, 15, 20, 25, 30) if d <= days]
    rows = []
    for d in pick:
        row = ["+%d" % d]
        for b, signs in configs:
            if not signs:
                continue
            vals = [dict((x[0], x[1]) for x in r["series"])[d] for r in results[(b, signs)]]
            formula = b * (1.0 - d / 30.0) + att_hash * d / 30.0
            row.append("%s / %s" % (span((100 * v for v in vals), "%.1f") + "%", pct(formula)))
        rows.append(row)
    out.append("Attacker weight share, simulated over the seeds / formula (attacker signs every checkpoint):")
    out.append("")
    out.append(md_table(["day after purchase"] + ["bought %s" % pct(b, 0) for b, s in configs if s], rows))
    out.append("")
    ev = []
    for b, signs in configs:
        rs = results[(b, signs)]
        ev.append([pct(b, 0), "signs" if signs else "silent", span((100 * r["got"] for r in rs), "%.1f") + "%",
                   span((100 * r["max_share"] for r in rs), "%.1f") + "%", span((100 * r["end_share"] for r in rs), "%.1f") + "%",
                   "never" if all(r["above13"] is None for r in rs) else "from day %s until day %s" % (span(r["above13"] for r in rs), span(r["last13"] for r in rs)),
                   span(r["stalls"] for r in rs), span(r["conflicts"] for r in rs)])
    out.append(md_table(["bought weight", "attacker", "bought, as picked", "peak share", "share at day %d" % days,
                         "holds at least 1/3 (can veto)", "stalled checkpoints in %d days" % days, "conflicting locks"], ev))
    out.append("")
    srows = []
    for b, signs in configs:
        if signs:
            continue
        rs = results[(b, signs)]
        for d in pick:
            vals = [dict((x[0], x[2]) for x in r["series"])[d] for r in rs]
            srows.append([pct(b, 0), "+%d" % d, span(vals)])
    out.append("Stalled checkpoints per day while the attacker withholds its votes (2,880 checkpoints a day):")
    out.append("")
    out.append(md_table(["bought weight", "day", "stalled that day"], srows))
    return "\n".join(out)


SCENARIOS = {"A": scenario_a, "B": scenario_b, "C": scenario_c, "D": scenario_d, "E": scenario_e, "F": scenario_f, "G": scenario_g,
             "H": scenario_h, "I": scenario_i, "J": scenario_j, "K": scenario_k}


def main(argv=None):
    ap = argparse.ArgumentParser(description="Igneum finality rule V2 simulation")
    ap.add_argument("--seed", type=int, default=7)
    ap.add_argument("--scenarios", default="A,B,C,D,E,F,G")
    ap.add_argument("--delay", type=float, default=2.0, help="one-way inter-region delay in seconds for the main runs")
    ap.add_argument("--grace", type=float, default=15.0)
    ap.add_argument("--quick", action="store_true", help="shortened runs for development")
    ap.add_argument("--seeds", default="", help="comma-separated seeds for scenarios H to K (default 7,11,13,17,19)")
    args = ap.parse_args(argv)
    q = args.quick
    args.days_a = 3 if q else 60
    args.days_delay = 1 if q else 3
    args.days_b = 4 if q else 35
    args.hours_c = 3 if q else 6
    args.hours_c_total = 3 if q else 72
    args.days_d35 = 1 if q else 3
    args.days_d50 = 1 if q else 12
    args.days_g = 3 if q else 25
    print("# finality_v2 output")
    print()
    p = P()
    print("seed %d, slot %.0f s, %d blocks per checkpoint, window %d h, presence %d checkpoints, dust %d, quorum 2/3, "
          "inter-region delay %.1f s (intra %.1f s, jitter sigma %.2f), grace %.0f s, uptime %.2f (mean outage %d slots), "
          "uptime %.3f for keys at or above %s of hashrate, honest keys %d Pareto %.1f, geography %s%s" % (
              args.seed, SLOT_S, BLOCKS_PER_CP, WINDOW_HOURS, p.presence, p.dust, args.delay, p.intra, p.jitter,
              args.grace, p.uptime, p.outage, p.uptime_big, pct(p.big_share, 0), N_HONEST, PARETO_SHAPE, "/".join(pct(g, 0) for g in GEOGRAPHY),
              ", QUICK" if q else ""))
    print()
    for s in args.scenarios.split(","):
        s = s.strip().upper()
        if s not in SCENARIOS:
            print("unknown scenario %s" % s, file=sys.stderr)
            return 2
        t = time.time()
        print(SCENARIOS[s](args))
        print()
        print("(%s took %.0f s)" % (s, time.time() - t), file=sys.stderr)
    return 0


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