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.
  igneum-v2  rule v2 (4 October 2026): the reference lane is the newest 600 blocks of the epoch (epoch lane only).
  igneum:key=value+key=value
             any variant of the Igneum class (ns, trig, trig_exit, soft=a/b, median3, clamp_window, ref_window,
             long_min, clamp_pct, ease_pct, cap, k0, ...).

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
  python3 sim.py --live records/live-2026-10-04.csv --hashrate records/live-2026-10-04-hashrate.csv \
                 --controllers igneum,igneum-v2   DAG replay of the live devnet against the record (4 October 2026)
"""
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,
                 trig_exit=None, soft=None, median3=False, clamp_window=None, ref_window=None):
        """Candidates of 4 October 2026 (live oscillation, docs/analysis/difficulty-2026-10-04-oscillation.md), all off by
        default so the rule of spec 2.3 is unchanged:
        trig_exit     hysteresis: the short lane, once engaged, stays engaged until it is within this fraction of the reference
        soft          (a, b): blend instead of switch; weight on the short lane rises linearly from 0 at |S/R - 1| = a to 1 at b
        median3       the short lane's rate is the median of its estimates at the block, its parent and its grandparent
        clamp_window  ms: the clamps also bind relative to the target of the newest chain block at least this far back on
                      the sanitised clock, so a burst of chain blocks inside one window moves the target once
        ref_window    blocks: the reference lane sees at most this many newest blocks of the epoch (the long lane's samples
                      and the epoch lane's steps alike), so a hash-rate step inside an epoch washes out of it in that many
                      blocks instead of lasting the whole epoch"""
        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.trig_exit, self.soft, self.median3 = trig_exit, soft, median3
        self.clamp_window, self.ref_window = clamp_window, ref_window
        self.lane = [False]     # per height: the short lane ruled the block's target (index 0 = genesis)
        self.flips = 0
        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
        eng = False
        if k < self.k_min:
            cand = parent
        else:
            ref_start = max(epoch_start, h - self.ref_window) if self.ref_window else epoch_start
            s_lo = max(ref_start, h - self.ns)
            r_s = self.rate_s(s_lo, h)
            if self.median3 and h - s_lo > 4:
                r_s = sorted([r_s, self.rate_s(s_lo, h - 1), self.rate_s(s_lo, h - 2)])[1]
            r_ref = None
            if k >= self.long_min:
                r_ref = self.rate_l(max(ref_start, h - self.window * self.rate))
            if r_ref is None:
                r_ref = self.rate_e(ref_start, h, h - ref_start)
            if self.literal:
                n = h - s_lo
                measured = max(self.pref_st[h] - self.pref_st[s_lo + 1], 1)
                eng = abs(T * (n - 1) / measured - 1) > self.trig
                r = r_s if eng else r_ref
            else:
                r, eng = choose_rate(r_s, r_ref, self.lane[-1], self.trig, self.trig_exit, self.soft)
            self.engaged += 1 if eng else 0
            cand = ONE / (r * T)
        base = None
        if self.clamp_window:
            j = h - 1
            while j > 0 and self.tss[h - 1] - self.tss[j] < self.clamp_window:
                j -= 1
            base = self.tg[j]
        if eng != self.lane[-1]:
            self.flips += 1
        self.lane.append(eng)
        return clamp_target(cand, parent, base, self.clamp_pct, self.ease_pct)


def choose_rate(r_s, r_ref, prev_engaged, trig, trig_exit=None, soft=None):
    """The trigger between the lanes. Returns (rate, short lane engaged).
    Spec 2.3 rule 3: the short lane rules while it differs from the reference by more than `trig`.
    trig_exit: hysteresis, an engaged short lane rules until it is within trig_exit of the reference.
    soft (a, b): the output is the reference plus a share of the difference that rises linearly from 0 at a to 1 at b."""
    if r_s is None:
        return r_ref, False
    if r_ref is None:
        return r_s, True
    dev = abs(r_s / r_ref - 1)
    if soft is not None:
        a, b = soft
        w = clamp((dev - a) / (b - a), 0.0, 1.0)
        return r_ref + w * (r_s - r_ref), w > 0
    thr = trig_exit if (prev_engaged and trig_exit is not None) else trig
    use_s = dev > thr
    return (r_s if use_s else r_ref), use_s


def clamp_target(cand, parent, base, clamp_pct, ease_pct):
    """Spec 2.3 rule 4: harden by at most clamp_pct and ease by at most ease_pct against the parent, then the floor and
    ceiling. `base` (the clamp-window candidate): the same bounds against the newest chain block at least one window back."""
    out = clamp(cand, parent / (1 + clamp_pct), parent * (1 + ease_pct))
    if base is not None:
        out = clamp(out, base / (1 + clamp_pct), base * (1 + 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)
    if name.startswith("igneum:"):
        # igneum:key=value+key=value with the Igneum constructor's names; soft=a/b; v2 expands to the adopted candidate
        opts = dict(kw)
        opts.update(parse_variant(name[len("igneum:"):]))
        return Igneum(g, **opts)
    if name == "igneum-v2":
        return Igneum(g, **dict(V2, **kw))
    raise ValueError(name)


# Rule v2, adopted 4 October 2026 (docs/analysis/difficulty-2026-10-04-oscillation.md): the reference lane sees the
# newest 600 blocks of the epoch through the epoch lane only (the sampled long lane is not consulted). `igneum-v2` runs it.
V2 = {"ref_window": 600, "long_min": 1 << 62}


def parse_variant(spec):
    out = {}
    for part in [p for p in spec.replace(",", "+").split("+") if p]:
        key, val = part.split("=")
        if key == "v2" or key == "soft" and "/" in val:
            if key == "v2":
                out.update(V2)
            else:
                a, b = val.split("/")
                out["soft"] = (float(a), float(b))
            continue
        if key in ("ns", "k_min", "k0", "long_min", "window", "rate", "lag_caps", "ref_window"):
            out[key] = int(val)
        elif key in ("cap", "clamp_window"):
            out[key] = float(val) * T if float(val) < 1000 else float(val)
        elif key in ("median3", "literal", "epoch_aware"):
            out[key] = val in ("1", "true", "yes")
        elif key == "trig_exit" and val in ("none", "None"):
            out[key] = None
        else:
            out[key] = float(val)
    return out


# ---------------------------------------------------------------- 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


# ---------------------------------------------------------------- DAG model (4 October 2026, the live oscillation)
# Several miners on several nodes. A miner fetches a template from its node (the node stamps it with its clock at the
# fetch and the miner never touches the stamp, as igneum-miner does), the template reaches the card `delay` later (the
# miner's job queue), and the miner refetches every `period` and right after its own block. A block carries the parents
# and the stamp of the template it was mined from, so blocks found inside one another's staleness are parallel and the
# next block merges them. GHOSTDAG with every merged block blue (k = 18 is far above the width seen here): selected
# parent = the parent with the most blue work, blue work = the selected parent's plus the work of its mergeset blues
# (the selected parent included), DAA score = the selected parent's plus the mergeset size. The rule runs exactly as
# `SampledDifficultyManager::igneum_difficulty_bits`: the chain walk from the virtual's selected parent with blue-work
# increments over sanitised-clock steps, and the sampled window (every 4th block of each mergeset, selected parent
# first then descending blue work, 661 samples) restricted to the epoch for the long lane.
class DagBlock:
    __slots__ = ("id", "parents", "sp", "ms", "work", "bw", "daa", "bs", "ts", "clock", "target", "past", "lane",
                 "found", "miner", "samples", "children", "arrival", "engaged_frac")

    def __init__(self, bid):
        self.id = bid
        self.parents = []
        self.sp = None
        self.ms = []          # mergeset without the selected parent, descending blue work
        self.children = []
        self.samples = []     # (ts, bw, daa) of the sampled blocks of this block's mergeset
        self.lane = False
        self.past = set()     # ids of recent ancestors (PAST_DEPTH ids back)


PAST_DEPTH = 400


class DagRule:
    """The Igneum rule on a DAG, parameters taken from an `Igneum` controller instance (so `igneum:...` variants apply)."""

    def __init__(self, ctrl):
        self.c = ctrl
        self.cache = {}

    def mergeset(self, parents, sp):
        out = []
        seen = set()
        stack = [p for p in parents if p is not sp]
        while stack:
            x = stack.pop()
            if x.id in seen or x is sp or x.id in sp.past or x.id < sp.id - PAST_DEPTH:
                continue
            seen.add(x.id)
            out.append(x)
            stack.extend(x.parents)
        out.sort(key=lambda b: (-b.bw, b.id))
        return out

    def virtual(self, tips):
        sp = max(tips, key=lambda b: (b.bw, -b.id))
        ms = self.mergeset(tips, sp)
        daa = sp.daa + 1 + len(ms)
        return sp, ms, daa

    @staticmethod
    def samples_of(sp, ms):
        """Kaspa's sampled_mergeset_iterator for a block whose selected parent is sp: index 1 = sp, then the mergeset in
        descending blue work; sampled where (daa(sp) + index) % rate == 0."""
        out = []
        for i, b in enumerate([sp] + ms, start=1):
            if (sp.daa + i) % 4 == 0:
                out.append((b.ts, b.bw, b.daa))
        return out

    def target_for(self, tips):
        """Target, DAA score and short-lane flag of the template built on `tips`."""
        key = tuple(sorted(b.id for b in tips))
        hit = self.cache.get(key)
        if hit is not None:
            return hit
        c = self.c
        sp, ms, daa = self.virtual(tips)
        parent = sp.target
        epoch_start = (daa // EPOCH) * EPOCH
        k = daa - epoch_start
        eng = False
        if k < c.k_min:
            cand = parent
        else:
            ref_start = max(epoch_start, daa - c.ref_window) if c.ref_window else epoch_start
            # long lane: the virtual's own sampled mergeset, then the chain's, newest first, 661 samples
            r_ref = None
            if k >= c.long_min:
                samples = self.samples_of(sp, ms)
                cur = sp
                while cur.id != 0 and len(samples) < c.window and cur.daa >= ref_start:
                    samples.extend(cur.samples)
                    cur = cur.sp
                samples = [s for s in samples[:c.window] if s[2] >= ref_start]
                if len(samples) >= c.long_min // c.rate:
                    mn = min(samples, key=lambda s: s[0])
                    mx = max(samples, key=lambda s: s[0])
                    if mx[1] > mn[1] and mx[0] > mn[0]:
                        r_ref = (mx[1] - mn[1]) / (mx[0] - mn[0])
            # chain walk: blue-work increments over sanitised-clock steps, oldest first
            limit = c.ns if r_ref is not None else max(c.long_min, c.ns)
            steps = []
            cur = sp
            while cur.id != 0 and cur.daa >= ref_start and len(steps) < limit:
                p = cur.sp
                steps.append((cur.bw - p.bw, min(cur.clock - p.clock, c.cap)))
                cur = p
            steps.reverse()
            r_s = self.short(steps, c.ns)
            if c.median3 and len(steps) > 4:
                r_s = sorted([r_s, self.short(steps[:-1], c.ns), self.short(steps[:-2], c.ns)])[1]
            if r_ref is None:
                work = sum(w for w, _ in steps)
                time = max(sum(s for _, s in steps), 1)
                kk = daa - ref_start
                r_ref = (kk * work / time + c.k0 * (ONE / parent) / T) / (kk + c.k0)
            r, eng = choose_rate(r_s, r_ref, sp.lane, c.trig, c.trig_exit, c.soft)
            cand = ONE / (r * T)
        base = None
        if c.clamp_window:
            b = sp
            while b.id != 0 and sp.clock - b.clock < c.clamp_window:
                b = b.sp
            base = b.target
        out = (clamp_target(cand, parent, base, c.clamp_pct, c.ease_pct), daa, eng)
        self.cache[key] = out
        return out

    @staticmethod
    def short(steps, ns):
        n = min(ns, len(steps))
        if n == 0:
            return None
        wsum = 0.0
        dsum = 0.0
        for i, (w, s) in enumerate(steps[len(steps) - n:], start=1):
            wsum += s * i
            dsum += w * i
        wsum = max(wsum, n * n * T / 20)
        return dsum / wsum


class DagMiner:
    def __init__(self, name, node, period_ms, delay_ms, rpc_ms=30):
        self.name, self.node, self.period, self.delay, self.rpc = name, node, period_ms, delay_ms, rpc_ms
        self.hr = 0.0          # hashes per ms
        self.template = None   # (tips, ts, target, daa, eng)
        self.pending = []      # (apply time, template) fetched but still in the job queue
        self.t_apply = None
        self.t_fetch = 0
        self.t_solve = None


def simulate_dag(miners, schedule, ctrl, t0, t_end, rng, genesis_d, genesis_ts, nodes=2, prop_ms=100, verify_ms=30):
    """Blocks of a DAG mined by `miners` under the hash-rate `schedule` (sorted (t_ms, miner name, hashes per ms)).
    Returns (blocks, chain ids, lane flips)."""
    rule = DagRule(ctrl)
    c = ctrl
    genesis = DagBlock(0)
    genesis.ts = genesis_ts
    genesis.clock = genesis_ts
    genesis.target = ONE / genesis_d
    genesis.work = genesis_d
    genesis.bw = 0.0
    genesis.daa = 0
    genesis.bs = 0
    genesis.sp = genesis
    genesis.found = t0
    genesis.arrival = [t0] * nodes
    genesis.miner = "genesis"
    blocks = [genesis]
    views = [[(t0, genesis)] for _ in range(nodes)]   # per node, (arrival, block) in arrival order
    by_name = {m.name: m for m in miners}
    for m in miners:
        m.t_fetch = t0 + rng.random() * m.period
    si = 0
    t = t0
    flips = 0
    last_lane = False

    def tips_at(node, tau):
        v = views[node]
        out = []
        for i in range(len(v) - 1, -1, -1):
            arr, b = v[i]
            if arr > tau:
                continue
            if all(ch.arrival[node] > tau for ch in b.children):
                out.append(b)
            if len(v) - i > 60 and out:
                break
        return out

    def redraw(m, now):
        if m.hr > 0 and m.template is not None:
            D = ONE / m.template[2]
            m.t_solve = now + rng.expovariate(1.0) * D / m.hr
        else:
            m.t_solve = None

    while True:
        # next event
        t_sched = schedule[si][0] if si < len(schedule) else None
        ev, tm, who = None, None, None
        for m in miners:
            for kind, tt in (("solve", m.t_solve), ("fetch", m.t_fetch), ("apply", m.t_apply)):
                if tt is not None and (tm is None or tt < tm):
                    ev, tm, who = kind, tt, m
        if t_sched is not None and (tm is None or t_sched <= tm):
            ev, tm, who = "sched", t_sched, None
        if tm is None or tm > t_end:
            break
        t = tm
        if ev == "sched":
            _, name, hr = schedule[si]
            si += 1
            m = by_name[name]
            m.hr = hr
            redraw(m, t)
        elif ev == "fetch":
            m = who
            tips = tips_at(m.node, t)
            target, daa, eng = rule.target_for(tips)
            m.pending.append((t + m.delay, (tips, t, target, daa, eng)))
            m.t_apply = m.pending[0][0]
            m.t_fetch = t + m.period
        elif ev == "apply":
            m = who
            m.template = m.pending.pop(0)[1]
            m.t_apply = m.pending[0][0] if m.pending else None
            redraw(m, t)
        elif ev == "solve":
            m = who
            tips, ts, target, daa, eng = m.template
            b = DagBlock(len(blocks))
            b.parents = tips
            b.sp = max(tips, key=lambda x: (x.bw, -x.id))
            b.ms = rule.mergeset(tips, b.sp)
            b.ts = ts
            b.target = target
            b.work = ONE / target
            b.bw = b.sp.bw + b.sp.work + sum(x.work for x in b.ms)
            b.daa = b.sp.daa + 1 + len(b.ms)
            b.bs = b.sp.bs + 1 + len(b.ms)
            b.clock = max(b.sp.clock + clamp(ts - b.sp.clock, -c.cap, c.cap), ts - c.lag_caps * c.cap)
            b.lane = eng
            b.found = t
            b.miner = m.name
            b.samples = rule.samples_of(b.sp, b.ms)
            past = set()
            for p in tips:
                past.add(p.id)
                past.update(p.past)
            b.past = {x for x in past if x >= b.id - PAST_DEPTH}
            b.arrival = [t + (verify_ms if n == m.node else prop_ms) for n in range(nodes)]
            for p in tips:
                p.children.append(b)
            blocks.append(b)
            for n in range(nodes):
                v = views[n]
                arr = b.arrival[n]
                i = len(v)
                while i > 0 and v[i - 1][0] > arr:
                    i -= 1
                v.insert(i, (arr, b))
            if eng != last_lane:
                flips += 1
                last_lane = eng
            m.t_fetch = min(m.t_fetch, t + m.rpc)
            redraw(m, t)
    # the selected chain of the final virtual
    tips = tips_at(0, t_end + 10 * prop_ms)
    chain = set()
    cur = max(tips, key=lambda x: (x.bw, -x.id)) if tips else genesis
    while True:
        chain.add(cur.id)
        if cur.id == 0:
            break
        cur = cur.sp
    return blocks, chain, flips


def window_stats(items, secs):
    """items: (ts_ms, is_chain, mergeset_blues, parents, D, step_ms or None), D in the node's convention (2^255 / target)
    i.e. half the expected hashes. Chain items in chain order."""
    ch = [x for x in items if x[1]]
    if len(ch) < 10:
        return {"blocks": len(items)}
    D = [x[4] for x in ch]
    ld = [math.log(x) for x in D]
    mean_ld = sum(ld) / len(ld)
    std_ld = math.sqrt(sum((x - mean_ld) ** 2 for x in ld) / len(ld))
    peaks = []
    trough = D[0]
    up = False
    peak = pi = None
    for i, x in enumerate(D):
        if not up:
            trough = min(trough, x)
            if x >= trough * 1.15:
                up, peak, pi = True, x, i
        else:
            if x > peak:
                peak, pi = x, i
            if x <= peak / 1.15:
                peaks.append((pi, peak / trough))
                up, trough = False, x
    mins = {}
    for x in items:
        mins[int(x[0] // 60000)] = mins.get(int(x[0] // 60000), 0) + 1
    counts = sorted(mins.values())[1:-1] if len(mins) > 2 else list(mins.values())
    steps = sorted(x[5] for x in ch if x[5] is not None)
    return {"blocks": len(items), "blocks_s": len(items) / secs, "chain_s": len(ch) / secs, "blues_per_chain": len(items) / len(ch),
            "merge2_pct": 100.0 * sum(1 for x in ch if x[2] >= 2) / len(ch), "parents2_pct": 100.0 * sum(1 for x in items if x[3] >= 2) / len(items),
            "step_median_ms": steps[len(steps) // 2] if steps else None,
            "D_mean_M": sum(D) / len(D) / 1e6, "D_min_M": min(D) / 1e6, "D_max_M": max(D) / 1e6, "std_logD": std_ld,
            "peaks15": len(peaks), "peak_ratio": (sum(p[1] for p in peaks) / len(peaks)) if peaks else None,
            "peak_spacing_blocks": ((peaks[-1][0] - peaks[0][0]) / (len(peaks) - 1)) if len(peaks) > 1 else None,
            "bpm_min": min(counts) if counts else None, "bpm_max": max(counts) if counts else None}


LIVE_COLS = [("blocks", "blocks"), ("blocks_s", "blocks/s"), ("chain_s", "chain/s"), ("blues_per_chain", "blues per chain block"),
             ("merge2_pct", "chain blocks merging 2+ %"), ("parents2_pct", "blocks with 2+ parents %"), ("step_median_ms", "chain step median ms"),
             ("D_mean_M", "difficulty mean M"), ("D_min_M", "min M"), ("D_max_M", "max M"), ("std_logD", "std log D"),
             ("peaks15", "peaks over 15%"), ("peak_ratio", "mean peak/trough"), ("peak_spacing_blocks", "peak spacing chain blocks"),
             ("bpm_min", "blocks/min min"), ("bpm_max", "blocks/min max"), ("flips", "lane flips"), ("engaged_pct", "short lane %")]


def record_items(rows, t_lo, t_hi):
    by = {r["hash"]: r for r in rows}
    out = []
    for r in sorted(rows, key=lambda r: (int(r["blue_score"]), int(r["daa_score"]))):
        ts = int(r["timestamp_ms"])
        if int(r["daa_score"]) == 0 or not (t_lo <= ts < t_hi):
            continue
        p = by.get(r.get("selected_parent"))
        step = (ts - int(p["timestamp_ms"])) if (p and r["is_chain_block"] == "1") else None
        out.append((ts, r["is_chain_block"] == "1", int(r.get("mergeset_blues") or 1), int(r["parents"]), float(r["difficulty"]), step))
    return out


def dag_items(blocks, chain, t_lo, t_hi):
    out = []
    for b in blocks[1:]:
        if not (t_lo <= b.ts < t_hi):
            continue
        is_chain = b.id in chain
        out.append((b.ts, is_chain, 1 + len(b.ms), len(b.parents), b.work / 2.0, (b.ts - b.sp.ts) if is_chain else None, b.lane))
    return out


def live_schedule(path, t0, metal_mhs):
    """Hash-rate events from the records CSV (STATUS lines: the last 30 s of every worker): (t_ms, worker, hashes/ms).
    A worker is off from 5 s after its last report when the next one is more than 45 s away. The Mac Metal miner has
    no upload; it mines `metal_mhs` from t0."""
    rows = list(csv.DictReader(open(path)))
    byw = {}
    for r in rows:
        byw.setdefault(r["worker"], []).append((float(r["t_unix"]) * 1000, float(r["mh_s_last_30s"]) * 1e3))
    ev = [(t0, "mac-metal", metal_mhs * 1e3)]
    nodes = {"mac-metal": 1}
    for w, pts in byw.items():
        nodes[w] = 0 if rows[[r["worker"] for r in rows].index(w)]["node"] == "pc-node" else 1
        pts.sort()
        prev = None
        for t, hr in pts:
            if prev is None or t - prev > 45000:
                if prev is not None:
                    ev.append((prev + 5000, w, 0.0))
                ev.append((t - 30000, w, hr))
            elif abs(hr - prev_hr) > 0.03 * max(prev_hr, 1):
                ev.append((t - 30000, w, hr))
            prev, prev_hr = t, hr
        ev.append((prev + 5000, w, 0.0))
    ev.sort()
    return ev, nodes


def run_live(args, names):
    import datetime
    rows = list(csv.DictReader(open(args.live)))
    rows.sort(key=lambda r: (int(r["blue_score"]), int(r["daa_score"])))
    genesis = rows[0]
    first = next(r for r in rows if int(r["daa_score"]) > 0 or int(r["blue_score"]) == 1)
    t0 = int(first["timestamp_ms"]) - 1000
    day = datetime.datetime.utcfromtimestamp(t0 / 1000).strftime("%Y-%m-%d")

    def at(hhmm):
        return int(datetime.datetime.fromisoformat("%sT%s:00+00:00" % (day, hhmm)).timestamp() * 1000)

    windows = []
    for w in args.live_windows.split(","):
        name, span = w.split("=")
        a, b = span.split("-")
        windows.append((name, at(a), at(b)))
    t_end = at(args.live_end) if args.live_end else max(int(r["timestamp_ms"]) for r in rows)
    genesis_d = ONE / bits_to_target(int(genesis["bits"], 16))
    schedule, node_of = live_schedule(args.hashrate, t0, args.metal_mhs)
    print("# Live replay of %s: %d headers, DAG model, seed %d, delay scale %.2f\n" % (args.live, len(rows), args.seed, args.dag_delay))
    print("Hash-rate schedule (from the log intake and the Mac's bench value), MH/s:\n")
    print("| time UTC | worker | node | MH/s |")
    print("|---|---|---|---|")
    for t, w, hr in schedule:
        print("| %s | %s | %s | %.1f |" % (datetime.datetime.utcfromtimestamp(t / 1000).strftime("%H:%M:%S"), w, "pc" if node_of[w] == 0 else "mac", hr / 1e3))
    results = {}
    for name, lo, hi in windows:
        results.setdefault("record", {})[name] = window_stats(record_items(rows, lo, hi), (hi - lo) / 1000)
    for cn in names:
        ctrl = make_controller(cn, ONE / genesis_d)
        if not isinstance(ctrl, Igneum):
            continue
        rng = random.Random(args.seed)
        miners = []
        for w, n in node_of.items():
            if w == "mac-metal":
                miners.append(DagMiner(w, n, 630, int(900 * args.dag_delay)))
            else:
                miners.append(DagMiner(w, n, 150, int(350 * args.dag_delay)))
        blocks, chain, flips = simulate_dag(miners, schedule, ctrl, t0, t_end, rng, genesis_d, int(genesis["timestamp_ms"]), prop_ms=int(100 * args.dag_delay) + 20)
        for name, lo, hi in windows:
            items = dag_items(blocks, chain, lo, hi)
            st = window_stats([x[:6] for x in items], (hi - lo) / 1000)
            seg = [b for b in blocks[1:] if lo <= b.ts < hi and b.id in chain]
            st["flips"] = sum(1 for a, b in zip(seg, seg[1:]) if a.lane != b.lane)
            st["engaged_pct"] = 100.0 * sum(1 for b in seg if b.lane) / len(seg) if seg else None
            results.setdefault(cn, {})[name] = st
        if args.live_dump:
            with open(args.live_dump % cn.replace(":", "_").replace("/", "_") if "%s" in args.live_dump else args.live_dump, "w") as f:
                f.write("id,found_ms,ts_ms,daa,blue_score,chain,parents,mergeset_blues,difficulty,miner,lane\n")
                for b in blocks[1:]:
                    f.write("%d,%d,%d,%d,%d,%d,%d,%d,%.0f,%s,%d\n" % (b.id, b.found, b.ts, b.daa, b.bs, b.id in chain, len(b.parents), 1 + len(b.ms), b.work / 2, b.miner, b.lane))
    for name, lo, hi in windows:
        print("\n### Window %s (%s to %s UTC)\n" % (name, datetime.datetime.utcfromtimestamp(lo / 1000).strftime("%H:%M"), datetime.datetime.utcfromtimestamp(hi / 1000).strftime("%H:%M")))
        cols = [(k, h) for k, h in LIVE_COLS if any(k in results[r].get(name, {}) for r in results)]
        print("| source | " + " | ".join(h for _, h in cols) + " |")
        print("|---|" + "---|" * len(cols))
        for r in results:
            st = results[r].get(name, {})
            print("| %s | " % r + " | ".join(fmt(st.get(k), 3 if k in ("std_logD", "blues_per_chain") else (2 if k in ("blocks_s", "chain_s", "peak_ratio") else 1)) for k, _ in cols) + " |")


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("--live", default="", help="DAG replay of a live record (records/live-2026-10-04.csv) against the controllers")
    ap.add_argument("--hashrate", default="", help="live replay: the hash-rate CSV from the log intake (records/live-2026-10-04-hashrate.csv)")
    ap.add_argument("--live-windows", default="solo=09:25-10:04,both=10:17-10:32,late=10:37-10:54", help="live replay: name=HH:MM-HH:MM UTC windows to compare")
    ap.add_argument("--live-end", default="", help="live replay: stop at HH:MM UTC (default: the record's last stamp)")
    ap.add_argument("--live-dump", default="", help="live replay: write every simulated block to this CSV (%s = controller name)")
    ap.add_argument("--dag-delay", type=float, default=1.0, help="live replay: scale of the miners' template staleness and propagation (1 = fitted)")
    ap.add_argument("--metal-mhs", type=float, default=26.7, help="live replay: the Mac Metal miner's hash rate (no upload; bench-log value)")
    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(",")
    if args.live:
        run_live(args, names)
        return
    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()