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