#!/usr/bin/env python3
"""Igneum economy simulator: GPU operators choosing between mining and proving.

Agent-based model of N operators (default 1,000) with heterogeneous card fleets and electricity
prices. Every operator runs its own client and maximises its own profit: every ten minutes (own
phase) it compares, per card class, the expected profit of MINING (lottery share of the 80%
emission), PROVING (internal shards from the 20% pool by sortition then open claiming, plus
external jobs in dollars settled in IGN with a burn) and HYBRID (mine, prove only the shards it is
assigned, program swap each way), and switches when the gain beats its own hysteresis. Nothing is
scheduled centrally.

Time is in ticks of TICK seconds (180 s). Everything that happens inside a tick (the 10-s sortition
window, shard times of 6 to 20 s, open-claim races) is resolved with closed-form rates; everything
across ticks (hash, difficulty, backlog, weights, modes, price) is state. Blocks, shards, jobs and
block attribution are sampled; shard allocation across operators is fluid (expected values).

All figures are approximate unless the source says measured (RTX 5090 229 MH/s is measured,
docs/bench-log.md). See docs/analysis/economy-2026-10-04.md for the assumptions table.

Usage
  python3 sim.py                          six scenarios, 5 seeds, markdown to stdout
  python3 sim.py --scenarios b --seeds 1  one run
  python3 sim.py --levers b --seeds 3     lever study on scenario b
  python3 sim.py --set pool=0.3,window=20 override parameters
  python3 sim.py --dump b --seeds 1       hourly trajectory CSV for scenario b
"""
import argparse
import math
import sys
import time

import numpy as np

# ---------------------------------------------------------------- parameters (all approximate)
P = dict(
    tick=180.0,            # s per tick
    days=30,
    n_ops=1000,
    emission=31.688,       # IGN per block, pre-halving, ramp complete (spec 2.5)
    pool=0.20,             # proving pool share of emission (spec 2.5, 5.3)
    window=10.0,           # sortition exclusive window, s (spec 7.2)
    assignees=8,           # provers drawn per shard (spec 7.2)
    burn=0.10,             # external job burn (spec 5.4)
    timeout=300.0,         # external job claim timeout, s (O-5.6, open; design 6 says 3,600 blocks to expire)
    shards_per_block=3.0,  # internal demand, shards (3060-20-s units) per block at launch traffic
    job_work=10.0,         # shard-equivalents of work per external job
    ext_usd_day=2000.0,    # external demand, dollars per day (customer brief: market low millions/yr)
    ext_mult=1.0,          # scenario multiplier on the dollar price per job
    price0=0.012,          # $ per IGN at t=0
    vol_day=0.05,          # price daily volatility
    swap=5.0,              # program swap, s each way (hybrid)
    aggregation=4.0,       # aggregation plus inclusion, s, added to every block proof
    waste=0.5,             # wasted duplicate attempts per open-claim shard (race losers)
    hyst_lo=0.05, hyst_hi=0.25,
    dwell_ticks=20,        # minimum ticks between switches of one (operator, class): 1 h
    decide_every=4,        # ticks between decisions: 12 min
    ema_ticks=20.0,        # observation window for realised rates: 1 h
    farm_share=0.20,       # operator 0 share of total hash
    backlog_rule=600.0,    # s of oldest age per halving of B_p (design 4.3)
)

# card classes: 5090 (measured hash), 3090, 3060, small (approximate)
CLS = ["5090", "3090", "3060", "small"]
HASH = np.array([229.0, 57.0, 46.0, 25.0])          # MH/s
PMINE = np.array([450.0, 320.0, 170.0, 120.0])      # W while hashing
PPROVE = np.array([500.0, 350.0, 170.0, 0.0])       # W while proving
PIDLE = np.array([60.0, 50.0, 30.0, 0.0])           # W proving-ready, idle
TPROVE = np.array([6.0, 12.0, 20.0, 1e9])           # s per shard (3060 = phase 2 gate target); small cannot prove
CANHYB = np.array([True, True, False, False])       # hybrid needs the dataset and the prover resident: 24 GB cards only
MIX = np.array([0.25, 0.25, 0.30, 0.20])            # fleet mix by card count
CANPROVE = np.array([True, True, True, False])
OFF, MINE, PROVE, HYB = 0, 1, 2, 3


def build_population(rng, p, extra=0):
    n = p["n_ops"] + extra
    size = np.maximum(1, np.round(rng.lognormal(math.log(5.0), 1.2, n))).astype(int)
    cards = np.zeros((n, 4), dtype=float)
    for i in range(n):
        cards[i] = rng.multinomial(size[i], MIX)
    # electricity: farms cheaper; lognormal around $0.10/kWh
    z = (np.log(size) - np.log(5.0)) / 1.2
    elec = np.exp(rng.normal(math.log(0.10) - 0.25 * z, 0.40, n))
    elec = np.clip(elec, 0.02, 0.40)
    # operator 0: the farm, all 5090s, sized to farm_share of total hash, cheap power
    cards[0] = 0
    rest = (cards[1:p["n_ops"]] * HASH).sum()
    farm_hash = rest * p["farm_share"] / (1 - p["farm_share"])
    cards[0, 0] = max(1, round(farm_hash / HASH[0]))
    elec[0] = 0.05
    hyst = rng.uniform(p["hyst_lo"], p["hyst_hi"], n)
    phase = rng.integers(0, p["decide_every"], n)
    return cards, elec, hyst, phase


