import numpy as np
from scipy.optimize import minimize
import scipy
 
 
 
m = 10 #m parameter from paper
 
def lamb(j, cs):
    return 1/m * sum(cs[i] for i in range(1, j))
 
 
#Computes f_j(alphas, alphat) in time O(m^2)
def fj(j, cs, l):
    part1 = 1-np.exp(-lamb(j, cs))
    part2 = 0
    for k in range(j, m+1):
        part2 += np.exp(-lamb(k, cs)) * (1-np.exp(-cs[k]/m)) * l/cs[k]
    return part1+part2
 
 
def evaluate_competitive_ratio(cs):
    for i in range(1, len(cs)):
        if cs[i]<cs[i-1]:
            raise Exception("Values are not increasing")
 
    competitive_ratio = 1-np.exp(-1/m * float(sum(cs)) )
    competitive_ratio = min(sum([np.exp(-lamb(k, cs)) * (1-np.exp(-cs[k]/m))/cs[k]  for k in range(1, m+1)]), competitive_ratio)
 
    for j in range(2, m+1): 
        alphat_bounds = [(cs[j-1],cs[j])]
        x0 = (cs[j-1]+cs[j])/2.0
 
        res = minimize(lambda l: fj(j, cs, l[0])/(1-np.exp(-l[0])), 
                       x0=x0, 
                       bounds=alphat_bounds)
        """As a sanity check, make sure res.fun <= a few values in the middle to make sure minimization worked"""
        for xx in np.linspace(alphat_bounds[0][0], alphat_bounds[0][1], 1000):
            assert res.fun <= fj(j, cs, xx)/(1-np.exp(-xx)), (alphat_bounds, xx, res)
 
        competitive_ratio = min(competitive_ratio, res.fun)
    return competitive_ratio
 
cs = [0.        , 0.07077646, 0.2268947 , 0.42146915, 0.60679691,
       0.8570195 , 1.17239753, 1.51036256, 1.9258193 , 2.88381902,
       3.97363258]
c = evaluate_competitive_ratio(cs)
print(c)
 

aW1wb3J0IG51bXB5IGFzIG5wCmZyb20gc2NpcHkub3B0aW1pemUgaW1wb3J0IG1pbmltaXplCmltcG9ydCBzY2lweQoKCgptID0gMTAgI20gcGFyYW1ldGVyIGZyb20gcGFwZXIKCmRlZiBsYW1iKGosIGNzKToKICAgIHJldHVybiAxL20gKiBzdW0oY3NbaV0gZm9yIGkgaW4gcmFuZ2UoMSwgaikpCgogICAgCiNDb21wdXRlcyBmX2ooYWxwaGFzLCBhbHBoYXQpIGluIHRpbWUgTyhtXjIpCmRlZiBmaihqLCBjcywgbCk6CiAgICBwYXJ0MSA9IDEtbnAuZXhwKC1sYW1iKGosIGNzKSkKICAgIHBhcnQyID0gMAogICAgZm9yIGsgaW4gcmFuZ2UoaiwgbSsxKToKICAgICAgICBwYXJ0MiArPSBucC5leHAoLWxhbWIoaywgY3MpKSAqICgxLW5wLmV4cCgtY3Nba10vbSkpICogbC9jc1trXQogICAgcmV0dXJuIHBhcnQxK3BhcnQyCgoKZGVmIGV2YWx1YXRlX2NvbXBldGl0aXZlX3JhdGlvKGNzKToKICAgIGZvciBpIGluIHJhbmdlKDEsIGxlbihjcykpOgogICAgICAgIGlmIGNzW2ldPGNzW2ktMV06CiAgICAgICAgICAgIHJhaXNlIEV4Y2VwdGlvbigiVmFsdWVzIGFyZSBub3QgaW5jcmVhc2luZyIpCiAgICAKICAgIGNvbXBldGl0aXZlX3JhdGlvID0gMS1ucC5leHAoLTEvbSAqIGZsb2F0KHN1bShjcykpICkKICAgIGNvbXBldGl0aXZlX3JhdGlvID0gbWluKHN1bShbbnAuZXhwKC1sYW1iKGssIGNzKSkgKiAoMS1ucC5leHAoLWNzW2tdL20pKS9jc1trXSAgZm9yIGsgaW4gcmFuZ2UoMSwgbSsxKV0pLCBjb21wZXRpdGl2ZV9yYXRpbykKCiAgICBmb3IgaiBpbiByYW5nZSgyLCBtKzEpOiAKICAgICAgICBhbHBoYXRfYm91bmRzID0gWyhjc1tqLTFdLGNzW2pdKV0KICAgICAgICB4MCA9IChjc1tqLTFdK2NzW2pdKS8yLjAKICAgICAgICAKICAgICAgICByZXMgPSBtaW5pbWl6ZShsYW1iZGEgbDogZmooaiwgY3MsIGxbMF0pLygxLW5wLmV4cCgtbFswXSkpLCAKICAgICAgICAgICAgICAgICAgICAgICB4MD14MCwgCiAgICAgICAgICAgICAgICAgICAgICAgYm91bmRzPWFscGhhdF9ib3VuZHMpCiAgICAgICAgIiIiQXMgYSBzYW5pdHkgY2hlY2ssIG1ha2Ugc3VyZSByZXMuZnVuIDw9IGEgZmV3IHZhbHVlcyBpbiB0aGUgbWlkZGxlIHRvIG1ha2Ugc3VyZSBtaW5pbWl6YXRpb24gd29ya2VkIiIiCiAgICAgICAgZm9yIHh4IGluIG5wLmxpbnNwYWNlKGFscGhhdF9ib3VuZHNbMF1bMF0sIGFscGhhdF9ib3VuZHNbMF1bMV0sIDEwMDApOgogICAgICAgICAgICBhc3NlcnQgcmVzLmZ1biA8PSBmaihqLCBjcywgeHgpLygxLW5wLmV4cCgteHgpKSwgKGFscGhhdF9ib3VuZHMsIHh4LCByZXMpCiAgICAgICAgCiAgICAgICAgY29tcGV0aXRpdmVfcmF0aW8gPSBtaW4oY29tcGV0aXRpdmVfcmF0aW8sIHJlcy5mdW4pCiAgICByZXR1cm4gY29tcGV0aXRpdmVfcmF0aW8KICAgIApjcyA9IFswLiAgICAgICAgLCAwLjA3MDc3NjQ2LCAwLjIyNjg5NDcgLCAwLjQyMTQ2OTE1LCAwLjYwNjc5NjkxLAogICAgICAgMC44NTcwMTk1ICwgMS4xNzIzOTc1MywgMS41MTAzNjI1NiwgMS45MjU4MTkzICwgMi44ODM4MTkwMiwKICAgICAgIDMuOTczNjMyNThdCmMgPSBldmFsdWF0ZV9jb21wZXRpdGl2ZV9yYXRpbyhjcykKcHJpbnQoYykK