igneum/sim/difficulty/sim.py

#!/usr/bin/env python3
"""Igneum difficulty controller simulator.

Replays block production from a hash-rate profile through candidate difficulty controllers and
reports how each one tracks the target rate. Chain model: one selected chain, DAA score = height,
1 block per second target, block times exponential with mean (expected hashes per block) / (hash rate).
Timestamps are the true block times in milliseconds (option --ts-noise adds miner clock jitter).

Controllers
  kaspa      Kaspa's sampled DAA (KIP-4) as the fork runs it: samples every 4th block, 661 samples,
             average target of the window without its earliest block, times measured over expected,
             fixed genesis bits for the first 600 blocks (150 samples).
  monero     Monero style: last 720 blocks, timestamps sorted, 60 cut from each end, lag 15.
  lwma60     Zawy LWMA-1, N = 60, solvetimes clamped to [-FTL, 6T], FTL = 132 s.
  lwma120    the same with N = 120.
  igneum     the Igneum rule (two lanes, epoch windows, 3% harden and 10% ease clamps, 20 T solvetime cap on a
             sanitised running clock, warm-up from block 8, target floor 2^128).
  igneum-literal
             the brief's literal trigger (fast lane engaged while the short-window RATE is more than
             25% off target) kept for the chatter measurement.

Profiles
  real       the devnet record (--record CSV) turned into a step profile of estimated hash rate.
  up50, down50, epoch30, walk10, hop3, hop10, warmup-hard, warmup-easy

Usage
  python3 sim.py                              every synthetic profile, every controller, markdown out
  python3 sim.py --record devnet-2026-10-03.csv   adds the real profile and the exact-bits replay check
  python3 sim.py --tune                       the parameter sweep for the Igneum candidate
"""
import argparse
import csv
import math
import random
import sys
from collections import deque

T = 1000                       # target time per block, ms
EPOCH = 3600                   # blocks per program epoch
FTL = 132 * 1000               # Kaspa's future timestamp tolerance, ms
MAX_TARGET = (1 << 255) - 1
ONE = float(1 << 255)          # genesis-independent scale: target = ONE / D, D = expected hashes
MIN_TARGET = ONE / float(1 << 128)   # the Igneum rule's floor: 2^128 expected hashes per block (spec 2.3 rule 4)


# ---------------------------------------------------------------- controllers
class Controller:
    name = "base"

    def __init__(self, genesis_target):
        self.g = genesis_target
        self.ts = []       # timestamps ms, index = height
        self.tg = []       # targets (float), index = height

    def push(self, ts, target):
        self.ts.append(ts)
        self.tg.append(target)

    def next_target(self):
        raise NotImplementedError


def clamp(x, lo, hi):
    return lo if x < lo else hi if x > hi else x


class Kaspa(Controller):
    """KIP-4 sampled DAA as in rusty-kaspa `SampledDifficultyManager::calculate_difficulty_bits`.
    Sampling: block at height h is a sample when (h_parent_daa + 1) % rate == 0, i.e. h % 4 == 3 on a chain;
    genesis is never a sample. Window keeps the 661 highest samples. Retarget needs 150 samples."""
    name = "kaspa"

    def __init__(self, g, window=661, rate=4, min_samples=150):
        super().__init__(g)
        self.window = window
        self.rate = rate
        self.min_samples = min_samples
        self.samples = deque()   # (ts, target) of sampled blocks, oldest first

    def push(self, ts, target):
        h = len(self.ts)
        super().push(ts, target)
        if h > 0 and h % self.rate == self.rate - 1:
            self.samples.append((ts, target))
            if len(self.samples) > self.window:
                self.samples.popleft()

    def next_target(self):
        n = len(self.samples)
        if n < self.min_samples:
            return self.tg[-1] if self.tg else self.g
        min_i = min(range(n), key=lambda i: self.samples[i][0])
        max_i = max(range(n), key=lambda i: self.samples[i][0])
        min_ts = self.samples[min_i][0]
        max_ts = self.samples[max_i][0]
        tsum = sum(t for i, (_, t) in enumerate(self.samples) if i != min_i)
        avg = tsum / (n - 1)
        measured = max(max_ts - min_ts, 1)
        expected = T * self.rate * (n - 1)
        return min(avg * measured / expected, ONE)


class Monero(Controller):
    """Monero next_difficulty: N 720, cut 60 each side on sorted timestamps, lag 15.
    Difficulty units, converted to target at the end."""
    name = "monero"

    def __init__(self, g, n=720, cut=60, lag=15):
        super().__init__(g)
        self.n, self.cut, self.lag = n, cut, lag

    def next_target(self):
        h = len(self.ts)
        end = h - self.lag
        if end < 2:
            return self.g
        start = max(0, end - self.n)
        ts = sorted(self.ts[start:end])
        length = len(ts)
        if length <= 2 * self.cut + 2:
            cb, ce = 0, length
        else:
            cb = (length - (self.n - 2 * self.cut) + 1) // 2 if length >= self.n else self.cut
            cb = clamp(cb, 0, self.cut)
            ce = length - cb
            if ce - cb < 2:
                cb, ce = 0, length
        span = max(ts[ce - 1] - ts[cb], 1)
        work = sum(ONE / t for t in self.tg[start + cb:start + ce])
        d = work * T / span
        return min(ONE / d, ONE)


