#! /usr/bin/python
import sys
import numpy as np
from scipy.stats import lognorm
import math

def normtest(x):
    x = np.matrix(x)
    s = np.matrix(np.cov(x,bias=1,rowvar=0))
    n,p = x.shape
    dift = x - np.mean(x,axis=0)
    dj = np.matrix(np.diag(dift*s.I*dift.T))
    y = x*s.I*x.T
    djk = - 2*y.T + np.diag(y.T).T*np.ones((1,n)) + np.ones((n,1))*np.diag(y.T)
    b = 1/(math.sqrt(2))*((2*p + 1)/4)**(1/(p + 4))*(n**(1/(p + 4)))
    
    if (np.rank(s) == p):
        hz = n * (1/(n**2) * np.sum(np.sum(np.exp( - (b**2)/2 * djk))) - 2 *
        ((1 + (b**2))**( - p/2)) * (1/n) * (np.sum(np.exp( - ((b**2)/(2 *
        (1 + (b**2)))) * dj))) + ((1 + (2 * (b**2)))**( - p/2)))
    else:
        hz = n*4

    wb = (1 + b**2)*(1 + 3*b**2)
    a = 1 + 2*b**2
    mu = 1 - a**(- p/2)*(1 + p*b**2/a + (p*(p + 2)*(b**4))/(2*a**2))
    si2 = 2*(1 + 4*b**2)**(- p/2) + 2*a**( - p)*(1 + (2*p*b**4)/a**2 + (3*p*
    (p + 2)*b**8)/(4*a**4)) - 4*wb**( - p/2)*(1 + (3*p*b**4)/(2*wb) + (p*
    (p + 2)*b**8)/(2*wb**2))

    pmu = math.log(math.sqrt(mu**4/(si2 + mu**2)))
    psi = math.sqrt(math.log((si2 + mu**2)/mu**2))
    p = 1 - lognorm.cdf(hz,psi,scale=math.exp(pmu))
    return p

if __name__ == "__main__":
    filename = sys.argv[1]
    x = np.loadtxt(filename,delimiter=",")
    x = x[:6000,7:]
    p = normtest(x)
    print p