#! /usr/bin/python
# Last edited on 2017-10-06 16:05:56 by jstolfi
# -*- coding: latin-1 -*-

import string;
import sys;
import random;
import optparse;
import copy;
from math import log, exp, sin, sqrt, hypot, pi;

## parse switches
usage = '''usage: %prog [options]

:: Simulates DAA algorithms for two competing chains'''
p = optparse.OptionParser(usage)
p.add_option('--someOption', dest='someVariable', action='store_true',
       help='Bla bla bla');
p.add_option('--anotherOption', dest='anotheVariable', action='store',
             help='Bla bla bla too');
p.set_defaults(someVariable=False, anotherVariable=100);

# In all functions below: 
#
#  {nc} is the number of chains, numbered {0..nc-1}.
#
#  {it} is the current clock time, in seconds since the start
#    of the simulation.
#
#  {dif[ic]} is the current difficulty of chain {ic}, defined as the 
#    expected amount of work (exahashes, EH) that it would take to solve a block
#    on that chain.
# 
#  {rev[ic]} is the current revenue that a miner would collect by solving the next
#    block in chain {ic}, measured in USD.  It includes
#    the block reward plus available fees, multiplied by the current
#    price of the coin.
#
#  {hashpw[ic]} is the hashpower (in EH/s) that 
#    is currently mining chain {ic}.
#
#  {bkclocks[ic][0..]} are the actual clock times (seconds) 
#    of the blocks solved so far on chain {ic}.
#
#  {bkstamps[ic][0..]} are the timestamps (seconds) 
#    of those blocks.
#
#  {bkdifs[ic][0..]} are the difficulties en force
#    when those blocks were solved.

def multichain_DAA_simulate(nc, dif, hashpw, nt, tstep):
  "Simulates DAA algorithms for two or more chains that" \
  " compete for the same set of miners.  The starting" \
  " difficulty of chain {ic} is {dif[ic]} and its" \
  " initial hashpower (EH/s) is {hashpw[ic]}.  The" \
  " simulation runs for {nt} steps of {tstep} seconds."
  
  # Initialize the revenue vector and lists of block times and timestamps:
  bkclocks = [None for ic in range(nc)];
  bkstamps = [None for ic in range(nc)];
  bkdifs = [None for ic in range(nc)];
  for ic in range(nc):
    bkclocks[ic] = copy.copy([]);
    bkstamps[ic] = copy.copy([]);
    bkdifs[ic] = copy.copy([]);

  for it in range(nt):
    # Simulate one more second of real time, from discrete time {it} to
    # time {it+1}.
    
    t = it*tstep + 0.5; # Time (seconds) at middle of step.
        
    # Choose the block revenue for each chain:
    rev = choose_revenue_per_block(nc, t);
    
    # Get the hashpower in each chain:
    adjust_hashpower(hashpw, nc, dif, rev, bkstamps, bkdifs, t, tstep);
    
    # Simulate mining:
    for ic in range(nc):
      # Determine whether chain {ic} solves a block or not:
      # Assumes that the probability of solving is small and 
      # of two or more blocks in the same second is zero. 
      # If the block rate gets close to 1 per {tstep}, should 
      # either reduce {tstep} or generate a block count with a 
      # Poisson distribution.
      prob = tstep*hashpw[ic]/dif[ic]; # Expected num of blocks in time step.
      # sys.stderr.write('%10.2f chain %d: hashpw = %12.6f prob = %12.8f\n' % (t,ic,hashpw[ic],prob));
      solved = (random.random() < prob);
      
      if solved:
        height = len(bkclocks[ic]);
        sys.stderr.write('t: %10.2f chain: %d blk: %8d' % (t,ic,height));
        rph = rev[ic]/dif[ic]; # Revenue rate (USD/EH).
        sys.stderr.write(' hpw: %12.5f dif: %12.4f rph: %12.3f' % (hashpw[ic],dif[ic],rph));
        # Append the block to the history of chain {ic}:
        bkdifs[ic].append(dif[ic]);
        bkclocks[ic].append(t);
        timestamp = choose_timestamp(t, ic, nc, dif, rev, bkstamps, bkdifs);
        bkstamps[ic].append(timestamp);
        # Adjust the difficulty of chain {ic}:
        dif[ic] = adjust_difficulty(ic, bkstamps[ic], bkdifs[ic]);
        sys.stderr.write(' next: %12.4f\n' % dif[ic]);

  for ic in range(nc):
    nb = len(bkstamps[ic]);
    sys.stderr.write('chain %d: %d blocks\n' % (ic,nb));
  # ------------------------------------------------------------
  