class Sim:
    def __init__(self, scenario, seed, p):
        self.p = dict(p)
        self.sc = scenario
        self.rng = np.random.default_rng(seed)
        p = self.p
        extra = 0
        if scenario == "f":
            extra = 200
        if scenario == "e":
            p["farm_share"] = 0.30
        self.cards, self.elec, self.hyst, self.phase = build_population(self.rng, p, extra)
        self.n = self.cards.shape[0]
        self.n_base = p["n_ops"]
        self.active = np.ones(self.n, bool)
        if extra:
            self.active[self.n_base:] = False
            # the incoming pool: 200 operators sized so their hash equals the base hash
            base_hash = (self.cards[: self.n_base] * HASH).sum()
            pool_hash = (self.cards[self.n_base:] * HASH).sum()
            self.cards[self.n_base:] *= base_hash / pool_hash
            self.cards[self.n_base:] = np.round(self.cards[self.n_base:])
            self.elec[self.n_base:] = np.clip(self.elec[self.n_base:], 0.02, 0.08)
        self.mode = np.full((self.n, 4), MINE, dtype=int)
        self.mode[self.cards == 0] = OFF
        # initial allocation: cards with cheap power start a fifth of their provable cards proving
        for i in range(self.n):
            for c in range(3):
                if self.cards[i, c] > 0 and self.rng.random() < 0.20:
                    self.mode[i, c] = PROVE
        if not self.active.all():
            self.mode[~self.active] = OFF
        self.fixed_mine = np.zeros(self.n, bool)
        self.withhold = np.zeros(self.n, bool)
        if scenario == "e":
            self.fixed_mine[0] = True
            self.withhold[0] = True
            self.mode[0] = np.where(self.cards[0] > 0, MINE, OFF)
        self.dead = np.zeros(self.n, bool)
        # weight: 30-day ring of mined blocks per operator, seeded from hash shares
        self.hist = np.zeros((self.n, 30))
        h0 = (self.cards * HASH * (self.mode != OFF)).sum(1)
        h0[~self.active] = 0
        self.hist[:] = (h0 / h0.sum() * 86400.0)[:, None]
        self.w = self.hist.sum(1)
        self.cards_active = self.cards[self.active].sum()
        self.day = 0
        self.last_switch = np.full((self.n, 4), -10**9, dtype=int)
        self.price = p["price0"]
        self.H_est = h0.sum() * 1e6
        self.D = self.H_est * 1.0   # expected hashes per block, 1 block/s
        self.q = 0.0                # backlog, shards
        self.age = 0.0              # oldest unproven block age, s
        self.duty = np.zeros((self.n, 4))   # hybrid proving duty
        # observed rates (EMA)
        self.R_hash = p["emission"] * (1 - p["pool"]) / self.H_est   # IGN per hash per s
        self.open_pc = np.zeros(4)   # $ per hour per PROVE card by class, open claims
        self.asg_pw = 0.0            # $ per hour per unit weight, assigned shards
        self.asg_pw_h = 0.0          # same for hybrid responders
        self.asg_ext_pw = 0.0        # $ per hour per unit weight, assigned external jobs
        self.t_open = p["window"] + TPROVE[0]
        self.eligible = CANPROVE.copy()
        self.ok60 = 1.0
        self.ok20 = 1.0
        self.hourly = []
        self.hour_acc = []
        self.cls_profit = np.zeros(4)
        self.cls_cards = np.zeros(4)
        self.cls_shards = np.zeros(4)
        self.ext_served = 0.0
        self.ext_demand = 0.0
        self.burned_ign = 0.0
        self.switches = 0

    # ------------------------------------------------------------ one tick
    def step(self, t):
        p = self.p
        T = p["tick"]
        day_f = t * T / 86400.0
        rng = self.rng
        # scenario events
        if self.sc == "d" and day_f >= 10 and not self.dead[0]:
            self.dead[0] = True
            self.mode[0] = OFF
        if self.sc == "f" and day_f >= 10 and not self.active[self.n_base]:
            self.active[self.n_base:] = True
            self.mode[self.n_base:] = np.where(self.cards[self.n_base:] > 0, MINE, OFF)
            self.cards_active = self.cards[self.active].sum()
        ext_mult = p["ext_mult"]
        if self.sc == "b" and day_f >= 7:
            ext_mult = 10.0
        if self.sc == "c":
            ext_mult = 0.0
        # price
        dt = T / 86400.0
        self.price *= math.exp(p["vol_day"] * math.sqrt(dt) * rng.standard_normal() - 0.5 * p["vol_day"] ** 2 * dt)
        if self.sc == "b" and 7 <= day_f < 14:
            self.price *= math.exp(math.log(0.30) / (7 * 86400.0 / T))
        price = self.price

        cards, mode = self.cards, self.mode
        is_m = mode == MINE
        is_p = mode == PROVE
        is_h = mode == HYB
        # hash
        hcards = cards * HASH * 1e6
        h_op = (hcards * (is_m + is_h * (1 - self.duty))).sum(1)
        H = h_op.sum()
        if H <= 0:
            H = 1.0
        # difficulty controller (dual-lane response approximated: 120-block estimate, 3%/10% clamps)
        lam = T * H / self.D
        blocks = rng.poisson(lam)
        self.H_est += (H - self.H_est) * min(1.0, blocks / 120.0)
        ratio = self.H_est / self.D
        ratio = min(max(ratio, 0.90 ** blocks), 1.03 ** blocks)
        self.D *= ratio
        E = p["emission"]
        # mined blocks per operator
        mined = rng.multinomial(blocks, h_op / H) if blocks > 0 else np.zeros(self.n, int)
        self.hist[:, self.day % 30] += mined
        self.w += mined
        w = self.w   # a dead operator's weight stays in the window until it ages out
        W = w.sum()
        mine_ign = mined * E * (1 - p["pool"])
        pool_ign = blocks * E * p["pool"]

        # internal shard demand with the backlog rule
        spb = p["shards_per_block"] * 0.5 ** math.floor(self.age / p["backlog_rule"])
        shards = rng.poisson(blocks * spb) if blocks > 0 else 0
        pool_per_shard = pool_ign / shards if shards > 0 else 0.0

        # capacities in shard-equivalents per tick
        capP = np.where(is_p, cards * T / TPROVE, 0.0)
        capP[:, 3] = 0
        capH = np.where(is_h, cards * T / (TPROVE + 2 * p["swap"]), 0.0)
        capH[:, 2:] = 0
        capP[self.dead | ~self.active] = 0
        capH[self.dead | ~self.active] = 0
        capP_op = capP.sum(1)
        capH_op = capH.sum(1)
        resp = ((capP_op + capH_op) >= 1.0) & ~self.withhold & ~self.dead
        a = (w * resp).sum() / W if W > 0 else 0.0
        p_resp = 1 - (1 - a) ** p["assignees"]
        shards_served = np.zeros((self.n, 4))
        # external jobs
        jobs_usd = p["ext_usd_day"] * ext_mult
        job_price = (p["ext_usd_day"] / (p["ext_usd_day"] / 50.0)) * ext_mult  # $50 per job baseline
        jobs = rng.poisson(p["ext_usd_day"] / 50.0 * T / 86400.0) if jobs_usd > 0 else 0
        ext_work = jobs * p["job_work"]
        ext_pay_per_shard = job_price / p["job_work"] * (1 - p["burn"])
        eligible = (p["job_work"] * TPROVE <= p["timeout"]) & CANPROVE
        ext_served = np.zeros((self.n, 4))
        ext_asg = [0.0]
        ext_open = np.zeros(4)
        self.ext_demand += ext_work
        internal_first = pool_per_shard * price >= ext_pay_per_shard

        def serve_external():
            # jobs are assigned by the same sortition as shards (design 6), claimed with a bond, so the
            # assignee keeps the job: the assigned part goes by weight, the rest is open and the fastest
            # eligible proving card claims it; what nobody can take before the deadline expires.
            nonlocal ext_work
            if ext_work <= 0:
                return
            capE_P = capP * eligible
            capE_H = capH * eligible
            capE_op = capE_P.sum(1) + capE_H.sum(1)
            respE = (capE_op >= 1.0) & ~self.withhold & ~self.dead
            aE = (w * respE).sum() / W if W > 0 else 0.0
            pE = 1 - (1 - aE) ** p["assignees"]
            s_asg = ext_work * pE
            over = 0.0
            wrE = w * respE
            WrE = wrE.sum()
            if s_asg > 0 and WrE > 0:
                x = s_asg * wrE / WrE
                got = np.minimum(x, capE_op)
                over = s_asg - got.sum()
                capE = capE_P + capE_H
                frac = capE / np.maximum(capE.sum(1, keepdims=True), 1e-9)
                alloc = got[:, None] * frac
                ext_served[:] += alloc
                ext_asg[0] += got.sum()
                usedP = alloc * (capE_P / np.maximum(capE, 1e-9))
                capP[:] -= usedP
                capH[:] -= alloc - usedP
            else:
                over = s_asg
            open_w = ext_work - s_asg + over
            for lat, cap, c in self.open_order(capP, capH):
                if not eligible[c] or open_w <= 0:
                    continue
                tot = cap[:, c].sum()
                if tot <= 0:
                    continue
                take = min(open_w, tot / (1 + p["waste"]))
                alloc = cap[:, c] * (take / tot)
                ext_served[:, c] += alloc
                ext_open[c] += take
                cap[:, c] -= alloc * (1 + p["waste"])
                open_w -= take
            ext_work = open_w

        if not internal_first:
            serve_external()
        # assignee race: class mix of responders
        wr = w * resp
        Wr = wr.sum()
        if Wr > 0:
            capPH = capP + capH
            share = capPH / np.maximum(capPH.sum(1, keepdims=True), 1e-9)
            mixP = (wr[:, None] * (capP / np.maximum(capPH.sum(1, keepdims=True), 1e-9))).sum(0) / Wr
            mixH = (wr[:, None] * (capH / np.maximum(capPH.sum(1, keepdims=True), 1e-9))).sum(0) / Wr
        else:
            mixP = np.zeros(4)
            mixH = np.zeros(4)
        t_open = p["window"] + np.inf
        for lat, cap, c in self.open_order(capP, capH):
            if cap[:, c].sum() > 0:
                t_open = p["window"] + lat
                break
        LA_P = TPROVE
        LA_H = TPROVE + p["swap"]
        winP = (LA_P <= t_open) * mixP
        winH = (LA_H <= t_open) * mixH
        p_win = winP[:3].sum() + winH[:3].sum()   # given p_resp
        s_asg = shards * p_resp * p_win
        # allocate assigned shards by weight, cap by capacity
        if s_asg > 0 and Wr > 0:
            x = s_asg * wr / Wr
            capPH_op = capP.sum(1) + capH.sum(1)
            got = np.minimum(x, capPH_op)
            overflow = s_asg - got.sum()
            # split within operator across classes by capacity
            capPH = capP + capH
            frac = capPH / np.maximum(capPH.sum(1, keepdims=True), 1e-9)
            alloc = got[:, None] * frac
            shards_served += alloc
            usedP = alloc * (capP / np.maximum(capPH, 1e-9))
            usedH = alloc - usedP
            capP -= usedP
            capH -= usedH
        else:
            overflow = s_asg
        open_new = shards - s_asg + overflow
        q_prev = self.q
        open_total = open_new + q_prev
        served_open = 0.0
        open_by_cls = np.zeros(4)
        open_lat = np.zeros(4)
        for lat, cap, c in self.open_order(capP, capH):
            tot = cap[:, c].sum()
            if tot <= 0 or open_total - served_open <= 0:
                continue
            take = min(open_total - served_open, tot / (1 + p["waste"]))
            alloc = cap[:, c] * (take / tot)
            shards_served[:, c] += alloc
            open_by_cls[c] += take
            open_lat[c] += take * lat
            cap[:, c] -= alloc * (1 + p["waste"])
            served_open += take
        self.q = max(0.0, open_total - served_open)
        if internal_first:
            serve_external()
        self.ext_served += ext_served.sum()
        # backlog age and proof latency
        arr_rate = max(shards / T, 1e-9)
        if self.q <= 0:
            self.age = 0.0
        elif served_open >= q_prev:
            # the old queue drained this tick; what is left arrived this tick
            self.age = min(T, self.q / arr_rate) + p["window"] + TPROVE[2]
        else:
            self.age += T
        wait = q_prev / max(served_open / T, 1e-9) if q_prev > 0 else 0.0
        thr60 = 60.0 - p["aggregation"]
        thr20 = 20.0 - p["aggregation"]
        p_q = min(1.0, max(0.0, (open_new - max(0.0, served_open - q_prev)) / max(shards, 1e-9))) if shards > 0 else 0.0
        # the assignee branch only counts where it wins; losers were served open
        pa = p_resp
        p_assigned_won = pa * p_win
        okA60 = pa * (((LA_P <= t_open) & (LA_P <= thr60)) * mixP)[:3].sum() + pa * (((LA_H <= t_open) & (LA_H <= thr60)) * mixH)[:3].sum()
        okA20 = pa * (((LA_P <= t_open) & (LA_P <= thr20)) * mixP)[:3].sum() + pa * (((LA_H <= t_open) & (LA_H <= thr20)) * mixH)[:3].sum()
        if open_by_cls[:3].sum() > 0:
            # open-served shards: latency = window + wait + class time, class by share served
            l_open = p["window"] + wait + open_lat[:3] / np.maximum(open_by_cls[:3], 1e-9)
            sh = open_by_cls[:3] / open_by_cls[:3].sum()
            okO60 = (sh * (l_open <= thr60)).sum()
            okO20 = (sh * (l_open <= thr20)).sum()
        else:
            okO60 = okO20 = 0.0
        p_open_served = max(0.0, 1 - p_assigned_won - p_q)
        ps60 = min(1.0, okA60 + p_open_served * okO60)
        ps20 = min(1.0, okA20 + p_open_served * okO20)
        lam_s = max(spb, 1e-9)
        def block_ok(ps):
            return (math.exp(-lam_s * (1 - ps)) - math.exp(-lam_s)) / (1 - math.exp(-lam_s))
        self.ok60 = block_ok(ps60)
        self.ok20 = block_ok(ps20)
        # income and power
        prove_ign = shards_served * pool_per_shard
        ext_usd = ext_served * ext_pay_per_shard
        self.burned_ign += ext_served.sum() * (job_price / p["job_work"]) * p["burn"] / max(price, 1e-9)
        busyP = np.where(is_p, np.minimum(1.0, (shards_served + ext_served) * TPROVE / np.maximum(cards * T, 1e-9)), 0.0)
        busyP[:, 3] = 0
        dutyH = np.where(is_h, np.minimum(1.0, (shards_served + ext_served) * (TPROVE + 2 * p["swap"]) / np.maximum(cards * T, 1e-9)), 0.0)
        dutyH[:, 3] = 0
        self.duty = dutyH
        kwh = cards * T / 3600.0 / 1000.0
        power_w = is_m * PMINE + is_p * (busyP * PPROVE + (1 - busyP) * PIDLE) + is_h * ((1 - dutyH) * PMINE + dutyH * PPROVE)
        cost = kwh * power_w * self.elec[:, None]
        # mining income per class within operator by hash share
        hsh = hcards * (is_m + is_h * (1 - self.duty))
        hfrac = hsh / np.maximum(hsh.sum(1, keepdims=True), 1e-9)
        mine_usd = (mine_ign * price)[:, None] * hfrac
        inc = mine_usd + prove_ign * price + ext_usd
        profit = inc - cost
        # observed rates
        al = 1.0 / p["ema_ticks"]
        self.R_hash += ((blocks * E * (1 - p["pool"]) / H / T) - self.R_hash) * al
        nPH_cls = (np.where(is_p | is_h, cards, 0.0)).sum(0)
        open_inc = open_by_cls * pool_per_shard * price + ext_open * ext_pay_per_shard   # $ by class from open claims
        open_pc = np.where(nPH_cls > 0, open_inc / np.maximum(nPH_cls, 1e-9), 0.0) * 3600.0 / T
        self.open_pc += (open_pc - self.open_pc) * al
        asg_usd = (s_asg - overflow) * pool_per_shard * price
        asg_pw = asg_usd / max(Wr, 1e-9) * 3600.0 / T
        self.asg_pw += (asg_pw - self.asg_pw) * al
        asg_ext_pw = ext_asg[0] * ext_pay_per_shard / max(W, 1e-9) * 3600.0 / T
        self.asg_ext_pw += (asg_ext_pw - self.asg_ext_pw) * al
        self.t_open = t_open
        self.eligible = eligible
        # decisions
        self.decide(t, price, w, Wr, p_resp, p_win, t_open)
        # accumulate hourly
        ca = max(self.cards_active, 1)
        n_p = (cards * is_p).sum()
        n_h = (cards * is_h).sum()
        n_m = (cards * is_m).sum()
        self.hour_acc.append((
            H / 1e6, (hcards * (mode != OFF)).sum() / 1e6, n_p / ca,
            n_h / ca, 1 - (n_p + n_h + n_m) / ca,
            blocks / T, self.D, self.q, self.age, self.ok60, self.ok20, price,
            ext_served.sum(), ext_work, shards,
        ))
        self.cls_profit += profit.sum(0)
        self.cls_cards += (cards * self.active[:, None]).sum(0)
        self.cls_shards += shards_served.sum(0)
        if len(self.hour_acc) * T >= 3600:
            arr = np.array(self.hour_acc)
            self.hourly.append(arr.mean(0).tolist() + [arr[:, 8].max(), arr[:, 0].min()])
            self.hour_acc = []
        if (t + 1) * T % 86400 == 0:
            self.day += 1
            self.w -= self.hist[:, self.day % 30]
            self.hist[:, self.day % 30] = 0

    def open_order(self, capP, capH):
        """Open-claim service order: the fastest responder wins the race. PROVE cards at their shard
        time, HYBRID cards at shard time plus one program swap."""
        p = self.p
        order = [(TPROVE[c], capP, c) for c in range(3)] + [(TPROVE[c] + p["swap"], capH, c) for c in range(2)]
        order.sort(key=lambda x: x[0])
        return order

    # ------------------------------------------------------------ decisions
    def decide(self, t, price, w, Wr, p_resp, p_win, t_open):
        p = self.p
        due = ((t + self.phase) % p["decide_every"] == 0) & self.active & ~self.dead & ~self.fixed_mine
        if not due.any():
            return
        idx = np.nonzero(due)[0]
        cards = self.cards[idx]
        elec = self.elec[idx][:, None]
        n_pcards = (self.cards[idx] * ((self.mode[idx] == PROVE) | (self.mode[idx] == HYB))).sum(1, keepdims=True)
        # $ per card-hour by class
        v_mine = HASH * 1e6 * self.R_hash * price * 3600.0 - PMINE * elec / 1000.0
        per_w = w[idx][:, None] / np.maximum(n_pcards + (n_pcards == 0) * cards, 1e-9)
        asg_int = self.asg_pw * per_w          # assigned internal shards, $ per card-hour at this operator's weight
        asg_ext = self.asg_ext_pw * per_w * self.eligible
        # a prover with no weight gets only open claims; one with weight gets assigned work shared over its proving cards
        v_prove = self.open_pc + asg_ext + asg_int * (TPROVE <= t_open) - PIDLE * elec / 1000.0 - 0.3 * (PPROVE - PIDLE) * elec / 1000.0
        v_hyb = v_mine * (1 - np.minimum(0.5, self.duty[idx])) + 0.9 * (asg_ext + asg_int * ((TPROVE + p["swap"]) <= t_open))
        v_prove[:, 3] = -1e9
        v_hyb[:, 2:] = -1e9
        v_off = np.zeros_like(v_mine)
        V = np.stack([v_off, v_mine, v_prove, v_hyb], 0)   # (4 modes, ops, classes)
        cur = self.mode[idx]
        v_cur = np.take_along_axis(V, cur[None], 0)[0]
        best = V.argmax(0)
        v_best = V.max(0)
        hy = self.hyst[idx][:, None]
        gain = v_best - v_cur
        ok = (gain > hy * np.abs(v_cur) + 1e-4) & (cards > 0) & ((t - self.last_switch[idx]) >= p["dwell_ticks"])
        ok &= best != cur
        if ok.any():
            newmode = np.where(ok, best, cur)
            self.mode[idx] = newmode
            ls = self.last_switch[idx]
            ls[ok] = t
            self.last_switch[idx] = ls
            self.switches += int(ok.sum())

    def run(self):
        nt = int(self.p["days"] * 86400 / self.p["tick"])
        for t in range(nt):
            self.step(t)
        return np.array(self.hourly)


