1 /***********************************************************************
2 TEMPLATE:
3 double weightedMean_<int|dbl>(ARGUMENTS_LIST)
4
5 ARGUMENTS_LIST:
6 X_C_TYPE *x, R_xlen_t nx, double *w, R_xlen_t *idxs, R_xlen_t nidxs, int narm, int refine
7
8 Copyright: Henrik Bengtsson, 2014
9 ***********************************************************************/
10 #include <R_ext/Constants.h>
11 #include "000.types.h"
12
13 /* Expand arguments:
14 X_TYPE => (X_C_TYPE, X_IN_C)
15 */
16 #include "000.templates-types.h"
17 #include <R_ext/Error.h>
18
19
CONCAT_MACROS(weightedMean,X_C_SIGNATURE)20 double CONCAT_MACROS(weightedMean, X_C_SIGNATURE)(X_C_TYPE *x, R_xlen_t nx, double *w, R_xlen_t *idxs, R_xlen_t nidxs, int narm, int refine) {
21 X_C_TYPE value;
22 double weight;
23 R_xlen_t i;
24 LDOUBLE sum = 0, wtotal = 0;
25 LDOUBLE avg = R_NaN;
26
27 for (i=0; i < nidxs; i++) {
28 weight = R_INDEX_GET(w, ((idxs == NULL) ? (i) : idxs[i]), NA_REAL);
29
30 /* Skip or early stopping? */
31 if (weight == 0) {
32 continue;
33 }
34
35 value = R_INDEX_GET(x, ((idxs == NULL) ? (i) : idxs[i]), X_NA);
36 #if X_TYPE == 'i'
37 if (X_ISNAN(value)) {
38 /* Skip or early stopping? */
39 if (narm) {
40 continue;
41 } else {
42 sum = R_NaReal;
43 break;
44 }
45 } else {
46 sum += (LDOUBLE)weight * (LDOUBLE)value;
47 wtotal += weight;
48 }
49 #elif X_TYPE == 'r'
50 if (!narm) {
51 sum += (LDOUBLE)weight * (LDOUBLE)value;
52 wtotal += weight;
53 /* Early stopping? Special for long LDOUBLE vectors */
54 if (i % 1048576 == 0 && ISNAN(sum)) break;
55 } else if (!X_ISNAN(value)) {
56 sum += (LDOUBLE)weight * (LDOUBLE)value;
57 wtotal += weight;
58 }
59 #endif
60 } /* for (i ...) */
61
62 if (wtotal > DOUBLE_XMAX || wtotal < -DOUBLE_XMAX) {
63 avg = R_NaN;
64 } else if (sum > DOUBLE_XMAX) {
65 avg = R_PosInf;
66 } else if (sum < -DOUBLE_XMAX) {
67 avg = R_NegInf;
68 } else {
69 avg = sum / wtotal;
70
71 #if X_TYPE == 'r'
72 /* Extra precision by summing over residuals? */
73 if (refine && R_FINITE(avg)) {
74 sum = 0;
75 for (i=0; i < nidxs; i++) {
76 weight = R_INDEX_GET(w, ((idxs == NULL) ? (i) : idxs[i]), NA_REAL);
77 /* Skip? */
78 if (weight == 0) {
79 continue;
80 }
81
82 value = R_INDEX_GET(x, ((idxs == NULL) ? (i) : idxs[i]), X_NA);
83 if (!narm) {
84 sum += (LDOUBLE)weight * (value - avg);
85 /* Early stopping? Special for long LDOUBLE vectors */
86 if (i % 1048576 == 0 && ISNAN(sum)) break;
87 } else if (!ISNAN(value)) {
88 sum += (LDOUBLE)weight * (value - avg);
89 }
90 }
91 avg += (sum / wtotal);
92 }
93 #endif
94 }
95
96 return (double)avg;
97 }
98
99
100 /***************************************************************************
101 HISTORY:
102 2015-06-07 [DJ]
103 o Supported subsetted computation.
104 2014-12-08 [HB]
105 o Created.
106 **************************************************************************/
107