medcouple: pymedcouple annotate

annotate pymedcouple @ 23:29a178b23219

Whitespace style fixes

author	Jordi Gutiérrez Hermoso <jordigh@octave.org>
date	Sat, 17 Jan 2015 20:04:46 -0500
parents	f5b4a2ab6204
children	cd3094a08e03

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
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	6 from itertools import tee, izip
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	7
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	8
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	9 def partsort(L, n):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	10 """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	11 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	12 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	13 purposes of this partial sort.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	14 """
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	15
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	16 beg = 0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	17 end = len(L)
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	18
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	19 while True:
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	20 pivot = random.randint(beg, end - 1)
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	21 pivotval = L[pivot][0]
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	22 L[pivot], L[end - 1] = L[end - 1], L[pivot]
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	23
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	24 idx = beg
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	25 for i in xrange(beg, end):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	26 if L[i][0] < pivotval:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	27 L[idx], L[i] = L[i], L[idx]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	28 idx += 1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	29
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	30 L[idx], L[end-1] = L[end-1], L[idx]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	31
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	32 if idx == n:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	33 return L[idx]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	34 elif idx < n:
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	35 beg = idx + 1
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	36 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	37 end = idx
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	38
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	39
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	40 def wmedian(A, W):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	41 """This computes the weighted median of array A with corresponding
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	42 weights W.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	43 """
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	44
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	45 AW = zip(A, W)
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	46 n = len(AW)
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 wtot = sum(W)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	49
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	50 beg = 0
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	51 end = n - 1
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	52
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	53 while True:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	54 mid = (beg + end)//2
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	55
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	56 partsort(AW, mid)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	57 trial = AW[mid][0]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	58
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	59 wleft = wright = 0
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	60 for (a, w) in AW:
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	61 if a < trial:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	62 wleft += w
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	63 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	64 # This also includes a == trial, i.e. the "middle"
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	65 # weight.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	66 wright += w
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	67
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	68 if 2*wleft > wtot:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	69 end = mid
21 d63e763f8ac1 wmedian: don't infinite loop if wleft == wright Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 20 diff changeset	70 elif 2*wright < wtot:
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	71 beg = mid
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	72 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	73 return trial
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	74
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	75
10 3d145bcc8694 Remove more debug code Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 9 diff changeset	76 def medcouple_1d(X, eps1 = 2-52, eps2 = 2-1022):
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	77 """Calculates the medcouple robust measure of skewness.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	78
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	79 Parameters
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	80 ----------
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	81 y : array-like, 1-d
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 Returns
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 mc : float
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	86 The medcouple statistic
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	87
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	88 .. [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	89 Skewness." Journal of Computational and Graphical Statistics, Vol.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	90 13, No. 4 (Dec., 2004), pp. 996- 1017
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 .. [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	93 in $X + Y$ and $X_1 + X_2 + \cdots X_m$". SIAM Journal of
9 e6bcaf38edaa Fix citation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 8 diff changeset	94 Computing, Vol. 7, No. 2 (May 1978), pp. 147-53.
4 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 """
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	97 # FIXME: Figure out what to do about NaNs.
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	98
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	99 n = len(X)
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	100 n2 = (n - 1)//2
4 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 < 3:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	103 return 0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	104
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	105 Z = sorted(X, reverse=True)
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 if n % 2 == 1:
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	108 Zmed = Z[n2]
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	109 else:
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	110 Zmed = (Z[n2] + Z[n2 + 1])/2
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	111
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	112 #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	113 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	114 return -1.0
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	115 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	116 return 1.0
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	117
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	118 # Centre Z wrt median, so that median(Z) = 0.
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	119 Z = [z - Zmed for z in Z]
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 # Scale inside [-0.5, 0.5], for greater numerical stability.
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	122 Zden = 2*max(Z[0], -Z[-1])
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	123 Z = [z/Zden for z in Z]
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	124 Zmed /= Zden
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	125
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	126 Zeps = eps1*(eps1 + abs(Zmed))
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	127
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	128 # 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	129 Zplus = [z for z in Z if z >= -Zeps]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	130 Zminus = [z for z in Z if Zeps >= z]
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 n_plus = len(Zplus)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	133 n_minus = len(Zminus)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	134
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	135
10 3d145bcc8694 Remove more debug code Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 9 diff changeset	136 def h_kern(i, j):
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	137 """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	138 Zplus and Zminus just defined above.
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 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	141 is the signum of their position.
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 """
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	144 a = Zplus[i]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	145 b = Zminus[j]
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	146
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	147 if abs(a - b) <= 2*eps2:
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	148 h = signum(n_plus - 1 - i - j)
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	149 else:
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	150 h = (a + b)/(a - b)
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	151
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	152 return h
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	153
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	154 # Init left and right borders
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	155 L = [0]*n_plus
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	156 R = [n_minus - 1]*n_plus
4 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 Ltot = 0
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	159 Rtot = n_minus*n_plus
22 f5b4a2ab6204 medcouple_1d: modify the definition of medc_idx to match R's definition Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 21 diff changeset	160 medc_idx = Rtot//2
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	161
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	162 # kth pair algorithm (Johnson & Mizoguchi)
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	163 while Rtot - Ltot > n_plus:
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	164
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	165 # First, compute the median inside the given bounds
6 e3b1dcc51e6a Fix tee usage Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 5 diff changeset	166 # (Be stingy, reuse same generator)
e3b1dcc51e6a Fix tee usage Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 5 diff changeset	167 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i])
e3b1dcc51e6a Fix tee usage Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 5 diff changeset	168
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	169 A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	170 W = [R[i] - L[i] + 1 for i in I2]
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	171 Am = wmedian(A, W)
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	172
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	173 Am_eps = eps1*(eps1 + abs(Am))
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	174
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	175 # Compute new left and right boundaries, based on the weighted
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	176 # median
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	177 P = []
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	178 Q = []
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	179
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	180 j = 0
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	181 for i in xrange(n_plus - 1, -1, -1):
20 8e02346789b2 Fix off by one error in P computation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13 diff changeset	182 while j < n_minus and h_kern(i, j) - Am > Am_eps:
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	183 j += 1
12 c81f6d263897 Spaces around plus and minus Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 11 diff changeset	184 P.append(j - 1)
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	185
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	186 P.reverse()
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	187
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	188 j = n_minus - 1
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	189 for i in xrange(0, n_plus):
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	190 while j >= 0 and h_kern(i, j) - Am < -Am_eps:
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	191 j -= 1
12 c81f6d263897 Spaces around plus and minus Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 11 diff changeset	192 Q.append(j + 1)
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	193
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	194 # 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	195 # the whole matrix may be.
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	196 sumP = sum(P) + len(P)
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	197 sumQ = sum(Q)
11 4c7d7ff28a0e Rename mid_idx to medc_idx Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 10 diff changeset	198
4c7d7ff28a0e Rename mid_idx to medc_idx Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 10 diff changeset	199 if medc_idx <= sumP - 1:
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	200 R = P
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	201 Rtot = sumP
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	202 else:
11 4c7d7ff28a0e Rename mid_idx to medc_idx Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 10 diff changeset	203 if medc_idx > sumQ - 1:
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	204 L = Q
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	205 Ltot = sumQ
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	206 else:
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	207 return Am
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	208
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	209 # Didn't find the median, but now we have a very small search
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	210 # space to find it in, just between the left and right boundaries.
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	211 # This space is of size Rtot - Ltot which is <= n_plus
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	212 A = []
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	213 for (i, (l, r)) in enumerate(izip(L, R)):
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	214 for j in xrange(l, r + 1):
10 3d145bcc8694 Remove more debug code Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 9 diff changeset	215 A.append(h_kern(i, j))
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	216 A.sort()
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	217 A.reverse()
6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	218
23 29a178b23219 Whitespace style fixes Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 22 diff changeset	219 Am = A[medc_idx - Ltot]
11 4c7d7ff28a0e Rename mid_idx to medc_idx Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 10 diff changeset	220
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	221 return Am
13 077261db7a58 Whitespace cleanup Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 12 diff changeset	222
4 74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	223
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	224 def signum(x):
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	225 if x > 0:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	226 return 1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	227 if x < 0:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	228 return -1
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	229 else:
74d0d08dbc95 Add Python implementation Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: diff changeset	230 return 0
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	231
8 b3e878bb793d Remove debug statements Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 7 diff changeset	232
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	233 def main():
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	234 import sys
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	235 fname = sys.argv[1]
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	236 with open(fname) as f:
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	237 data = [float(x) for x in f.readlines() if x.strip() != ""]
8 b3e878bb793d Remove debug statements Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 7 diff changeset	238
7 6c7c7dc9d8ef hairball Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 6 diff changeset	239 print "%.16g" % medcouple_1d(data)
5 cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	240
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	241 if __name__ == "__main__":
cbe17f888c79 Make it executable Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 4 diff changeset	242 main()

Mercurial > hg > medcouple

annotate pymedcouple @ 23:29a178b23219