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 bits/double, just as a sanity check
  12.  
  13.  
  14. def stable_qtk(x):
  15.     #The function (1-e^(-x))/x is highly unstable for small x, so we will 
  16.     # Lower bound it using the summation in Equation 13 in the paper
  17.     ans = 0
  18.     for beta in range(30):
  19.         ans += mp.exp(-x) * x**beta / mp.factorial(beta) * 1/(beta+1)
  20.     return ans
  21.  
  22.  
  23.  
  24. #Computes f_j(alphas, alphat) in time O(m^2)
  25. def fj(j, alphas, alphat):
  26.     part1 = 0
  27.     for k in range(1, j): #Goes from 1 to j-1 as in paper
  28.         part1 += 1/m * (1-alphas[k])
  29.  
  30.     #alphas_hat[nu]=alphas[nu] if nu<=j-1 and alphat if nu==j
  31.     alphas_hat = [alphas[nu] for nu in range(j)] + [alphat]    
  32.     part2 = 0
  33.     for k in range(j, m+1): #Goes from j to m as in paper
  34.         product = 1
  35.         for nu in range(1, k): #Goes from 1 to k-1
  36.             product *= (alphas[nu]**(1/m))
  37.  
  38.         wk = 0
  39.         s_nu = 0
  40.         for nu in range(j): #from 0 to j-1
  41.             r_nu = (m-(k-1)+nu)/m * mp.log(alphas_hat[nu]/alphas_hat[nu+1])
  42.             wk += mp.exp(-s_nu)*(1-mp.exp(-r_nu)) * 1/(m-(k-1)+nu)
  43.             s_nu += r_nu
  44.  
  45.         q_t_k = stable_qtk( 1/m * mp.log(alphat/alphas[k]) )
  46.         part2 += product * wk * q_t_k
  47.  
  48.     return part1 + part2
  49.  
  50.  
  51. def evaluate_competitive_ratio(alphas):
  52.     assert np.isclose(np.float64(alphas[0]), 1) #first should be 1
  53.     assert np.isclose(np.float64(alphas[-1]), 0) #Last should be 0
  54.     assert len(alphas)==(m+2)
  55.  
  56.     competitive_ratio = 1
  57.  
  58.     for j in range(1, m+2): #Goes from 1 to m+1 as in paper
  59.         #Avoid precision errors  when alphat~~1, subtract 1e-8
  60.         alphat_bounds = [(np.float64(alphas[j]),np.float64(min(alphas[j-1], 1-1e-8)))] 
  61.  
  62.         x0 = [np.float64((alphas[j]+alphas[j-1])/2)]
  63.         res = minimize(lambda alphat: fj(j, alphas, alphat[0])/(1-alphat[0]), 
  64.                        x0=x0, 
  65.                        bounds=alphat_bounds)
  66.         """
  67.         As a sanity check, we will evaluate fj(alphas, x)/(1-x) for x in alphat_bounds
  68.         and assert that res.fun (the minimum value we got) is <= fj(alphas, x)/(1-x). 
  69.         This is just a sanity check to increase the confidence that the minimizer 
  70.         actually got the right minimum
  71.         """
  72.         trials = np.linspace(alphat_bounds[0][0], alphat_bounds[0][1], 30) #30 breaks 
  73.         min_in_trials = min([ fj(j, alphas, x)/(1-x) for x in trials ])
  74.         assert res.fun <= min_in_trials
  75.         """
  76.         End of sanity check
  77.         """
  78.         competitive_ratio = min(competitive_ratio, res.fun)
  79.  
  80.     return competitive_ratio
  81.  
  82.  
  83. alphas = [mpf('1.0'), mpf('0.66758603836404173'), mpf('0.62053145929311715'), 
  84. mpf('0.57324846512425975'),
  85.  mpf('0.52577742556626594'), mpf('0.47816906417879007'), mpf('0.43049233470891257'),
  86.  mpf('0.38283722646593055'), mpf('0.33533950489086961'), mpf('0.28831226925828957'),
  87.  mpf('0.23273108361807243'), mpf('0.19315610994691487'), mpf('0.16547915613363387'),
  88.  mpf('0.13558301500280728'), mpf('0.10412501367635961'), mpf('0.071479537771643828'),
  89.  mpf('0.036291830527618585'), mpf('0.0')]
  90. c = evaluate_competitive_ratio(alphas)
  91. print(c)
  92.  
Runtime error #stdin #stdout #stderr 0.36s 63848KB
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/ibg3eH
language:
Python 3 (python  3.12)
created:
1 year 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!