COLS = ["hash_MHs", "hash_pot", "frac_prove", "frac_hyb", "frac_off", "bps", "D", "q", "age", "ok60", "ok20",
        "price", "ext_served", "ext_lost", "shards", "age_max", "hash_min"]


def metrics(hourly, sim):
    h = hourly
    d = {}
    hrs = h.shape[0]
    nd = hrs // 24
    def day_mean(col, d0, d1):
        d0, d1 = min(d0, nd - 1), min(d1, nd)
        return h[d0 * 24:d1 * 24, COLS.index(col)].mean()
    def day_max(col, dd):
        dd = min(dd, nd)
        return h[(dd - 1) * 24:dd * 24, COLS.index(col)].max()
    pre_hash = day_mean("hash_MHs", 1, 7) if nd >= 7 else h[:, COLS.index("hash_MHs")].mean()
    d["hash_pre"] = pre_hash
    d["hash_min_ratio"] = h[24:, COLS.index("hash_min")].min() / pre_hash if nd >= 2 else 1.0
    d["hash_d30_ratio"] = day_mean("hash_MHs", 29, 30) / pre_hash
    # hours below 50% of pre-event hash
    d["hours_hash_lt50"] = int((h[24:, COLS.index("hash_MHs")] < 0.5 * pre_hash).sum())
    d["mining_hash_share"] = (h[:, COLS.index("hash_MHs")] / np.maximum(h[:, COLS.index("hash_pot")], 1e-9)).mean()
    for dd in (1, 7, 10, 15, 20, 30):
        d[f"prove_d{dd}"] = day_mean("frac_prove", dd - 1, dd)
        d[f"hyb_d{dd}"] = day_mean("frac_hyb", dd - 1, dd)
        d[f"off_d{dd}"] = day_mean("frac_off", dd - 1, dd)
    d["age_max"] = h[:, COLS.index("age_max")].max()
    d["q_max"] = h[:, COLS.index("q")].max()
    d["age_d10"] = day_max("age_max", 10)
    d["age_d20"] = day_max("age_max", 20)
    d["age_d30"] = day_max("age_max", 30)
    daily_age = h[:, COLS.index("age_max")].reshape(-1, 24).max(1)
    x = np.arange(10)
    slope = np.polyfit(x, daily_age[-10:], 1)[0] if nd >= 10 else 0.0
    d["backlog_growing"] = bool(slope > 0 and daily_age[-1] > 60)
    d["backlog_600"] = bool(d["age_max"] > 600)
    d["ok60_mean"] = h[:, COLS.index("ok60")].mean()
    d["ok60_min_day"] = h[:, COLS.index("ok60")].reshape(-1, 24).mean(1).min()
    d["ok20_mean"] = h[:, COLS.index("ok20")].mean()
    d["bps_mean"] = h[:, COLS.index("bps")].mean()
    d["bps_max"] = h[:, COLS.index("bps")].max()
    d["bps_min"] = h[:, COLS.index("bps")].min()
    d["D_end_over_start"] = h[-1, COLS.index("D")] / h[min(24, hrs - 1), COLS.index("D")]
    d["price_end"] = h[-1, COLS.index("price")]
    es, el = h[:, COLS.index("ext_served")].sum(), h[:, COLS.index("ext_lost")].sum()
    d["ext_delivered"] = es / (es + el) if es + el > 0 else float("nan")
    # oscillation: zero crossings of the proving fraction around its 24-h mean, per day, last 10 days
    fp = h[:, COLS.index("frac_prove")]
    ma = np.convolve(fp, np.ones(24) / 24, mode="same")
    dev = (fp - ma)[-240:]
    d["osc_cross_per_day"] = float((np.diff(np.sign(dev)) != 0).sum() / 10.0)
    d["osc_amp_pts"] = float((np.percentile(fp[-240:], 90) - np.percentile(fp[-240:], 10)) * 100)
    d["osc_std_pts"] = float(fp[-240:].std() * 100)
    days = sim.p["days"]
    d["profit_card_day"] = sim.cls_profit / np.maximum(sim.cls_cards / (days * 86400 / sim.p["tick"]), 1e-9) / days
    d["shard_share"] = sim.cls_shards / max(sim.cls_shards.sum(), 1e-9)
    d["switches"] = sim.switches
    d["miner_shortage"] = bool(d["hours_hash_lt50"] >= 1)
    d["window_miss"] = bool(d["ok60_min_day"] < 0.90)
    d["oscillation"] = bool(d["osc_amp_pts"] > 10 and d["osc_cross_per_day"] >= 1)
    return d


