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