def adjust_hashpower(hashpw, nc, dif, rev, bkstamps, bkdifs, t, tstep):
  "Modifies the hashpower {hashpw[0..nc-1]} mining each chain" \
  " at about time {t}, given the" \
  " next-block difficulties {dif[0..nc-1]} and expected block revenues {rev[0..nc-1]}."
  
  nb_steady = 20; # Preserve the hashpower until we have this many blocks.
  
  # Find the index {ic_max} of the most profitable chain and its
  # expected revenue rate {rr_max} (USD per exahash):
  ic_max = -1; 
  rr_max = 0.0;
  for ic in range(nc):
    rr = rev[ic]/dif[ic]; 
    if (ic_max == -1) or (rr > rr_max):
      ic_max = ic; 
      rr_max = rr;
  assert rr_max > 0.0;
  
  # Compute the new hashpower for the losing chains:
  thash = 0.0; # Current total hashpower.
  tloser = 0.0; # New total hashpower of losing chains.
  for ic in range(nc):
    thash = thash + hashpw[ic];
    if (ic != ic_max):
      # Loser chain:
      nb = len(bkstamps[ic]);
      # Do not change hashpower for the first {nb_steady} blocks.
      if (nb >= nb_steady):
        # Revenue per EH in chain {ic}
        rr = rev[ic]/dif[ic];
        hashpw[ic] = hashpw[ic]*(1 - fleeing_hashpower_fraction(rr, rr_max, tstep)); 
      # Accumulate adjusted hashpower of the losing chains:
      tloser = tloser + hashpw[ic];
    
  # Recompute hashpower of winning chain:
  hashpw[ic_max] = thash - tloser;
  
  # Here should simulate variation in total hashpower
  # from time {t-tstep} to time {t}.  Assume constant for now.
  # -----------------------------------------------------------------
  
def fleeing_hashpower_fraction(rr, rr_max, tstep):
  "Guesses the fraction of the hashrate of a less-profitable chain that will" \
  " switch to the most profitable chain in {tstep} seconds, given the revenue" \
  " rates {rr} and {rr_max} (in USD/exahash) of the two chains."
  
  # The version below is rather pessimistc and assumes that 
  # a 10% difference between {rr} and {rr_max} will motivate
  # a fraction {mfrac} of the current miners to switch to
  # the most profitable chain, in each second.
  
  assert rr_max >= 1.0e-3;
  mfrac = 0.01;
  ratio = rr/rr_max;
  if ratio <= 1.0e-6:
    sys.stderr.write('%10.2f chain %d: !! revenue ratio is too small = %12.8f\n' % (t, ic, rr));
    ratio = 1.0e-6;
  elif ratio >= 1.0e+6:
    sys.stderr.write('%10.2f chain %d: !! revenue ratio is too large = %12.8f\n' % (t, ic, rr));
    ratio = 1.0e+6;

  # Compute the hashpower reduction factor {hred} for this chain:
  z = log(ratio)/log(0.9);
  flee_frac = exp(z*tstep*log(mfrac));
  
  return flee_frac;
  # ----------------------------------------------------------------------

def choose_timestamp(t, ic, nc, dif, rev, bkstamps, bkdifs):
  "Defines the timestamp of a new block of chain {ic} that was solved" \
  " at time {t}.\n\n" \
  "Miners who wish to game the system by manipulating the timestamps" \
  " can use the data that they know about all chains: the" \
  " current difficulties {dif[0..nc-1]} and next-block revenues {rev[0..nc-1]} and" \
  " the difficulties {bkdifs[0..nc-1][0..]}" \
  " and the timestamps {bkstamps[0..nc-1][0..]} of all previously solved" \
  " blocks (excluding the one just solved).\n\n" \
  "Note that {dif[ic]} is the" \
  " difficulty of the just-solved block, while {dif[jc]} for {jc!=ic} is" \
  " the difficulty of solving the next block in chain {jc}."
        
  # This implementation assumes that miners have accurate clocks
  # and do not try to game thetimestamps:
  return t;
  # ----------------------------------------------------------------------
  
