Mercurial > hg > medcouple

annotate talk/code/show_mc.py @ 67:1aa38f4846e0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
showalgo: show the rank
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Wed, 18 May 2016 17:43:45 -0400
parents 916349600c4b
children 5510ac95a3d3
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")
67
1aa38f4846e0 showalgo: show the rank
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 66
diff changeset
98 medc_idxfile = open("medc_idx", "w")
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
99 jackpot = open("jackpot", "w")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
100
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101 # 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
102
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 n = len(X)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 n2 = (n - 1)//2
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 if n < 3:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107 return 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109 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
110
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 if n % 2 == 1:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 Zmed = Z[n2]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
113 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114 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
115
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116 #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
117 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
118 return -1.0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 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
120 return 1.0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122 # 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
123 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
124
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 # 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
126 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
127 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
128 Zmed /= Zden
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130 Zeps = eps1*(eps1 + abs(Zmed))
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
131
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
132 # 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
133 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
134 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
135
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
136 n_plus = len(Zplus)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
137 n_minus = len(Zminus)
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
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
140 def h_kern(i, j):
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
141 """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
142 Zplus and Zminus just defined above.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
143
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
144 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
145 is the signum of their position.
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 """
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
148 a = Zplus[i]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
149 b = Zminus[j]
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
150
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
151 if abs(a - b) <= 2*eps2:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
152 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
153 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
154 h = (a + b)/(a - b)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
155
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
156 return h
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
157
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
158 for i, zp in enumerate(Zplus):
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
159 for j, zm in enumerate(Zminus):
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
160 Hfile.write("%s " % h_kern(i, j)),
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
161 Hfile.write("\n")
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
162
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
163 # Init left and right borders
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
164 L = [0]*n_plus
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
165 R = [n_minus - 1]*n_plus
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
166
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
167 Ltot = 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
168 Rtot = n_minus*n_plus
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
169 medc_idx = Rtot//2
67
1aa38f4846e0 showalgo: show the rank
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 66
diff changeset
170 medc_idxfile.write("%s\n" % medc_idx)
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
171
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
172 Lfile.write("%s\n" % L)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
173 Rfile.write("%s\n" % R)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
174 RemFile.write("%s\n" % Rtot)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
175
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
176 # kth pair algorithm (Johnson & Mizoguchi)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
177 while Rtot - Ltot > n_plus:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
178
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
179 # 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
180 # (Be stingy, reuse same generator)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
181 [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
182
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
183 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
184 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
185 Am = wmedian(A, W)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
186
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
187 Am_eps = eps1*(eps1 + abs(Am))
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
188
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
189 Amfile.write("%s\n" % Am)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
190
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
191 # 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
192 # median
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
193 P = []
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
194 Q = []
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
195
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
196 j = 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
197 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
198 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
199 j += 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
200 P.append(j - 1)
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 P.reverse()
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
203
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
204 j = n_minus - 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
205 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
206 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
207 j -= 1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
208 Q.append(j + 1)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
209
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
210 # 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
211 # the whole matrix may be.
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
212 sumP = sum(P) + len(P)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
213 sumQ = sum(Q)
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
214
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
215 Ptotfile.write("%s\n" % sumP)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
216 Qtotfile.write("%s\n" % sumQ)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
217 Pfile.write("%s\n" % P)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
218 Qfile.write("%s\n" % Q)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
219
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
220 if medc_idx <= sumP - 1:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
221 R = P
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
222 Rtot = sumP
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
223 Lfile.write("%s\n" % L)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
224 Rfile.write("%s\n" % R)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
225 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
226 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
227 if medc_idx > sumQ - 1:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
228 L = Q
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
229 Ltot = sumQ
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
230 Lfile.write("%s\n" % L)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
231 Rfile.write("%s\n" % R)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
232 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
233 else:
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
234 jackpot.write("1")
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
235 return Am
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
236
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
237 # 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
238 # 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
239 # 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
240 jackpot.write("0")
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
241 A = []
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
242 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
243 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
244 A.append(h_kern(i, j))
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
245
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
246 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
247 A.reverse()
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 Am = A[medc_idx - Ltot]
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 return Am
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
252
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
253
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
254 def signum(x):
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 if x < 0:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
258 return -1
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
259 else:
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
260 return 0
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
261
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
262
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
263 def main():
67
1aa38f4846e0 showalgo: show the rank
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 66
diff changeset
264 data = np.random.uniform(size=26)
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
265
66
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
266 mc = medcouple_1d(data)
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
267 with open("mc", "w") as f:
916349600c4b algorithm-visualisation code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 65
diff changeset
268 f.write("%s\n" % mc)
65
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
269
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
270 if __name__ == "__main__":
364a0c9df823 show_mc.py: copy from medcouple.py
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
271 main()