pi_compare.py raw

   1  #!/usr/bin/env python3
   2  """
   3  I vs PI comparison - clean, controlled, same seeds.
   4  
   5  Naming corrected from the swapped 2021 lineage:
   6    P (proportional) = gain on the instantaneous error (kp)
   7    I (integral)     = gain on the EMA accumulator (ki)
   8    D (derivative)   = gain on the band-passed error difference (kd)
   9  
  10  The claim under test: the proportional kick (P) added to the integral
  11  accumulator (I) tracks hashpower changes faster. (Previously reported
  12  as "P vs PI" under the swapped labels - the physics is unchanged, the
  13  labels are corrected.)
  14  
  15  Four configs, same seeds, D=60s delay floor, monotonic stamps, direct e:
  16  
  17    sym-I    kp=0,  ki=10M                    (integral only)
  18    sym-PI   kp=10M, ki=10M                   (symmetric P+I)
  19    asym-I   kp=0,  ki_up=5M, ki_down=50M     (asymmetric integral only)
  20    asym-PI  kp=10M, ki_up=5M, ki_down=50M    (asymmetric P+I - final design)
  21  
  22  Metrics:
  23    steady: mean interval, d wander (log-std of d - oscillation amplitude)
  24    step 3x at 20d: t50 (hours to reach d=2), overshoot (max d)
  25    collapse -90%: t50-down (hours for d to fall halfway to 0.1)
  26  
  27  Run:  python3 pi_compare.py
  28  """
  29  import math
  30  import os
  31  import random
  32  import sys
  33  
  34  sys.path.insert(0, os.path.dirname(__file__))
  35  import mtp_sim as M
  36  
  37  
  38  class Cfg:
  39      def __init__(self, kp, ki_up, ki_down, seed=7, blocks=9000):
  40          self.window = 101
  41          self.future = 60
  42          self.kp = kp
  43          self.ki = ki_up
  44          self.kd = 0
  45          self.ki_up = ki_up
  46          self.ki_down = ki_down
  47          self.alpha = 0.05
  48          self.lp = 0.15
  49          self.fill = 6600
  50          self.ebasis = "direct"
  51          self.monotonic = True
  52          self.delay = 60.0
  53          self.att_hw = 1.0
  54          self.prop_base = 0.0
  55          self.bandwidth = 100e6
  56          self.honest_miners = 1
  57          self.scenario = "steady"
  58          self.seed = seed
  59          self.blocks = blocks
  60          self.drop = 0.9
  61          self.ramp_x = 10.0
  62          self.flood_x = 50
  63          self.flood_min = 30
  64          self.stall_frac = 0.001
  65          self.stall_days = 4
  66          self.attack_share = 0.6
  67  
  68  
  69  def simulate(cfg, step_t=None, collapse_t=None):
  70      """Run the chain. Returns (chain, step_t50, step_overshoot, coll_t50)."""
  71      rng = random.Random(cfg.seed)
  72      world = M.World()
  73      world.const(0, 1.0)
  74      if step_t is not None:
  75          world.const(step_t, 3.0)
  76      if collapse_t is not None:
  77          world.const(collapse_t, 1.0 - cfg.drop)
  78      chain = M.Chain(cfg)
  79      chain.add_block(0.0, False, seed=True)
  80      step_t50 = None
  81      step_over = 1.0
  82      coll_t50 = None
  83      coll_target = 0.5 + (1.0 - cfg.drop) / 2.0
  84      for _ in range(cfg.blocks):
  85          h, a, mode = world.rates_at(chain.t)
  86          total = h + a
  87          if total <= 0:
  88              break
  89          dt = rng.expovariate(total / (M.T * chain.d))
  90          arr = chain.t + dt
  91          nxt = world.next_boundary(chain.t)
  92          if nxt is not None and arr >= nxt:
  93              chain.advance_time(nxt)
  94              continue
  95          hw = chain.att_hw if (a > 0 and rng.random() < a / total) else 1.0
  96          attacker = hw != 1.0
  97          min_arr = chain.last_block_t + chain.delay * hw
  98          if arr < min_arr:
  99              arr = min_arr
 100          if nxt is not None and arr >= nxt:
 101              chain.advance_time(nxt)
 102              continue
 103          chain.add_block(arr, attacker)
 104          if step_t is not None and chain.t > step_t:
 105              step_over = max(step_over, chain.d)
 106              if step_t50 is None and chain.d >= 2.0:
 107                  step_t50 = (chain.t - step_t) / 3600
 108          if collapse_t is not None and chain.t > collapse_t and coll_t50 is None:
 109              if chain.d <= coll_target:
 110                  coll_t50 = (chain.t - collapse_t) / 3600
 111      return chain, step_t50, step_over, coll_t50
 112  
 113  
 114  def steady_metrics(cfg):
 115      chain, _, _, _ = simulate(cfg)
 116      ints = sorted(b[2] for b in chain.blocks if b[2] > 0)
 117      n = max(1, len(ints))
 118      dlogs = [math.log(b[3]) for b in chain.blocks[cfg.blocks // 3:]]
 119      mean = sum(dlogs) / len(dlogs)
 120      var = sum((x - mean) ** 2 for x in dlogs) / len(dlogs)
 121      return sum(ints) / n, math.sqrt(var), chain.d
 122  
 123  
 124  def report(name, kp, ki_up, ki_down, seeds=(7, 8)):
 125      rows = []
 126      for seed in seeds:
 127          cfg = Cfg(kp, ki_up, ki_down, seed=seed)
 128          mean_int, wander, d_fin = steady_metrics(cfg)
 129          cfg2 = Cfg(kp, ki_up, ki_down, seed=seed)
 130          _, t50, over, _ = simulate(cfg2, step_t=20 * 86400)
 131          cfg3 = Cfg(kp, ki_up, ki_down, seed=seed, blocks=12000)
 132          _, _, _, coll = simulate(cfg3, collapse_t=100 * M.T)
 133          rows.append((seed, mean_int, wander, d_fin, t50, over, coll))
 134      print(f"{name:10s} | mean_int  wander  d_fin | step_t50  overshoot | coll_t50")
 135      for seed, mean_int, wander, d_fin, t50, over, coll in rows:
 136          t50s = f"{t50:7.1f}h" if t50 else "    -"
 137          overs = f"{over:9.2f}x" if over else "     -"
 138          colls = f"{coll:7.1f}h" if coll else "    -"
 139          print(f"  seed {seed}  | {mean_int:7.0f}s {wander:7.3f} {d_fin:6.2f} | "
 140                f"{t50s} {overs} | {colls}")
 141      print()
 142  
 143  
 144  def main():
 145      print("I vs PI, same seeds, D=60s delay floor, monotonic, direct e")
 146      print("(labels corrected: P = instantaneous error gain, I = EMA gain)\n")
 147      report("sym-I", 0, 10e6, 10e6)
 148      report("sym-PI", 10e6, 10e6, 10e6)
 149      report("asym-I", 0, 5e6, 50e6)
 150      report("asym-PI", 10e6, 5e6, 50e6)
 151  
 152  
 153  if __name__ == "__main__":
 154      main()
 155