Mercurial > hg > medcouple

annotate medcouple.d @ 74:305b7361a5bd default tip @

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
showalgo: save a snapshot instead of waiting for keyboard input
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Sun, 29 May 2016 19:05:01 -0400
parents 18e7cd6ff057
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
38
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
1 // medcouple.d ---
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
2
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3 // Copyright © 2016 Jordi Gutiérrez Hermoso <jordigh@octave.org>
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5 // Author: Jordi Gutiérrez Hermoso
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7 // This program is free software; you can redistribute it and/or
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8 // modify it under the terms of the GNU General Public License
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9 // as published by the Free Software Foundation; either version 3
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10 // of the License, or (at your option) any later version.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12 // This program is distributed in the hope that it will be useful,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15 // GNU General Public License for more details.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 // You should have received a copy of the GNU General Public License
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18 // along with this program. If not, see <http://www.gnu.org/licenses/>.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 import std.stdio;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 import std.array;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22 import std.conv: to;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23 import std.range: iota, zip;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
24 import std.string: strip;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25 import std.typecons: Tuple;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26 import std.algorithm: filter, map, max, reduce, sort, topN;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27 import std.math: abs;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 // Could also just be reduce!"a+b", they call this "string lambda"
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30 alias sum = reduce!((a,b) => a + b);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32 int signum(T)(T val) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 return (cast(T)(0) < val ) - (val < cast(T)(0));
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 double wmedian(double[] A, long[] W) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37 auto AW = zip(A,W).array;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 long n = AW.length;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 auto wtot = sum(W);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41 typeof(n)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42 beg = 0,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43 end = n - 1;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45 while(true) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46 auto mid = (beg+end)/2;
51
18e7cd6ff057 medcouple.d: narrow computation of topN
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 38
diff changeset
47 topN!((x, y) => x[0] > y[0])(AW[beg..end+1], mid);
38
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 auto trial = AW[mid][0];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51 typeof(n) wleft = 0, wright = 0;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
53 foreach(aw; AW) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 auto a = aw[0];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55 auto w = aw[1];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56 if (a > trial)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57 wleft += w;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
58 else
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
59 wright += w;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62 if (2*wleft > wtot)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63 end = mid;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64 else if (2*wright < wtot)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65 beg = mid;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66 else
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67 return trial;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
69 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
70
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
71 double medcouple(const double[] X,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
72 double eps1 = double.epsilon,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
73 double eps2 = double.min)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 long n = X.length, n2 = (n-1)/2;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
76 if (n < 3)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77 return 0;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 auto Z = X.dup;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80 sort!((a,b) => a > b)(Z);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82 auto Zmed = n % 2 == 1 ? Z[n2] : (Z[n2] + Z[n2 + 1])/2;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84 // Check if the median is at the edges up to relative epsilon
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 if (abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)))
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
86 return -1.0;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
87 if (abs(Z[n - 1] - Zmed) < eps1*(eps1 + abs(Zmed)))
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
88 return 1.0;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
89
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
90 // Centre Z wrt median, so that median(Z) = 0.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
91 Z[] -= Zmed;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
92
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
93 // Scale inside [-0.5, 0.5], for greater numerical stability.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
94 auto Zden = 2*max(Z[0], -Z[n - 1]);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
95 Z[] /= Zden;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
96 Zmed /= Zden;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
97
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
98
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
99 auto
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 Zeps = eps1*(eps1 + abs(Zmed)),
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102 // These overlap on the entries that are tied with the median
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 Zplus = filter!(z => z >= -Zeps)(Z).array,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 Zminus = filter!(z => Zeps >= z)(Z).array;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 long
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107 n_plus = Zplus.length,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 n_minus = Zminus.length;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 /*
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 Kernel function h for the medcouple, closing over the values of
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
113 Zplus and Zminus just defined above.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
115 In case a and be are within epsilon of the median, the kernel
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116 is the signum of their position.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
117 */
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
118 auto h_kern(long i, long j){
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 auto
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
120 a = Zplus[i],
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121 b = Zminus[j];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
123 if (abs(a - b) <= 2*eps2)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
124 return signum(n_plus - 1 - i - j);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 else
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
126 return (a+b)/(a-b);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
127 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
128
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130 auto
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
131 // Init left and right borders
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
132 L = replicate([0L], n_plus),
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
133 R = replicate([n_minus - 1], n_plus),
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
134
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
135 Ltot = 0L,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
136 Rtot = n_minus*n_plus,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
137 medc_idx = Rtot/2;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
138
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
139 // kth pair algorithm (Johnson & Mizoguchi)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
140 while (Rtot - Ltot > n_plus) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
141
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
142 // First, compute the weighted median of row medians inside the
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
143 // given bounds
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
144 auto
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
145 I = filter!(i => L[i] <= R[i])(iota(0, n_plus-1)),
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
146 A = map!(i => h_kern(i, (L[i] + R[i])/2))(I).array,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
147 W = map!(i => R[i] - L[i] + 1)(I).array,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
148
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
149 Am = wmedian(A, W),
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
150 Am_eps = eps1*(eps1 + abs(Am));
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
151
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
152 // Compute new left and right boundaries, based on the weighted
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
153 // median
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
154 auto P = new long[n_plus], Q = new long[n_plus];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
155 {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
156 auto j = 0L;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
157 for (auto i = n_plus - 1; i >= 0; i--) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
158 while (j < n_minus && h_kern(i, j) - Am > Am_eps)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
159 j++;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
160 P[i] = j - 1;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
161 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
162 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
163
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
164 {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
165 auto j = n_minus - 1;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
166 for (auto i = 0; i < n_plus; i++) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
167 while (j >= 0 && h_kern(i, j) - Am < -Am_eps)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
168 j--;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
169 Q[i] = j + 1;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
170 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
171 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
172
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
173 // Discard left or right entries, depending on which side the
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
174 // medcouple is.
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
175 auto
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
176 sumP = sum(P) + n_plus,
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
177 sumQ = sum(Q);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
178
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
179 if (medc_idx <= sumP - 1) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
180 R = P;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
181 Rtot = sumP;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
182 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
183 else {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
184 if (medc_idx > sumQ - 1) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
185 L = Q;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
186 Ltot = sumQ;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
187 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
188 else
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
189 return Am;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
190 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
191 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
192
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
193 // Didn't find the median, but now we have a very small search space
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
194 // to find it in, just between the left and right boundaries. This
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
195 // space is of size Rtot - Ltot which is <= n_plus
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
196 auto
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
197 A = new double[Rtot - Ltot],
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
198 idx = 0L;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
199
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
200 for (auto i = 0L; i < n_plus; i++) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
201 for (auto j = L[i]; j <= R[i]; j++){
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
202 A[idx] = h_kern(i, j);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
203 ++idx;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
204 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
205 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
206
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
207 topN!((x,y) => x > y)(A, medc_idx - Ltot);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
208 return A[medc_idx - Ltot];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
209
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
210 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
211
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
212 int main(string[] argv) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
213 if (argv.length < 2) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
214 stderr.writeln("Please provide input file.");
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
215 return 1;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
216 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
217
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
218 auto fname = argv[1];
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
219 double[] data;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
220
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
221 foreach(line; File(fname).byLine) {
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
222 if (! line.empty)
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
223 data ~= to!(double)(line.strip);
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
224 }
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
225
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
226 writefln("%0.16f", medcouple(data));
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
227
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
228 return 0;
45fb52ba3057 Add D implementation of the medcouple
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
229 }