fork download 
  1. import numpy as np
  2. from scipy.optimize import minimize
  3. import scipy
  4. import math
  5. np.random.seed(0)
  6.  
  7. k = 3
  8. m = 420
  9. p = 6
  10. comp_ratio = 0.8901 #The competitive ratio we claim.
  11.  
  12. def verify_competitive_ratio(cs, eps):
  13.     assert m%p == 0
  14.  
  15.     C = [[cs[outer + inner] for inner in range(p-1, -1, -1) for _ in range(m//p)] for outer in range(0, len(cs)-1,p)]
  16.     C = np.array(C)
  17.  
  18.     Prob = [[[0 for i in range(m+5)] for j in range(k+5)] for b in range(2)]
  19.  
  20.     def cost(l):
  21.         l = l[0]
  22.         for i in range(m+1):
  23.             for j in range(k+1):
  24.                 for b in range(2):
  25.                     if j>=k or i>=m:
  26.                         Prob[b][j][i]=b
  27.  
  28.         for i in range(m+1, -1, -1):
  29.             for j in range(k+1, -1, -1):
  30.                 for b in range(2):
  31.                     if j>=k or i>=m:
  32.                         continue 
  33.                     if l<=C[j, i]:
  34.                         Prob[b][j][i] = np.exp(-C[j,i]/m)*Prob[b][j][i+1] + (1-np.exp(-C[j,i]/m))*(l/C[j, i] * Prob[1][j+1][i+1] + (C[j, i]-l)/C[j, i] * Prob[0][j+1][i+1])
  35.                     else:
  36.                         Prob[b][j][i] = np.exp(-C[j,i]/m)*Prob[b][j][i+1] + (1-np.exp(-C[j,i]/m))*Prob[1][j+1][i+1]
  37.  
  38.         return Prob[0][0][0] - comp_ratio*(1-np.exp(-l))
  39.  
  40.     defeciency = 1 #This is the minimum of Prob[0][0][0] - comp_ratio*(1-e^(-l))
  41.                     #Needs to be >=0 at end of execution to that claimed ratio is True
  42.  
  43.     defeciency = min(defeciency, 1-np.exp(-sum(C[0])/m) - comp_ratio) #Case of l'>max(C)
  44.  
  45.     """Run global optimization on cost in the range (0, max(C)]"""
  46.     res = scipy.optimize.shgo(cost, bounds=[(0,C.max())], iters=10,  
  47.                               options={'disp':False, 'f_tol':1e-9})
  48.     defeciency = min(defeciency, res.fun)
  49.  
  50.  
  51.     "As a sanity check, make sure that global minimizer succeeded"
  52.     ls = np.linspace(0.0, C.max(), math.ceil(C.max()/eps))
  53.     for l in ls:
  54.         assert res.fun <= cost([l])
  55.     """End of sanity check"""
  56.  
  57.     if defeciency>=0:
  58.         print(f"Claimed bound of {comp_ratio} is True")
  59.     else:
  60.         print(f"Claimed bound of {comp_ratio} is False")
  61.  
  62.  
  63. def MonteCarlo(cs):
  64.     C = [[cs[outer + inner] for inner in range(p-1, -1, -1) for _ in range(m//p)] for outer in range(0, len(cs)-1,p)]
  65.     C = np.array(C)
  66.  
  67.     N = 4000 #Large n, bound converges for n->Infinity
  68.     epochs = 10000 
  69.     prophet = 0
  70.     alg = 0
  71.     for _ in range(epochs):
  72.         X = np.random.uniform(0, 1, (N, 2)) #First dim=time, second dim=Xi~U(0, 1)
  73.         X = X[X[:, 0].argsort()] #Sort by time of arrival 
  74.         r = 0
  75.         clock = 0
  76.         i_star = None
  77.         for i in range(len(X)):
  78.             ti, vi = X[i][0], X[i][1]
  79.             if r<k and vi>=1- C[r][math.floor(ti*m)]/N  and ti>=clock:
  80.                 r = r + 1
  81.                 i_star = i
  82.                 clock = math.ceil(ti*m)/m
  83.         prophet += X[:, 1].max()
  84.         if i_star:
  85.             alg += X[i_star][1]
  86.  
  87.     print(f"The competitive ratio on n={N} IID U(0, 1) random variables is {alg/prophet} using {epochs} epochs.")
  88.  
  89.  
  90. cs0 = [3.64589394e+00, 3.58116098e+00, 2.03323633e+00, 1.93319241e+00,
  91.        1.15603731e+00, 9.92652855e-01, 6.10147568e-01, 3.94833386e-01,
  92.        2.41093283e-01, 1.36659577e-01, 4.80563875e-02, 2.83455285e-02,
  93.        8.39298670e-02, 1.91858842e-02, 0.00133218127, 1.33218127e-03,
  94.        1.05769060e-03, 1.05769044e-03]
  95.  
  96. verify_competitive_ratio(cs0, eps=0.0001) #Takes  3-4 minutes
  97. MonteCarlo(cs0) #Takes ~1 min
Time limit exceeded #stdin #stdout 5s 63080KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/6G7hou
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!