Mercurial > hg > medcouple

annotate medcouple.py @ 4:74d0d08dbc95

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Python implementation
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Tue, 13 Jan 2015 16:55:48 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
1 import random
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
2 import numpy as np
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4 from itertools import tee
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6 def partsort(L, n):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7 """This is a partial sort of the weighted list L = [(a,w) for a in A,
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8 w in W], so that the element at position n is in the correct spot.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9 Also known as quickselect. The elements of W are ignored for the
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10 purposes of this partial sort.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11 """
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 beg = 0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14 end = len(L)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16 while True:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 pivot = random.randint(beg,end-1)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18 pivotval = L[pivot][0]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19 L[pivot], L[end-1] = L[end-1], L[pivot]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 idx = beg
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22 for i in xrange(beg, end):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23 if L[i][0] < pivotval:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
24 L[idx], L[i] = L[i], L[idx]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25 idx += 1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27 L[idx], L[end-1] = L[end-1], L[idx]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 if idx == n:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30 return L[idx]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31 elif idx < n:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32 beg = idx+1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 else:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 end = idx
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 def wmedian(A, W):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37 """This computes the weighted median of array A with corresponding
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 weights W.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 """
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41 AW = zip(A,W)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42 n = len(AW)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 wtot = sum(W)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46 beg = 0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
47 end = n
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 while True:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50 mid = (beg + end)//2
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52 partsort(AW, mid)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
53 trial = AW[mid][0]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55 wleft = wright = 0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56 for (a,w) in AW:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57 if a < trial:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
58 wleft += w
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
59 else:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60 # This also includes a == trial, i.e. the "middle"
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61 # weight.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62 wright += w
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64 if 2*wleft > wtot:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65 end = mid
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66 elif 2*wright <= wtot:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67 beg = mid
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68 else:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
69 return trial
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
70
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
71 def medcouple_1d(X, eps1 = 2**-52, eps2 = 2**-1022 ):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
72 """Calculates the medcouple robust measure of skewness.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
73
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 Parameters
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 ----------
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
76 y : array-like, 1-d
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78 Returns
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 -------
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80 mc : float
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81 The medcouple statistic
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83 .. [1] G. Brys, M. Hubert, and A. Struyf "A Robust Measure of
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84 Skewness." Journal of Computational and Graphical Statistics, Vol.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 13, No. 4 (Dec., 2004), pp. 996- 1017
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
86
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
87 .. [2] D. B. Johnson and T. Mizoguchi "Selecting the Kth Element
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
88 in $X + Y$ and $X_1 + X_2 + \cdots X_m$". SIAM Journal of
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
89 Computing, 7, 147-53.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
90
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
91 """
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
92 # FIXME: Figure out what to do about NaNs.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
93
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
94 n = len(X)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
95 n2 = (n-1)//2
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
96
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
97 if n < 3:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
98 return 0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
99
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 Z = sorted(X, reverse=True)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102 if n % 2 == 1:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 Zmed = Z[n2 + 1]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 else:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105 Zmed = (Z[n2] + Z[n2 + 1])/2
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107 #Check if the median is at the edges up to relative epsilon
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109 return -1.0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 return 1.0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
113 # Centre Z wrt median, so that median(Z) = 0.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114 Z = [z - Zmed for z in Z]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
115
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116 # Scale inside [-0.5, 0.5], for greater numerical stability.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
117 Zden = -2*max(-Z[0], Z[-1])
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
118 Z = [z/zden for z in Zden]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 Zmed /= Zden
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
120
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121 Zeps = eps1*(eps1 + abs(Zmed))
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
123 # These overlap on the entries that are tied with the median
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
124 Zplus = [z for z in Z if z >= -Zeps]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 Zminus = [z for z in Z if Zeps >= z]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
126
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
127 n_plus = len(Zplus)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
128 n_minus = len(Zminus)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130 def h_kern(i, j):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
131 """Kernel function h for the medcouple, closing over the values of
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
132 Zplus and Zminus just defined above.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
133
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
134 In case a and be are within epsilon of the median, the kernel
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
135 is the signum of their position.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
136
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
137 """
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
138 a = Zplus[i]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
139 b = Zminus[j]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
140
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
141 if abs(a - b) <= 2*eps2:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
142 return signum(n_plus - 1 - i - j)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
143
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
144 return (a+b)/(a-b)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
145
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
146 # Init left and right borders
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
147 L = [0]*n_plus
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
148 R = [n_plus-1]*n_plus
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
149
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
150 Ltot = 0
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
151 Rtot = n_plus*n_minus
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
152
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
153 mid_idx = Rtot//2
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
154
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
155 # kth pair algorithm (Johnson & Mizoguchi)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
156 while Rtot - Ltot > n_plus:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
157
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
158 # First, compute the median inside the given bounds
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
159 I1 = (i for i in xrange(0, n_plus) if L[i] <= R[i])
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
160 I2 = tee(I1) # Be stingy, reuse same generator
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
161 A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
162 W = [R[i] - L[i] for i in I2]
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
163 Am = wmedian(A,W)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
164
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
165 # Compute new left and right boundaries, based on the weighted
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
166 # median
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
167 P = Q = []
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
168
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
169 j = -1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
170 for i in xrange(n_plus - 1, -1, -1):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
171 while j < n_plus - 1 and h_kern(i, j) > Am:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
172 j += 1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
173 P.append(j)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
174 P.reverse()
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
175
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
176 j = n_plus
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
177 for i in xrange(0, n_plus):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
178 while j > -1 and h_kern(i, j) < Am:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
179 j -= 1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
180 Q.append(j)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
181
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
182 # Check on which side of those bounds the desired median of
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
183 # the whole matrix may be.
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
184 sumP = sum(P)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
185 sumQ = sum(Q)
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
186 if mid_idx <= sumP:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
187 R = P
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
188 Rtot = sumP
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
189 elif mid_idx > sumQ:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
190 L = Q
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
191 Ltot = sumQ
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
192 else:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
193 return Am
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
194
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
195 # Uh, implement this part...
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
196 import ipdb; ipdb.set_trace()
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
197
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
198
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
199 def signum(x):
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
200 if x > 0:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
201 return 1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
202 if x < 0:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
203 return -1
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
204 else:
74d0d08dbc95 Add Python implementation
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
205 return 0