Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
comment
[simgrid.git] / contrib / network_model / regress.py
1 #!/usr/bin/env python
2
3 import sys
4 from math import sqrt
5
6
7 if len(sys.argv) != 2 and len(sys.argv) != 4:
8         print("Usage : %s datafile", sys.argv[0])
9         print("or    : %s datafile p1 p2", sys.argv[0])
10         print("where : p1 < p2 belongs to sizes in datafiles")
11         sys.exit(-1)
12
13 if len(sys.argv) == 4:
14         p1=int(sys.argv[2])
15         p2=int(sys.argv[3])
16
17 ##-----------------------------------------
18 ## avg : return average of a list of values
19 ## param l list of values
20 ##-----------------------------------------
21 def avg (l):
22         sum=0
23         for e in l:
24                 sum=sum+e;
25         return sum/len(l)
26
27 ##-------------------------------------------------
28 ## cov : covariance 
29 ## param X first data vector (..x_i..)
30 ## param Y second data vector (..x_i..)
31 ## = 1/n \Sum_{i=1}^n (x_i - avg(x)) * (y_i - avg(y))
32 ##--------------------------------------------------
33 def cov(X,Y):
34
35         n=len(X)   #  n=len(X)=len(Y)
36         avg_X = avg( X )
37         avg_Y = avg( Y )
38         S_XY=0
39         for i in range(n):
40                 S_XY = S_XY + ((X[i]-avg_X)*(Y[i]-avg_Y))
41
42         return (S_XY/n) 
43
44
45 ##----------------------------------
46 ## variance : variance 
47 ## param X data vector ( ..x_i.. )
48 ## (S_X)^2 = (Sum ( x_i - avg(x) )^2 ) / n
49 ##----------------------------------
50 def variance( X ):
51         S_X2 = 0
52         n = len( X )
53         avg_X = avg ( X )
54         for i in range(n):
55                 S_X2 = S_X2 + ((X[i] - avg_X)**2)
56
57         return (S_X2/n)
58
59 ##-----------------------------------------------------------------------------------------------
60 ## correl_split_weighted : compute regression on each segment and
61 ## return the weigthed sum of correlation coefficients
62 ## param X first data vector (..x_i..)
63 ## param Y second data vector (..x_i..)
64 ## param segments list of pairs (i,j) where i refers to the ith value in X, and  jth value in X
65 ## return (C,[(i1,j1,X[i1],X[j1]), (i2,j2,X[i2],X[j2]), ....]
66 ##    where i1,j1 is the first segment,  c1 the correlation coef on this segment, n1 the number of values
67 ##          i2,j2 is the second segment, c2 the correlation coef on this segment, n2 the number of values
68 ##          ...
69 ##    and C=c1/n1+c2/n2+...
70 ##-----------------------------------------------------------------------------------------------
71 def correl_split_weighted( X , Y , segments ):
72         # expects segments = [(0,i1-1),(i1-1,i2-1),(i2,len-1)] 
73         correl = list();
74         interv = list();  # regr. line coeffs and range
75         glob_corr=0
76         sum_nb_val=0
77         for (start,stop) in segments:
78                 sum_nb_val = sum_nb_val + stop - start;
79                 #if start==stop :
80                 #       return 0
81                 S_XY= cov( X [start:stop+1], Y [start:stop+1] )
82                 S_X2 = variance(  X [start:stop+1] ) 
83                 S_Y2 = variance(  Y [start:stop+1] )  # to compute correlation
84                 if S_X2*S_Y2 == 0:
85                         return (0,[])
86                 c = S_XY/(sqrt(S_X2)*sqrt(S_Y2)) 
87                 a = S_XY/S_X2                           # regr line coeffs
88                 b= avg ( Y[start:stop+1] ) - a * avg( X[start:stop+1] )
89                 print("   range [%d,%d] corr=%f, coeff det=%f [a=%f, b=%f]" % (X[start],X[stop],c,c**2,a, b)) 
90                 correl.append( (c, stop-start) );   # store correl. coef + number of values (segment length)
91                 interv.append( (a,b, X[start],X[stop]) );
92         
93         for (c,l) in correl:
94                 glob_corr = glob_corr + (l/sum_nb_val)*c  # weighted product of correlation 
95                 print('-- %f * %f' % (c,l/sum_nb_val))
96                 
97         print("-> glob_corr=%f\n" % glob_corr)
98         return (glob_corr,interv);
99         
100
101
102
103 ##-----------------------------------------------------------------------------------------------
104 ## correl_split : compute regression on each segment and
105 ## return the product of correlation coefficient
106 ## param X first data vector (..x_i..)
107 ## param Y second data vector (..x_i..)
108 ## param segments list of pairs (i,j) where i refers to the ith value in X, and  jth value in X
109 ## return (C,[(i1,j1,X[i1],X[j1]), (i2,j2,X[i2],X[j2]), ....]
110 ##    where i1,j1 is the first segment,  c1 the correlation coef on this segment,
111 ##          i2,j2 is the second segment, c2 the correlation coef on this segment,
112 ##          ...
113 ##    and C=c1*c2*...
114 ##-----------------------------------------------------------------------------------------------
115 def correl_split( X , Y , segments ):
116         # expects segments = [(0,i1-1),(i1-1,i2-1),(i2,len-1)] 
117         correl = list();
118         interv = list();  # regr. line coeffs and range
119         glob_corr=1
120         for (start,stop) in segments:
121                 #if start==stop :
122                 #       return 0
123                 S_XY= cov( X [start:stop+1], Y [start:stop+1] )
124                 S_X2 = variance(  X [start:stop+1] ) 
125                 S_Y2 = variance(  Y [start:stop+1] )  # to compute correlation
126                 if S_X2*S_Y2 == 0:
127                         return (0,[])
128                 c = S_XY/(sqrt(S_X2)*sqrt(S_Y2)) 
129                 a = S_XY/S_X2                           # regr line coeffs
130                 b= avg ( Y[start:stop+1] ) - a * avg( X[start:stop+1] )
131                 print("   range [%d,%d] corr=%f, coeff det=%f [a=%f, b=%f]" % (X[start],X[stop],c,c**2,a, b)) 
132                 correl.append( (c, stop-start) );   # store correl. coef + number of values (segment length)
133                 interv.append( (a,b, X[start],X[stop]) );
134         
135         for (c,l) in correl:
136                 glob_corr = glob_corr * c  # product of correlation coeffs
137         print("-> glob_corr=%f\n" % glob_corr)
138         return (glob_corr,interv);
139         
140
141
142 ##-----------------------------------------------------------------------------------------------
143 ## main
144 ##-----------------------------------------------------------------------------------------------
145 sum=0
146 nblines=0
147 skampidat = open(sys.argv[1], "r")
148
149 timings = []
150 sizes = []
151 for line in skampidat:
152         l = line.split();
153         if line[0] != '#' and len(l)>=3:   # is it a comment ?
154
155 ## expected format
156 ## ---------------
157 #count= 8388608  8388608  144916.1       7.6       32  144916.1  143262.0
158 #("%s %d %d %f %f %d %f %f\n" % (countlbl, count, countn, time, stddev, iter, mini, maxi)
159                 timings.append (float(l[3]))
160                 sizes.append(int(l[1]))
161                 nblines=nblines+1
162
163
164
165 ##----------------------- search for best break points-----------------
166 ## example
167 ## p1=2048  -> p1inx=11  delta=3 -> [8;14]
168 ##                                              8 : segments[(0,7),(8,13),(13,..)]
169 ##                                              ....                                            
170 ## p2=65536 -> p2inx=16  delta=3 -> [13;19]
171
172 if len(sys.argv) == 4:
173
174         p1inx = sizes.index( p1 );
175         p2inx = sizes.index( p2 );
176         max_glob_corr = 0;
177         max_p1inx = p1inx
178         max_p2inx = p2inx
179
180         ## tweak parameters here to extend/reduce search
181         search_p1 = 30          # number of values to search +/- around p1
182         search_p2 = 45          # number of values to search +/- around p2
183         min_seg_size = 3
184
185         lb1 = max(1, p1inx-search_p1)           
186         ub1 = min(p1inx+search_p1,search_p1, p2inx);
187         lb2 = max(p1inx,p2inx-search_p2)        # breakpoint +/- delta
188         ub2 = min(p2inx+search_p2,len(sizes)-1);    
189
190         print("** evaluating over \n");
191         print("interv1:\t %d <--- %d ---> %d" % (sizes[lb1],p1,sizes[ub1]))
192         print("rank:   \t (%d)<---(%d)--->(%d)\n" % (lb1,p1inx,ub1))
193         print("interv2:\t\t %d <--- %d ---> %d" % (sizes[lb2],p2,sizes[ub2]))
194         print("rank:   \t\t(%d)<---(%d)--->(%d)\n" % (lb2,p2inx,ub2))
195         for i in range(lb1,ub1+1):
196                 for j in range(lb2,ub2+1):
197                         if i<j:         # segments must not overlap
198                                 if i+1 >=min_seg_size and j-i+1 >= min_seg_size and len(sizes)-1-j >= min_seg_size : # not too small segments
199                                         print("** i=%d,j=%d" % (i,j))
200                                         segments = [(0,i),(i,j),(j,len(sizes)-1)]
201                                         (glob_cor, interv) = correl_split( sizes, timings, segments)
202                                         if ( glob_cor > max_glob_corr):
203                                                 max_glob_corr = glob_cor
204                                                 max_interv = interv
205
206         for (a,b,i,j) in max_interv:
207                 print("** OPT: [%d .. %d]" % (i,j))
208         print("** Product of correl coefs = %f" % (max_glob_corr))
209
210         print("#-------------------- cut here the gnuplot code -----------------------------------------------------------\n");
211         preamble='set output "regr.eps"\n\
212 set terminal postscript eps color\n\
213 set key left\n\
214 set xlabel "Each message size in bytes"\n\
215 set ylabel "Time in us"\n\
216 set logscale x\n\
217 set logscale y\n\
218 set grid'
219
220         print(preamble);
221         print('plot "%s" u 3:4:($5) with errorbars title "skampi traces %s",\\' % (sys.argv[1],sys.argv[1]));
222         for (a,b,i,j) in max_interv:
223                 print('"%s" u (%d<=$3 && $3<=%d? $3:0/0):(%f*($3)+%f) w linespoints title "regress. %s-%s bytes",\\' % (sys.argv[1],i,j,a,b,i,j))
224         
225         print("#-------------------- /cut here the gnuplot code -----------------------------------------------------------\n");
226
227
228 else:
229         print('\n** Linear regression on %d values **\n' % (nblines))
230         print('\n   sizes=',sizes,'\n\n')
231         avg_sizes = avg( sizes )
232         avg_timings = avg( timings )
233         print("avg_timings=%f, avg_sizes=%f, nblines=%d\n" % (avg_timings,avg_sizes,nblines))
234
235         S_XY= cov( sizes, timings )
236         S_X2 = variance(  sizes ) 
237         S_Y2 = variance(  timings )  # to compute correlation
238
239         a = S_XY/S_X2;
240         correl = S_XY/(sqrt(S_X2)*sqrt(S_Y2))  # corealation coeff (Bravais-Pearson)
241
242
243         b= avg_timings - a * avg_sizes
244         print("[S_XY=%f, S_X2=%f]\n[correlation=%f, coeff det=%f]\n[a=%f, b=%f]\n" % (S_XY, S_X2, correl,correl**2,a, b)) 
245         
246