1  import math, random
  2  
  3  hashpower = [float(x) for x in open('hashpower.csv').readlines()]
  4  
  5  # Target block time
  6  TARGET = 12
  7  # Should be 86400, but can reduce for a quicker sim
  8  SECONDS_IN_DAY = 86400
  9  # Look at the 1/x day exponential moving average
 10  EMA_FACTOR = 0.01
 11  # Damping factor for simple difficulty adjustment
 12  SIMPLE_ADJUST_DAMPING_FACTOR = 20
 13  # Maximum per-block diff adjustment (as fraction of current diff)
 14  SIMPLE_ADJUST_MAX = 0.5
 15  # Damping factor for quadratic difficulty adjustment
 16  QUADRATIC_ADJUST_DAMPING_FACTOR = 3
 17  # Maximum per-block diff adjustment (as fraction of current diff)
 18  QUADRATIC_ADJUST_MAX = 0.5
 19  # Threshold for bounded adjustor
 20  BOUNDED_ADJUST_THRESHOLD = 1.3
 21  # Bounded adjustment factor
 22  BOUNDED_ADJUST_FACTOR = 0.01
 23  # How many blocks back to look
 24  BLKS_BACK = 10
 25  # Naive difficulty adjustment factor
 26  NAIVE_ADJUST_FACTOR = 1/1024.
 27  
 28  
 29  # Produces a value according to the exponential distribution; used
 30  # to determine the time until the next block given an average block
 31  # time of t
 32  def expdiff(t):
 33      return -math.log(random.random()) * t
 34  
 35  
 36  # abs_sqr(3) = 9, abs_sqr(-7) = -49, etc
 37  def abs_sqr(x):
 38      return -(x**2) if x < 0 else x**2
 39  
 40  
 41  # Given an array of the most recent timestamps, and the most recent
 42  # difficulties, compute the next difficulty
 43  def simple_adjust(timestamps, diffs):
 44      if len(timestamps) < BLKS_BACK + 2:
 45          return diffs[-1]
 46      # Total interval between previous block and block a bit further back
 47      delta = timestamps[-2] - timestamps[-2-BLKS_BACK] + 0.0
 48      # Expected interval
 49      expected = TARGET * BLKS_BACK
 50      # Compute adjustment factor
 51      fac = 1 - (delta / expected - 1) / SIMPLE_ADJUST_DAMPING_FACTOR
 52      fac = max(min(fac, 1 + SIMPLE_ADJUST_MAX), 1 - SIMPLE_ADJUST_MAX)
 53      return diffs[-1] * fac
 54  
 55  
 56  # Alternative adjustment algorithm
 57  def quadratic_adjust(timestamps, diffs):
 58      if len(timestamps) < BLKS_BACK + 2:
 59          return diffs[-1]
 60      # Total interval between previous block and block a bit further back
 61      delta = timestamps[-2] - timestamps[-2-BLKS_BACK] + 0.0
 62      # Expected interval
 63      expected = TARGET * BLKS_BACK
 64      # Compute adjustment factor
 65      fac = 1 - abs_sqr(delta / expected - 1) / QUADRATIC_ADJUST_DAMPING_FACTOR
 66      fac = max(min(fac, 1 + QUADRATIC_ADJUST_MAX), 1 - QUADRATIC_ADJUST_MAX)
 67      return diffs[-1] * fac
 68  
 69  
 70  # Alternative adjustment algorithm
 71  def bounded_adjust(timestamps, diffs):
 72      if len(timestamps) < BLKS_BACK + 2:
 73          return diffs[-1]
 74      # Total interval between previous block and block a bit further back
 75      delta = timestamps[-2] - timestamps[-2-BLKS_BACK] + 0.0
 76      # Expected interval
 77      expected = TARGET * BLKS_BACK
 78      if delta / expected > BOUNDED_ADJUST_THRESHOLD:
 79          fac = (1 - BOUNDED_ADJUST_FACTOR)
 80      elif delta / expected < 1 / BOUNDED_ADJUST_THRESHOLD:
 81          fac = (1 + BOUNDED_ADJUST_FACTOR) ** (delta / expected)
 82      else:
 83          fac = 1
 84      return diffs[-1] * fac
 85  
 86  
 87  # Old Ethereum algorithm
 88  def old_adjust(timestamps, diffs):
 89      if len(timestamps) < 2:
 90          return diffs[-1]
 91      delta = timestamps[-1] - timestamps[-2]
 92      expected = TARGET * 0.693
 93      if delta > expected:
 94          fac = 1 - NAIVE_ADJUST_FACTOR
 95      else:
 96          fac = 1 + NAIVE_ADJUST_FACTOR
 97      return diffs[-1] * fac
 98  
 99  
100  def test(source, adjust):
101      # Variables to keep track of for stats purposes
102      ema = maxema = minema = TARGET
103      lthalf, gtdouble, lttq, gtft = 0, 0, 0, 0
104      count = 0
105      # Block times
106      times = [0]
107      # Block difficulty values
108      diffs = [source[0]]
109      # Next time to print status update
110      nextprint = 10**6
111      # Main loop
112      while times[-1] < len(source) * SECONDS_IN_DAY:
113          # Print status update every 10**6 seconds
114          if times[-1] > nextprint:
115              print '%d out of %d processed, ema %f' % \
116                  (times[-1], len(source) * SECONDS_IN_DAY, ema)
117              nextprint += 10**6
118          # Grab hashpower from data source
119          hashpower = source[int(times[-1] // SECONDS_IN_DAY)]
120          # Calculate new difficulty
121          diffs.append(adjust(times, diffs))
122          # Calculate next block time
123          times.append(times[-1] + expdiff(diffs[-1] / hashpower))
124          # Calculate min and max ema
125          ema = ema * (1 - EMA_FACTOR) + (times[-1] - times[-2]) * EMA_FACTOR
126          minema = min(minema, ema)
127          maxema = max(maxema, ema)
128          count += 1
129          # Keep track of number of blocks we are below 75/50% or above
130          # 133/200% of target
131          if ema < TARGET * 0.75:
132              lttq += 1
133              if ema < TARGET * 0.5:
134                  lthalf += 1
135          elif ema > TARGET * 1.33333:
136              gtft += 1
137              if ema > TARGET * 2:
138                  gtdouble += 1
139          # Pop items to save memory
140          if len(times) > 2000:
141              times.pop(0)
142              diffs.pop(0)
143      print 'min', minema, 'max', maxema, 'avg', times[-1] / count, \
144          'ema < half', lthalf * 1.0 / count, \
145          'ema > double', gtdouble * 1.0 / count, \
146          'ema < 3/4', lttq * 1.0 / count, \
147          'ema > 4/3', gtft * 1.0 / count
148  
149  # Example usage
150  # blkdiff.test(blkdiff.hashpower, blkdiff.simple_adjust)