igneum/sim/difficulty/attacks/attacks.py

#!/usr/bin/env python3
"""Adversarial runs against the Igneum difficulty rule (spec 2.3) and Kaspa's sampled DAA.

Builds on ../sim.py (the controllers are imported unchanged). Adds what an attacker needs that the base
simulator does not model: several miners with their own on/off strategies, attribution of each block to the
miner that found it (drawn by hash share at the moment of the solve), hashes spent per miner, and timestamp
forging that stays inside the consensus rules (header timestamp at most FTL = 132 s ahead of the clock, and
strictly above the sampled past median time: 27 samples at 10-block intervals, mean of the central 11, as
`SampledPastMedianTimeManager::calc_past_median_time` does it). Honest miners stamp real time, raised to the
past median plus one when the rule forces it.

Scenarios (README.md): hop (pool hopping for profit), pulse (50x bursts for 10 min every hour, and once),
ts (timestamp stretching), osc (short-lane oscillation), epoch (hold dodger and hold flooder), pollute (leave
as the long lane takes over). Every scenario runs for each rule named: igneum, kaspa, and the candidates
igneum-san (sanitised running clock, README.md), igneum-san10 (the same with a 10 T cap), igneum-sanflat
(the same with flat weights on the short lane).

Blocks per hash of a miner is the expected-gap estimator sum(1 / H_b) / sum(D_b / H_b) over the blocks it was
on for (H_b total hash rate, D_b expected hashes per block), which has no block-count noise; the realised
counts are kept where they matter (pulse block share).

    python3 attacks.py                                   every scenario, 3 seeds, igneum and kaspa
    python3 attacks.py --scenario ts --rules igneum,igneum-san,kaspa
    python3 attacks.py --base --rules igneum,igneum-san  the base simulator's profiles through the candidates
"""
import argparse
import math
import os
import random
import sys
from multiprocessing import Pool

sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
import sim as S  # noqa: E402

T = S.T
EPOCH = S.EPOCH
ONE = S.ONE
SCALE = S.SCALE
H1 = S.H1
HOUR = 3600 * 1000
FTL_MS = 132 * 1000
PMT_SAMPLES = 27
PMT_INTERVAL = 10
PMT_CENTRE = 11
clamp = S.clamp


# ---------------------------------------------------------------- timestamp rules
def past_median_time(ts, h):
    """Past median time for the block at height h from the chain timestamps ts[0..h-1]:
    27 samples at 10-block intervals ending at the parent, sorted, mean of the central 11."""
    samp = [ts[j] for j in range(h - 1, -1, -PMT_INTERVAL)][:PMT_SAMPLES]
    samp.sort()
    n = len(samp)
    if n <= PMT_CENTRE:
        return sum(samp) // n
    lo = (n - PMT_CENTRE) // 2
    return sum(samp[lo:lo + PMT_CENTRE]) // PMT_CENTRE


class Rules:
    """Timestamp rules the forger must respect. Kaspa's: 132 s ahead of the clock, above the sampled past
    median. The proposed tight rules (README.md): 10 s ahead of the clock, above the past median and at least the
    selected parent's timestamp minus 10 s."""

    def __init__(self, ftl_ms=FTL_MS, back_ms=None):
        self.ftl_ms = ftl_ms
        self.back_ms = back_ms

    def bounds(self, ts, h, now):
        lo = past_median_time(ts, h) + 1
        if self.back_ms is not None:
            lo = max(lo, ts[h - 1] - self.back_ms)
        hi = now + self.ftl_ms
        return lo, max(lo, hi)


TS_RULES = {"kaspa": lambda: Rules(), "tight": lambda: Rules(ftl_ms=10 * 1000, back_ms=10 * 1000)}


