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
stdin
Standard input is empty
stdout
Standard output is empty