class Lwma(Controller):
    """Zawy LWMA-1 (2018 reference): linear weights, solvetime in [-FTL, 6T], average target over N.
    Grows from 4 blocks instead of waiting for N."""
    name = "lwma"

    def __init__(self, g, n=120, lo=-FTL, hi=6 * T, min_blocks=4):
        super().__init__(g)
        self.n, self.lo, self.hi, self.min_blocks = n, lo, hi, min_blocks
        self.name = "lwma%d" % n

    def next_target(self):
        h = len(self.ts)
        if h < self.min_blocks + 1:
            return self.g
        n = min(self.n, h - 1)
        ts = self.ts
        wsum = 0.0
        tsum = 0.0
        for i in range(1, n + 1):
            j = h - n - 1 + i
            st = clamp(ts[j] - ts[j - 1], self.lo, self.hi)
            wsum += st * i
            tsum += self.tg[j]
        wsum = max(wsum, n * n * T / 20)
        avg = tsum / n
        return min(avg * wsum / (T * n * (n + 1) / 2), ONE)


class Igneum(Controller):
    """Igneum dual-lane controller (docs/spec/02-consensus.md section 2.3).

    Every lane estimates the hash rate as work over time (what Kaspa's estimateNetworkHashesPerSecond does),
    never as average target times average solvetime, which is biased while targets ramp.
    k = blocks of the current epoch already in the past (epoch = height // 3600).
    Every chain step is measured on a sanitised running clock (4 October 2026, after the timestamp-stretching
    attack of sim/difficulty/attacks/): c(genesis) = t(genesis); c(b) = max(c(p) + clamp(t(b) - c(p), -CAP, +CAP),
    t(b) - LAG_CAPS x CAP); step(b) = min(c(b) - c(p), CAP). The steps of a lane sum to the clock's span over it, so
    a forged stamp is paid back by the honest blocks after it instead of being cancelled to zero measured time.
    S (short): linear-weighted work over the clock steps of the last min(NS, k) chain blocks.
    E (epoch): work over the clock steps of the whole epoch so far, shrunk toward the parent's implied rate by
               k / (k + K0); used as the reference while k < LONG_MIN.
    L (long):  blue work between the earliest and latest sample of Kaspa's sampled window restricted to the
               epoch, over their time span; the reference once k >= LONG_MIN.
    Candidate = S when S and the reference disagree by more than TRIG, else the reference. The output may
    fall (harder) by at most CLAMP and rise (easier) by at most EASE per block, then is bounded to
    [MIN_TARGET, ONE]. The parent's target is held for the first K_MIN blocks of an epoch. Genesis bits only for
    blocks 1 to K_MIN."""
    name = "igneum"

    def __init__(self, g, ns=120, k_min=8, k0=16, long_min=600, trig=0.25, clamp_pct=0.03, ease_pct=0.10, cap=20 * T,
                 window=661, rate=4, literal=False, epoch_aware=True, kaspa_formula=False, lag_caps=3):
        super().__init__(g)
        self.ns, self.k_min, self.k0, self.long_min = ns, k_min, k0, long_min
        self.trig, self.clamp_pct, self.cap = trig, clamp_pct, cap
        self.lag_caps = lag_caps
        self.tss = []           # sanitised clock per height
        self.sts = []           # clock step per height (the solvetime the lanes see)
        self.ease_pct = clamp_pct if ease_pct is None else ease_pct
        self.window, self.rate = window, rate
        self.literal = literal
        self.epoch_aware = epoch_aware
        self.kaspa_formula = kaspa_formula
        self.samples = deque()  # (height, ts, blue_work)
        self.pref_d = [0.0]     # prefix sums of D (work), index h+1 = work through block h = blue work
        self.pref_t = [0.0]     # prefix sums of targets (only for the kaspa_formula variant)
        self.pref_st = [0.0]    # prefix sums of capped solvetimes
        self.engaged = 0

    def push(self, ts, target):
        h = len(self.ts)
        if h == 0:
            st, tss = 0, ts
        else:
            # the sanitised clock: one capped step from the parent's clock, never more than LAG_CAPS caps behind
            # the raw stamp (an honest idle gap is not paid back as cap-sized steps for many blocks)
            tss = self.tss[-1] + clamp(ts - self.tss[-1], -self.cap, self.cap)
            tss = max(tss, ts - self.lag_caps * self.cap)
            st = min(tss - self.tss[-1], self.cap)
        self.tss.append(tss)
        self.sts.append(st)
        super().push(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]))   # the long lane keeps the raw sampled span
            if len(self.samples) > self.window:
                self.samples.popleft()

    def rate_s(self, lo, hi):
        """Linear-weighted hash rate over chain steps lo+1..hi-1 (work of block j over the clock step of block j)."""
        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
            wsum += self.sts[j] * i
            dsum += (self.pref_d[j + 1] - self.pref_d[j]) * i
        wsum = max(wsum, n * n * T / 20)
        return dsum / wsum

    def rate_e(self, lo, hi, k):
        n = hi - lo - 1
        if n < 1:
            return None
        work = self.pref_d[hi] - self.pref_d[lo + 1]
        measured = max(self.pref_st[hi] - self.pref_st[lo + 1], 1)
        r = work / measured
        prior = (self.pref_d[hi] - self.pref_d[hi - 1]) / T   # the parent's implied rate
        return (k * r + self.k0 * prior) / (k + self.k0)

    def rate_l(self, lo):
        samp = [s for s in self.samples if s[0] >= lo]
        n = len(samp)
        if n < self.long_min // self.rate:
            return None
        mn = min(samp, key=lambda s: s[1])
        mx = max(samp, key=lambda s: s[1])
        measured = max(mx[1] - mn[1], 1)
        if self.kaspa_formula:
            tsum = sum(self.tg[s[0]] for s in samp if s is not mn)
            return 1.0 / ((tsum / (n - 1)) * measured / (T * self.rate * (n - 1))) / ONE * ONE * (1.0 / T) * T if False else None
        return (mx[2] - mn[2]) / measured

    def next_target(self):
        h = len(self.ts)
        if h == 0:
            return self.g
        parent = self.tg[-1]
        epoch_start = (h // EPOCH) * EPOCH if self.epoch_aware else 0
        k = h - epoch_start
        if k < self.k_min:
            cand = parent
        else:
            s_lo = max(epoch_start, h - self.ns)
            r_s = self.rate_s(s_lo, h)
            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, h, k)
            if self.literal:
                n = h - s_lo
                measured = max(self.pref_st[h] - self.pref_st[s_lo + 1], 1)
                use_s = abs(T * (n - 1) / measured - 1) > self.trig
            else:
                use_s = abs(r_s / r_ref - 1) > self.trig
            r = r_s if use_s else r_ref
            self.engaged += 1 if use_s else 0
            cand = ONE / (r * T)
        out = clamp(cand, parent / (1 + self.clamp_pct), parent * (1 + self.ease_pct))
        return clamp(out, MIN_TARGET, ONE)


