Mercurial > hg > medcouple

annotate wgt_himed_templ.h @ 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 e374cdd5e1e7
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
1 /*------ Definition of a template for whimed(_i) : *
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
2 * -------- ~~~~~~
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3 * i.e., included several times from ./wgt_himed.c
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4 * ~~~~~~~~~~~~~
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8 #ifdef _d_whimed_
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10 # define _WHIMED_ whimed
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11 # define _WGT_TYPE_ double
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12 # define _WGT_SUM_TYPE_ double
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 # undef _d_whimed_
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15 #elif defined (_i_whimed_)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 # define _WHIMED_ whimed_i
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18 # define _WGT_TYPE_ int
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19 # define _WGT_SUM_TYPE_ int64_t
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 # undef _i_whimed_
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22 #else
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23 # error "must define correct whimed_ macro !"
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
24 #endif
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26 #include <algorithm>
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 double _WHIMED_(double *a, _WGT_TYPE_ *w, int n,
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30 double* a_cand, double *a_srt, _WGT_TYPE_* w_cand)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31 {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 /*
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 Algorithm to compute the weighted high median in O(n) time.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 The whimed is defined as the smallest a[j] such that the sum
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37 of the weights of all a[i] <= a[j] is strictly greater than
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 half of the total weight.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40 Arguments:
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42 a: double array containing the observations
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43 n: number of observations
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 w: array of (int/double) weights of the observations.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
47 int n2, i, kcand;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48 /* sum of weights: `int' do overflow when n ~>= 1e5 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 _WGT_SUM_TYPE_ wleft, wmid, wright, w_tot, wrest;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51 double trial;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
53 w_tot = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 for (i = 0; i < n; ++i)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55 w_tot += w[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56 wrest = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
58 /* REPEAT : */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
59 do {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60 for (i = 0; i < n; ++i)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61 a_srt[i] = a[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62 n2 = n/2;/* =^= n/2 +1 with 0-indexing */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63 std::nth_element(a_srt, a_srt + n2, a_srt + n);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64 trial = a_srt[n2];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66 wleft = 0; wmid = 0; wright= 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67 for (i = 0; i < n; ++i) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68 if (a[i] < trial)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
69 wleft += w[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
70 else if (a[i] > trial)
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
71 wright += w[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
72 else
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
73 wmid += w[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 /* wleft = sum_{i; a[i] < trial} w[i]
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
76 * wmid = sum_{i; a[i] == trial} w[i] at least one 'i' since trial is one a[]!
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77 * wright= sum_{i; a[i] > trial} w[i]
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 kcand = 0;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80 if (2 * (wrest + wleft) > w_tot) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81 for (i = 0; i < n; ++i) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82 if (a[i] < trial) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83 a_cand[kcand] = a[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84 w_cand[kcand] = w[i]; ++kcand;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
86 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
87 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
88 else if (2 * (wrest + wleft + wmid) <= w_tot) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
89 for (i = 0; i < n; ++i) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
90 if (a[i] > trial) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
91 a_cand[kcand] = a[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
92 w_cand[kcand] = w[i]; ++kcand;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
93 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
94 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
95 wrest += wleft + wmid;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
96 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
97 else {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
98 return trial;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
99 /*==========*/
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101 n = kcand;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102 for (i = 0; i < n; ++i) {
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 a[i] = a_cand[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 w[i] = w_cand[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 } while(1);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 } /* _WHIMED_ */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110 #undef _WHIMED_
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 #undef _WGT_TYPE_
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 #undef _WGT_SUM_TYPE_