SCEN = {
    "a": "baseline",
    "b": "external pays 10x from day 7, coin price falls 70% over days 7 to 14",
    "c": "no external demand",
    "d": "the 20% operator disappears at day 10",
    "e": "a 30% operator never fulfils its assignments",
    "f": "a pool with hash equal to the network's arrives at day 10 (2x hash)",
}


def fmt(v, nd=2):
    if isinstance(v, (bool, np.bool_)):
        return "yes" if v else "no"
    if isinstance(v, (int, np.integer)):
        return f"{v:,}"
    if isinstance(v, float):
        if abs(v) >= 1000:
            return f"{v:,.0f}"
        return f"{v:.{nd}f}"
    return str(v)


def agg(ms, key):
    vals = np.array([m[key] for m in ms], dtype=float)
    return vals.mean(), vals.min(), vals.max()


def run_scenarios(scen, seeds, p, quiet=False):
    out = {}
    for s in scen:
        ms = []
        for seed in seeds:
            t0 = time.time()
            sim = Sim(s, seed, p)
            h = sim.run()
            ms.append(metrics(h, sim))
            if not quiet:
                print(f"<!-- scenario {s} seed {seed}: {time.time() - t0:.1f} s -->", file=sys.stderr)
        out[s] = ms
    return out


def table_scenarios(out):
    keys = [
        ("mining_hash_share", "hash mining share (mean)"),
        ("prove_d1", "cards proving, day 1"), ("prove_d10", "cards proving, day 10"), ("prove_d20", "cards proving, day 20"), ("prove_d30", "cards proving, day 30"),
        ("hyb_d30", "cards hybrid, day 30"), ("off_d30", "cards off, day 30"),
        ("hash_min_ratio", "hash min / pre-event"), ("hash_d30_ratio", "hash day 30 / pre-event"), ("hours_hash_lt50", "hours hash under 50%"),
        ("q_max", "backlog max, shards"), ("age_max", "oldest unproven age max, s"),
        ("age_d10", "age max day 10, s"), ("age_d20", "age max day 20, s"), ("age_d30", "age max day 30, s"),
        ("ok60_mean", "blocks proven within 60 s"), ("ok60_min_day", "worst day within 60 s"), ("ok20_mean", "blocks proven within 20 s"),
        ("bps_mean", "blocks/s mean"), ("bps_max", "blocks/s hourly max"), ("bps_min", "blocks/s hourly min"), ("D_end_over_start", "difficulty end / start"),
        ("ext_delivered", "external jobs delivered"), ("price_end", "price at day 30, $"),
        ("osc_cross_per_day", "proving-share crossings per day"), ("osc_amp_pts", "proving-share 10-90 pct range, points"),
        ("switches", "mode switches (operator x class)"),
    ]
    scen = list(out)
    lines = ["| Metric | " + " | ".join(f"{s} | {s} min | {s} max" for s in scen) + " |",
             "|---|" + "---|" * (3 * len(scen))]
    for k, label in keys:
        row = [label]
        for s in scen:
            m, lo, hi = agg(out[s], k)
            nd = 4 if k == "price_end" else 2
            row += [fmt(float(m), nd), fmt(float(lo), nd), fmt(float(hi), nd)]
        lines.append("| " + " | ".join(row) + " |")
    print("\n".join(lines))
    print()
    print("| Flag | " + " | ".join(scen) + " |")
    print("|---|" + "---|" * len(scen))
    for k, label in [("miner_shortage", "miner shortage (hash under 50% of pre-event for 1 h or more)"),
                     ("backlog_600", "backlog over 600 s (design's backlog-rule trigger)"),
                     ("backlog_growing", "growing backlog (positive 10-day trend and over 60 s at day 30)"),
                     ("window_miss", "window miss (a day under 90% of blocks within 60 s)"),
                     ("oscillation", "oscillation (proving-share 10-90 range over 10 points in the last 10 days)")]:
        row = [label]
        for s in scen:
            n = sum(1 for m in out[s] if m[k])
            row.append(f"{n}/{len(out[s])}")
        print("| " + " | ".join(row) + " |")
    print()
    print("Operator profit by card class, $ per card-day (mean over seeds, all modes including off):")
    print()
    print("| Scenario | " + " | ".join(CLS) + " | shard share " + " | shard share ".join(CLS[:3]) + " |")
    print("|---|" + "---|" * 7)
    for s in scen:
        pr = np.mean([m["profit_card_day"] for m in out[s]], 0)
        sh = np.mean([m["shard_share"] for m in out[s]], 0)
        print(f"| {s} | " + " | ".join(f"{v:.3f}" for v in pr) + " | " + " | ".join(f"{v:.2f}" for v in sh[:3]) + " |")
    print()


