fork download
  1. import numpy as np
  2. from scipy.optimize import minimize
  3. import mpmath as mp
  4. from mpmath import mpf
  5. import scipy
  6.  
  7.  
  8.  
  9. m = 16 #m parameter from paper
  10. mp.dps = 500 #This will force mpmath to use a precision of
  11. #500 decimal places, just as a sanity check
  12.  
  13. comp_ratio = mpf('0.6724') #This is the competitive ratio we claim
  14.  
  15. def stable_qtk(x):
  16. #The function (1-e^(-x))/x is unstable for small x, so we will
  17. # Lower bound it using the summation in Equation 13 in the paper
  18. ans = mpf('0')
  19. for beta in range(30):
  20. ans += mp.exp(-x) * x**beta / mp.factorial(beta) * 1/(beta+1)
  21. return ans
  22.  
  23.  
  24.  
  25. #Computes f_j(alphas, eta) in time O(m^2)
  26. def fj(j, alphas, eta):
  27. part1 = mpf('0')
  28. for k in range(1, j): #Goes from 1 to j-1 as in paper
  29. part1 += mpf('1')*1/m * (1-alphas[k])
  30.  
  31. #alphas_hat[nu]=alphas[nu] if nu<=j-1 and eta if nu==j
  32. alphas_hat = [alphas[nu] for nu in range(j)] + [eta]
  33. part2 = mpf('0')
  34. for k in range(j, m+1): #Goes from j to m as in paper
  35. product = mpf('1')
  36. for nu in range(1, k): #Goes from 1 to k-1
  37. product *= (alphas[nu]**(1/m))
  38.  
  39. wk = mpf('0')
  40. s_nu = mpf('0')
  41. for nu in range(j): #from 0 to j-1
  42. r_nu = (m-(k-1)+nu)/m * mp.log(alphas_hat[nu]/alphas_hat[nu+1])
  43. wk += mp.exp(-s_nu)*(1-mp.exp(-r_nu)) * 1/(m-(k-1)+nu)
  44. s_nu += r_nu
  45.  
  46. q_t_k = stable_qtk( 1/m * mp.log(eta/alphas[k]) )
  47. part2 += product * wk * q_t_k
  48.  
  49. return part1 + part2
  50.  
  51.  
  52. def verify_competitive_ratio(alphas):
  53. assert np.isclose(np.float64(alphas[0]), 1) #first should be 1
  54. assert np.isclose(np.float64(alphas[-1]), 0) #Last should be 0
  55. assert len(alphas)==(m+2)
  56.  
  57. defeciency = mpf('1') #This quantity is the mninimum of
  58. # f_j(alpha1, ..., alpha m, eta)- comp_ratio*(1-eta)
  59. #This needs to be >=0 at the end of the code
  60.  
  61.  
  62. for j in range(1, m+2): #Goes from 1 to m+1 as in paper
  63. eta_bounds = [(np.float64(alphas[j]),np.float64(alphas[j-1]))]
  64.  
  65. x0 = [np.float64((alphas[j]+alphas[j-1])/2)]
  66. res = minimize(lambda alphat:
  67. fj(j, alphas, alphat[0]) - comp_ratio*(1-alphat[0]),
  68. x0=x0,
  69. bounds=eta_bounds)
  70. """
  71. As a sanity check, we will evaluate fj(alphas, x) - comp_ratio*(1-x) for x in
  72. eta_bounds and assert that res.fun (the minimum value we got) is <= that.
  73. This is just a sanity check to increase the confidence that the minimizer
  74. actually got the right minimum
  75. """
  76. opt = fj(j, alphas, res.x[0]) - comp_ratio*(1-res.x[0])
  77. trials = np.linspace(eta_bounds[0][0], eta_bounds[0][1], 200) #200 breaks
  78. min_in_trials = min([ fj(j, alphas, x) - comp_ratio*(1-x) for x in trials ])
  79. assert opt <= min_in_trials
  80. """
  81. End of sanity check
  82. """
  83. defeciency = min(defeciency, opt)
  84.  
  85. if defeciency>=mpf('0'):
  86. print("Claimed bound is true")
  87. else:
  88. print("Claimed bound is False")
  89.  
  90.  
  91. alphas = [mpf('1.0'), mpf('0.66758603836404173'), mpf('0.62053145929311715'),
  92. mpf('0.57324846512425975'),
  93. mpf('0.52577742556626594'), mpf('0.47816906417879007'), mpf('0.43049233470891257'),
  94. mpf('0.38283722646593055'), mpf('0.33533950489086961'), mpf('0.28831226925828957'),
  95. mpf('0.23273108361807243'), mpf('0.19315610994691487'), mpf('0.16547915613363387'),
  96. mpf('0.13558301500280728'), mpf('0.10412501367635961'), mpf('0.071479537771643828'),
  97. mpf('0.036291830527618585'), mpf('0.0')]
  98.  
  99. verify_competitive_ratio(alphas) #Takes roughly 20 seconds
Runtime error #stdin #stdout #stderr 0.41s 63204KB
stdin
Standard input is empty
stdout
Standard output is empty
stderr
Traceback (most recent call last):
  File "./prog.py", line 3, in <module>
ModuleNotFoundError: No module named 'mpmath'