Mercurial > hg > medcouple

comparison pymedcouple @ 8:b3e878bb793d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove debug statements
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Thu, 15 Jan 2015 16:25:44 -0500
parents 6c7c7dc9d8ef
children e6bcaf38edaa
comparison
equal deleted inserted replaced
7:6c7c7dc9d8ef 8:b3e878bb793d
110 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)): 110 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)):
111 return -1.0 111 return -1.0
112 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)): 112 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
113 return 1.0 113 return 1.0
114 114
115 if debug:
116 print "Zmed = ", Zmed
117
118 # Centre Z wrt median, so that median(Z) = 0. 115 # Centre Z wrt median, so that median(Z) = 0.
119 Z = [z - Zmed for z in Z] 116 Z = [z - Zmed for z in Z]
120 117
121 # Scale inside [-0.5, 0.5], for greater numerical stability. 118 # Scale inside [-0.5, 0.5], for greater numerical stability.
122 Zden = 2*max(Z[0], -Z[-1]) 119 Zden = 2*max(Z[0], -Z[-1])
123 if debug:
124 print "Zden =", Zden
125
126 Z = [z/Zden for z in Z] 120 Z = [z/Zden for z in Z]
127 Zmed /= Zden 121 Zmed /= Zden
128 122
129 Zeps = eps1*(eps1 + abs(Zmed)) 123 Zeps = eps1*(eps1 + abs(Zmed))
130 124
131 # These overlap on the entries that are tied with the median 125 # These overlap on the entries that are tied with the median
132 Zplus = [z for z in Z if z >= -Zeps] 126 Zplus = [z for z in Z if z >= -Zeps]
133 Zminus = [z for z in Z if Zeps >= z] 127 Zminus = [z for z in Z if Zeps >= z]
134 128
135 if debug:
136 for zp in Zplus:
137 print "Zplus : %g" % zp
138 for zm in Zminus:
139 print "Zminus: %g" % zm
140
141 n_plus = len(Zplus) 129 n_plus = len(Zplus)
142 n_minus = len(Zminus) 130 n_minus = len(Zminus)
143
144 if debug:
145 print "n_plus =", n_plus
146 print "n_minus = ", n_minus
147 131
148 def h_kern(i, j, debug=False): 132 def h_kern(i, j, debug=False):
149 """Kernel function h for the medcouple, closing over the values of 133 """Kernel function h for the medcouple, closing over the values of
150 Zplus and Zminus just defined above. 134 Zplus and Zminus just defined above.
151 135
159 if abs(a - b) <= 2*eps2: 143 if abs(a - b) <= 2*eps2:
160 h = signum(n_plus - 1 - i - j) 144 h = signum(n_plus - 1 - i - j)
161 else: 145 else:
162 h = (a+b)/(a-b) 146 h = (a+b)/(a-b)
163 147
164 if False:
165 print "a = {a}, b = {b}, h({i},{j}) = {h}".format(**locals())
166
167 return h 148 return h
168 149
169 # Init left and right borders 150 # Init left and right borders
170 L = [0]*n_plus 151 L = [0]*n_plus
171 R = [n_minus - 1]*n_plus 152 R = [n_minus - 1]*n_plus
172 153
173 Ltot = 0 154 Ltot = 0
174 Rtot = n_minus*n_plus 155 Rtot = n_minus*n_plus
175 mid_idx = (Rtot-1)//2 156 mid_idx = (Rtot-1)//2
176 157
177 if debug:
178 print "mid_idx = ", mid_idx
179
180 # kth pair algorithm (Johnson & Mizoguchi) 158 # kth pair algorithm (Johnson & Mizoguchi)
181 while Rtot - Ltot > n_plus: 159 while Rtot - Ltot > n_plus:
182 if debug: 160
183 print "L = ", L
184 print "R = ", R
185
186 # First, compute the median inside the given bounds 161 # First, compute the median inside the given bounds
187 # (Be stingy, reuse same generator) 162 # (Be stingy, reuse same generator)
188 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i]) 163 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i])
189 164
190 A = [h_kern(i, (L[i] + R[i])//2) for i in I1] 165 A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
191 W = [R[i] - L[i] + 1 for i in I2] 166 W = [R[i] - L[i] + 1 for i in I2]
192 Am = wmedian(A,W) 167 Am = wmedian(A,W)
193 if debug:
194 print "Am = ", Am
195 168
196 Am_eps = eps1*(eps1 + abs(Am)) 169 Am_eps = eps1*(eps1 + abs(Am))
197 170
198 # Compute new left and right boundaries, based on the weighted 171 # Compute new left and right boundaries, based on the weighted
199 # median 172 # median
216 189
217 # Check on which side of those bounds the desired median of 190 # Check on which side of those bounds the desired median of
218 # the whole matrix may be. 191 # the whole matrix may be.
219 sumP = sum(P) + len(P) 192 sumP = sum(P) + len(P)
220 sumQ = sum(Q) 193 sumQ = sum(Q)
221 if debug:
222 print "P: ", sumP, P
223 print "Q: ", sumQ, Q
224 print
225 194
226 if mid_idx <= sumP - 1: 195 if mid_idx <= sumP - 1:
227 R = P 196 R = P
228 Rtot = sumP 197 Rtot = sumP
229 else: 198 else:
234 return Am 203 return Am
235 204
236 # Didn't find the median, but now we have a very small search 205 # Didn't find the median, but now we have a very small search
237 # space to find it in, just between the left and right boundaries. 206 # space to find it in, just between the left and right boundaries.
238 # This space is of size Rtot - Ltot which is <= n_plus 207 # This space is of size Rtot - Ltot which is <= n_plus
239 if debug:
240 print "L =", L
241 print "R =", R
242 A = [] 208 A = []
243 for (i, (l, r)) in enumerate(izip(L, R)): 209 for (i, (l, r)) in enumerate(izip(L, R)):
244 for j in xrange(l, r + 1): 210 for j in xrange(l, r + 1):
245 A.append(h_kern(i, j, debug=debug)) 211 A.append(h_kern(i, j, debug=debug))
246 A.sort() 212 A.sort()
247 A.reverse() 213 A.reverse()
248 214
249 if debug:
250 print "len(A) = ", len(A)
251 Am = A[mid_idx - Ltot] 215 Am = A[mid_idx - Ltot]
252 return Am 216 return Am
253 217
254 218
255 def signum(x): 219 def signum(x):
258 if x < 0: 222 if x < 0:
259 return -1 223 return -1
260 else: 224 else:
261 return 0 225 return 0
262 226
227
263 def main(): 228 def main():
264 import sys 229 import sys
265 fname = sys.argv[1] 230 fname = sys.argv[1]
266 with open(fname) as f: 231 with open(fname) as f:
267 data = [float(x) for x in f.readlines() if x.strip() != ""] 232 data = [float(x) for x in f.readlines() if x.strip() != ""]
233
268 print "%.16g" % medcouple_1d(data) 234 print "%.16g" % medcouple_1d(data)
269 235
270 if __name__ == "__main__": 236 if __name__ == "__main__":
271 main() 237 main()