#!/usr/bin/env python3
"""
Limenka MTP-cadence simulation - core discrete-event model.

Models the settled single-lane 600s fork design:

  e       = nTime(i) - nTime(i-1), direct stamp delta (MTP-101 is used
            for the activation gate only, not for e)
  reward  = R_full * e / 600          (uncapped, floor zero)
  payload = S_max * e / 600           (uncapped; header+coinbase subtracted)
  DAA     = PI controller, error = 600 - e, per-block update.
            P (proportional) = kp on the instantaneous error, symmetric.
            I (integral)     = ki on the EMA accumulator, asymmetric
                               raise-slow/lower-fast for flood resistance.
            D (derivative)   = kd on the band-passed error difference -
                               measured dead in this system, kd = 0.
  time    = fork blocks stamped <= now + F (default 60s), honest NTP
            assumed: no stamp may exceed arrival + F

Attackers may stamp blocks up to arrival + F (reward inflation via time
fakery bounded by F) and self-pad blocks to the full payload allowance
with zero-cost self-pay txs, so realized attacker block size equals the
allowance, not the mempool backlog.

Scenarios: steady, ramp_slow, ramp_fast, collapse, flood, monopoly,
faketime, stall (see scenarios()).

Metrics: mean/median/p90 interval, issuance vs parity, difficulty
max/final, stall window after flood exit, reward split, backlog stats,
realized max block size.

Run:  python3 mtp_sim.py --scenario flood --blocks 30000
"""
import argparse
import bisect
import random

T = 600                 # target interval, seconds
S_MAX = 4_000_000       # standard payload bytes (witness discount removed)
R_FULL = 1.0            # normalized full subsidy: 1 per 600s

PID_SCALE = 1_000_000
DEN_MIN = 100_000       # per-block target move clamp: 0.1x
DEN_MAX = 10_000_000    # per-block target move clamp: 10x


class World:
    """Piecewise-linear hashrate schedule. const(t, h, a), ramp(t0, t1, h0, h1).

    Segments carry the attacker's stamp mode: 'honest' (stamp at arrival),
    'future' (stamp at arrival+F, max allowed), 'past' (stamp at MTP floor,
    the earliest valid - can be ~window/2*T behind wall clock)."""

    def __init__(self):
        self.segs = []

    def const(self, start, h, a=0.0, mode="honest"):
        self.segs.append((start, None, h, h, a, mode))

    def ramp(self, start, end, h0, h1, a=0.0, mode="honest"):
        self.segs.append((start, end, h0, h1, a, mode))

    def rates_at(self, t):
        h = a = 0.0
        mode = "honest"
        for start, end, h0, h1, aa, md in self.segs:
            if t >= start:
                h, a, mode = h0, aa, md
                if end is None:
                    continue
                if t < end:
                    frac = (t - start) / (end - start)
                    h = h0 + (h1 - h0) * frac
                else:
                    h = h1
        return h, a, mode

    def next_boundary(self, t):
        """Next schedule change strictly after t, or None."""
        nxt = None
        for start, end, _, _, _, _ in self.segs:
            for cand in (start, end):
                if cand is None:
                    continue
                if cand > t:
                    nxt = cand if nxt is None else min(nxt, cand)
        return nxt


