import numpy as np
from scipy.optimize import minimize
import mpmath as mp
from mpmath import mpf
import scipy
m = 16 #m parameter from paper
mp.dps = 500 #This will force mpmath to use a precision of
#500 decimal places, just as a sanity check
comp_ratio = mpf('0.6724') #This is the competitive ratio we claim
def stable_qtk(x):
#The function (1-e^(-x))/x is unstable for small x, so we will
# Lower bound it using the summation in Equation 13 in the paper
ans = mpf('0')
for beta in range(30):
ans += mp.exp(-x) * x**beta / mp.factorial(beta) * 1/(beta+1)
return ans
#Computes f_j(alphas, eta) in time O(m^2)
def fj(j, alphas, eta):
part1 = mpf('0')
for k in range(1, j): #Goes from 1 to j-1 as in paper
part1 += mpf('1')*1/m * (1-alphas[k])
#alphas_hat[nu]=alphas[nu] if nu<=j-1 and eta if nu==j
alphas_hat = [alphas[nu] for nu in range(j)] + [eta]
part2 = mpf('0')
for k in range(j, m+1): #Goes from j to m as in paper
product = mpf('1')
for nu in range(1, k): #Goes from 1 to k-1
product *= (alphas[nu]**(1/m))
wk = mpf('0')
s_nu = mpf('0')
for nu in range(j): #from 0 to j-1
r_nu = (m-(k-1)+nu)/m * mp.log(alphas_hat[nu]/alphas_hat[nu+1])
wk += mp.exp(-s_nu)*(1-mp.exp(-r_nu)) * 1/(m-(k-1)+nu)
s_nu += r_nu
q_t_k = stable_qtk( 1/m * mp.log(eta/alphas[k]) )
part2 += product * wk * q_t_k
return part1 + part2
def verify_competitive_ratio(alphas):
assert np.isclose(np.float64(alphas[0]), 1) #first should be 1
assert np.isclose(np.float64(alphas[-1]), 0) #Last should be 0
assert len(alphas)==(m+2)
defeciency = mpf('1') #This quantity is the mninimum of
# f_j(alpha1, ..., alpha m, eta)- comp_ratio*(1-eta)
#This needs to be >=0 at the end of the code
for j in range(1, m+2): #Goes from 1 to m+1 as in paper
eta_bounds = [(np.float64(alphas[j]),np.float64(alphas[j-1]))]
x0 = [np.float64((alphas[j]+alphas[j-1])/2)]
res = minimize(lambda alphat:
fj(j, alphas, alphat[0]) - comp_ratio*(1-alphat[0]),
x0=x0,
bounds=eta_bounds)
"""
As a sanity check, we will evaluate fj(alphas, x) - comp_ratio*(1-x) for x in
eta_bounds and assert that res.fun (the minimum value we got) is <= that.
This is just a sanity check to increase the confidence that the minimizer
actually got the right minimum
"""
opt = fj(j, alphas, res.x[0]) - comp_ratio*(1-res.x[0])
trials = np.linspace(eta_bounds[0][0], eta_bounds[0][1], 200) #200 breaks
min_in_trials = min([ fj(j, alphas, x) - comp_ratio*(1-x) for x in trials ])
assert opt <= min_in_trials
"""
End of sanity check
"""
defeciency = min(defeciency, opt)
if defeciency>=mpf('0'):
print("Claimed bound is true")
else:
print("Claimed bound is False")
alphas = [mpf('1.0'), mpf('0.66758603836404173'), mpf('0.62053145929311715'),
mpf('0.57324846512425975'),
mpf('0.52577742556626594'), mpf('0.47816906417879007'), mpf('0.43049233470891257'),
mpf('0.38283722646593055'), mpf('0.33533950489086961'), mpf('0.28831226925828957'),
mpf('0.23273108361807243'), mpf('0.19315610994691487'), mpf('0.16547915613363387'),
mpf('0.13558301500280728'), mpf('0.10412501367635961'), mpf('0.071479537771643828'),
mpf('0.036291830527618585'), mpf('0.0')]
verify_competitive_ratio(alphas) #Takes roughly 20 seconds
aW1wb3J0IG51bXB5IGFzIG5wCmZyb20gc2NpcHkub3B0aW1pemUgaW1wb3J0IG1pbmltaXplCmltcG9ydCBtcG1hdGggYXMgbXAKZnJvbSBtcG1hdGggaW1wb3J0IG1wZgppbXBvcnQgc2NpcHkKCgoKbSA9IDE2ICNtIHBhcmFtZXRlciBmcm9tIHBhcGVyCm1wLmRwcyA9IDUwMCAjVGhpcyB3aWxsIGZvcmNlIG1wbWF0aCB0byB1c2UgYSBwcmVjaXNpb24gb2YgCiAgICAgICAgICAgICAjNTAwIGRlY2ltYWwgcGxhY2VzLCBqdXN0IGFzIGEgc2FuaXR5IGNoZWNrCgpjb21wX3JhdGlvID0gbXBmKCcwLjY3MjQnKSAjVGhpcyBpcyB0aGUgY29tcGV0aXRpdmUgcmF0aW8gd2UgY2xhaW0KCmRlZiBzdGFibGVfcXRrKHgpOgogICAgI1RoZSBmdW5jdGlvbiAoMS1lXigteCkpL3ggaXMgdW5zdGFibGUgZm9yIHNtYWxsIHgsIHNvIHdlIHdpbGwgCiAgICAjIExvd2VyIGJvdW5kIGl0IHVzaW5nIHRoZSBzdW1tYXRpb24gaW4gRXF1YXRpb24gMTMgaW4gdGhlIHBhcGVyCiAgICBhbnMgPSBtcGYoJzAnKQogICAgZm9yIGJldGEgaW4gcmFuZ2UoMzApOgogICAgICAgIGFucyArPSBtcC5leHAoLXgpICogeCoqYmV0YSAvIG1wLmZhY3RvcmlhbChiZXRhKSAqIDEvKGJldGErMSkKICAgIHJldHVybiBhbnMKCiAgICAKICAgIAojQ29tcHV0ZXMgZl9qKGFscGhhcywgZXRhKSBpbiB0aW1lIE8obV4yKQpkZWYgZmooaiwgYWxwaGFzLCBldGEpOgogICAgcGFydDEgPSBtcGYoJzAnKQogICAgZm9yIGsgaW4gcmFuZ2UoMSwgaik6ICNHb2VzIGZyb20gMSB0byBqLTEgYXMgaW4gcGFwZXIKICAgICAgICBwYXJ0MSArPSBtcGYoJzEnKSoxL20gKiAoMS1hbHBoYXNba10pCgogICAgI2FscGhhc19oYXRbbnVdPWFscGhhc1tudV0gaWYgbnU8PWotMSBhbmQgZXRhIGlmIG51PT1qCiAgICBhbHBoYXNfaGF0ID0gW2FscGhhc1tudV0gZm9yIG51IGluIHJhbmdlKGopXSArIFtldGFdICAgIAogICAgcGFydDIgPSBtcGYoJzAnKQogICAgZm9yIGsgaW4gcmFuZ2UoaiwgbSsxKTogI0dvZXMgZnJvbSBqIHRvIG0gYXMgaW4gcGFwZXIKICAgICAgICBwcm9kdWN0ID0gbXBmKCcxJykKICAgICAgICBmb3IgbnUgaW4gcmFuZ2UoMSwgayk6ICNHb2VzIGZyb20gMSB0byBrLTEKICAgICAgICAgICAgcHJvZHVjdCAqPSAoYWxwaGFzW251XSoqKDEvbSkpCiAgICAgICAgCiAgICAgICAgd2sgPSBtcGYoJzAnKQogICAgICAgIHNfbnUgPSBtcGYoJzAnKQogICAgICAgIGZvciBudSBpbiByYW5nZShqKTogI2Zyb20gMCB0byBqLTEKICAgICAgICAgICAgcl9udSA9IChtLShrLTEpK251KS9tICogbXAubG9nKGFscGhhc19oYXRbbnVdL2FscGhhc19oYXRbbnUrMV0pCiAgICAgICAgICAgIHdrICs9IG1wLmV4cCgtc19udSkqKDEtbXAuZXhwKC1yX251KSkgKiAxLyhtLShrLTEpK251KQogICAgICAgICAgICBzX251ICs9IHJfbnUKICAgICAgICAgICAgCiAgICAgICAgcV90X2sgPSBzdGFibGVfcXRrKCAxL20gKiBtcC5sb2coZXRhL2FscGhhc1trXSkgKQogICAgICAgIHBhcnQyICs9IHByb2R1Y3QgKiB3ayAqIHFfdF9rCiAgICAKICAgIHJldHVybiBwYXJ0MSArIHBhcnQyCiAgICAgICAgCgpkZWYgdmVyaWZ5X2NvbXBldGl0aXZlX3JhdGlvKGFscGhhcyk6CiAgICBhc3NlcnQgbnAuaXNjbG9zZShucC5mbG9hdDY0KGFscGhhc1swXSksIDEpICNmaXJzdCBzaG91bGQgYmUgMQogICAgYXNzZXJ0IG5wLmlzY2xvc2UobnAuZmxvYXQ2NChhbHBoYXNbLTFdKSwgMCkgI0xhc3Qgc2hvdWxkIGJlIDAKICAgIGFzc2VydCBsZW4oYWxwaGFzKT09KG0rMikKICAgIAogICAgZGVmZWNpZW5jeSA9IG1wZignMScpICAjVGhpcyBxdWFudGl0eSBpcyB0aGUgbW5pbmltdW0gb2YKICAgICAgICAgICAgICAgICAgICAjIGZfaihhbHBoYTEsIC4uLiwgYWxwaGEgbSwgZXRhKS0gY29tcF9yYXRpbyooMS1ldGEpCiAgICAgICAgICAgICAgICAgICAgI1RoaXMgbmVlZHMgdG8gYmUgPj0wIGF0IHRoZSBlbmQgb2YgdGhlIGNvZGUKICAgICAgICAKICAgIAogICAgZm9yIGogaW4gcmFuZ2UoMSwgbSsyKTogI0dvZXMgZnJvbSAxIHRvIG0rMSBhcyBpbiBwYXBlcgogICAgICAgIGV0YV9ib3VuZHMgPSBbKG5wLmZsb2F0NjQoYWxwaGFzW2pdKSxucC5mbG9hdDY0KGFscGhhc1tqLTFdKSldIAogICAgICAgIAogICAgICAgIHgwID0gW25wLmZsb2F0NjQoKGFscGhhc1tqXSthbHBoYXNbai0xXSkvMildCiAgICAgICAgcmVzID0gbWluaW1pemUobGFtYmRhIGFscGhhdDogCiAgICAgICAgICAgICAgICAgICAgICAgICAgIGZqKGosIGFscGhhcywgYWxwaGF0WzBdKSAtIGNvbXBfcmF0aW8qKDEtYWxwaGF0WzBdKSwgCiAgICAgICAgICAgICAgICAgICAgICAgeDA9eDAsIAogICAgICAgICAgICAgICAgICAgICAgIGJvdW5kcz1ldGFfYm91bmRzKQogICAgICAgICIiIgogICAgICAgIEFzIGEgc2FuaXR5IGNoZWNrLCB3ZSB3aWxsIGV2YWx1YXRlIGZqKGFscGhhcywgeCkgLSBjb21wX3JhdGlvKigxLXgpIGZvciB4IGluIAogICAgICAgIGV0YV9ib3VuZHMgYW5kIGFzc2VydCB0aGF0IHJlcy5mdW4gKHRoZSBtaW5pbXVtIHZhbHVlIHdlIGdvdCkgaXMgPD0gdGhhdC4gCiAgICAgICAgVGhpcyBpcyBqdXN0IGEgc2FuaXR5IGNoZWNrIHRvIGluY3JlYXNlIHRoZSBjb25maWRlbmNlIHRoYXQgdGhlIG1pbmltaXplciAKICAgICAgICBhY3R1YWxseSBnb3QgdGhlIHJpZ2h0IG1pbmltdW0KICAgICAgICAiIiIKICAgICAgICBvcHQgPSBmaihqLCBhbHBoYXMsIHJlcy54WzBdKSAtIGNvbXBfcmF0aW8qKDEtcmVzLnhbMF0pCiAgICAgICAgdHJpYWxzID0gbnAubGluc3BhY2UoZXRhX2JvdW5kc1swXVswXSwgZXRhX2JvdW5kc1swXVsxXSwgMjAwKSAjMjAwIGJyZWFrcyAKICAgICAgICBtaW5faW5fdHJpYWxzID0gbWluKFsgZmooaiwgYWxwaGFzLCB4KSAtIGNvbXBfcmF0aW8qKDEteCkgZm9yIHggaW4gdHJpYWxzIF0pCiAgICAgICAgYXNzZXJ0IG9wdCA8PSBtaW5faW5fdHJpYWxzCiAgICAgICAgIiIiCiAgICAgICAgRW5kIG9mIHNhbml0eSBjaGVjawogICAgICAgICIiIgogICAgICAgIGRlZmVjaWVuY3kgPSBtaW4oZGVmZWNpZW5jeSwgb3B0KQogICAgICAgICAgICAKICAgIGlmIGRlZmVjaWVuY3k+PW1wZignMCcpOgogICAgICAgIHByaW50KCJDbGFpbWVkIGJvdW5kIGlzIHRydWUiKQogICAgZWxzZToKICAgICAgICBwcmludCgiQ2xhaW1lZCBib3VuZCBpcyBGYWxzZSIpCiAgICAgICAgICAgIAogICAgCmFscGhhcyA9IFttcGYoJzEuMCcpLCBtcGYoJzAuNjY3NTg2MDM4MzY0MDQxNzMnKSwgbXBmKCcwLjYyMDUzMTQ1OTI5MzExNzE1JyksIAptcGYoJzAuNTczMjQ4NDY1MTI0MjU5NzUnKSwKIG1wZignMC41MjU3Nzc0MjU1NjYyNjU5NCcpLCBtcGYoJzAuNDc4MTY5MDY0MTc4NzkwMDcnKSwgbXBmKCcwLjQzMDQ5MjMzNDcwODkxMjU3JyksCiBtcGYoJzAuMzgyODM3MjI2NDY1OTMwNTUnKSwgbXBmKCcwLjMzNTMzOTUwNDg5MDg2OTYxJyksIG1wZignMC4yODgzMTIyNjkyNTgyODk1NycpLAogbXBmKCcwLjIzMjczMTA4MzYxODA3MjQzJyksIG1wZignMC4xOTMxNTYxMDk5NDY5MTQ4NycpLCBtcGYoJzAuMTY1NDc5MTU2MTMzNjMzODcnKSwKIG1wZignMC4xMzU1ODMwMTUwMDI4MDcyOCcpLCBtcGYoJzAuMTA0MTI1MDEzNjc2MzU5NjEnKSwgbXBmKCcwLjA3MTQ3OTUzNzc3MTY0MzgyOCcpLAogbXBmKCcwLjAzNjI5MTgzMDUyNzYxODU4NScpLCBtcGYoJzAuMCcpXQoKdmVyaWZ5X2NvbXBldGl0aXZlX3JhdGlvKGFscGhhcykgI1Rha2VzIHJvdWdobHkgMjAgc2Vjb25kcw==