import math
import numpy

def read_file(filepos,ratio=1):
    out=[]
    for line in file(filepos):
        data=line.split()
        out.append(ratio*float(data[1]))
    return out

def read_correlation(filepos):

    text=file(filepos).read()
    text=text[:-1].replace('\n',';')
    out=numpy.matrix(text)
    return out


class Chi_carre:
    
    def __init__(self,expected,variance,correlation,ratio=1):
        self.expected=read_file(expected,ratio)
        self.variance=read_file(variance,math.sqrt(ratio))
        self.correlation=read_correlation(correlation)
        self.len=len(self.expected)
        
    def compute_from_value(self,dico,mode=1):
        total=0
        if 1==mode: # standard chi-carre
            for key in range(0,self.len):
                total+=((dico[key]-self.expected[key])/self.variance[key])**2        
            return total/len(self.expected)
        if 2==mode: #test ofcoherence
            for key in range(0,self.len):
                prov=0
                fact1=(dico[key]-self.expected[key])/self.variance[key]
                for key2 in range(0,self.len):
                    prov+=fact1*(dico[key2]-self.expected[key2])/(self.variance[key2]*self.correlation[key,key2])
                total+=prov
            total=math.fabs(total)/(len(self.expected)**2)
            return total
        if 3==mode:
            inv_mat=self.correlation.I
            total=0
            for key in range(0,self.len):
                prov_num=0
                prov_den=0
                for key2 in range(0,self.len):
                    prov_num+=inv_mat[key,key2]*(dico[key]-self.expected[key])
                    prov_den=inv_mat[key,key2]*self.expected[key2]
                total+=(prov_num/prov_den)**2
            return total/float(self.len**2)
        if 4==mode:
            inv_mat=self.correlation.I
            total=0
            for key in range(0,self.len):
                prov_num=0
                prov_den=self.variance[key]
                for key2 in range(0,self.len):
                    prov_num+=inv_mat[key,key2]*(dico[key]-self.expected[key])
                total+=(prov_num/prov_den)**2
            return total/float(self.len**2)
        if 5==mode:
            cov=numpy.matrix(self.correlation)
            for key in range(0,self.len):
                for key2 in range(0,self.len):
                    cov[key,key2]=cov[key,key2]*self.variance[key]*self.variance[key2]
            inv_mat=cov.I
            total=0
            for key in range(0,self.len):
                prov=0
                for key2 in range(0,self.len):
                    prov+=(dico[key]-self.expected[key])*inv_mat[key,key2]*(dico[key2]-self.expected[key2])
                    total+=(dico[key]-self.expected[key])*inv_mat[key,key2]*(dico[key2]-self.expected[key2])
                print prov
            return total/float(self.len)

    def compute_from_file(self,filepos,mode=1):
        dico=read_file(filepos)
        return self.compute_from_value(dico,mode=mode)

class p_value:

    def __init__(self):
        self.p_value=[value/10.0 for value in range(0,26)]
        self.ntot=0
        self.value={}
        for value in self.p_value:
            self.value[value]=0

    def read_file(self,filename):

        for line in file(filename):
            self.add_value(float(line))

    def add_value(self,value):
            self.ntot+=1
            for pvalue in self.p_value:
                if value>pvalue:
                    self.value[pvalue]+=1
                else: 
                    break

    def compute_result(self):
        self.result()
        for value in self.p_value:
            self.value[value]/=float(self.ntot)
            self.value[value]*=100

    def result(self):
        text=''
        for pvalue in self.p_value:
            text+= str(pvalue)+' '+str(self.value[pvalue])+'\n'
        return text


def compute_pvalue():
    stat_obj=[]
    range_test=1
    for i in range(0,range_test): stat_obj.append(p_value())
#    obj1=Chi_carre('value','error','correlation')
    obj2=Chi_carre('dist_value_200','dist_error_200','dist_correlation_200')

    ff=open('chi_carre','w')
    for i in range(0,10000):
        for j in range(0,range_test):
#            chi_carre=obj1.compute_from_file('10bin_'+str(i),mode=5)
#            stat_obj[0].add_value(chi_carre)
            chi_carre=obj2.compute_from_file('dist_'+str(i),mode=5)
            stat_obj[0].add_value(chi_carre)
            ff.writelines(str(chi_carre)+'\n')
 
    for j in range(0,range_test):
        stat_obj[j].compute_result()
        text=stat_obj[j].result()
        print text
    ff.close()
    ff=open('pvalue','w')
    ff.writelines(text)
    ff.close()

if '__main__'==__name__:
    obj2=Chi_carre('MW_mean','MW_error','MW_correlation',1000)
    chi_carre=obj2.compute_from_file('MW_SM_event',mode=5)
    print chi_carre