fork download
copy 
  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
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
Traceback (most recent call last):
  File "./prog.py", line 3, in <module>
ModuleNotFoundError: No module named 'mpmath'
https://ideone.com/AuSVbj
language:
Python 3 (python  3.9.5)
created:
170 days ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!