medcouple: pymedcouple annotate

annotate pymedcouple @ 5:cbe17f888c79

Make it executable

author	Jordi Gutiérrez Hermoso <jordigh@octave.org>
date	Tue, 13 Jan 2015 17:00:49 -0500
parents	medcouple.py@74d0d08dbc95
children	e3b1dcc51e6a

rev	line source
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	1 #!/usr/bin/python
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	2
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	3 import random
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	4 import numpy as np
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 from itertools import tee
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	7
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	8 def partsort(L, n):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	9 """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	10 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	11 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	12 purposes of this partial sort.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	13 """
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	14
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	15 beg = 0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	16 end = len(L)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	17
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	18 while True:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	19 pivot = random.randint(beg,end-1)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	20 pivotval = L[pivot][0]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	21 L[pivot], L[end-1] = L[end-1], L[pivot]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	22
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	23 idx = beg
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	24 for i in xrange(beg, end):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	25 if L[i][0] < pivotval:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	26 L[idx], L[i] = L[i], L[idx]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	27 idx += 1
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 L[idx], L[end-1] = L[end-1], L[idx]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	30
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	31 if idx == n:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	32 return L[idx]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	33 elif idx < n:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	34 beg = idx+1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	35 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	36 end = idx
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	37
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	38 def wmedian(A, W):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	39 """This computes the weighted median of array A with corresponding
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	40 weights W.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	41 """
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	42
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	43 AW = zip(A,W)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	44 n = len(AW)
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 wtot = sum(W)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	47
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	48 beg = 0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	49 end = n
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	50
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	51 while True:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	52 mid = (beg + end)//2
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	53
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	54 partsort(AW, mid)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	55 trial = AW[mid][0]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	56
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	57 wleft = wright = 0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	58 for (a,w) in AW:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	59 if a < trial:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	60 wleft += w
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	61 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	62 # This also includes a == trial, i.e. the "middle"
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	63 # weight.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	64 wright += w
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	65
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	66 if 2*wleft > wtot:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	67 end = mid
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	68 elif 2*wright <= wtot:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	69 beg = mid
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	70 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	71 return trial
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	72
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	73 def medcouple_1d(X, eps1 = 2-52, eps2 = 2-1022 ):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	74 """Calculates the medcouple robust measure of skewness.
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 Parameters
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 y : array-like, 1-d
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 Returns
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	81 -------
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	82 mc : float
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	83 The medcouple statistic
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	84
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	85 .. [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	86 Skewness." Journal of Computational and Graphical Statistics, Vol.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	87 13, No. 4 (Dec., 2004), pp. 996- 1017
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	88
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	89 .. [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	90 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	91 Computing, 7, 147-53.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	92
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 # FIXME: Figure out what to do about NaNs.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	95
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	96 n = len(X)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	97 n2 = (n-1)//2
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	98
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	99 if n < 3:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	100 return 0
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 Z = sorted(X, reverse=True)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	103
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	104 if n % 2 == 1:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	105 Zmed = Z[n2 + 1]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	106 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	107 Zmed = (Z[n2] + Z[n2 + 1])/2
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	108
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	109 #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	110 if abs(Z[0] - 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 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	113 return 1.0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	114
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	115 # Centre Z wrt median, so that median(Z) = 0.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	116 Z = [z - Zmed for z in Z]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	117
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	118 # Scale inside [-0.5, 0.5], for greater numerical stability.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	119 Zden = -2*max(-Z[0], Z[-1])
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	120 Z = [z/Zden for z in Z]
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	121 Zmed /= Zden
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 Zeps = eps1*(eps1 + abs(Zmed))
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	124
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	125 # 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	126 Zplus = [z for z in Z if z >= -Zeps]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	127 Zminus = [z for z in Z if Zeps >= z]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	128
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	129 n_plus = len(Zplus)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	130 n_minus = len(Zminus)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	131
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	132 def h_kern(i, j):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	133 """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	134 Zplus and Zminus just defined above.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	135
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	136 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	137 is the signum of their position.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	138
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	139 """
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	140 a = Zplus[i]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	141 b = Zminus[j]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	142
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	143 if abs(a - b) <= 2*eps2:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	144 return signum(n_plus - 1 - i - j)
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 return (a+b)/(a-b)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	147
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	148 # Init left and right borders
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	149 L = [0]*n_plus
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	150 R = [n_plus-1]*n_plus
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	151
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	152 Ltot = 0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	153 Rtot = n_plus*n_minus
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 mid_idx = Rtot//2
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	156
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	157 # kth pair algorithm (Johnson & Mizoguchi)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	158 while Rtot - Ltot > n_plus:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	159
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	160 # First, compute the median inside the given bounds
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	161 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	162 I2 = tee(I1) # Be stingy, reuse same generator
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	163 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	164 W = [R[i] - L[i] for i in I2]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	165 Am = wmedian(A,W)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	166
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	167 # Compute new left and right boundaries, based on the weighted
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	168 # median
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	169 P = Q = []
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	170
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	171 j = -1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	172 for i in xrange(n_plus - 1, -1, -1):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	173 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	174 j += 1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	175 P.append(j)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	176 P.reverse()
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	177
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	178 j = n_plus
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	179 for i in xrange(0, n_plus):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	180 while j > -1 and h_kern(i, j) < Am:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	181 j -= 1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	182 Q.append(j)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	183
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	184 # 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	185 # the whole matrix may be.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	186 sumP = sum(P)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	187 sumQ = sum(Q)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	188 if mid_idx <= sumP:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	189 R = P
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	190 Rtot = sumP
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	191 elif mid_idx > sumQ:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	192 L = Q
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	193 Ltot = sumQ
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	194 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	195 return Am
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	196
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	197 # Uh, implement this part...
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	198 import ipdb; ipdb.set_trace()
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	199
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	200
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	201 def signum(x):
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 if x < 0:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	205 return -1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	206 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	207 return 0
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	208
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	209 def main():
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	210 import sys
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	211 fname = sys.argv[1]
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	212 with open(fname) as f:
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	213 data = [float(x) for x in f.readlines() if x.strip() != ""]
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	214 print medcouple_1d(data)
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	215
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	216 if __name__ == "__main__":
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	217 main()

Mercurial > hg > medcouple

annotate pymedcouple @ 5:cbe17f888c79