import numpy as np
from scipy.optimize import minimize
import scipy
import math
np.random.seed(0)
k = 3
m = 420
p = 6
comp_ratio = 0.8901 #The competitive ratio we claim.
def verify_competitive_ratio(cs, eps):
assert m%p == 0
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)]
C = np.array(C)
Prob = [[[0 for i in range(m+5)] for j in range(k+5)] for b in range(2)]
def cost(l):
l = l[0]
for i in range(m+1):
for j in range(k+1):
for b in range(2):
if j>=k or i>=m:
Prob[b][j][i]=b
for i in range(m+1, -1, -1):
for j in range(k+1, -1, -1):
for b in range(2):
if j>=k or i>=m:
continue
if l<=C[j, i]:
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])
else:
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]
return Prob[0][0][0] - comp_ratio*(1-np.exp(-l))
defeciency = 1 #This is the minimum of Prob[0][0][0] - comp_ratio*(1-e^(-l))
#Needs to be >=0 at end of execution to that claimed ratio is True
defeciency = min(defeciency, 1-np.exp(-sum(C[0])/m) - comp_ratio) #Case of l'>max(C)
"""Run global optimization on cost in the range (0, max(C)]"""
res = scipy.optimize.shgo(cost, bounds=[(0,C.max())], iters=10,
options={'disp':False, 'f_tol':1e-9})
defeciency = min(defeciency, res.fun)
"As a sanity check, make sure that global minimizer succeeded"
ls = np.linspace(0.0, C.max(), math.ceil(C.max()/eps))
for l in ls:
assert res.fun <= cost([l])
"""End of sanity check"""
if defeciency>=0:
print(f"Claimed bound of {comp_ratio} is True")
else:
print(f"Claimed bound of {comp_ratio} is False")
def MonteCarlo(cs):
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)]
C = np.array(C)
N = 4000 #Large n, bound converges for n->Infinity
epochs = 10000
prophet = 0
alg = 0
for _ in range(epochs):
X = np.random.uniform(0, 1, (N, 2)) #First dim=time, second dim=Xi~U(0, 1)
X = X[X[:, 0].argsort()] #Sort by time of arrival
r = 0
clock = 0
i_star = None
for i in range(len(X)):
ti, vi = X[i][0], X[i][1]
if r<k and vi>=1- C[r][math.floor(ti*m)]/N and ti>=clock:
r = r + 1
i_star = i
clock = math.ceil(ti*m)/m
prophet += X[:, 1].max()
if i_star:
alg += X[i_star][1]
print(f"The competitive ratio on n={N} IID U(0, 1) random variables is {alg/prophet} using {epochs} epochs.")
cs0 = [3.64589394e+00, 3.58116098e+00, 2.03323633e+00, 1.93319241e+00,
1.15603731e+00, 9.92652855e-01, 6.10147568e-01, 3.94833386e-01,
2.41093283e-01, 1.36659577e-01, 4.80563875e-02, 2.83455285e-02,
8.39298670e-02, 1.91858842e-02, 0.00133218127, 1.33218127e-03,
1.05769060e-03, 1.05769044e-03]
verify_competitive_ratio(cs0, eps=0.0001) #Takes 3-4 minutes
MonteCarlo(cs0) #Takes ~1 min
aW1wb3J0IG51bXB5IGFzIG5wCmZyb20gc2NpcHkub3B0aW1pemUgaW1wb3J0IG1pbmltaXplCmltcG9ydCBzY2lweQppbXBvcnQgbWF0aApucC5yYW5kb20uc2VlZCgwKQoKayA9IDMKbSA9IDQyMApwID0gNgpjb21wX3JhdGlvID0gMC44OTAxICNUaGUgY29tcGV0aXRpdmUgcmF0aW8gd2UgY2xhaW0uCgpkZWYgdmVyaWZ5X2NvbXBldGl0aXZlX3JhdGlvKGNzLCBlcHMpOgogICAgYXNzZXJ0IG0lcCA9PSAwCiAgICAKICAgIEMgPSBbW2NzW291dGVyICsgaW5uZXJdIGZvciBpbm5lciBpbiByYW5nZShwLTEsIC0xLCAtMSkgZm9yIF8gaW4gcmFuZ2UobS8vcCldIGZvciBvdXRlciBpbiByYW5nZSgwLCBsZW4oY3MpLTEscCldCiAgICBDID0gbnAuYXJyYXkoQykKICAgIAogICAgUHJvYiA9IFtbWzAgZm9yIGkgaW4gcmFuZ2UobSs1KV0gZm9yIGogaW4gcmFuZ2Uoays1KV0gZm9yIGIgaW4gcmFuZ2UoMildCiAgICAKICAgIGRlZiBjb3N0KGwpOgogICAgICAgIGwgPSBsWzBdCiAgICAgICAgZm9yIGkgaW4gcmFuZ2UobSsxKToKICAgICAgICAgICAgZm9yIGogaW4gcmFuZ2UoaysxKToKICAgICAgICAgICAgICAgIGZvciBiIGluIHJhbmdlKDIpOgogICAgICAgICAgICAgICAgICAgIGlmIGo+PWsgb3IgaT49bToKICAgICAgICAgICAgICAgICAgICAgICAgUHJvYltiXVtqXVtpXT1iCiAgICAgICAgCiAgICAgICAgZm9yIGkgaW4gcmFuZ2UobSsxLCAtMSwgLTEpOgogICAgICAgICAgICBmb3IgaiBpbiByYW5nZShrKzEsIC0xLCAtMSk6CiAgICAgICAgICAgICAgICBmb3IgYiBpbiByYW5nZSgyKToKICAgICAgICAgICAgICAgICAgICBpZiBqPj1rIG9yIGk+PW06CiAgICAgICAgICAgICAgICAgICAgICAgIGNvbnRpbnVlIAogICAgICAgICAgICAgICAgICAgIGlmIGw8PUNbaiwgaV06CiAgICAgICAgICAgICAgICAgICAgICAgIFByb2JbYl1bal1baV0gPSBucC5leHAoLUNbaixpXS9tKSpQcm9iW2JdW2pdW2krMV0gKyAoMS1ucC5leHAoLUNbaixpXS9tKSkqKGwvQ1tqLCBpXSAqIFByb2JbMV1baisxXVtpKzFdICsgKENbaiwgaV0tbCkvQ1tqLCBpXSAqIFByb2JbMF1baisxXVtpKzFdKQogICAgICAgICAgICAgICAgICAgIGVsc2U6CiAgICAgICAgICAgICAgICAgICAgICAgIFByb2JbYl1bal1baV0gPSBucC5leHAoLUNbaixpXS9tKSpQcm9iW2JdW2pdW2krMV0gKyAoMS1ucC5leHAoLUNbaixpXS9tKSkqUHJvYlsxXVtqKzFdW2krMV0KICAgICAgICAKICAgICAgICByZXR1cm4gUHJvYlswXVswXVswXSAtIGNvbXBfcmF0aW8qKDEtbnAuZXhwKC1sKSkKICAgIAogICAgZGVmZWNpZW5jeSA9IDEgI1RoaXMgaXMgdGhlIG1pbmltdW0gb2YgUHJvYlswXVswXVswXSAtIGNvbXBfcmF0aW8qKDEtZV4oLWwpKQogICAgICAgICAgICAgICAgICAgICNOZWVkcyB0byBiZSA+PTAgYXQgZW5kIG9mIGV4ZWN1dGlvbiB0byB0aGF0IGNsYWltZWQgcmF0aW8gaXMgVHJ1ZQogICAgICAgIAogICAgZGVmZWNpZW5jeSA9IG1pbihkZWZlY2llbmN5LCAxLW5wLmV4cCgtc3VtKENbMF0pL20pIC0gY29tcF9yYXRpbykgI0Nhc2Ugb2YgbCc+bWF4KEMpCiAgICAKICAgICIiIlJ1biBnbG9iYWwgb3B0aW1pemF0aW9uIG9uIGNvc3QgaW4gdGhlIHJhbmdlICgwLCBtYXgoQyldIiIiCiAgICByZXMgPSBzY2lweS5vcHRpbWl6ZS5zaGdvKGNvc3QsIGJvdW5kcz1bKDAsQy5tYXgoKSldLCBpdGVycz0xMCwgIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBvcHRpb25zPXsnZGlzcCc6RmFsc2UsICdmX3RvbCc6MWUtOX0pCiAgICBkZWZlY2llbmN5ID0gbWluKGRlZmVjaWVuY3ksIHJlcy5mdW4pCiAgICAKICAgIAogICAgIkFzIGEgc2FuaXR5IGNoZWNrLCBtYWtlIHN1cmUgdGhhdCBnbG9iYWwgbWluaW1pemVyIHN1Y2NlZWRlZCIKICAgIGxzID0gbnAubGluc3BhY2UoMC4wLCBDLm1heCgpLCBtYXRoLmNlaWwoQy5tYXgoKS9lcHMpKQogICAgZm9yIGwgaW4gbHM6CiAgICAgICAgYXNzZXJ0IHJlcy5mdW4gPD0gY29zdChbbF0pCiAgICAiIiJFbmQgb2Ygc2FuaXR5IGNoZWNrIiIiCiAgICAKICAgIGlmIGRlZmVjaWVuY3k+PTA6CiAgICAgICAgcHJpbnQoZiJDbGFpbWVkIGJvdW5kIG9mIHtjb21wX3JhdGlvfSBpcyBUcnVlIikKICAgIGVsc2U6CiAgICAgICAgcHJpbnQoZiJDbGFpbWVkIGJvdW5kIG9mIHtjb21wX3JhdGlvfSBpcyBGYWxzZSIpCgoKZGVmIE1vbnRlQ2FybG8oY3MpOgogICAgQyA9IFtbY3Nbb3V0ZXIgKyBpbm5lcl0gZm9yIGlubmVyIGluIHJhbmdlKHAtMSwgLTEsIC0xKSBmb3IgXyBpbiByYW5nZShtLy9wKV0gZm9yIG91dGVyIGluIHJhbmdlKDAsIGxlbihjcyktMSxwKV0KICAgIEMgPSBucC5hcnJheShDKQoKICAgIE4gPSA0MDAwICNMYXJnZSBuLCBib3VuZCBjb252ZXJnZXMgZm9yIG4tPkluZmluaXR5CiAgICBlcG9jaHMgPSAxMDAwMCAKICAgIHByb3BoZXQgPSAwCiAgICBhbGcgPSAwCiAgICBmb3IgXyBpbiByYW5nZShlcG9jaHMpOgogICAgICAgIFggPSBucC5yYW5kb20udW5pZm9ybSgwLCAxLCAoTiwgMikpICNGaXJzdCBkaW09dGltZSwgc2Vjb25kIGRpbT1YaX5VKDAsIDEpCiAgICAgICAgWCA9IFhbWFs6LCAwXS5hcmdzb3J0KCldICNTb3J0IGJ5IHRpbWUgb2YgYXJyaXZhbCAKICAgICAgICByID0gMAogICAgICAgIGNsb2NrID0gMAogICAgICAgIGlfc3RhciA9IE5vbmUKICAgICAgICBmb3IgaSBpbiByYW5nZShsZW4oWCkpOgogICAgICAgICAgICB0aSwgdmkgPSBYW2ldWzBdLCBYW2ldWzFdCiAgICAgICAgICAgIGlmIHI8ayBhbmQgdmk+PTEtIENbcl1bbWF0aC5mbG9vcih0aSptKV0vTiAgYW5kIHRpPj1jbG9jazoKICAgICAgICAgICAgICAgIHIgPSByICsgMQogICAgICAgICAgICAgICAgaV9zdGFyID0gaQogICAgICAgICAgICAgICAgY2xvY2sgPSBtYXRoLmNlaWwodGkqbSkvbQogICAgICAgIHByb3BoZXQgKz0gWFs6LCAxXS5tYXgoKQogICAgICAgIGlmIGlfc3RhcjoKICAgICAgICAgICAgYWxnICs9IFhbaV9zdGFyXVsxXQogICAgCiAgICBwcmludChmIlRoZSBjb21wZXRpdGl2ZSByYXRpbyBvbiBuPXtOfSBJSUQgVSgwLCAxKSByYW5kb20gdmFyaWFibGVzIGlzIHthbGcvcHJvcGhldH0gdXNpbmcge2Vwb2Noc30gZXBvY2hzLiIpCiAgICAgICAgICAgIAoKY3MwID0gWzMuNjQ1ODkzOTRlKzAwLCAzLjU4MTE2MDk4ZSswMCwgMi4wMzMyMzYzM2UrMDAsIDEuOTMzMTkyNDFlKzAwLAogICAgICAgMS4xNTYwMzczMWUrMDAsIDkuOTI2NTI4NTVlLTAxLCA2LjEwMTQ3NTY4ZS0wMSwgMy45NDgzMzM4NmUtMDEsCiAgICAgICAyLjQxMDkzMjgzZS0wMSwgMS4zNjY1OTU3N2UtMDEsIDQuODA1NjM4NzVlLTAyLCAyLjgzNDU1Mjg1ZS0wMiwKICAgICAgIDguMzkyOTg2NzBlLTAyLCAxLjkxODU4ODQyZS0wMiwgMC4wMDEzMzIxODEyNywgMS4zMzIxODEyN2UtMDMsCiAgICAgICAxLjA1NzY5MDYwZS0wMywgMS4wNTc2OTA0NGUtMDNdCiAgCnZlcmlmeV9jb21wZXRpdGl2ZV9yYXRpbyhjczAsIGVwcz0wLjAwMDEpICNUYWtlcyAgMy00IG1pbnV0ZXMKTW9udGVDYXJsbyhjczApICNUYWtlcyB+MSBtaW4=