# ---------------------------------------------------------------- candidate controllers
class IgneumSan(S.Igneum):
    """Spec 2.3 with one change (README.md): every chain step's solvetime is measured on a sanitised running
    clock. ts'(genesis) = ts(genesis); ts'(b) = ts'(parent) + clamp(ts(b) - ts'(parent), -cap, +cap). The short
    and epoch lanes sum sanitised steps, which telescope to ts'(newest) - ts'(oldest), so a forged stamp is paid
    back exactly by the honest blocks after it instead of being cancelled to zero measured time. The long lane
    keeps the raw sampled span. `flat` drops the linear weights of the short lane; `cap` as in the spec."""
    name = "igneum-san"

    def __init__(self, g, flat=False, lag=0, snap_ms=60 * 1000, **kw):
        super().__init__(g, **kw)
        self.flat = flat
        self.lag = lag
        self.snap_ms = snap_ms
        self.tss = []
        self.sts = []

    def push(self, ts, target):
        h = len(self.ts)
        if h == 0:
            st, tss = 0, ts
        else:
            tss = self.tss[-1] + clamp(ts - self.tss[-1], -self.cap, self.cap)
            if self.snap_ms is not None:
                # the clock never lags a raw stamp by more than snap_ms (an honest idle gap must not be paid back
                # as cap-sized steps for many blocks); the jump itself is not counted as solvetime. A backward
                # forger would need a run of snap_ms / 10 s blocks under the tight rules to reach it.
                tss = max(tss, ts - self.snap_ms)
            st = min(tss - self.tss[-1], self.cap)
        self.tss.append(tss)
        self.sts.append(st)
        S.Controller.push(self, ts, target)
        self.pref_d.append(self.pref_d[-1] + ONE / target)
        self.pref_t.append(self.pref_t[-1] + target)
        self.pref_st.append(self.pref_st[-1] + st)
        if h > 0 and h % self.rate == self.rate - 1:
            self.samples.append((h, ts, self.pref_d[-1]))
            if len(self.samples) > self.window:
                self.samples.popleft()

    def rate_s(self, lo, hi):
        n = hi - lo - 1
        if n < 1:
            return None
        wsum = 0.0
        dsum = 0.0
        for i in range(1, n + 1):
            j = lo + i
            w = 1 if self.flat else i
            wsum += self.sts[j] * w
            dsum += (self.pref_d[j + 1] - self.pref_d[j]) * w
        wsum = max(wsum, (n * T / 10) if self.flat else (n * n * T / 20))
        return dsum / wsum

    def next_target(self):
        """As the spec, with the newest `lag` chain steps left out of the short and epoch lanes."""
        if self.lag == 0:
            return super().next_target()
        h = len(self.ts)
        parent = self.tg[-1]
        epoch_start = (h // EPOCH) * EPOCH
        k = h - epoch_start
        if k < self.k_min:
            return parent
        hh = max(h - self.lag, epoch_start + 2)
        s_lo = max(epoch_start, hh - self.ns)
        r_s = self.rate_s(s_lo, hh)
        r_ref = None
        if k >= self.long_min:
            r_ref = self.rate_l(max(epoch_start, h - self.window * self.rate))
        if r_ref is None:
            r_ref = self.rate_e(epoch_start, hh, hh - epoch_start)
        if r_s is None:
            cand = parent
        else:
            use_s = abs(r_s / r_ref - 1) > self.trig
            self.engaged += 1 if use_s else 0
            cand = ONE / ((r_s if use_s else r_ref) * T)
        out = clamp(cand, parent / (1 + self.clamp_pct), parent * (1 + self.ease_pct))
        return min(out, ONE)


_orig_make = S.make_controller


def make_controller(name, g, **kw):
    if name == "igneum-san":
        return IgneumSan(g, **kw)
    if name == "igneum-san-nosnap":
        return IgneumSan(g, snap_ms=None, **kw)
    if name.startswith("igneum-ease"):        # igneum-ease<percent>: the spec rule with another ease clamp
        return S.Igneum(g, ease_pct=int(name[len("igneum-ease"):]) / 100.0, **kw)
    if name == "igneum-san-j120":
        return IgneumSan(g, snap_ms=120 * 1000, **kw)
    if name == "igneum-san10":
        return IgneumSan(g, cap=10 * T, **kw)
    if name == "igneum-sanflat":
        return IgneumSan(g, flat=True, **kw)
    if name.startswith("igneum-san-"):      # igneum-san-<flat|lwma>-lag<n>-cap<n>
        parts = name.split("-")[2:]
        flat = parts[0] == "flat"
        lag = int(parts[1][3:]) if len(parts) > 1 else 0
        cap = int(parts[2][3:]) * T if len(parts) > 2 else 20 * T
        return IgneumSan(g, flat=flat, lag=lag, cap=cap, **kw)
    return _orig_make(name, g, **kw)


S.make_controller = make_controller


def make(rule, genesis_d):
    return make_controller(rule, ONE / genesis_d)


# ---------------------------------------------------------------- miners
class Miner:
    """hash in units of H1 (the hash rate that gives 1 block/s at D = SCALE). policy(state) -> on for the next
    block (evaluated at block boundaries); gate(t) -> on at time t and schedule(t) -> next gate change, for
    time-keyed bursts; stamp: honest | future | past | alt (future and past on alternate blocks)."""

    def __init__(self, name, hash_h1, policy=None, stamp="honest", start_ms=0):
        self.name = name
        self.hash = hash_h1 * H1
        self.policy = policy or (lambda st: st["t"] >= self.start_ms)
        self.gate = None
        self.schedule = None
        self.stamp = stamp
        self.start_ms = start_ms
        self.hashes = 0.0
        self.blocks = 0
        self.exp_blocks = 0.0
        self.on = False
        self.last_switch = -1e18
        self.flip = 0

    def live(self, t):
        return self.on and t >= self.start_ms and (self.gate is None or self.gate(t))

    def stamp_for(self, t, lo, hi):
        if self.stamp == "future":
            return hi
        if self.stamp == "past":
            return lo
        if self.stamp == "alt":
            self.flip ^= 1
            return hi if self.flip else lo
        return max(int(t), lo)


class AvgD:
    """Trailing time-weighted mean of D over the last 6 hours (the hopper's long-run average)."""

    def __init__(self, span_ms=6 * HOUR):
        self.span = span_ms
        self.q = []
        self.sum = 0.0
        self.wsum = 0.0
        self.i = 0

    def add(self, t, D, dt):
        self.q.append((t, D, dt))
        self.sum += D * dt
        self.wsum += dt
        while self.i < len(self.q) and self.q[self.i][0] < t - self.span:
            _, d, w = self.q[self.i]
            self.sum -= d * w
            self.wsum -= w
            self.i += 1

    def value(self, default):
        return self.sum / self.wsum if self.wsum > 0 else default


def run(miners, ctrl, rng, duration_ms, genesis_d, rules=None, log=None):
    """Mine duration_ms of chain. Rows: (t, h, D, H, r, finder, active names)."""
    rules = rules or Rules()
    t = 0.0
    target = ONE / genesis_d
    ctrl.push(0, target)
    h = 1
    rows = []
    avg = AvgD()
    st = {"t": 0.0, "h": 1, "D": genesis_d, "avgD": genesis_d, "k": 1}
    while t < duration_ms:
        target = ctrl.next_target()
        D = ONE / target
        st.update(t=t, h=h, D=D, k=h % EPOCH, avgD=avg.value(genesis_d))
        for m in miners:
            want = bool(m.policy(st))
            if want != m.on:
                m.on = want
                m.last_switch = t
        work = rng.expovariate(1.0) * D
        tt = t
        while True:
            live = [m for m in miners if m.live(tt)]
            hr = sum(m.hash for m in live)
            nxt = None
            for m in miners:
                for x in ([m.start_ms] if m.start_ms > tt else []) + ([m.schedule(tt)] if m.schedule else []):
                    if x is not None and x > tt and (nxt is None or x < nxt):
                        nxt = x
            if hr > 0 and (nxt is None or work <= hr * (nxt - tt)):
                dt = work / hr
                for m in live:
                    m.hashes += m.hash * dt
                tt += dt
                break
            if nxt is None:
                tt = duration_ms + 1
                break
            dt = nxt - tt
            for m in live:
                m.hashes += m.hash * dt
            work -= hr * dt
            tt = nxt
        if tt > duration_ms:
            break
        gap = max(tt - t, 1.0)
        t += gap
        live = [m for m in miners if m.live(t)]
        tot = sum(m.hash for m in live)
        u = rng.random() * tot
        finder = live[-1]
        acc = 0.0
        for m in live:
            acc += m.hash
            if u <= acc:
                finder = m
                break
        for m in live:
            m.exp_blocks += m.hash / tot
        finder.blocks += 1
        lo, hi = rules.bounds(ctrl.ts, h, int(t))
        ts = min(max(finder.stamp_for(t, lo, hi), lo), hi)
        ctrl.push(ts, target)
        avg.add(t, D, gap)
        rows.append((t, h, D, tot, tot * T / D, finder.name, tuple(m.name for m in live)))
        if log is not None:
            log.append((t, h, D, tot, finder.name, ts, int(t)))
        h += 1
    return rows


# ---------------------------------------------------------------- metrics
def bph(rows, name, t_from=0):
    """Expected blocks per hash of a miner over the blocks it was on for after t_from: sum(1/H) / sum(D/H)."""
    num = den = 0.0
    for r in rows:
        if r[0] >= t_from and name in r[6]:
            num += 1.0 / r[3]
            den += r[2] / r[3]
    return num / den if den > 0 else float("nan")


def on_frac(rows, name, t_from):
    act = tot = 0.0
    prev = None
    for r in rows:
        if prev is not None and r[0] >= t_from:
            g = r[0] - prev
            tot += g
            if name in r[6]:
                act += g
        prev = r[0]
    return act / tot if tot else 0.0


def rate_stats(rows, t_from, t_to):
    seg = [r for r in rows if t_from <= r[0] < t_to]
    if len(seg) < 10:
        return None, None, None, None
    rs = [r[4] for r in seg]
    mean = sum(rs) / len(rs)
    std = math.sqrt(sum((x - mean) ** 2 for x in rs) / len(rs))
    gaps = [seg[i][0] - seg[i - 1][0] for i in range(1, len(seg))]
    return mean, std, max(gaps) / 1000, len(seg) / ((t_to - t_from) / 1000)


def max_blocks_in(rows, t0, bin_ms):
    bins = {}
    for r in rows:
        if r[0] >= t0:
            b = int(r[0] // bin_ms)
            bins[b] = bins.get(b, 0) + 1
    return max(bins.values()) if bins else 0


# ---------------------------------------------------------------- scenarios
def scen_hop(rule, seed, x_pct, dwell_s=0, hours=24):
    """Pool hopper with x% of the honest base, on only while D is below its trailing 6-hour mean (1% band),
    base constant. Starts at hour 2; measured from hour 3."""
    rng = random.Random(seed)
    base = Miner("base", 1.0)
    band = 0.01

    def policy(st):
        m = hopper
        if st["t"] < m.start_ms:
            return False
        if st["t"] - m.last_switch < dwell_s * 1000:
            return m.on
        if m.on:
            return st["D"] <= st["avgD"] * (1 + band)
        return st["D"] < st["avgD"] * (1 - band)

    hopper = Miner("hopper", x_pct / 100.0, policy=policy, start_ms=2 * HOUR)
    rows = run([base, hopper], make(rule, SCALE), rng, hours * HOUR, SCALE)
    t0 = 3 * HOUR
    switches = 0
    on = False
    for r in rows:
        if r[0] >= t0 and ("hopper" in r[6]) != on:
            switches += 1
            on = not on
    _, std, gap, bps = rate_stats(rows, t0, rows[-1][0])
    return {"rule": rule, "seed": seed, "x": x_pct, "dwell": dwell_s, "adv": bph(rows, "hopper", t0) / bph(rows, "base", t0) - 1,
            "hop_on_frac": on_frac(rows, "hopper", t0), "switches": switches, "std": std, "worst_gap": gap, "blocks_per_s": bps}


def scen_pulse(rule, seed, mult=50.0, burst_s=600, period_s=3600, hours=8):
    """A rented burst of mult x the base for burst_s every period_s from hour 2 (period_s >= hours*3600: once)."""
    rng = random.Random(seed)
    base = Miner("base", 1.0)
    pulser = Miner("pulser", mult, start_ms=2 * HOUR)
    P = period_s * 1000
    B = burst_s * 1000
    s0 = pulser.start_ms
    pulser.gate = lambda t: t >= s0 and ((t - s0) % P) < B
    pulser.schedule = lambda t: s0 if t < s0 else s0 + ((t - s0) // P) * P + (B if ((t - s0) % P) < B else P)
    rows = run([base, pulser], make(rule, SCALE), rng, hours * HOUR, SCALE)
    t0 = 2 * HOUR
    seg = [r for r in rows if r[0] >= t0]
    b_base = bph(rows, "base", t0)
    b_pulse = bph(rows, "pulser", t0)
    steady = bph([r for r in rows if HOUR <= r[0] < 2 * HOUR], "base")
    share_blocks = sum(1 for r in seg if r[5] == "pulser") / len(seg)
    share_hash = pulser.hashes / (pulser.hashes + sum(base.hash * (seg[i][0] - seg[i - 1][0]) for i in range(1, len(seg))))
    _, std, gap, bps = rate_stats(rows, t0, rows[-1][0])
    return {"rule": rule, "seed": seed, "period_s": period_s, "pulse_vs_base": b_pulse / b_base - 1, "pulse_vs_steady": b_pulse / steady - 1,
            "base_vs_steady": b_base / steady - 1, "weight_per_hash": share_blocks / share_hash, "pulse_block_share": share_blocks,
            "pulse_hash_share": share_hash, "peak_blocks_per_min": max_blocks_in(rows, t0, 60000), "worst_gap": gap, "blocks_per_s": bps}


def scen_ts(rule, seed, share_pct, stamp, hours=6, rules="kaspa"):
    """A miner with share_pct of the total hash rate stamping every block at the latest allowed future offset
    (future), the earliest allowed (past) or alternating (alt) from hour 2; honest for the first 2 hours. Base
    honest and constant; genesis sized for both."""
    rng = random.Random(seed)
    a = share_pct / (100.0 - share_pct)
    base = Miner("base", 1.0)
    att = Miner("forger", a)
    gd = SCALE * (1 + a)
    orig = att.stamp_for

    def stamp_for(t, lo, hi):
        if t < 2 * HOUR:
            return max(int(t), lo)
        att.stamp = stamp
        return orig(t, lo, hi)

    att.stamp_for = stamp_for
    log = []
    rows = run([base, att], make(rule, gd), rng, hours * HOUR, gd, TS_RULES[rules](), log)
    pre = rate_stats(rows, 1 * HOUR, 2 * HOUR)
    post = rate_stats(rows, 3 * HOUR, hours * HOUR)
    tail = rate_stats(rows, (hours - 1) * HOUR, hours * HOUR)
    seg = [r for r in rows if r[0] >= 3 * HOUR]
    mean_D = sum(r[2] for r in seg) / len(seg) if seg else float("nan")
    offs = [l[5] - l[6] for l in log if l[0] >= 3 * HOUR and l[4] == "forger"]
    return {"rule": rule, "seed": seed, "ts_rules": rules, "share": share_pct, "stamp": stamp, "forger_offset_s": (sum(offs) / len(offs) / 1000) if offs else 0.0,
            "rate_pre": pre[3], "rate_post": post[3], "rate_last_hour": tail[3], "drift": post[3] - 1.0, "D_ratio": mean_D / gd,
            "std": post[1], "worst_gap": post[2], "forger_adv": bph(rows, "forger", 3 * HOUR) / bph(rows, "base", 3 * HOUR) - 1}


def scen_osc(rule, seed, share_pct=25, period_blocks=120, hours=6):
    """A miner of share_pct of the full hash rate on and off every period_blocks from hour 2."""
    rng = random.Random(seed)
    a = share_pct / (100.0 - share_pct)
    base = Miner("base", 1.0)
    osc = Miner("osc", a, start_ms=2 * HOUR)
    osc.policy = lambda st: st["t"] >= 2 * HOUR and ((st["h"] // period_blocks) % 2 == 0)
    rows = run([base, osc], make(rule, SCALE), rng, hours * HOUR, SCALE)
    pre = rate_stats(rows, 1 * HOUR, 2 * HOUR)
    post = rate_stats(rows, 3 * HOUR, hours * HOUR)
    t0 = 3 * HOUR
    seg = [r for r in rows if r[0] >= t0]
    hrs = [r[3] for r in seg]
    mh = sum(hrs) / len(hrs)
    floor = math.sqrt(sum((x / mh - 1) ** 2 for x in hrs) / len(hrs))
    return {"rule": rule, "seed": seed, "share": share_pct, "period": period_blocks, "std_pre": pre[1], "std_post": post[1],
            "std_ratio": post[1] / pre[1], "std_floor": floor, "worst_gap_pre": pre[2], "worst_gap_post": post[2],
            "blocks_per_s": post[3], "osc_adv": bph(rows, "osc", t0) / bph(rows, "base", t0) - 1}


def scen_epoch(rule, seed, kind, share_pct=30, hours=8):
    """Epoch-boundary games from hour 2. dodge: a share_pct miner is off for the 8 hold blocks of each epoch.
    flood: a 10x miner is on only for the hold blocks. flood30: a share_pct miner on only for the hold blocks."""
    rng = random.Random(seed)
    a = share_pct / (100.0 - share_pct)
    base = Miner("base", 1.0)
    if kind == "dodge":
        g = Miner("gamer", a, start_ms=2 * HOUR)
        g.policy = lambda st: st["t"] >= 2 * HOUR and st["k"] >= 8
        gd = SCALE * (1 + a)
    elif kind == "flood":
        g = Miner("gamer", 10.0, start_ms=2 * HOUR)
        g.policy = lambda st: st["t"] >= 2 * HOUR and st["k"] < 8
        gd = SCALE
    elif kind == "flood30":
        g = Miner("gamer", a, start_ms=2 * HOUR)
        g.policy = lambda st: st["t"] >= 2 * HOUR and st["k"] < 8
        gd = SCALE
    else:
        raise ValueError(kind)
    rows = run([base, g], make(rule, gd), rng, hours * HOUR, gd)
    t0 = 3 * HOUR
    post = rate_stats(rows, t0, hours * HOUR)
    return {"rule": rule, "seed": seed, "kind": kind, "gamer_adv": bph(rows, "gamer", t0) / bph(rows, "base", t0) - 1,
            "gamer_on_frac": on_frac(rows, "gamer", t0), "std": post[1], "worst_gap": post[2], "blocks_per_s": post[3]}


def scen_pollute(rule, seed, mult=10.0, join_k=0, leave_k=600, hours=4):
    """A mult x miner joins at block join_k of epoch 2 and leaves at block leave_k of the same epoch (600 is where
    the reference lane switches from the epoch window to the sampled long lane)."""
    rng = random.Random(seed)
    base = Miner("base", 1.0)
    big = Miner("big", mult)
    e = 2 * EPOCH
    big.policy = lambda st: e + join_k <= st["h"] < e + leave_k
    rows = run([base, big], make(rule, SCALE), rng, hours * HOUR, SCALE)
    t_join = next((r[0] for r in rows if r[1] >= e + join_k), None)
    t_leave = next((r[0] for r in rows if r[1] >= e + leave_k), None)
    blocks = [r[:5] for r in rows]
    m_leave = S.metrics(blocks, None, t_step=t_leave) if t_leave else {}
    m_join = S.metrics(blocks, None, t_step=t_join, t_end=t_leave) if t_join else {}
    return {"rule": rule, "seed": seed, "mult": mult, "join_k": join_k, "leave_k": leave_k,
            "join_first10_s": m_join.get("first10_s"), "join_settled_s": m_join.get("recover_s"),
            "leave_first10_s": m_leave.get("first10_s"), "leave_settled_s": m_leave.get("recover_s"),
            "leave_worst_gap_s": m_leave.get("worst_gap_s"), "leave_below05": m_leave.get("below05"), "leave_above2x": m_leave.get("above2x")}


def base_profiles(rule, seed):
    """The base simulator's synthetic profiles through one controller (its own metrics)."""
    out = []
    for p in S.synthetic_profiles(seed):
        if p.name in ("up50", "down50", "epoch30", "walk10", "hop10", "polluted", "warmup-hard"):
            r = S.run_profile(p, [rule], seed, 0)[0]
            out.append({"rule": rule, "seed": seed, "profile": p.name, "first10_s": r.get("first10_s"), "settled_s": r.get("recover_s"),
                        "worst_gap_s": r.get("worst_gap_s") or r.get("ss_worst_gap_s"), "steady_std": r.get("ss_std"),
                        "peak_rate": r.get("peak_rate"), "not_recovered": r.get("recover_fail")})
    return out


# ---------------------------------------------------------------- driver
def job(args):
    fn, kw = args
    return globals()[fn](**kw)


def fmt(x):
    if x is None:
        return "none"
    if isinstance(x, float):
        if math.isnan(x):
            return "nan"
        return ("%.3f" % x) if abs(x) < 10 else ("%.1f" % x)
    return str(x)


def table(title, rows, cols):
    print("\n### " + title + "\n")
    print("| " + " | ".join(cols) + " |")
    print("|" + "---|" * len(cols))
    for r in rows:
        print("| " + " | ".join(fmt(r.get(c)) for c in cols) + " |")


def mean_rows(rows, keys, worst=()):
    groups = {}
    for r in rows:
        groups.setdefault(tuple(r[x] for x in keys), []).append(r)
    out = []
    for k, rs in groups.items():
        m = dict(rs[0])
        m["seeds"] = len(rs)
        for f in rs[0]:
            if f in keys or f == "seed":
                continue
            vals = [r[f] for r in rs if isinstance(r[f], (int, float)) and not (isinstance(r[f], float) and math.isnan(r[f]))]
            if len(vals) == len(rs):
                m[f] = sum(vals) / len(vals)
                if f in worst:
                    m[f + "_worst"] = max(vals, key=abs)
            else:
                m[f] = None if not vals else sum(vals) / len(vals)
                m[f + "_none"] = len(rs) - len(vals)
        out.append(m)
    return out


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--scenario", default="all")
    ap.add_argument("--seeds", type=int, default=3)
    ap.add_argument("--rules", default="igneum,kaspa")
    ap.add_argument("--base", action="store_true", help="also run the base simulator's profiles through each rule")
    ap.add_argument("--jobs", type=int, default=4)
    ap.add_argument("--hop-hours", type=int, default=24)
    ap.add_argument("--ts-rules", default="kaspa", help="timestamp rules the forger obeys: kaspa or tight (README.md)")
    args = ap.parse_args()
    seeds = list(range(7, 7 + args.seeds))
    rules = args.rules.split(",")
    want = args.scenario.split(",") if args.scenario != "all" else ["hop", "pulse", "ts", "osc", "epoch", "pollute"]
    jobs = []
    for rule in rules:
        for seed in seeds:
            if "hop" in want:
                for x in (10, 30, 50, 100):
                    for dwell in (0, 60):
                        jobs.append(("scen_hop", dict(rule=rule, seed=seed, x_pct=x, dwell_s=dwell, hours=args.hop_hours)))
            if "pulse" in want:
                jobs.append(("scen_pulse", dict(rule=rule, seed=seed)))
                jobs.append(("scen_pulse", dict(rule=rule, seed=seed, period_s=10 ** 6, hours=4)))
            if "ts" in want:
                for share in (30, 50):
                    for stamp in ("future", "past", "alt"):
                        jobs.append(("scen_ts", dict(rule=rule, seed=seed, share_pct=share, stamp=stamp, rules=args.ts_rules)))
            if "osc" in want:
                jobs.append(("scen_osc", dict(rule=rule, seed=seed)))
            if "epoch" in want:
                for kind in ("dodge", "flood", "flood30"):
                    jobs.append(("scen_epoch", dict(rule=rule, seed=seed, kind=kind)))
            if "pollute" in want:
                for join_k, leave_k in ((0, 600), (480, 600), (480, 720), (0, 1200)):
                    jobs.append(("scen_pollute", dict(rule=rule, seed=seed, join_k=join_k, leave_k=leave_k)))
            if args.base:
                jobs.append(("base_profiles", dict(rule=rule, seed=seed)))
    try:
        os.nice(19)
    except OSError:
        pass
    with Pool(args.jobs) as pool:
        results = pool.map(job, jobs)
    by = {}
    for (fn, _), r in zip(jobs, results):
        by.setdefault(fn, []).extend(r if isinstance(r, list) else [r])
    print("# Difficulty attacks, seeds %s, rules %s\n" % (seeds, rules))
    if "scen_hop" in by:
        table("1. Pool hopping: hopper's blocks per hash against the always-on base (adv = ratio - 1), %d h" % args.hop_hours,
              mean_rows(by["scen_hop"], ["rule", "x", "dwell"], worst=("adv",)),
              ["rule", "x", "dwell", "adv", "adv_worst", "hop_on_frac", "switches", "std", "worst_gap", "blocks_per_s", "seeds"])
    if "scen_pulse" in by:
        table("2. Pulsed hash rate: 50x for 10 min, every hour (period 3600) and once (period 1000000); weight_per_hash = block share / hash share",
              mean_rows(by["scen_pulse"], ["rule", "period_s"], worst=("pulse_vs_base", "weight_per_hash")),
              ["rule", "period_s", "pulse_vs_base", "pulse_vs_base_worst", "weight_per_hash", "weight_per_hash_worst", "pulse_vs_steady", "base_vs_steady",
               "pulse_block_share", "pulse_hash_share", "peak_blocks_per_min", "worst_gap", "blocks_per_s", "seeds"])
    if "scen_ts" in by:
        table("3. Timestamp stretching inside the rules: block rate after 1 h of forging (drift = rate - 1), 6 h runs",
              mean_rows(by["scen_ts"], ["rule", "ts_rules", "share", "stamp"], worst=("drift",)),
              ["rule", "ts_rules", "share", "stamp", "forger_offset_s", "rate_pre", "rate_post", "rate_last_hour", "drift", "drift_worst", "D_ratio", "std", "worst_gap", "forger_adv", "seeds"])
    if "scen_osc" in by:
        table("4. Short-lane oscillation: 25% miner on and off every 120 blocks (std_floor = the hash square wave at a fixed D)",
              mean_rows(by["scen_osc"], ["rule", "share", "period"], worst=("std_ratio",)),
              ["rule", "share", "period", "std_pre", "std_post", "std_ratio", "std_ratio_worst", "std_floor", "worst_gap_pre", "worst_gap_post", "blocks_per_s", "osc_adv", "seeds"])
    if "scen_epoch" in by:
        table("5. Epoch-boundary games: hold dodger (30%, off for blocks 0 to 7), hold flooder (10x, on for blocks 0 to 7), 30% flooder",
              mean_rows(by["scen_epoch"], ["rule", "kind"], worst=("gamer_adv",)),
              ["rule", "kind", "gamer_adv", "gamer_adv_worst", "gamer_on_frac", "std", "worst_gap", "blocks_per_s", "seeds"])
    if "scen_pollute" in by:
        table("6. Polluted window under adversarial timing: 10x miner joins at block join_k of an epoch and leaves at leave_k (600 = long lane takes over)",
              mean_rows(by["scen_pollute"], ["rule", "join_k", "leave_k"]),
              ["rule", "join_k", "leave_k", "join_first10_s", "join_settled_s", "leave_first10_s", "leave_settled_s", "leave_worst_gap_s", "leave_below05", "leave_above2x", "seeds"])
    if "base_profiles" in by:
        table("Base profiles (sim.py metrics) through each rule",
              mean_rows(by["base_profiles"], ["rule", "profile"]),
              ["rule", "profile", "first10_s", "settled_s", "settled_s_none", "worst_gap_s", "steady_std", "peak_rate", "not_recovered", "seeds"])


if __name__ == "__main__":
    main()