import numpy as np
import scipy.stats as st 
import scipy.interpolate as intrp
import math

def getCriticalValuesKS(n, alpha):
	#GETCRITICALVALUEKS compute the critical value for the two-sided KS test
	#statistic

	# References:
	#   Massey, F.J., (1951) "The Kolmogorov-Smirnov Test for Goodness of Fit",
	#         Journal of the American Statistical Association, 46(253):68-78.
	#   Miller, L.H., (1956) "Table of Percentage Points of Kolmogorov Statistics",
	#         Journal of the American Statistical Association, 51(273):111-121.
	#   Marsaglia, G., W.W. Tsang, and J. Wang (2003), "Evaluating Kolmogorov`s
	#         Distribution", Journal of Statistical Software, vol. 8, issue 18.

	# The critical value table used below is expressed in reference to a
	# 1-sided significance level.  Need to halve the significance level for
	# a basic two-sided test.

	alpha1 = alpha / 2.

	if n <= 20:  # Small sample exact values.
		# Exact K-S test critical values are solutions of an nth order polynomial.
		# Miller's approximation is excellent for sample sizes n > 20. For n <= 20,
		# Miller tabularized the exact critical values by solving the nth order
		# polynomial. These exact values for n <= 20 are hard-coded into the matrix
		# 'exact' shown below. Rows 1:20 correspond to sample sizes n = 1:20.
		a1 = np.array([0.00500,  0.01000,  0.02500,  0.05000,  0.10000])  # 1-sided significance level

		exact = np.array([[0.99500,  0.99000,  0.97500,  0.95000,  0.90000],
        				[0.92929,  0.90000,  0.84189,  0.77639,  0.68377],
        				[0.82900,  0.78456,  0.70760,  0.63604,  0.56481],
				        [0.73424,  0.68887,  0.62394,  0.56522,  0.49265],
        				[0.66853,  0.62718,  0.56328,  0.50945,  0.44698],
        				[0.61661,  0.57741,  0.51926,  0.46799,  0.41037],
        				[0.57581,  0.53844,  0.48342,  0.43607,  0.38148],
        				[0.54179,  0.50654,  0.45427,  0.40962,  0.35831],
        				[0.51332,  0.47960,  0.43001,  0.38746,  0.33910],
	    				[0.48893,  0.45662,  0.40925,  0.36866,  0.32260],
	  			        [0.46770,  0.43670,  0.39122,  0.35242,  0.30829],
						[0.44905,  0.41918,  0.37543,  0.33815,  0.29577],
    				    [0.43247,  0.40362,  0.36143,  0.32549,  0.28470],
    				    [0.41762,  0.38970,  0.34890,  0.31417,  0.27481],
    				    [0.40420,  0.37713,  0.33760,  0.30397,  0.26588],
   					    [0.39201,  0.36571,  0.32733,  0.29472,  0.25778],
    				    [0.38086,  0.35528,  0.31796,  0.28627,  0.25039],
    				    [0.37062,  0.34569,  0.30936,  0.27851,  0.24360],
    				    [0.36117,  0.33685,  0.30143,  0.27136,  0.23735],
    				    [0.35241,  0.32866,  0.29408,  0.26473,  0.23156]])

		criticalValue  =  intrp.interpolate1d(a1 , exact[n-1,:], kind = 'cubic')(alpha1)
	else: # Large sample approximate values.
		#  alpha is a 1-sided significance level
		A = 0.09037 * (-math.log(alpha1, 10))**1.5 + 0.01515 * math.log(alpha1,10)**2 - 0.08467 * alpha1 - 0.11143
		asymptoticStat =  np.sqrt(-0.5*np.log(alpha1)/n)
		criticalValue  =  asymptoticStat - 0.16693 / n - A / n**1.5
    	
		criticalValue  =  np.min([criticalValue, 1-alpha1])
	
	return criticalValue


def cdf(X, bnds = [-5, 5]):
	"""
	Empirical CDF
	"""
	n = 1000
	x = np.linspace(bnds[0], bnds[1], n)
	Nx = np.array([(X < x_i).sum() for x_i in x]) / float(X.shape[0]) # Computes empirical CDF
	xx = np.hstack([x.reshape(n,1), x.reshape(n,1)]).flatten()
	Nxx = np.hstack([ np.hstack([Nx.reshape(n,1), Nx.reshape(n,1)]).flatten()[1:2*n], 1])
	return xx, Nxx

def Burr12(x, c = 2, d = 4):
	return st.burr12(c,d).cdf(x)

def KSTest(X, alpha = 0.1):
	xgrid = np.linspace(0, 1, 10001)
	n = X.shape[0]
	Fhat = [(X<=x0).sum()/float(n) for x0 in xgrid]
	Dn = np.max(np.abs(Fhat - Burr12(xgrid)))
	#[Dn1, p] = st.kstest(X, F) # Comparison with the scipy.stats' builtin function
	print( 'KS Test statistic: ' + str(Dn) )

	val = getCriticalValuesKS(n, alpha)
	true = 1 * (Dn > val)
	rej = ['cannot be', 'is']
	stat = {'Statistic': Dn, 'Quantile': val, 'Significance': alpha}
	message = 'KS test: the null hypothesis H0 ' + rej[true] + ' rejected at level alpha = ' + str(alpha) 
	return message, stat