def lever_study(scenario, seeds, p, with_sens=True, sens_only=False):
    levers = {
        "pool": [0.10, 0.20, 0.30, 0.40],
        "window": [5.0, 10.0, 20.0, 30.0],
        "burn": [0.0, 0.10, 0.25, 0.50],
        "timeout": [60.0, 120.0, 300.0, 600.0],
        "shards_per_block": [3.0, 30.0, 100.0, 300.0],
        "ema_ticks": [7.0, 20.0, 60.0, 480.0],
    }
    if sens_only:
        levers = {k: v for k, v in levers.items() if k in ("shards_per_block", "ema_ticks")}
    elif not with_sens:
        levers = {k: v for k, v in levers.items() if k not in ("shards_per_block", "ema_ticks")}
    print(f"Lever study on scenario {scenario} ({SCEN[scenario]}), seeds {list(seeds)}, traffic {p['shards_per_block']:.0f} shards per block. One parameter at a time, the rest at the design values.")
    print()
    print("Rows pool, window, burn, timeout are protocol levers; shards_per_block (chain traffic) and ema_ticks (operators' observation window, ticks of 180 s) are model sensitivities.")
    print()
    print("| Lever | Value | hash min / pre | hash d30 / pre | hours hash under 50% | age max s | within 60 s (mean) | worst day within 60 s | cards proving d30 | cards off d30 | external delivered | 3060 shard share | balance score |")
    print("|---|---|---|---|---|---|---|---|---|---|---|---|---|")
    rows = []
    for lever, vals in levers.items():
        for v in vals:
            pp = dict(p)
            pp[lever] = v
            out = run_scenarios([scenario], seeds, pp, quiet=True)[scenario]
            m = {k: agg(out, k)[0] for k in ["hash_min_ratio", "hash_d30_ratio", "hours_hash_lt50", "age_max", "ok60_mean", "ok60_min_day", "prove_d30", "off_d30", "ext_delivered"]}
            sh = np.mean([x["shard_share"] for x in out], 0)[2]
            # balance score: mean of (hash d30 ratio capped at 1), worst-day proof share, 1 - age_max/600 capped
            score = (min(1.0, m["hash_d30_ratio"]) + m["ok60_min_day"] + max(0.0, 1 - m["age_max"] / 600.0)) / 3.0
            star = " (design)" if abs(v - P[lever]) < 1e-9 and lever in ("pool", "window", "burn", "timeout") else (" (default)" if abs(v - P[lever]) < 1e-9 else "")
            print(f"| {lever} | {v}{star} | {m['hash_min_ratio']:.2f} | {m['hash_d30_ratio']:.2f} | {m['hours_hash_lt50']:.0f} | {m['age_max']:.0f} | {m['ok60_mean']:.3f} | {m['ok60_min_day']:.3f} | {m['prove_d30']:.3f} | {m['off_d30']:.3f} | {m['ext_delivered']:.2f} | {sh:.2f} | {score:.3f} |")
            rows.append((lever, v, score))
            sys.stdout.flush()
    print()
    return rows


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--scenarios", default="a,b,c,d,e,f")
    ap.add_argument("--seeds", type=int, default=5)
    ap.add_argument("--seed0", type=int, default=1)
    ap.add_argument("--levers", default=None, help="scenario letter for the lever study")
    ap.add_argument("--sens", action="store_true", help="with --levers: only the two sensitivities (traffic, observation window)")
    ap.add_argument("--nosens", action="store_true", help="with --levers: protocol levers only")
    ap.add_argument("--set", default="", help="k=v,k=v parameter overrides")
    ap.add_argument("--dump", default=None, help="scenario letter: print the hourly CSV for the first seed")
    ap.add_argument("--days", type=int, default=None)
    args = ap.parse_args()
    p = dict(P)
    for kv in filter(None, args.set.split(",")):
        k, v = kv.split("=")
        p[k] = float(v)
    if args.days:
        p["days"] = args.days
    seeds = range(args.seed0, args.seed0 + args.seeds)
    if args.dump:
        sim = Sim(args.dump, args.seed0, p)
        h = sim.run()
        print("hour," + ",".join(COLS))
        for i, row in enumerate(h):
            print(f"{i}," + ",".join(f"{v:.6g}" for v in row))
        m = metrics(h, sim)
        for k, v in m.items():
            print(f"# {k}: {v}", file=sys.stderr)
        return
    if args.levers:
        lever_study(args.levers, seeds, p, with_sens=not args.nosens, sens_only=args.sens)
        return
    scen = args.scenarios.split(",")
    print(f"# Igneum economy simulation, {p['n_ops']} operators, {p['days']} days, seeds {list(seeds)}, tick {p['tick']:.0f} s")
    print()
    print("Parameters: " + ", ".join(f"{k}={v}" for k, v in p.items()))
    print()
    for s in scen:
        print(f"- {s}: {SCEN[s]}")
    print()
    t0 = time.time()
    out = run_scenarios(scen, seeds, p)
    table_scenarios(out)
    print(f"<!-- total {time.time() - t0:.0f} s -->")


if __name__ == "__main__":
    main()