DifferentialWeight: correlation.py.txt

File correlation.py.txt, 10.5 KB (added by trac, 7 years ago)
Line 
1import numpy,math
2
3#1 ###########################################################################################################################
4class histo:                                                                         #########################################
5#1 ###########################################################################################################################
6
7    #2 ####################################################################################       
8    def __init__(self):
9        """ initialisation of variable """
10        self.nb_point=0
11        self.value={}
12        self.bin=[]
13        self.len=0
14       
15    #2 ####################################################################################       
16    def init_dict(self,data):
17        """ define 0 entry for each bin  """
18
19        for bin,value in data.items():
20            self.bin.append('%f' % float(bin))
21            self.value['%f' % float(bin)]=0
22        self.len=len(data)
23
24    #2 ####################################################################################                             
25    def add(self,value):
26   
27        if self.nb_point==0:
28                self.init_dict(value)
29        self.nb_point+=1               
30        for bin in self.bin:
31            self.value[bin]+=float(value[bin])
32
33    #2 ####################################################################################                             
34    def write(self,filename):
35   
36            ff=open(filename,'w')
37            for i in range(0,self.len):
38                    ff.writelines(str(i)+'\t'+str(self.value[i])+'\n')
39                   
40    #2 ####################################################################################                             
41    def out(self):
42        return self.value
43       
44    #2 ####################################################################################       
45    def charge_file(self,name,auto=1):
46        """ read file and return a dict of type {bin:value} """
47
48        out={}
49        for line in file(name):
50            data=line.split()
51            out['%f' % float(data[0])] = float(data[1])
52       
53        if auto==1:
54            self.add(out)
55           
56        return out
57
58
59#1 ###########################################################################################################################
60class Error_in_histo(histo):                                                         #########################################
61#1 ###########################################################################################################################
62
63    #2 ####################################################################################
64    def __init__(self):
65        """ initialisation of variable """
66        self.nb_point=0
67        self.sum={}
68        self.sum_square={}
69        self.quad_term={}
70       
71    #2 ####################################################################################       
72    def init_dict(self,data):
73        """ define 0 entry for each bin  """
74
75        if isinstance(data,dict):
76            data=data.keys()
77
78        for bin in data:
79            self.sum['%f' % float(bin)]=0
80            self.sum_square['%f' % float(bin)]=0
81            self.quad_term['%f' % float(bin)]={}
82            list_bin2=[ bin2 for bin2 in data if float(bin2)<float(bin)]
83            for bin2 in list_bin2:
84                self.quad_term['%f'% float(bin)]['%f' % float(bin2)]=0
85       
86    #2 ####################################################################################       
87    def charge_file(self,name,auto=1):
88        """ read file and return a dict of type {bin:value} """
89       
90        data=histo.charge_file(self,name,auto=0)
91       
92        if auto==1:
93            self.update_value(data)
94               
95        return data
96       
97    #2 ####################################################################################       
98    def update_value(self,data):
99        """ update self.sum,self.sum_square,self.quad_term,...
100            data is a dict of type {bin:value}
101            doesn't return anything
102        """
103        if self.nb_point==0 and self.sum=={}:
104            self.init_dict(data)
105
106        self.nb_point+=1
107        prov= data.keys()
108        prov.sort()
109
110        for bin,value in data.items():
111            self.sum['%f' % float(bin)] += float(value)
112            self.sum_square['%f' % float(bin)] += float(value)**2
113            list_bin2=[ bin2 for bin2 in data.keys() if float(bin2)<float(bin)]
114            for bin2 in list_bin2:
115                self.quad_term['%f' % float(bin)]['%f' % float(bin2)]+=float(value)*float(data[bin2])
116
117    #2 ####################################################################################       
118    def compute_result(self):
119        """ compute the actual status of the average,variance and correlation
120            return self.average(dict),self.variance(dict),self.correlation(dict^2)
121        """
122        average={}
123        variance={}
124        covariance={}
125        correlation={}
126
127        key_list=[float(value) for value in self.sum.keys()]
128        key_list.sort()
129        key_list=['%f' % value for value in key_list]
130        for bin in key_list:
131            #average
132            average[bin]= self.sum[bin]/self.nb_point
133            #variance
134            try:
135                variance[bin]=math.sqrt((self.sum_square[bin]-self.nb_point*average[bin]**2)/(self.nb_point))#add -1 to have the other estimator
136            except ValueError:
137                print (self.sum_square[bin]-self.nb_point*average[bin]),'/',(self.nb_point-1)
138                variance[bin]=0.0
139            #covariance
140            covariance[bin]={}
141            correlation[bin]={}
142            list_bin2=[bin2 for bin2 in self.sum.keys() if float(bin2)<float(bin)]
143            for bin2 in list_bin2:
144                covariance[bin][bin2]=self.quad_term[bin][bin2]/self.nb_point-average[bin]*average[bin2]
145               
146                if variance[bin] and variance[bin2]: correlation[bin][bin2]=covariance[bin][bin2]/(variance[bin]*variance[bin2])
147                else:                                correlation[bin][bin2]=float('nan')
148
149        self.average=average
150        self.variance=variance
151        #self.covariance=covariance
152        self.correlation=correlation
153        return average,variance,correlation
154       
155    #2 ####################################################################################       
156    def write_result(self,filename1,filename2,filename3,n=1):
157        """ write result in a file """
158
159        key_list=[ float(value) for value in self.sum.keys()]
160        key_list.sort()
161        key_list=['%f' % value for value in key_list]
162
163        ff=open(filename1,'w')
164        for bin in key_list:
165            ff.writelines(str(bin)+'\t'+str(n*self.average[bin])+'\n')
166        ff.close()
167
168        ff=open(filename2,'w')
169        for bin in key_list:
170            ff.writelines(str(bin)+'\t'+str(math.sqrt(n)*self.variance[bin])+'\n')
171        ff.close()
172
173
174        ff=open(filename3,'w')
175        #ff.writelines('\t'.join([' ']+[str(name) for name in key_list])+'\n')
176        for i in range(0,len(key_list)):
177                line=''
178                for j in range(0,len(key_list)):
179                        if j<i:    line+=str(self.correlation[key_list[i]][key_list[j]])+'\t'
180                        elif j==i: line+='1.'+'0'*12+'\t'
181                        else:      line+=str(self.correlation[key_list[j]][key_list[i]])+'\t'
182                ff.writelines(line+'\n')
183        ff.close()
184
185
186#1 ###########################################################################################################################
187def MW_correlation(nb_event):                                                        #########################################
188#1 ###########################################################################################################################
189    """ compute sigma and correlation from a sample of event """
190
191    stat_tool=Error_in_histo()
192    j=0
193    while j<nb_event:
194        if j%250==0: print 'status',j
195        step=numpy.random.poisson(1)
196        hist=histo()
197        for k in range(j,j+step):
198            hist.charge_file('card_1event_'+str(j)+'.txt')
199        stat_tool.update_value(hist.out())
200        j+=step
201    stat_tool.compute_result()
202    stat_tool.write_result('../SM_PGS_mean_5','../SM_PGS_error_5','../SM_PGS_correlation_5')
203
204
205#1 ###########################################################################################################################
206def MW_correlation_mult(nb_event_by_dir, cur_var, bin):          #########################################
207#1 ###########################################################################################################################
208    """ compute sigma and correlation from a sample of event """
209
210    #initialize data
211    stat_tool = Error_in_histo()
212    #data = stat_tool.charge_file(dir_pos(0)+'/0%d_card_1event_0.txt' % (cur_var), auto=0)
213    stat_tool.init_dict(bin)
214
215    cur_event = 0
216    while 1:
217        if cur_event%100 == 0: print 'status', cur_var, cur_event
218        step=numpy.random.poisson(1)
219        if cur_event + step > nb_event_by_dir:
220            break
221        hist=histo()
222        for k in range(cur_event, cur_event+step):
223            hist.charge_file(init_pos+'/0%d_card_1event_%s.txt' % (cur_var, k))
224        stat_tool.update_value(hist.out())
225        cur_event += step
226                             
227    stat_tool.compute_result()
228    stat_tool.write_result('%s/mean_var_%s' % (output_dir, cur_var),'%s/error_var_%s' % (output_dir, cur_var),'%s/correlation_var_%s' % (output_dir, cur_var))
229
230   
231
232#1 ########################################################################################################################### 
233class MW_chi_square:                                                                 ######################################### 
234#1 ###########################################################################################################################   
235
236    def __init__(self,total_event=380000):
237        self.read_event=0
238        self.total_event=total_event
239
240    def next_sample(self,N):
241        step=numpy.random.poisson(N)
242        if not self.read_event+step<self.total_event:
243            print 'end of event'
244            return {}
245        hist=histo()
246        for nb_data in range(self.read_event,self.read_event+step):
247            print 'charge all/event_'+str(nb_data)
248            hist.charge_file('all/event_'+str(nb_data))
249
250        return hist.out()
251
252
253
254bin = lambda min,gap: [ '%f' % float(min+i*gap) for i in range(0,50)]
255
256if '__main__'==__name__:
257
258    init_dir = 'where find the data ' #./Events/MYNAME/Graph/all_event
259    output_dir 'where to write the results'
260
261
262
263    MW_correlation_mult(5000, 1, bin(112, 24))
264    MW_correlation_mult(5000, 2, bin(5, 10))
265    MW_correlation_mult(5000, 3, bin(-0.98,0.04))
266