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