Mercurial > hg > medcouple

annotate talk/code/show_mc.py @ 66:916349600c4b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
algorithm-visualisation code
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Wed, 18 May 2016 17:03:43 -0400
parents 364a0c9df823
children 1aa38f4846e0
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
1 #!/usr/bin/python
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
2 # -*- coding: utf-8 -*-
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4 # medcouple.py ---
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6 # Copyright © 2015 Jordi Gutiérrez Hermoso <jordigh@octave.org>
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8 # Author: Jordi Gutiérrez Hermoso
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10 # This program is free software; you can redistribute it and/or
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11 # modify it under the terms of the GNU General Public License
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12 # as published by the Free Software Foundation; either version 3
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 # of the License, or (at your option) any later version.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15 # This program is distributed in the hope that it will be useful,
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18 # GNU General Public License for more details.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 # You should have received a copy of the GNU General Public License
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 # along with this program. If not, see <http://www.gnu.org/licenses/>.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
24 import numpy as np
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26 from itertools import tee, izip
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28 def wmedian(A, W):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 """This computes the weighted median of array A with corresponding
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30 weights W.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31 """
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 AW = zip(A, W)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 n = len(AW)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 wtot = sum(W)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 beg = 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 end = n - 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41 while True:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42 mid = (beg + end)//2
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 AW = sorted(AW, key = lambda x: x[0]) # A partial sort would suffice here
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46 trial = AW[mid][0]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
47
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48 wleft = wright = 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 for (a, w) in AW:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50 if a < trial:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51 wleft += w
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
53 # This also includes a == trial, i.e. the "middle"
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 # weight.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55 wright += w
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57 if 2*wleft > wtot:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
58 end = mid
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
59 elif 2*wright < wtot:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60 beg = mid
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62 return trial
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65 def medcouple_1d(X, eps1 = 2**-52, eps2 = 2**-1022):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66 """Calculates the medcouple robust measure of skewness.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68 Parameters
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
69 ----------
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
70 y : array-like, 1-d
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
71
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
72 Returns
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
73 -------
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 mc : float
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 The medcouple statistic
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
76
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77 .. [1] G. Brys, M. Hubert, and A. Struyf "A Robust Measure of
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78 Skewness." Journal of Computational and Graphical Statistics, Vol.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 13, No. 4 (Dec., 2004), pp. 996- 1017
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81 .. [2] D. B. Johnson and T. Mizoguchi "Selecting the Kth Element
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82 in $X + Y$ and $X_1 + X_2 + \cdots X_m$". SIAM Journal of
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83 Computing, Vol. 7, No. 2 (May 1978), pp. 147-53.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 """
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
86
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
87 # Yeah, yeah, leaky file handles; it'll be handled by the time the
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
88 # program closes.
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
89 Hfile = open("H", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
90 Amfile = open("Am", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
91 Rfile = open("R", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
92 Lfile = open("L", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
93 Pfile = open("P", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
94 Qfile = open("Q", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
95 RemFile = open("remaining", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
96 Ptotfile = open("Ptotal", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
97 Qtotfile = open("Qtotal", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
98 jackpot = open("jackpot", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
99
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 # FIXME: Figure out what to do about NaNs.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102 n = len(X)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 n2 = (n - 1)//2
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105 if n < 3:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 return 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 Z = sorted([float(x) for x in X], reverse=True)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110 if n % 2 == 1:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 Zmed = Z[n2]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
113 Zmed = (Z[n2] + Z[n2 + 1])/2
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
115 #Check if the median is at the edges up to relative epsilon
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
117 return -1.0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
118 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 return 1.0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
120
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121 # Centre Z wrt median, so that median(Z) = 0.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122 Z = [z - Zmed for z in Z]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
123
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
124 # Scale inside [-0.5, 0.5], for greater numerical stability.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 Zden = 2*max(Z[0], -Z[-1])
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
126 Z = [z/Zden for z in Z]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
127 Zmed /= Zden
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
128
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129 Zeps = eps1*(eps1 + abs(Zmed))
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
131 # These overlap on the entries that are tied with the median
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
132 Zplus = [z for z in Z if z >= -Zeps]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
133 Zminus = [z for z in Z if Zeps >= z]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
134
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
135 n_plus = len(Zplus)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
136 n_minus = len(Zminus)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
137
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
138
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
139 def h_kern(i, j):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
140 """Kernel function h for the medcouple, closing over the values of
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
141 Zplus and Zminus just defined above.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
142
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
143 In case a and be are within epsilon of the median, the kernel
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
144 is the signum of their position.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
145
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
146 """
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
147 a = Zplus[i]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
148 b = Zminus[j]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
149
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
150 if abs(a - b) <= 2*eps2:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
151 h = signum(n_plus - 1 - i - j)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
152 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
153 h = (a + b)/(a - b)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
154
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
155 return h
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
156
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
157 for i, zp in enumerate(Zplus):
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
158 for j, zm in enumerate(Zminus):
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
159 Hfile.write("%s " % h_kern(i, j)),
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
160 Hfile.write("\n")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
161
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
162 # Init left and right borders
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
163 L = [0]*n_plus
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
164 R = [n_minus - 1]*n_plus
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
165
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
166 Ltot = 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
167 Rtot = n_minus*n_plus
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
168 medc_idx = Rtot//2
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
169
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
170 Lfile.write("%s\n" % L)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
171 Rfile.write("%s\n" % R)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
172 RemFile.write("%s\n" % Rtot)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
173
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
174 # kth pair algorithm (Johnson & Mizoguchi)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
175 while Rtot - Ltot > n_plus:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
176
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
177 # First, compute the median inside the given bounds
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
178 # (Be stingy, reuse same generator)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
179 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i])
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
180
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
181 A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
182 W = [R[i] - L[i] + 1 for i in I2]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
183 Am = wmedian(A, W)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
184
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
185 Am_eps = eps1*(eps1 + abs(Am))
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
186
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
187 Amfile.write("%s\n" % Am)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
188
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
189 # Compute new left and right boundaries, based on the weighted
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
190 # median
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
191 P = []
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
192 Q = []
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
193
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
194 j = 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
195 for i in xrange(n_plus - 1, -1, -1):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
196 while j < n_minus and h_kern(i, j) - Am > Am_eps:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
197 j += 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
198 P.append(j - 1)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
199
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
200 P.reverse()
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
201
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
202 j = n_minus - 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
203 for i in xrange(0, n_plus):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
204 while j >= 0 and h_kern(i, j) - Am < -Am_eps:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
205 j -= 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
206 Q.append(j + 1)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
207
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
208 # Check on which side of those bounds the desired median of
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
209 # the whole matrix may be.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
210 sumP = sum(P) + len(P)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
211 sumQ = sum(Q)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
212
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
213 Ptotfile.write("%s\n" % sumP)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
214 Qtotfile.write("%s\n" % sumQ)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
215 Pfile.write("%s\n" % P)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
216 Qfile.write("%s\n" % Q)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
217
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
218 if medc_idx <= sumP - 1:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
219 R = P
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
220 Rtot = sumP
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
221 Lfile.write("%s\n" % L)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
222 Rfile.write("%s\n" % R)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
223 RemFile.write("%s\n" % (Rtot - Ltot))
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
224 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
225 if medc_idx > sumQ - 1:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
226 L = Q
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
227 Ltot = sumQ
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
228 Lfile.write("%s\n" % L)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
229 Rfile.write("%s\n" % R)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
230 RemFile.write("%s\n" % (Rtot - Ltot))
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
231 else:
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
232 jackpot.write("1")
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
233 return Am
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
234
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
235 # Didn't find the median, but now we have a very small search
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
236 # space to find it in, just between the left and right boundaries.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
237 # This space is of size Rtot - Ltot which is <= n_plus
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
238 jackpot.write("0")
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
239 A = []
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
240 for (i, (l, r)) in enumerate(izip(L, R)):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
241 for j in xrange(l, r + 1):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
242 A.append(h_kern(i, j))
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
243
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
244 A.sort() # A partial sort would suffice here
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
245 A.reverse()
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
246
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
247 Am = A[medc_idx - Ltot]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
248
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
249 return Am
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
250
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
251
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
252 def signum(x):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
253 if x > 0:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
254 return 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
255 if x < 0:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
256 return -1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
257 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
258 return 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
259
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
260
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
261 def main():
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
262 data = np.random.uniform(size=18)
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
263
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
264 mc = medcouple_1d(data)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
265 with open("mc", "w") as f:
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
266 f.write("%s\n" % mc)
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
267
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
268 if __name__ == "__main__":
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
269 main()