1 /* included 2 x from ./rowMedians.c
2 ~~~~~~~~~~~~
3 ***********************************************************************
4 TEMPLATE:
5 SEXP rowMedians_<Integer|Real>(...)
6
7 GENERATES:
8 SEXP rowMedians_Integer(SEXP x, int nrow, int ncol, int narm, int hasna, int byrow)
9 SEXP rowMedians_Real(SEXP x, int nrow, int ncol, int narm, int hasna, int byrow)
10
11 Arguments:
12 The following macros ("arguments") should be defined for the
13 template to work as intended.
14
15 - METHOD: the name of the resulting function
16 - X_TYPE: 'i' or 'r'
17
18 Authors:
19 Adopted from rowQuantiles.c by R. Gentleman.
20 Template by Henrik Bengtsson.
21
22 Copyright: Martin Maechler 2014-2021; Henrik Bengtsson, 2007-2013.
23 ***********************************************************************/
24 #include <R.h>
25 #include <Rinternals.h>
26 #include <Rmath.h>
27
28 /* Expand arguments:
29 X_TYPE => (X_C_TYPE, X_IN_C, X_ISNAN, [METHOD_NAME])
30 */
31 #include "templates-types.h"
32
33 #if X_TYPE == 'i'
34 #define PSORT iPsort
35 #elif X_TYPE == 'r'
36 #define PSORT rPsort
37 #endif
38
METHOD_NAME(SEXP x,int nrow,int ncol,int narm,int hasna,int byrow)39 SEXP METHOD_NAME(SEXP x, int nrow, int ncol, int narm, int hasna, int byrow)
40 {
41 /* R allocate a double vector of length 'nrow'
42 Note that 'nrow' means 'ncol' if byrow=FALSE. */
43 SEXP ans = PROTECT(allocVector(REALSXP, nrow));
44
45 C_METHOD_NAME(X_IN_C(x), REAL(ans), nrow, ncol, narm, hasna, byrow);
46
47 UNPROTECT(1);
48 return(ans);
49 }
50
51 // "C-only" method [SEXP free, can be called from "pure C" :
C_METHOD_NAME(X_C_TYPE * x,double * res,int nrow,int ncol,int narm,int hasna,int byrow)52 void C_METHOD_NAME(X_C_TYPE *x, double *res,
53 int nrow, int ncol, int narm, int hasna, int byrow)
54 {
55 Rboolean isOdd;
56 int ii, jj, kk, qq;
57 int *colOffset;
58 X_C_TYPE value,
59 /* R allocate memory for the 'rowData'. This will be
60 taken care of by the R garbage collector later on. */
61 *rowData = (X_C_TYPE *) R_alloc(ncol, sizeof(X_C_TYPE));
62
63 if (!hasna) // If there are no missing values, don't try to remove them
64 narm = FALSE;
65
66 /* When narm is false, isOdd and qq are the same for all rows */
67 if (!narm) {
68 isOdd = (ncol % 2 == 1);
69 qq = (int)(ncol/2) - 1;
70 } else {
71 isOdd = FALSE;
72 qq = 0;
73 }
74
75 value = 0;
76
77 /* Pre-calculate the column offsets */
78 colOffset = (int *) R_alloc(ncol, sizeof(int));
79
80 // HJ begin
81 if (byrow) {
82 for(jj=0; jj < ncol; jj++)
83 colOffset[jj] = (int)jj*nrow;
84 } else {
85 for(jj=0; jj < ncol; jj++)
86 colOffset[jj] = jj;
87 }
88 // HJ end
89
90 if (hasna) {
91 for(ii=0; ii < nrow; ii++) {
92 if(ii % 1000 == 0)
93 R_CheckUserInterrupt();
94
95 int rowIdx = byrow ? ii : ncol*ii; //HJ
96
97 kk = 0; /* The index of the last non-NA value detected */
98 for(jj=0; jj < ncol; jj++) {
99 value = x[rowIdx+colOffset[jj]]; //HJ
100
101 if (X_ISNAN(value)) {
102 if (!narm) {
103 kk = -1;
104 break;
105 }
106 } else {
107 rowData[kk] = value;
108 kk = kk + 1;
109 }
110 }
111
112 if (kk == 0) {
113 res[ii] = R_NaN;
114 } else if (kk == -1) {
115 res[ii] = R_NaReal;
116 } else {
117 /* When narm is true, isOdd and qq may change with row */
118 if (narm) {
119 isOdd = (kk % 2 == 1);
120 qq = (int)(kk/2) - 1;
121 }
122
123 /* Permute x[0:kk-1] so that x[qq] is in the correct
124 place with smaller values to the left, ... */
125 PSORT(rowData, kk, qq+1);
126 value = rowData[qq+1];
127
128 if (isOdd) {
129 res[ii] = (double)value;
130 } else {
131 if (narm || !X_ISNAN(value)) {
132 /* Permute x[0:qq-2] so that x[qq-1] is in the correct
133 place with smaller values to the left, ... */
134 PSORT(rowData, qq+1, qq);
135 if (X_ISNAN(rowData[qq]))
136 res[ii] = R_NaReal;
137 else
138 res[ii] = ((double)(rowData[qq] + value))/2;
139 } else {
140 res[ii] = (double)value;
141 }
142 }
143 }
144 } // for(..)
145 } else { // no NAs
146 for(ii=0; ii < nrow; ii++) {
147 if(ii % 1000 == 0)
148 R_CheckUserInterrupt();
149
150 int rowIdx = byrow ? ii : ncol*ii; //HJ
151
152 for(jj=0; jj < ncol; jj++)
153 rowData[jj] = x[rowIdx+colOffset[jj]]; //HJ
154
155 /* Permute x[0:ncol-1] so that x[qq] is in the correct
156 place with smaller values to the left, ... */
157 PSORT(rowData, ncol, qq+1);
158 value = rowData[qq+1];
159
160 if (isOdd) {
161 res[ii] = (double)value;
162 } else {
163 /* Permute x[0:qq-2] so that x[qq-1] is in the correct
164 place with smaller values to the left, ... */
165 PSORT(rowData, qq+1, qq);
166 res[ii] = (double)((rowData[qq] + value))/2;
167 }
168 } // for(..)
169 } /* if (hasna ...) */
170 }
171
172 /* Undo template macros */
173 #undef PSORT
174 #include "templates-types_undef.h"
175
176
177 /***************************************************************************
178 HISTORY:
179 2014-12-09 [MMaechler]
180 o do not use '== TRUE' '== FALSE' -- as we have no NA here
181 o resolve REAL(ans) outside for(ii ..)
182 o add "SEXP-free" C routines, others can call: C_rowMedians_(Real|Integer)
183
184 2013-04-23 [HB]
185 o BUG FIX: The integer template of rowMedians_<Integer|Real>() would
186 not handle ties properly. This was because ties were calculated as
187 '(double)((rowData[qq] + value)/2)' instead of
188 '((double)(rowData[qq] + value))/2'.
189 2013-01-13 [HB]
190 o Merged rowMedians_Integer() and rowMedians_Read() into template
191 rowMedians_<Integer|Real>().
192 2013-01-13 [HB]
193 o Using internal arguments 'by_row' instead of 'by_column'.
194 2011-12-11 [HB]
195 o BUG FIX: rowMediansReal(..., na.rm=TRUE) did not handle NaN:s, only NA:s.
196 Note that NaN:s does not exist for integers.
197 2011-10-12 [HJ]
198 o Added colMedians().
199 o Now rowMediansInteger/Real() can operate also by columns, cf. argument
200 'by_column'.
201 2007-08-14 [HB]
202 o Added checks for user interrupts every 1000 line.
203 o Added argument 'hasNA' to rowMedians().
204 2005-12-07 [HB]
205 o BUG FIX: When calculating the median of an even number (non-NA) values,
206 the length of the second sort was one element too short, which made the
207 method to freeze, i.e. rPsort(rowData, qq, qq) is now (...qq+1, qq).
208 2005-11-24 [HB]
209 o By implementing a special version for integers, there is no need to
210 coerce to double in R, which would take up twice the amount of memory.
211 o rowMedians() now handles NAs too.
212 o Adopted from rowQuantiles.c in Biobase of Bioconductor.
213 **************************************************************************/
214