Mercurial > hg > medcouple

annotate jmedcouple.c++ @ 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 4c8b4c1e851e
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
37
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
1 // jmedcouple.c++ ---
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
2
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
3 // Copyright © 2015 Jordi Gutiérrez Hermoso <jordigh@octave.org>
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
4
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
5 // Author: Jordi Gutiérrez Hermoso
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
6
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
7 // This program is free software; you can redistribute it and/or
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
8 // modify it under the terms of the GNU General Public License
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
9 // as published by the Free Software Foundation; either version 3
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
10 // of the License, or (at your option) any later version.
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
11
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
12 // This program is distributed in the hope that it will be useful,
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
15 // GNU General Public License for more details.
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
16
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
17 // You should have received a copy of the GNU General Public License
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
18 // along with this program. If not, see <http://www.gnu.org/licenses/>.
4fb3b87b8610 Add GPL'ed header
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 35
diff changeset
19
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 #include <vector>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 #include <limits>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22 #include <iostream>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23 #include <string>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
24 #include <fstream>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25 #include <iomanip>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26 #include <algorithm>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27 #include <utility>
39
4c8b4c1e851e jmedcouple: add missing include
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 37
diff changeset
28 #include <tuple>
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30 using namespace std;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
32 template <typename T>
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 int signum(T val) {
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
34 return (T(0) < val) - (val < T(0));
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37 template <typename Container, typename Out = typename Container::value_type>
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 Out sum(const Container& C) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 return accumulate(C.begin(), C.end(), Out());
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
40 }
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42 /*
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43 This computes the weighted median of array A with corresponding
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 weights W.
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45 */
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46 double wmedian(const vector<double>& A,
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
47 const vector<long>& W)
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48 {
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
49 typedef pair<double, long> aw_t;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
51 long n = A.size();
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52 vector<aw_t> AW(n);
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
53 for (long i = 0; i < n; i ++)
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 AW[i] = make_pair(A[i], W[i]);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56 long wtot = sum(W);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
58 long beg = 0;
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
59 long end = n - 1;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61 while (true) {
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
62 long mid = (beg + end)/2;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63 nth_element(AW.begin(), AW.begin() + mid, AW.end(),
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64 [](const aw_t& l, const aw_t& r){return l.first > r.first;});
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65 double trial = AW[mid].first;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66 long wleft = 0, wright = 0;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68 for (aw_t& aw: AW ) {
35
6c264f2053bc jmedcouple: use std::tie and std::move for tuple unpacking
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 30
diff changeset
69 double a;
6c264f2053bc jmedcouple: use std::tie and std::move for tuple unpacking
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 30
diff changeset
70 long w;
6c264f2053bc jmedcouple: use std::tie and std::move for tuple unpacking
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 30
diff changeset
71 tie(a, w) = move(aw);
6c264f2053bc jmedcouple: use std::tie and std::move for tuple unpacking
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 30
diff changeset
72 if (a > trial)
6c264f2053bc jmedcouple: use std::tie and std::move for tuple unpacking
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 30
diff changeset
73 wleft += w;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 else
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 // This also includes a == trial, i.e. the "middle" weight.
35
6c264f2053bc jmedcouple: use std::tie and std::move for tuple unpacking
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 30
diff changeset
76 wright += w;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 if (2*wleft > wtot)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80 end = mid;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81 else if (2*wright < wtot)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82 beg = mid;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83 else
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84 return trial;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 }
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
86
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
87 return 0;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
88 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
89
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
90 double medcouple(const std::vector<double>& X,
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
91 double eps1 = numeric_limits<double>::epsilon(),
28
65f6fad65f49 jmedcouple: get rid of unused trace_lev parameter
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 27
diff changeset
92 double eps2 = numeric_limits<double>::min())
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
93 {
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
94 long n = X.size(), n2 = (n - 1)/2;
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
95
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
96 if (n < 3)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
97 return 0;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
98
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
99 vector<double> Z = X;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 sort(Z.begin(), Z.end(), std::greater<double>());
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102 double Zmed;
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
103 if (n % 2 == 1)
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 Zmed = Z[n2];
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105 else
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 Zmed = (Z[n2] + Z[n2 + 1])/2;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 // Check if the median is at the edges up to relative epsilon
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109 if (abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)))
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110 return -1.0;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 if (abs(Z[n - 1] - Zmed) < eps1*(eps1 + abs(Zmed)))
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 return 1.0;
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
113
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114 // Centre Z wrt median, so that median(Z) = 0.
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
115 for_each(Z.begin(), Z.end(), [&](double& z){ z -= Zmed;});
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
117 // Scale inside [-0.5, 0.5], for greater numerical stability.
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
118 double Zden = 2*max(Z[0], -Z[n - 1]);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 for_each(Z.begin(), Z.end(), [&](double& z){z /= Zden;});
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
120 Zmed /= Zden;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122 double Zeps = eps1*(eps1 + abs(Zmed));
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
123
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
124 // These overlap on the entries that are tied with the median
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 vector<double> Zplus, Zminus;
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
126 copy_if(Z.begin(), Z.end(), back_inserter(Zplus),
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
127 [=](double z){return z >= -Zeps;});
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
128 copy_if(Z.begin(), Z.end(), back_inserter(Zminus),
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129 [=](double z){return Zeps >= z;});
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
131 long
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
132 n_plus = Zplus.size(),
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
133 n_minus = Zminus.size();
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
134
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
135 /*
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
136 Kernel function h for the medcouple, closing over the values of
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
137 Zplus and Zminus just defined above.
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
138
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
139 In case a and be are within epsilon of the median, the kernel
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
140 is the signum of their position.
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
141 */
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
142 auto h_kern = [&](long i, long j) {
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
143 double a = Zplus[i];
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
144 double b = Zminus[j];
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
145
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
146 double h;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
147 if (abs(a - b) <= 2*eps2)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
148 h = signum(n_plus - 1 - i - j);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
149 else
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
150 h = (a + b)/(a - b);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
151
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
152 return h;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
153 };
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
154
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
155 // Init left and right borders
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
156 vector<long> L(n_plus, 0);
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
157 vector<long> R(n_plus, n_minus - 1);
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
158
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
159 long Ltot = 0;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
160 long Rtot = n_minus*n_plus;
27
40cbfb64e952 jmedcouple: change type of medc_idx to long
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 25
diff changeset
161 long medc_idx = Rtot/2;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
162
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
163 // kth pair algorithm (Johnson & Mizoguchi)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
164 while (Rtot - Ltot > n_plus) {
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
165
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
166 // First, compute the median inside the given bounds
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
167 vector<double> A;
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
168 vector<long> W;
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
169 for (long i = 0; i < n_plus; i++) {
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
170 if (L[i] <= R[i]) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
171 A.push_back(h_kern(i, (L[i] + R[i])/2));
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
172 W.push_back(R[i] - L[i] + 1);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
173 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
174 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
175 double Am = wmedian(A, W);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
176
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
177 double Am_eps = eps1*(eps1 + abs(Am));
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
178
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
179 // Compute new left and right boundaries, based on the weighted
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
180 // median
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
181 vector<long> P(n_plus), Q(n_plus);
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
182
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
183 {
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
184 long j = 0;
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
185 for (long i = n_plus - 1; i >= 0; i--) {
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
186 while (j < n_minus and h_kern(i, j) - Am > Am_eps)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
187 j++;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
188 P[i] = j - 1;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
189 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
190 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
191
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
192 {
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
193 long j = n_minus - 1;
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
194 for (long i = 0; i < n_plus; i++) {
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
195 while (j >= 0 and h_kern(i, j) - Am < -Am_eps)
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
196 j--;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
197 Q[i] = j + 1;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
198 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
199 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
200
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
201 long
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
202 sumP = sum(P) + n_plus,
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
203 sumQ = sum(Q);
30
2e277129e81b jmedcouple: whitespace fixes
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 29
diff changeset
204
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
205 if (medc_idx <= sumP - 1) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
206 R = P;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
207 Rtot = sumP;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
208 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
209 else {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
210 if (medc_idx > sumQ - 1) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
211 L = Q;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
212 Ltot = sumQ;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
213 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
214 else
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
215 return Am;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
216 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
217 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
218
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
219 // Didn't find the median, but now we have a very small search space
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
220 // to find it in, just between the left and right boundaries. This
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
221 // space is of size Rtot - Ltot which is <= n_plus
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
222 vector<double> A;
29
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
223 for (long i = 0; i < n_plus; i++) {
cf8c8c825c6d jmedcouple: change all ints to long, to avoid overflow
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 28
diff changeset
224 for (long j = L[i]; j <= R[i]; j++)
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
225 A.push_back(h_kern(i, j));
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
226 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
227 nth_element(A.begin(), A.begin() + medc_idx - Ltot, A.end(),
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
228 [](double x, double y) {return x > y;});
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
229
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
230 double Am = A[medc_idx - Ltot];
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
231 return Am;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
232 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
233
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
234
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
235 int main(int argc, char** argv) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
236 double
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
237 eps1 = numeric_limits<double>::epsilon(),
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
238 eps2 = numeric_limits<double>::min();
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
239
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
240 string fname;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
241 if (argc > 1) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
242 fname = argv[1];
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
243 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
244 else {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
245 cerr << "Please provide input file." << endl;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
246 return 1;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
247 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
248
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
249 ifstream ifs(fname);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
250 if (not ifs) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
251 cerr << "File " << fname << " not found." << endl;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
252 return 1;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
253 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
254
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
255 vector<double> data;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
256 double datum;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
257 while (ifs >> datum) {
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
258 data.push_back(datum);
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
259 }
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
260
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
261 cout << setprecision(16)
28
65f6fad65f49 jmedcouple: get rid of unused trace_lev parameter
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 27
diff changeset
262 << medcouple(data, eps1, eps2) << endl;
25
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
263
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
264 return 0;
518ef922033e Add C++ translation of Python code
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
265 }