class Chain:
    def __init__(self, cfg):
        self.cfg = cfg
        self.d = 1.0                    # difficulty multiplier (1 = parent)
        self.t = 0.0                    # wall clock, seconds
        self.mtp_sorted = []            # sorted stamps of last W blocks
        self.mtp_deque = []
        self.avg_error = 0.0            # EMA of error (P accumulator)
        self.err_lp = 0.0              # low-passed error (D state)
        self.err_lp_prev = 0.0         # previous low-passed error (D state)
        self.last_block_t = 0.0         # completion time of last block
        self.height = 0
        self.issuance = 0.0
        self.backlog = 0.0              # mempool backlog, bytes
        self.max_backlog = 0.0
        self.max_block_size = 0.0
        self.honest_reward = 0.0
        self.attacker_reward = 0.0
        self.blocks = []                # (arrival, attacker, e, d, size)
        self.stamp_violations = 0
        self.orphans = 0                # blocks lost to propagation races
        # delay floor, asymmetric gains, propagation - optional cfg attrs
        # Naming (corrected from the swapped 2021 lineage):
        #   P (proportional) = gain on the instantaneous error
        #   I (integral)     = gain on the EMA accumulator (avg_error)
        #   D (derivative)   = gain on the band-passed error difference
        self.delay = getattr(cfg, "delay", 0.0)
        self.att_hw = getattr(cfg, "att_hw", 1.0)
        self.prop_base = getattr(cfg, "prop_base", 0.0)
        self.bandwidth = getattr(cfg, "bandwidth", 100e6)
        self.honest_miners = getattr(cfg, "honest_miners", 1)
        # DAA setpoint, separate from the reward time basis T (600s).
        # Calibrated above 600 to cancel the asymmetric-integral
        # rectification bias (the ~17s steady dilation).
        self.target = getattr(cfg, "target", T)
        self.kp = cfg.kp
        self.ki_up = getattr(cfg, "ki_up", cfg.ki)
        self.ki_down = getattr(cfg, "ki_down", cfg.ki)

    def mtp(self):
        n = len(self.mtp_sorted)
        return self.mtp_sorted[(n - 1) // 2] if n else 0

    def advance_time(self, t_new):
        dt = t_new - self.t
        if dt > 0:
            self.backlog += self.cfg.fill * dt
            self.max_backlog = max(self.max_backlog, self.backlog)
            self.t = t_new

    def push_stamp(self, s):
        bisect.insort(self.mtp_sorted, s)
        self.mtp_deque.append(s)
        if len(self.mtp_deque) > self.cfg.window:
            old = self.mtp_deque.pop(0)
            del self.mtp_sorted[bisect.bisect_left(self.mtp_sorted, old)]

    def add_block(self, arrival, attacker, seed=False, mode="honest"):
        mtp_prev = self.mtp()
        last_stamp = self.mtp_deque[-1] if self.mtp_deque else arrival
        cap = arrival + self.cfg.future
        if self.cfg.monotonic:
            # strict stamp monotonicity: nTime > prev nTime (fork rule)
            floor = last_stamp + 1
        else:
            # parent-style rule: nTime > MTP(prev) only
            floor = mtp_prev + 1
        if mode == "past":
            # stamp at the earliest consensus-valid time, possibly far in
            # the past relative to wall clock (MTP floor only)
            stamp = floor
        elif mode == "future":
            # stamp as far ahead as the future limit allows
            stamp = min(cap, max(floor, arrival))
        else:
            # honest: stamp at real arrival
            stamp = min(cap, max(arrival, floor))
        if cap < floor:
            # both rules cannot be satisfied at `arrival`: the miner waits
            # for the clock to catch up, then stamps at the floor
            self.stamp_violations += 1
            arrival = floor - self.cfg.future
            cap = arrival + self.cfg.future
            stamp = floor

        if seed:
            # first block: continue parent difficulty, no DAA update
            self.advance_time(arrival)
            self.last_block_t = arrival
            self.push_stamp(stamp)
            self.height += 1
            return

        self.advance_time(arrival)
        self.last_block_t = arrival
        self.push_stamp(stamp)
        if self.cfg.ebasis == "mtp":
            e = max(1, self.mtp() - mtp_prev)
        else:  # "direct": stamp-to-stamp delta of consecutive blocks
            e = max(1, stamp - last_stamp)

        reward = R_FULL * e / T
        allowance = S_MAX * e / T
        size = allowance if attacker else min(allowance, self.backlog)
        self.backlog = max(0.0, self.backlog - size)
        self.max_block_size = max(self.max_block_size, size)

        self.issuance += reward
        if attacker:
            self.attacker_reward += reward
        else:
            self.honest_reward += reward

        error = self.target - e
        self.avg_error = self.avg_error * (1.0 - self.cfg.alpha) + error * self.cfg.alpha
        # band-pass derivative: difference of the low-passed error.
        # The low-pass kills high-frequency noise; the differencing
        # removes the DC component - a band-pass in the standard PID
        # filtered-derivative sense. (Measured dead in this system - the
        # single-sample exponential measurement has SNR 1, the derivative
        # carries ~3-5x more noise than signal. kd stays 0.)
        self.err_lp_prev = self.err_lp
        self.err_lp = self.err_lp * (1.0 - self.cfg.lp) + error * self.cfg.lp
        deriv = self.err_lp - self.err_lp_prev
        # P: proportional on the instantaneous error (symmetric).
        # I: integral on the EMA accumulator, asymmetric raise-slow /
        # lower-fast for flood resistance.
        ki = self.ki_up if error > 0 else self.ki_down
        correction = (ki * self.avg_error + self.kp * error +
                      self.cfg.kd * deriv) / PID_SCALE
        denom = max(DEN_MIN, min(DEN_MAX, PID_SCALE + correction))
        self.d *= denom / PID_SCALE

        self.height += 1
        self.blocks.append((arrival, attacker, e, self.d, size))


def scenarios(cfg):
    w = World()
    if cfg.scenario == "steady":
        w.const(0, 1.0)
    elif cfg.scenario == "ramp_slow":
        w.ramp(0, 30 * 86400, 1.0, 2.0)
    elif cfg.scenario in ("ramp_fast", "ramp100", "ramp1000", "ramp1e6"):
        w.ramp(0, 86400, 1.0, cfg.ramp_x)
    elif cfg.scenario in ("collapse", "collapse99"):
        w.const(0, 1.0)
        w.const(100 * T, 1.0 - cfg.drop)
    elif cfg.scenario == "flood":
        w.const(0, 1.0)
        w.const(100 * T, 1.0, cfg.flood_x)
        w.const(100 * T + cfg.flood_min * 60, 1.0)
    elif cfg.scenario == "monopoly":
        w.const(0, 0.1, 0.9)
    elif cfg.scenario in ("faketime", "future_attack", "future_minority",
                          "future_majority"):
        share = cfg.attack_share if cfg.scenario != "faketime" else 0.6
        w.const(0, 1.0 - share, share, mode="future")
    elif cfg.scenario in ("past_attack", "past_minority", "past_majority"):
        share = cfg.attack_share
        w.const(0, 1.0 - share, share, mode="past")
    elif cfg.scenario == "stall":
        w.const(0, 1.0)
        w.const(100 * T, cfg.stall_frac)
        w.const(100 * T + cfg.stall_days * 86400, 1.0)
    else:
        raise SystemExit(f"unknown scenario {cfg.scenario}")
    return w


def flood_exit(cfg):
    if cfg.scenario == "flood":
        return 100 * T + cfg.flood_min * 60
    return None


def run(cfg):
    """Run the simulation. Returns a dict of metrics."""
    rng = random.Random(cfg.seed)
    world = scenarios(cfg)
    chain = Chain(cfg)
    exit_t = flood_exit(cfg)

    chain.add_block(0.0, False, seed=True)
    for _ in range(cfg.blocks):
        h, a, mode = world.rates_at(chain.t)
        total = h + a
        if total <= 0:
            break
        dt = rng.expovariate(total / (T * chain.d))
        arr = chain.t + dt
        nxt = world.next_boundary(chain.t)
        if nxt is not None and arr >= nxt:
            chain.advance_time(nxt)  # rate changes before any block arrives
            continue
        attacker = rng.random() < a / total
        hw = chain.att_hw if attacker else 1.0
        # sequential delay floor: block can't complete before prev + D*hw
        min_arr = chain.last_block_t + chain.delay * hw
        if arr < min_arr:
            arr = min_arr
        if nxt is not None and arr >= nxt:
            chain.advance_time(nxt)
            continue
        # propagation race: competing miners finishing within our
        # propagation window orphan our block. The delay floor serializes
        # block production when D > propagation, so competitors cannot
        # exist inside the window (compet_window -> 0).
        if chain.prop_base > 0:
            est_e = max(1.0, arr - chain.last_block_t)
            est_size = (S_MAX * est_e / T) if attacker else min(S_MAX * est_e / T, chain.backlog)
            prop = chain.prop_base + est_size / chain.bandwidth
            hw_other = chain.att_hw if not attacker else 1.0
            # competing hashrate: the rest of the honest population (split
            # into N miners) plus the attacker side
            if attacker:
                other_rate = h / (T * chain.d)
            else:
                other_rate = (h * (chain.honest_miners - 1) / chain.honest_miners + a) / (T * chain.d)
            # competitor completion spacing: PoW mean, floored by their
            # delay (delay-bound competitors serialize at the floor)
            spacing = max((1.0 / other_rate) if other_rate > 0 else float("inf"),
                          chain.delay * hw_other)
            p_orphan = min(1.0, 2.0 * prop / spacing) * 0.5
            if rng.random() < p_orphan:
                chain.orphans += 1
                continue
        chain.add_block(arr, attacker, mode=mode if attacker else "honest")

    ints = sorted(b[2] for b in chain.blocks if b[2] > 0)
    n = max(1, len(ints))
    parity = R_FULL * chain.t / T
    res = {
        "scenario": cfg.scenario,
        "blocks": chain.height,
        "sim_hours": chain.t / 3600,
        "mean_int": sum(ints) / n,
        "med_int": ints[n // 2],
        "p90_int": ints[int(n * 0.9)],
        "issuance_ratio": chain.issuance / parity if parity > 0 else 0.0,
        "d_max": max(b[3] for b in chain.blocks) if chain.blocks else 1.0,
        "d_final": chain.d,
        "max_backlog": chain.max_backlog,
        "max_block_size": chain.max_block_size,
        "honest_reward": chain.honest_reward,
        "attacker_reward": chain.attacker_reward,
        "stamp_violations": chain.stamp_violations,
        "stall_first_blk": float("nan"),
        "stall_restore": float("nan"),
    }

    if exit_t is not None:
        post = [b for b in chain.blocks if b[0] >= exit_t and not b[1]]
        if post:
            res["stall_first_blk"] = post[0][0] - exit_t
            good = 0
            for b in post:
                if abs(b[2] - T) <= 0.5 * T:
                    good += 1
                else:
                    good = 0
                if good >= 3:
                    res["stall_restore"] = b[0] - exit_t
                    break
    return res


def main():
    ap = argparse.ArgumentParser(description="Limenka MTP-cadence sim")
    ap.add_argument("--scenario", default="steady")
    ap.add_argument("--blocks", type=int, default=30000)
    ap.add_argument("--seed", type=int, default=42)
    ap.add_argument("--window", type=int, default=101)
    ap.add_argument("--future", type=int, default=60)
    ap.add_argument("--kp", type=float, default=50e6)
    ap.add_argument("--ki", type=float, default=10e6)
    ap.add_argument("--kd", type=float, default=16.7e6)
    ap.add_argument("--alpha", type=float, default=0.10)
    ap.add_argument("--lp", type=float, default=0.15)
    ap.add_argument("--fill", type=float, default=6600)
    ap.add_argument("--ebasis", choices=["mtp", "direct"], default="mtp")
    ap.add_argument("--drop", type=float, default=0.9)
    ap.add_argument("--flood-x", type=float, default=50)
    ap.add_argument("--flood-min", type=float, default=30)
    ap.add_argument("--stall-frac", type=float, default=0.001)
    ap.add_argument("--stall-days", type=float, default=4)
    ap.add_argument("--ramp-x", type=float, default=10)
    ap.add_argument("--attack-share", type=float, default=0.2)
    ap.add_argument("--monotonic", action="store_true")
    cfg = ap.parse_args()

    res = run(cfg)
    for k, v in res.items():
        if isinstance(v, float):
            print(f"{k:18s} {v:,.3f}")
        else:
            print(f"{k:18s} {v}")


if __name__ == "__main__":
    main()