#! /usr/bn/python3

def simulate_halving():
  H = 90;       # Total active hashpower, in EH/s.
  R0 = 1.0/600; # Target bblock rate, blk/s.
  D = H/R0;     # Difficulty at halving, in EH/blk.
  P = 7000;     # Price of bitcoin, in USD/BTC.
  C = 2400;     # Block capacity, in tx/blk.
  fee = initial_fee_distr(C*R); # Distr of fees among incoming txs.
  mem = initial_mempool_state();    # Distr of txs per fee in mempool.
  mst = initial_miner_state(D, 12.5, P, 0.9); # Distrib of miner effs.
  t = 0;        # Time of last block with 12.5 BTC reward, in secs.
  t0 = t;       # Time at last diffculty adjustment.
  epr = 0.05;   # Prob of solving a block in next integration step, in blk.
  for b in range(14400):
    # Simulate situation after {d} blocks mined since halving:
    if b > 0 and (b %% 2100) == 0:
      # Adjust difficulty
      Ract = 2100/(t - t0);  
      adj = max(0.25, min(4.0, R0/Ract)); 
      D = D*adj
      t0 = t
    # Simulate mining:
    X = 0; # Total miners cost for the next block.
    nt = int(log(1.0e-)/log(1.0 - epr) + 1); # Number of integration steps.
    for k in range(nt):
      # Compute block rate{R} at current active hashpower:
      R = H/D;
      # Choose a time step {dt} with about {epr} probability of block solving:
      dt = epr/R;
      # Advance time:
      t = t + dt;
      # Update incoming fee distribution and traffic {T} in that time step:
      T = update_incoming_fee_distr(fee,mem,C)
      # Update mempool state and compute tot fees of next block:
      F = update_mempool_state(mem,fee,dt,C)
      # Update miner state and active hashpower:
      H = update_miners_state(mst,D,P,F,dt)
      if boolrand(epr):
        # Pretend that block was solved
        break
    # Reward miners, distributing the reward among them:
    reward_miners(mst,P,F)
    # Remove the solved block from the mempool:
    remove_solved_block(mem,C)
  return
  # ----------------------------------------------------------------------
  
def update_incoming_fee_distr(fee,mem,C,R):
  # Updates the distribution {fee} of transaction 
  # fees in the incoming traffic, simulating the 
  # users' reaction to the mempool state {mem}.
  # Returns as result the short-term average incoming
  # traffic {T} in  tx/s.  Assumes that the block capacity
  # is {C} (in txs), and the recent block rate is {R}
  # (in blk/s).
  #
  # The parameter {fee} is a dictionary object.
  # {fee.hist} is a list; {fee.hist[k]} is the 
  # incoming rate (in tx/s) of transactions that pay
  # a fee of {class_fee(k)} USD of fee.
  #
  # This procedure assumes that, depending on the 
  # size of the backlog and the recent block rate, 
  # some users will raise their fees in an attempt 
  # to get through, while others
  # will just give up.  The latter are adjusted 
  # so that the backlog will not grow much beyond
  # some maximum size.
  #
  
  # Compute the total backlog size:
  M = mempool_size(mem);

  # Estimate time {tc} for backlog to clear if traffic 
  # were to drop to 80% of capacity {C/R}immediately:
  tc = ceil(M/(0.2*C) + 0.001)/R

  # Estimate the total transaction rate based on capacity and 
  # current backlog size:
  Mmax = 500*C;  # Max mempool size, in tx. 
  T = ???

  # Compute the fee distrbution based on mempool fee distribution
  # for an arbitrary incoming traffic:
  
  hist = fee.hist;
  nh = len(hist)
  k = 0; # Next entry of histogram.
  while True:
    f = class_fee(k); # Fee paid by transactions in this group.
    # Estimate how many tx/s will pay fee {f}:
    t = ??? f  M  C R
    if t < 1.0e-6: t = 0; # Truncate distrib where too small.
    if k >= nh and t == 0:
      break
    if k < len(h):
      hist[k] = t
    else
      hist.append(t); nh = nh + 1
    k = k + 1
  
  # Scale ditribution to total traffic {T}:
  sum = 1.0e-100; for t in hist: sum = sum + t;
  scale = T/sum;
  for k in range(hist): hist[k] = hist[k]*scale
  return T
  # ----------------------------------------------------------------------
  
def update_mempool_state(mem,fee,dt,C):
  # Updates the mempool state {mem} to account for the
  # traffic with the fee distribution {fee} arriving
  # for {dt} seconds.
  #
  # The mempool state is a dictionary object.
  # The element {mem.hist} is a list.
  # {mem.hist[k]} is how many transactions are waiting to 
  # be processed that pay fee {class_fee(k)}.  The 
  # number may be fractional, to compensate
  # for the quantization of the fee.
  #
  # Returns the total fees {F} that would be paid by the 
  # transactions in the next block, if it was solved 
  # at the end of that time interval, assuming that its 
  # capacity is {C}.
  
  nm = len(mem.hist);
  nf = len(fee.hist);
  for k in range(max(nm,nf)):
    f = class_fee(k);
    m = mem.hist[k] if k < nm else 0; # Mempool txs that pay {f}.
    t = fee.hist[k] if k < nf else 0; # Incoming tx/s that pay {f}.
    m = m + t*dt;
    if k < nm:
      mem.hist[k] = m
    else:
      mem.hist,append(m); nm = nm + 1
  
  # Simulate the filling of the next block:
  tB = 0; # Transactions put in block so far (maybe fractional).
  k = nm-1; # Next fee class to include in block.
  F = 0; # Total fees in block.
  while k >= 0 and tB + mem.hist[k] <= C:
    tk = mem.hist[k]
    tB = tB + tk
    F = F + tk*class_fee(k)
    k = k - 1
  if k >= 0:
    # Top up block with some of fee class {k}:
    tk = C - tB
    tB = C
    F = F + tk*class_fee(k)
  return F
  # ----------------------------------------------------------------------
  
def class_fee(k):
  # The mean transaction fee (in USD/tx) of entry {k} in 
  # mempool and fee distribution histograms.
  
  return 0.005 + 0.10*k
  # ----------------------------------------------------------------------