def adjust_difficulty(ic, bkstamps_ic, bkdifs_ic):
  "The difficulty adjustment algorithm of chain {ic}.\n\n" \
  "Assumes that a new block has just been solved on chain {ic}," \
  " and the timestamps and difficulties of all" \
  " the blocks solved so far on that chain are {bkstamps[_ic0..]} and {bkdifs_ic[0..]." 
    
  # Insert here the proper algorithms, e.g. Satoshi's DAA
  # for chain 0 and the proposed Bitcoin Cash DAA for chain 1.
  
  if ic == 0:
    return adjust_difficulty_BTC(bkstamps_ic, bkdifs_ic);
  elif ic == 1:
    return adjust_difficulty_BCH(bkstamps_ic, bkdifs_ic);
  else:
    assert False;
    
def adjust_difficulty_BTC(bkstamps_ic, bkdifs_ic):
  "The difficulty adjustment algorithm of the BTC chain.\n\n" \
  "Receives the difficulties {bkdifs_ic[0..]} and" \
  " the timestamps {bkstamps[0..]} of the previous blocks.\n\n" \
  "Should be Satoshi's original algorithm that runs every" \
  " 2016 blocks.  Assumes that the most recent difficulty" \
  " adjustment happened just before the start of the simulation."
  
  # !!! PLACEHOLDER: NO ADJUSTMENT EVER !!!
  nb = len(bkstamps_ic);
  assert nb >= 1;
  dif_last = bkdifs_ic[nb-1]; # Difficulty of last block.
  return dif_last;
  # ----------------------------------------------------------------------
  
def adjust_difficulty_BCH(bkstamps_ic, bkdifs_ic):
  "The (proposed) difficulty adjustment algorithm of the BCH chain.\n\n" \
  "Receives the difficulties {bkdifs_ic[0..]} and" \
  " the timestamps {bkstamps_ic[0..]} of the previous blocks."
  
  # !!! PLACEHOLDER: NAIVE ADJUSTMENT !!!
  
  target_time = 600; # Target time between blocks.
  nw = 6;            # Number of blocks in window.
  slow = 0.25;       # Mismatch reduction factor.
  
  nb = len(bkstamps_ic);
  assert nb >= 1;
  dif_last = bkdifs_ic[nb-1]; # Difficulty of last block.
  
  if nb < nw:
    # Preserve the difficulty until we have enough blocks:
    return bkdifs_ic[nb-1];
    
  # Compute the total expected work (exahashes) that was expended to generate the
  # last {nw-1} blocks:
  tot_work = 0;
  for ib in range(nw-1):
    # The workexpected to generate a block is 
    tot_work = tot_work + bkdifs_ic[nb-1-ib];

  # Get the total time {tot_time} that was apparently needed
  # to create those blocks:
  tot_time = bkstamps_ic[nb-1] - bkstamps_ic[nb-nw];
  if tot_time <= 0.0:
    # Probably gamed timestamps:
    tot_time = 0.0;
  # Fudge {tot_time} to be away from 0:
  tot_time = hypot(1.0, tot_time);

  # Assuming that the work {tot_work} produced the 
  # last {nw-1} blocks in time {tot_time}, computes the
  # difficulty (work per block) {dif_ideal}
  # that would give the desired {target_time}:
  dif_ideal = tot_work*target_time/tot_time;

  # Adjust {dif} exponentially towards {dif_ideal}:
  dif_new = dif_last*exp(slow*log(dif_ideal/dif_last))
  
  return dif_new;
  # ----------------------------------------------------------------------
  
def choose_revenue_per_block(nc, t):
  "Returns the expected revenue that a miner will get if he" \
  " solves a block on each chain {ic} in {0..nc-1} around time {t}."
  
  # Amplitude and period of price fluctuations of chain 1:
  amp1 = 200;
  per1 = 3600;
  
  rev = [0 for ic in range(nc)];
  rev[0] = 50000.0; # Assume that price is constant on chain 0.
  rev[1] = 4000 + amp1*sin(2*pi*t/per1);
  return rev;
  # ----------------------------------------------------------------------

# MAIN

# Get arguments:
(opts, args) = p.parse_args()

# Initial conditions:
nc = 2;
hashpw = [9.0, 1.0];   # Hashpower on each chain (EH/s).
dif = [600*hashpw[0], 600*hashpw[1]];  # Expected time to solve 1 block with 1 EH/s.

# Time interval:
tstep = 1.0;   # Time step (seconds)
nt = 100000;   # Number of steps.

multichain_DAA_simulate(nc, dif, hashpw, nt, tstep); 

sys.stderr.write('done.')