def make_controller(name, g, **kw):
    if name == "kaspa":
        return Kaspa(g)
    if name == "monero":
        return Monero(g)
    if name == "lwma60":
        return Lwma(g, n=60)
    if name == "lwma120":
        return Lwma(g, n=120)
    if name == "igneum":
        return Igneum(g, **kw)
    if name == "igneum-literal":
        return Igneum(g, literal=True, **kw)
    if name == "igneum-noepoch":
        return Igneum(g, epoch_aware=False, **kw)
    raise ValueError(name)


# ---------------------------------------------------------------- profiles
class Profile:
    """Hash rate in hashes per ms as a function of (time ms, height). `steps` lists (t_ms, factor) events
    for the recovery metrics; `epoch_steps` are height-keyed factors applied at epoch boundaries."""

    def __init__(self, name, base, duration_ms, genesis_d=None, seed=7):
        self.name = name
        self.base = base
        self.duration = duration_ms
        self.time_steps = []    # (t_ms, factor)
        self.epoch_steps = {}   # height -> factor
        self.walk = None        # (sigma per minute)
        self.genesis_d = genesis_d if genesis_d is not None else base * T
        self.rng = random.Random(seed)
        self._walk_cache = {}

    def rate(self, t, h):
        f = 1.0
        for ts, fac in self.time_steps:
            if t >= ts:
                f *= fac
        for hh, fac in self.epoch_steps.items():
            if h >= hh:
                f *= fac
        if self.walk:
            m = int(t // 60000)
            if m not in self._walk_cache:
                prev = self._walk_cache.get(m - 1, 0.0) if m > 0 else 0.0
                self._walk_cache[m] = prev + self.rng.gauss(0, self.walk)
            f *= math.exp(self._walk_cache[m])
        return self.base * f

    def events(self):
        """Step events as (t_ms or None, height or None, factor)."""
        out = [(t, None, f) for t, f in self.time_steps]
        out += [(None, h, f) for h, f in self.epoch_steps.items()]
        return out


SCALE = float(1 << 28)   # synthetic difficulty unit: 2^28 hashes per block, the devnet genesis
H1 = SCALE / T           # hash rate that gives one block per second at D = SCALE


def synthetic_profiles(seed):
    P = []
    hour = 3600 * 1000
    p = Profile("up50", H1, 5 * hour, seed=seed); p.time_steps = [(3 * hour, 50.0)]; P.append(p)
    p = Profile("down50", 50 * H1, 7 * hour, seed=seed); p.time_steps = [(3 * hour, 1 / 50.0)]; P.append(p)
    p = Profile("epoch30", H1, 8 * hour, seed=seed)
    rng = random.Random(seed)
    fac = 1.0
    for e in range(2, 8):
        s = rng.choice([1.3, 1 / 1.3])
        p.epoch_steps[e * EPOCH] = s
    P.append(p)
    p = Profile("walk10", H1, 6 * hour, seed=seed); p.walk = 0.10 / math.sqrt(60); P.append(p)
    p = Profile("hop3", H1, 6 * hour, seed=seed)
    p.time_steps = [(3 * hour + i * 15 * 60000, 3.0 if i % 2 == 0 else 1 / 3.0) for i in range(12)]
    P.append(p)
    p = Profile("hop10", H1, 6 * hour, seed=seed)
    p.time_steps = [(3 * hour + i * 15 * 60000, 10.0 if i % 2 == 0 else 1 / 10.0) for i in range(12)]
    P.append(p)
    # the devnet of 3 Oct 2026: genesis sized for 0.6 blocks/s of GPU, a 0.11 blocks/s Metal worker for 21 min,
    # 13 min of nothing, then the PC; Kaspa's first retarget at block 600 sees the whole slow history
    p = Profile("polluted", 0.6 * H1, 5 * hour, genesis_d=1.0 * SCALE, seed=seed)
    p.time_steps = [(0, 0.11 / 0.6), (21 * 60000, 1e-9), (34 * 60000, 0.6 / 0.11 / 1e-9)]
    P.append(p)
    p = Profile("warmup-hard", H1, 2 * hour, genesis_d=10.0 * SCALE, seed=seed); p.time_steps = [(0, 1.0)]; P.append(p)
    p = Profile("warmup-easy", H1, 2 * hour, genesis_d=0.1 * SCALE, seed=seed); p.time_steps = [(0, 1.0)]; P.append(p)
    return P


# ---------------------------------------------------------------- record (real devnet)
def load_record(path):
    """Returns (chain, all): chain blocks (one per blue score, for the exact replay) and every block
    (for the hash-rate profile), as (blue_score, daa_score, timestamp_ms, bits)."""
    rows = list(csv.DictReader(open(path)))
    rows.sort(key=lambda r: (int(r["blue_score"]), int(r["daa_score"]), int(r["timestamp_ms"])))
    every = [(int(r["blue_score"]), int(r["daa_score"]), int(r["timestamp_ms"]), int(r["bits"], 16)) for r in rows]
    chain = []
    seen = set()
    for r in rows:
        if r["is_chain_block"] != "1" and r["blue_score"] != "0":
            continue
        bs = int(r["blue_score"])
        if bs in seen:
            continue
        seen.add(bs)
        chain.append((bs, int(r["daa_score"]), int(r["timestamp_ms"]), int(r["bits"], 16)))
    return chain, every


def bits_to_target(bits):
    exp = bits >> 24
    mant = bits & 0x7FFFFF
    if exp <= 3:
        return mant >> (8 * (3 - exp))
    return mant << (8 * (exp - 3))


def target_to_bits(t):
    """Uint256::compact_target_bits: size in bytes, mantissa top 3 bytes, sign-bit carry."""
    if t == 0:
        return 0
    size = (t.bit_length() + 7) // 8
    if size <= 3:
        mant = t << (8 * (3 - size))
    else:
        mant = t >> (8 * (size - 3))
    if mant & 0x800000:
        mant >>= 8
        size += 1
    return (size << 24) | mant


def kaspa_bits_exact(record, window=661, rate=4, min_samples=150):
    """Replay the record's own timestamps and bits through rusty-kaspa's integer arithmetic; return
    (checked, mismatches). Sample rule as in `sampled_mergeset_iterator` on a chain."""
    samples = deque()
    checked = mism = 0
    first_bad = None
    exact_before_first = 0
    for i, (bs, daa, ts, bits) in enumerate(record):
        if i > 0:
            n = len(samples)
            if n >= min_samples:
                min_i = min(range(n), key=lambda j: samples[j][0])
                max_i = max(range(n), key=lambda j: samples[j][0])
                tsum = sum(t for j, (_, t) in enumerate(samples) if j != min_i)
                avg = tsum // (n - 1)
                measured = max(samples[max_i][0] - samples[min_i][0], 1)
                expected = T * rate * (n - 1)
                new = min(avg * measured // expected, MAX_TARGET)
                want = target_to_bits(new)
                checked += 1
                if want != bits:
                    mism += 1
                    if first_bad is None:
                        first_bad = (daa, hex(bits), hex(want))
                elif first_bad is None:
                    exact_before_first += 1
        parent_daa = record[i - 1][1] if i > 0 else None
        if i > 0 and (parent_daa + 1) % rate == 0:
            samples.append((ts, bits_to_target(bits)))
            if len(samples) > window:
                samples.popleft()
    return checked, mism, first_bad, exact_before_first


def record_profile(record, every, events, seed):
    """Piecewise-constant hash rate estimated from the record between the given event times.
    events: sorted list of (t_ms, label). Rate in a segment = sum of expected hashes / duration."""
    t0 = record[1][2]
    tend = record[-1][2]
    bounds = [t0] + [t for t, _ in events if t0 < t < tend] + [tend]
    segs = []
    for a, b in zip(bounds, bounds[1:]):
        hashes = sum(2 * ONE / bits_to_target(bits) for (_, _, ts, bits) in every[1:] if a < ts <= b)
        dur = max(b - a, 1)
        segs.append((a, b, hashes / dur))
    p = Profile("real", segs[0][2], tend - t0, genesis_d=ONE / bits_to_target(record[0][3]), seed=seed)
    for (a, b, r), (a2, b2, r2) in zip(segs, segs[1:]):
        p.time_steps.append((a2 - t0, r2 / r if r > 0 else 1.0))
    p.segments = [(a - t0, b - t0, r) for a, b, r in segs]
    p.min_rate = 1e-12 * H1
    return p


# ---------------------------------------------------------------- simulation
def simulate(profile, ctrl, rng, ts_noise=0, max_blocks=400000):
    """Returns list of (t_ms, height, D, hashrate, r) per block, r = expected block rate / target."""
    t = 0
    out = []
    target = ONE / profile.genesis_d
    ctrl.push(0, target)                   # genesis at t = 0
    h = 1
    while t < profile.duration and h < max_blocks:
        target = ctrl.next_target()
        D = ONE / target
        work = rng.expovariate(1.0) * D      # hashes needed for this block
        t_start = t
        tt = t
        while True:
            hr = profile.rate(tt, h)
            nxt = min((st for st, _ in profile.time_steps if st > tt), default=None)
            if hr > 0 and (nxt is None or work <= hr * (nxt - tt)):
                tt += work / hr
                break
            if nxt is None:                      # zero hash rate for ever: end the run
                tt = profile.duration + 1
                break
            work -= hr * (nxt - tt)
            tt = nxt
        gap = tt - t_start
        hr = profile.rate(tt, h)
        gap = max(gap, 1.0)
        t += gap
        ts = int(t) + (rng.randint(-ts_noise, ts_noise) if ts_noise else 0)
        ctrl.push(ts, target)
        r = hr * T / D
        out.append((t, h, D, hr, r))
        h += 1
    return out


def smoothed(xs, half):
    """Centred moving mean over 2*half+1 values (shorter at the ends)."""
    n = len(xs)
    pref = [0.0]
    for x in xs:
        pref.append(pref[-1] + x)
    out = []
    for i in range(n):
        lo = max(0, i - half)
        hi = min(n, i + half + 1)
        out.append((pref[hi] - pref[lo]) / (hi - lo))
    return out


def metrics(blocks, profile, t_step=None, h_step=None, settle_blocks=100, t_end=None):
    """Recovery after one step, overshoot, steady state and gap statistics."""
    res = {}
    if t_step is not None or h_step is not None:
        if t_step is None:
            t_step = next((b[0] for b in blocks if b[1] >= h_step), None)
        post = [b for b in blocks if b[0] >= t_step and (t_end is None or b[0] < t_end)] if t_step is not None else []
        if not post:
            res["recover_s"] = None
            res["overshoot"] = None
        else:
            rec = None
            sm = smoothed([b[4] for b in post], 60)
            first = next((i for i, x in enumerate(sm) if abs(x - 1) <= 0.10), None)
            res["first10_s"] = None if first is None else (post[first][0] - t_step) / 1000
            for i in range(len(post)):
                seg = sm[i:i + settle_blocks]
                if len(seg) == settle_blocks and all(abs(x - 1) <= 0.10 for x in seg):
                    rec = post[i][0] - t_step
                    break
            res["recover_s"] = None if rec is None else rec / 1000
            res["recover_blocks"] = None if rec is None else sum(1 for b in post if b[0] - t_step <= rec)
            r0 = post[0][4]
            side = 1 if r0 > 1 else -1
            cross = next((i for i, x in enumerate(sm) if (x - 1) * side <= 0), None)
            if cross is None:
                res["overshoot"] = 0.0
            else:
                tail = sm[cross:cross + 3600]
                res["overshoot"] = max(0.0, max((1 - x) * side for x in tail))
            res["above2x"] = sum(1 for b in post if b[4] > 2)
            res["below05"] = sum(1 for b in post if b[4] < 0.5)
            res["worst_gap_s"] = max((post[i][0] - post[i - 1][0]) for i in range(1, len(post))) / 1000 if len(post) > 1 else None
    return res


def steady(blocks, t_from, t_to):
    seg = [b for b in blocks if t_from <= b[0] < t_to]
    if len(seg) < 10:
        return None, None, None
    rs = [b[4] for b in seg]
    mean = sum(rs) / len(rs)
    std = math.sqrt(sum((x - mean) ** 2 for x in rs) / len(rs))
    # blocks per minute
    mins = {}
    for b in seg:
        mins[int(b[0] // 60000)] = mins.get(int(b[0] // 60000), 0) + 1
    counts = list(mins.values())[1:-1]
    cm = sum(counts) / len(counts) if counts else 0
    cs = math.sqrt(sum((c - cm) ** 2 for c in counts) / len(counts)) if counts else 0
    gaps = [seg[i][0] - seg[i - 1][0] for i in range(1, len(seg))]
    return std, (cs / cm if cm else None), max(gaps) / 1000


def fmt(x, nd=1):
    if x is None:
        return "none"
    if isinstance(x, float):
        return ("%." + str(nd) + "f") % x
    return str(x)


def run_profile(profile, names, seed, ts_noise, kw=None):
    rows = []
    hour = 3600 * 1000
    for name in names:
        rng = random.Random(seed)
        ctrl = make_controller(name, ONE / profile.genesis_d, **(kw or {}))
        blocks = simulate(profile, ctrl, rng, ts_noise=ts_noise)
        row = {"ctrl": name, "blocks": len(blocks)}
        ev = profile.events()
        if profile.name.startswith("warmup"):
            m = metrics(blocks, profile, t_step=0)
            row.update(m)
            s = steady(blocks, 1 * hour, 2 * hour)
        elif profile.name == "epoch30":
            recs, ovs, a2, b5, wg = [], [], 0, 0, 0
            firsts = []
            for h, f in sorted(profile.epoch_steps.items()):
                t_end = next((b[0] for b in blocks if b[1] >= h + EPOCH), None)
                m = metrics(blocks, profile, h_step=h, t_end=t_end)
                recs.append(m.get("recover_s"))
                firsts.append(m.get("first10_s"))
                ovs.append(m.get("overshoot") or 0)
                a2 += m.get("above2x", 0)
                b5 += m.get("below05", 0)
            m = metrics(blocks, profile, h_step=2 * EPOCH)
            good = [r for r in recs if r is not None]
            gf = [r for r in firsts if r is not None]
            row["first10_s"] = (sum(gf) / len(gf)) if gf else None
            row["recover_s"] = (sum(good) / len(good)) if good else None
            row["recover_max_s"] = max(good) if good else None
            row["recover_fail"] = len(recs) - len(good)
            row["overshoot"] = max(ovs)
            row["above2x"] = sum(1 for b in blocks if b[4] > 2 and b[1] >= 2 * EPOCH)
            row["below05"] = sum(1 for b in blocks if b[4] < 0.5 and b[1] >= 2 * EPOCH)
            row["worst_gap_s"] = max(blocks[i][0] - blocks[i - 1][0] for i in range(1, len(blocks))) / 1000
            s = steady(blocks, 1 * hour, 2 * hour)
        elif profile.name.startswith("hop"):
            recs = []
            ovs = []
            firsts = []
            for i, (t, f) in enumerate(profile.time_steps):
                t_end = profile.time_steps[i + 1][0] if i + 1 < len(profile.time_steps) else None
                m = metrics(blocks, profile, t_step=t, t_end=t_end)
                recs.append(m.get("recover_s"))
                firsts.append(m.get("first10_s"))
                ovs.append(m.get("overshoot") or 0)
            good = [r for r in recs if r is not None]
            gf = [r for r in firsts if r is not None]
            row["first10_s"] = (sum(gf) / len(gf)) if gf else None
            row["recover_s"] = (sum(good) / len(good)) if good else None
            row["recover_max_s"] = max(good) if good else None
            row["recover_fail"] = len(recs) - len(good)
            row["overshoot"] = max(ovs)
            t0 = profile.time_steps[0][0]
            row["above2x"] = sum(1 for b in blocks if b[4] > 2 and b[0] >= t0)
            row["below05"] = sum(1 for b in blocks if b[4] < 0.5 and b[0] >= t0)
            row["worst_gap_s"] = max(blocks[i][0] - blocks[i - 1][0] for i in range(1, len(blocks)) if blocks[i][0] >= t0) / 1000
            s = steady(blocks, 1 * hour, 3 * hour)
        elif profile.name == "walk10":
            row["above2x"] = sum(1 for b in blocks if b[4] > 2 and b[0] >= hour)
            row["below05"] = sum(1 for b in blocks if b[4] < 0.5 and b[0] >= hour)
            row["worst_gap_s"] = max(blocks[i][0] - blocks[i - 1][0] for i in range(1, len(blocks)) if blocks[i][0] >= hour) / 1000
            s = steady(blocks, 1 * hour, 6 * hour)
        elif profile.name == "polluted":
            t_pc = profile.time_steps[2][0]
            m = metrics(blocks, profile, t_step=t_pc)
            row.update(m)
            post = [b for b in blocks if b[0] >= t_pc]
            row["peak_rate"] = max(b[4] for b in post)
            row["trough_rate"] = min(b[4] for b in post if b[0] >= t_pc + 600000) if any(b[0] >= t_pc + 600000 for b in post) else None
            s = steady(blocks, 3 * hour, 5 * hour)
        elif getattr(profile, "custom", False):
            recs, firsts, ovs = [], [], []
            for i, (t, f) in enumerate(profile.time_steps):
                t_end = profile.time_steps[i + 1][0] if i + 1 < len(profile.time_steps) else None
                m = metrics(blocks, profile, t_step=t, t_end=t_end)
                row["step%d_first" % (i + 1)] = m.get("first10_s")
                row["step%d_settled" % (i + 1)] = m.get("recover_s")
                row["step%d_gap" % (i + 1)] = m.get("worst_gap_s")
            row["blocks_within10"] = sum(1 for b in blocks if abs(b[4] - 1) <= 0.1)
            s = (None, None, None)
        elif profile.name == "real":
            segs = profile.segments
            # the biggest step in the record
            big = max(profile.time_steps, key=lambda e: abs(math.log(e[1])))
            m = metrics(blocks, profile, t_step=big[0])
            row.update(m)
            row["step"] = big[1]
            s = steady(blocks, segs[-1][0] + 600000, segs[-1][1]) if segs[-1][1] - segs[-1][0] > 1200000 else (None, None, None)
            row["blocks_within10"] = sum(1 for b in blocks if abs(b[4] - 1) <= 0.1)
        else:
            t, f = profile.time_steps[0]
            m = metrics(blocks, profile, t_step=t)
            row.update(m)
            s = steady(blocks, 1 * hour, 3 * hour)
        row["ss_std"], row["ss_bpm_cv"], row["ss_worst_gap_s"] = s
        if isinstance(ctrl, Igneum):
            row["engaged"] = ctrl.engaged
        rows.append(row)
    return rows


COLS = [("ctrl", "controller"), ("first10_s", "first within 10% s"), ("recover_s", "settled s"), ("recover_blocks", "recovery blocks"),
        ("recover_max_s", "worst recovery s"), ("recover_fail", "not recovered"), ("overshoot", "overshoot"),
        ("ss_std", "steady std of rate"), ("ss_bpm_cv", "steady blocks/min CV"), ("worst_gap_s", "worst gap s"),
        ("above2x", "blocks above 2x"), ("below05", "blocks below 0.5x"), ("blocks", "blocks"), ("engaged", "fast-lane blocks"),
        ("blocks_within10", "blocks within 10%"), ("step", "step factor"), ("peak_rate", "peak rate"), ("trough_rate", "trough rate after 10 min")]


def print_table(title, rows):
    cols = [(k, h) for k, h in COLS if any(k in r for r in rows)]
    print("\n### " + title + "\n")
    print("| " + " | ".join(h for _, h in cols) + " |")
    print("|" + "---|" * len(cols))
    for r in rows:
        print("| " + " | ".join(fmt(r.get(k), 3 if k in ("ss_std", "ss_bpm_cv", "overshoot") else 1) for k, _ in cols) + " |")


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--record")
    ap.add_argument("--seed", type=int, default=7)
    ap.add_argument("--ts-noise", type=int, default=0, help="uniform timestamp jitter, ms")
    ap.add_argument("--profiles", default="all")
    ap.add_argument("--controllers", default="kaspa,monero,lwma60,lwma120,igneum,igneum-literal")
    ap.add_argument("--tune", action="store_true")
    ap.add_argument("--events", default="", help="real record: comma list of ISO times of miner events")
    ap.add_argument("--dump", default="", help="profile:controller, write per-block trajectory to stdout as CSV and exit")
    ap.add_argument("--custom", default="", help="extra profile: name=genesis_bits_hex,base_MHs,duration_s,t_s:factor,t_s:factor,...")
    ap.add_argument("--quick", action="store_true",
                    help="smoke run for CI: profiles up50 and warmup-hard, controllers kaspa and igneum (about a minute); --profiles and --controllers still override")
    args = ap.parse_args()
    if args.quick:
        if args.profiles == "all":
            args.profiles = "up50,warmup-hard"
        if args.controllers == "kaspa,monero,lwma60,lwma120,igneum,igneum-literal":
            args.controllers = "kaspa,igneum"
    names = args.controllers.split(",")
    profiles = synthetic_profiles(args.seed)
    if args.profiles != "all":
        want = args.profiles.split(",")
        profiles = [p for p in profiles if p.name in want]
    if args.custom:
        name, rest = args.custom.split("=")
        parts = rest.split(",")
        gbits, mhs, dur = int(parts[0], 16), float(parts[1]), float(parts[2])
        p = Profile(name, mhs * 1e6 / 1000, dur * 1000, genesis_d=ONE / bits_to_target(gbits), seed=args.seed)
        p.time_steps = [(float(x.split(":")[0]) * 1000, float(x.split(":")[1])) for x in parts[3:]]
        p.custom = True
        profiles.append(p)
    if args.dump:
        pn, cn = args.dump.split(":")
        p = next(q for q in synthetic_profiles(args.seed) if q.name == pn)
        ctrl = make_controller(cn, ONE / p.genesis_d)
        for t, h, D, hr, r in simulate(p, ctrl, random.Random(args.seed)):
            print("%d,%d,%.6g,%.6g,%.4f" % (t, h, D, hr, r))
        return
    print("# Difficulty controller simulation, seed %d, ts noise %d ms\n" % (args.seed, args.ts_noise))
    if args.tune:
        p_names = ["up50", "down50", "epoch30", "hop10", "walk10", "polluted"]
        base = {}
        grid = [("ns", [60, 120, 180]), ("trig", [0.15, 0.25, 0.35]), ("clamp_pct", [0.02, 0.03, 0.05]),
                ("ease_pct", [0.03, 0.06, 0.10, 0.15]),
                ("k0", [8, 16, 32]), ("cap", [6 * T, 20 * T, FTL])]
        for key, vals in grid:
            rows = []
            for v in vals:
                kw = dict(base)
                kw[key] = v
                agg = {"ctrl": "%s=%s" % (key, v if key != "cap" else v // 1000)}
                for p in [q for q in synthetic_profiles(args.seed) if q.name in p_names]:
                    r = run_profile(p, ["igneum"], args.seed, args.ts_noise, kw)[0]
                    agg[p.name + " rec"] = r.get("recover_s")
                    agg[p.name + " first"] = r.get("first10_s")
                    agg[p.name + " std"] = r.get("ss_std")
                    if p.name == "down50":
                        agg["down50 gap"] = r.get("worst_gap_s")
                    if p.name == "polluted":
                        agg["polluted peak"] = r.get("peak_rate")
                    if p.name == "walk10":
                        agg["walk10 >2x/<0.5x"] = "%d/%d" % (r.get("above2x", 0), r.get("below05", 0))
                rows.append(agg)
            keys = list(rows[0].keys())
            print("\n### Tune %s\n" % key)
            print("| " + " | ".join(keys) + " |")
            print("|" + "---|" * len(keys))
            for r in rows:
                print("| " + " | ".join(fmt(r[k], 3 if k.endswith("std") else 1) for k in keys) + " |")
        return
    for p in profiles:
        rows = run_profile(p, names, args.seed, args.ts_noise)
        if getattr(p, "custom", False):
            continue
        desc = {"up50": "hash rate x50 at 3 h", "down50": "hash rate /50 at 3 h",
                "epoch30": "x1.3 or /1.3 at each epoch boundary from epoch 2", "walk10": "10% per hour log random walk",
                "hop3": "3x pool hops in and out every 15 min from 3 h", "hop10": "10x pool hops in and out every 15 min from 3 h",
                "polluted": "window polluted by a pre-step slow period (21 min at 0.11 blocks/s, 13 min idle, then 0.6 blocks/s at genesis difficulty)", "warmup-hard": "genesis difficulty 10x too hard", "warmup-easy": "genesis difficulty 10x too easy"}[p.name]
        print_table("%s: %s" % (p.name, desc), rows)
    for p in profiles:
        if getattr(p, "custom", False):
            rows = run_profile(p, names, args.seed, args.ts_noise)
            keys = [k for k in rows[0].keys() if k != "ss_std" and k != "ss_bpm_cv" and k != "ss_worst_gap_s"]
            print("\n### %s: custom profile\n" % p.name)
            print("| " + " | ".join(keys) + " |")
            print("|" + "---|" * len(keys))
            for r in rows:
                print("| " + " | ".join(fmt(r.get(k)) for k in keys) + " |")
    if args.record:
        rec, every = load_record(args.record)
        checked, mism, first_bad, exact = kaspa_bits_exact(rec)
        print("\n### Exact replay of Kaspa's rule on the record\n")
        print("Blocks %d, chain blocks %d, retargets checked %d, exact before the first mismatch %d, mismatches %d%s" % (
            len(every), len(rec), checked, exact, mism, "" if first_bad is None else ", first at daa %d: have %s want %s" % first_bad))
        print("(A chain-only replay cannot see the merged blocks that Kaspa's sampled window includes once the DAG widens.)")
        import datetime
        ev = []
        for s in [x for x in args.events.split(",") if x]:
            ev.append((int(datetime.datetime.fromisoformat(s).timestamp() * 1000), s))
        p = record_profile(rec, every, sorted(ev), args.seed)
        print("\nReal profile segments (hash rate estimated as expected hashes over wall time):\n")
        print("| from s | to s | blocks/s at genesis difficulty | MH/s |")
        print("|---|---|---|---|")
        g = p.genesis_d
        for a, b, r in p.segments:
            print("| %d | %d | %.3f | %.1f |" % (a / 1000, b / 1000, r * T / (2 * g), r * 1000 / 1e6))
        rows = run_profile(p, names, args.seed, args.ts_noise)
        print_table("real: the devnet record replayed as a hash-rate profile", rows)


if __name__ == "__main__":
    main()