Mercurial > hg > medcouple

comparison jmedcouple.c++ @ 30:2e277129e81b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
jmedcouple: whitespace fixes
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Sun, 18 Jan 2015 20:42:25 -0500
parents cf8c8c825c6d
children 6c264f2053bc
comparison
equal deleted inserted replaced
29:cf8c8c825c6d 30:2e277129e81b
7 #include <algorithm> 7 #include <algorithm>
8 #include <utility> 8 #include <utility>
9 9
10 using namespace std; 10 using namespace std;
11 11
12 template <typename T> 12 template <typename T>
13 int signum(T val) { 13 int signum(T val) {
14 return (T(0) < val) - (val < T(0)); 14 return (T(0) < val) - (val < T(0));
15 } 15 }
16 16
17 template <typename Container, typename Out = typename Container::value_type> 17 template <typename Container, typename Out = typename Container::value_type>
18 Out sum(const Container& C) { 18 Out sum(const Container& C) {
19 return accumulate(C.begin(), C.end(), Out()); 19 return accumulate(C.begin(), C.end(), Out());
20 } 20 }
21 21
22 /* 22 /*
23 This computes the weighted median of array A with corresponding 23 This computes the weighted median of array A with corresponding
24 weights W. 24 weights W.
25 */ 25 */
58 else if (2*wright < wtot) 58 else if (2*wright < wtot)
59 beg = mid; 59 beg = mid;
60 else 60 else
61 return trial; 61 return trial;
62 } 62 }
63 63
64 return 0; 64 return 0;
65 } 65 }
66 66
67 double medcouple(const std::vector<double>& X, 67 double medcouple(const std::vector<double>& X,
68 double eps1 = numeric_limits<double>::epsilon(), 68 double eps1 = numeric_limits<double>::epsilon(),
69 double eps2 = numeric_limits<double>::min()) 69 double eps2 = numeric_limits<double>::min())
70 { 70 {
71 long n = X.size(), n2 = (n - 1)/2; 71 long n = X.size(), n2 = (n - 1)/2;
72 72
73 if (n < 3) 73 if (n < 3)
74 return 0; 74 return 0;
75 75
76 vector<double> Z = X; 76 vector<double> Z = X;
77 sort(Z.begin(), Z.end(), std::greater<double>()); 77 sort(Z.begin(), Z.end(), std::greater<double>());
78 78
79 double Zmed; 79 double Zmed;
80 if (n % 2 == 1) 80 if (n % 2 == 1)
81 Zmed = Z[n2]; 81 Zmed = Z[n2];
82 else 82 else
83 Zmed = (Z[n2] + Z[n2 + 1])/2; 83 Zmed = (Z[n2] + Z[n2 + 1])/2;
84 84
85 // Check if the median is at the edges up to relative epsilon 85 // Check if the median is at the edges up to relative epsilon
86 if (abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed))) 86 if (abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)))
87 return -1.0; 87 return -1.0;
88 if (abs(Z[n - 1] - Zmed) < eps1*(eps1 + abs(Zmed))) 88 if (abs(Z[n - 1] - Zmed) < eps1*(eps1 + abs(Zmed)))
89 return 1.0; 89 return 1.0;
90 90
91 // Centre Z wrt median, so that median(Z) = 0. 91 // Centre Z wrt median, so that median(Z) = 0.
92 for_each(Z.begin(), Z.end(), [&](double& z){ z -= Zmed;}); 92 for_each(Z.begin(), Z.end(), [&](double& z){ z -= Zmed;});
93 93
94 // Scale inside [-0.5, 0.5], for greater numerical stability. 94 // Scale inside [-0.5, 0.5], for greater numerical stability.
95 double Zden = 2*max(Z[0], -Z[n - 1]); 95 double Zden = 2*max(Z[0], -Z[n - 1]);
98 98
99 double Zeps = eps1*(eps1 + abs(Zmed)); 99 double Zeps = eps1*(eps1 + abs(Zmed));
100 100
101 // These overlap on the entries that are tied with the median 101 // These overlap on the entries that are tied with the median
102 vector<double> Zplus, Zminus; 102 vector<double> Zplus, Zminus;
103 copy_if(Z.begin(), Z.end(), back_inserter(Zplus), 103 copy_if(Z.begin(), Z.end(), back_inserter(Zplus),
104 [=](double z){return z >= -Zeps;}); 104 [=](double z){return z >= -Zeps;});
105 copy_if(Z.begin(), Z.end(), back_inserter(Zminus), 105 copy_if(Z.begin(), Z.end(), back_inserter(Zminus),
106 [=](double z){return Zeps >= z;}); 106 [=](double z){return Zeps >= z;});
107 107
108 long 108 long
109 n_plus = Zplus.size(), 109 n_plus = Zplus.size(),
110 n_minus = Zminus.size(); 110 n_minus = Zminus.size();
111 111
112 /* 112 /*
113 Kernel function h for the medcouple, closing over the values of 113 Kernel function h for the medcouple, closing over the values of
114 Zplus and Zminus just defined above. 114 Zplus and Zminus just defined above.
126 else 126 else
127 h = (a + b)/(a - b); 127 h = (a + b)/(a - b);
128 128
129 return h; 129 return h;
130 }; 130 };
131 131
132 // Init left and right borders 132 // Init left and right borders
133 vector<long> L(n_plus, 0); 133 vector<long> L(n_plus, 0);
134 vector<long> R(n_plus, n_minus - 1); 134 vector<long> R(n_plus, n_minus - 1);
135 135
136 long Ltot = 0; 136 long Ltot = 0;
151 } 151 }
152 double Am = wmedian(A, W); 152 double Am = wmedian(A, W);
153 153
154 double Am_eps = eps1*(eps1 + abs(Am)); 154 double Am_eps = eps1*(eps1 + abs(Am));
155 155
156 //Compute new left and right boundaries, based on the weighted 156 // Compute new left and right boundaries, based on the weighted
157 //median 157 // median
158 vector<long> P(n_plus), Q(n_plus); 158 vector<long> P(n_plus), Q(n_plus);
159 159
160 { 160 {
161 long j = 0; 161 long j = 0;
162 for (long i = n_plus - 1; i >= 0; i--) { 162 for (long i = n_plus - 1; i >= 0; i--) {
176 } 176 }
177 177
178 long 178 long
179 sumP = sum(P) + n_plus, 179 sumP = sum(P) + n_plus,
180 sumQ = sum(Q); 180 sumQ = sum(Q);
181 181
182 if (medc_idx <= sumP - 1) { 182 if (medc_idx <= sumP - 1) {
183 R = P; 183 R = P;
184 Rtot = sumP; 184 Rtot = sumP;
185 } 185 }
186 else { 186 else {
207 double Am = A[medc_idx - Ltot]; 207 double Am = A[medc_idx - Ltot];
208 return Am; 208 return Am;
209 } 209 }
210 210
211 211
212
213 int main(int argc, char** argv) { 212 int main(int argc, char** argv) {
214 double 213 double
215 eps1 = numeric_limits<double>::epsilon(), 214 eps1 = numeric_limits<double>::epsilon(),
216 eps2 = numeric_limits<double>::min(); 215 eps2 = numeric_limits<double>::min();
217 216
239 cout << setprecision(16) 238 cout << setprecision(16)
240 << medcouple(data, eps1, eps2) << endl; 239 << medcouple(data, eps1, eps2) << endl;
241 240
242 return 0; 241 return 0;
243 } 242 }
244