Empirical bounds checking on approximations related to the performance of a parallel queue structure

queue_math.py

#!/usr/bin/python
#
# queue_math.py
#
# Checking some calculations on the memory efficiency of a
# hypothetical queue data structure.
#

import random
import sys

# Parameters

num_trials = 1000

# Functions

def lower(n, k):
  prod = 1.0
  for i in range(k):
    prod *= (n - i)
  return prod

def exp_m_for_k_is_1(n):
  sum = 0
  for i in range(1, n + 1):
    sum += lower(n, i) / float(n ** i)
  return sum

def monte_carlo_m_for_k_is_1(n):
  global num_trials
  m_sum = 0
  for trial in xrange(num_trials):
    array = [0] * n
    for i in xrange(n + 1):
      j = random.randrange(n)
      if array[j]:
        m_sum += i
        break
      array[j] = 1
  return m_sum / float(num_trials)

def p_sub_ell(n, k, ell):
  prod = 1.0
  for i in range(ell):
    prod *= (1 - (i / float(n)) ** k)
  return prod

def exp_m(n, k):
  sum = 0
  for ell in range(1, n + 1):
    sum += p_sub_ell(n, k, ell)
  return sum

def monte_carlo_m(n, k):
  global num_trials
  m_sum = 0
  for trial in xrange(num_trials):
    array = [0] * n
    for i in xrange(n + 1):
      for j in xrange(k):
        h = random.randrange(n)
        if array[h]: continue
        array[h] = 1
        break
      else:  # Only runs if break never happens.
        m_sum += i
        break
  return m_sum / float(num_trials)

def lower_bound_m(n, k):
  a = n ** (k / float(k + 1))
  b = 1 - 1 / float((k + 2) * (k + 1))
  return a * b

def lower_bound_m2(n, k):
  if k + 1 > n: return 0, 0  # This bound is only useful for k + 1 <= n.
  lam = (k + 1) ** (1 / float(k + 1))
  alpha = k / float(k + 1)
  t = int(lam * (n ** alpha))
  val = t - (t ** (k + 2)) / float((k + 2) * (k + 1) * (n ** k))
  return val, t

def print_table(max_n, k):
  print('')
  print('k=%d' % k)
  print('n\tappr1\tappr2\tsum\tmonte carlo')
  print('----------------------------------------')
  for n in range(1, max_n + 1):
    a1 = lower_bound_m(n, k)
    a2, t = lower_bound_m2(n, k)
    a2 = '%.3f' % a2 if a2 else '-'
    s = exp_m(n, k)
    m = monte_carlo_m(n, k)
    print('%d\t%.3f\t%s\t%.3f\t%.3f' % (n, a1, a2, s, m))
  print('')

# Main

if __name__ == '__main__':

  k = 2
  if len(sys.argv) > 1:
    k = int(sys.argv[1])

  print_table(30, k)