Mercurial > hg > medcouple

comparison pymedcouple @ 31:1a2bb52722c2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
pymedcouple: replace slow Python partsort with faster built-in full sort
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Sun, 18 Jan 2015 20:45:35 -0500
parents cd3094a08e03
children
comparison
equal deleted inserted replaced
30:2e277129e81b 31:1a2bb52722c2
1 #!/usr/bin/python 1 #!/usr/bin/python
2 2
3 import random 3 import random
4 4
5 from itertools import tee, izip 5 from itertools import tee, izip
6
7
8 def partsort(L, n):
9 """This is a partial sort of the weighted list L = [(a,w) for a in A,
10 w in W], so that the element at position n is in the correct spot.
11 Also known as quickselect. The elements of W are ignored for the
12 purposes of this partial sort.
13 """
14
15 beg = 0
16 end = len(L)
17
18 while True:
19 pivot = random.randint(beg, end - 1)
20 pivotval = L[pivot][0]
21 L[pivot], L[end - 1] = L[end - 1], L[pivot]
22
23 idx = beg
24 for i in xrange(beg, end):
25 if L[i][0] < pivotval:
26 L[idx], L[i] = L[i], L[idx]
27 idx += 1
28
29 L[idx], L[end-1] = L[end-1], L[idx]
30
31 if idx == n:
32 return L[idx]
33 elif idx < n:
34 beg = idx + 1
35 else:
36 end = idx
37
38 6
39 def wmedian(A, W): 7 def wmedian(A, W):
40 """This computes the weighted median of array A with corresponding 8 """This computes the weighted median of array A with corresponding
41 weights W. 9 weights W.
42 """ 10 """
50 end = n - 1 18 end = n - 1
51 19
52 while True: 20 while True:
53 mid = (beg + end)//2 21 mid = (beg + end)//2
54 22
55 partsort(AW, mid) 23 AW = sorted(AW, key = lambda x: x[0]) # A partial sort would suffice here
24
56 trial = AW[mid][0] 25 trial = AW[mid][0]
57 26
58 wleft = wright = 0 27 wleft = wright = 0
59 for (a, w) in AW: 28 for (a, w) in AW:
60 if a < trial: 29 if a < trial:
210 # This space is of size Rtot - Ltot which is <= n_plus 179 # This space is of size Rtot - Ltot which is <= n_plus
211 A = [] 180 A = []
212 for (i, (l, r)) in enumerate(izip(L, R)): 181 for (i, (l, r)) in enumerate(izip(L, R)):
213 for j in xrange(l, r + 1): 182 for j in xrange(l, r + 1):
214 A.append(h_kern(i, j)) 183 A.append(h_kern(i, j))
215 A.sort() 184
185 A.sort() # A partial sort would suffice here
216 A.reverse() 186 A.reverse()
217 187
218 Am = A[medc_idx - Ltot] 188 Am = A[medc_idx - Ltot]
219 189
220 return Am 190 return Am