Rcpp Version 1.0.9
rowSums.h
Go to the documentation of this file.
1 // -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
2 //
3 // rowSums.h: Rcpp R/C++ interface class library -- rowSums, colSums, rowMeans, colMeans
4 //
5 // Copyright (C) 2016 Nathan Russell
6 //
7 // This file is part of Rcpp.
8 //
9 // Rcpp is free software: you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License as published by
11 // the Free Software Foundation, either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // Rcpp is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 // GNU General Public License for more details.
18 //
19 // You should have received a copy of the GNU General Public License
20 // along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
21 
22 #ifndef Rcpp__sugar__rowSums_h
23 #define Rcpp__sugar__rowSums_h
24 
25 namespace Rcpp {
26 namespace sugar {
27 namespace detail {
28 
29 
30 inline bool check_na(double x) {
31  return ISNAN(x);
32 }
33 
34 inline bool check_na(int x) {
35  return x == NA_INTEGER;
36 }
37 
38 inline bool check_na(Rboolean x) {
39  return x == NA_LOGICAL;
40 }
41 
42 inline bool check_na(SEXP x) {
43  return x == NA_STRING;
44 }
45 
46 inline bool check_na(Rcomplex x) {
47  return ISNAN(x.r) || ISNAN(x.i);
48 }
49 
50 
51 inline void incr(double* lhs, double rhs) {
52  *lhs += rhs;
53 }
54 
55 inline void incr(int* lhs, int rhs) {
56  *lhs += rhs;
57 }
58 
59 inline void incr(Rcomplex* lhs, const Rcomplex& rhs) {
60  lhs->r += rhs.r;
61  lhs->i += rhs.i;
62 }
63 
64 
65 inline void div(double* lhs, R_xlen_t rhs) {
66  *lhs /= static_cast<double>(rhs);
67 }
68 
69 inline void div(Rcomplex* lhs, R_xlen_t rhs) {
70  lhs->r /= static_cast<double>(rhs);
71  lhs->i /= static_cast<double>(rhs);
72 }
73 
74 
75 inline void set_nan(double* x) {
76  *x = R_NaN;
77 }
78 
79 inline void set_nan(Rcomplex* x) {
80  x->r = R_NaN;
81  x->i = R_NaN;
82 }
83 
84 
85 template <int RTYPE>
86 struct RowSumsReturn {
87  typedef Vector<RTYPE> type;
88  enum { rtype = RTYPE };
89 };
90 
91 template <>
92 struct RowSumsReturn<LGLSXP> {
93  typedef Vector<INTSXP> type;
94  enum { rtype = INTSXP };
95 };
96 
97 template <int RTYPE>
98 struct ColSumsReturn
99  : public RowSumsReturn<RTYPE> {};
100 
101 
102 template <int RTYPE>
103 struct RowMeansReturn {
104  typedef Vector<REALSXP> type;
105  enum { rtype = REALSXP };
106 };
107 
108 template <>
109 struct RowMeansReturn<CPLXSXP> {
110  typedef Vector<CPLXSXP> type;
111  enum { rtype = CPLXSXP };
112 };
113 
114 template <int RTYPE>
115 struct ColMeansReturn
116  : public RowMeansReturn<RTYPE> {};
117 
118 
119 } // detail
120 
121 
122 // RowSums
123 // na.rm = FALSE
124 // default input
125 // default output
126 //
127 template <int RTYPE, bool NA, typename T, bool NA_RM = false>
128 class RowSumsImpl :
129  public Lazy<typename detail::RowSumsReturn<RTYPE>::type, RowSumsImpl<RTYPE, NA, T, NA_RM> > {
130 private:
131  const MatrixBase<RTYPE, NA, T>& ref;
132 
133  typedef detail::RowSumsReturn<RTYPE> return_traits;
134  typedef typename return_traits::type return_vector;
135  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
136 
137 public:
138  RowSumsImpl(const MatrixBase<RTYPE, NA, T>& ref_)
139  : ref(ref_)
140  {}
141 
142  return_vector get() const {
143  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
144  return_vector res(nr);
145 
146  for (j = 0; j < nc; j++) {
147  for (i = 0; i < nr; i++) {
148  detail::incr(&res[i], ref(i, j));
149  }
150  }
151 
152  return res;
153  }
154 };
155 
156 // RowSums
157 // na.rm = FALSE
158 // LGLSXP / INTSXP input
159 // INTSXP output
160 //
161 // int + NA_LOGICAL (NA_INTEGER) != NA_INTEGER, as is the
162 // case with NA_REAL, so we specialize for these two SEXPTYPES
163 // and do explicit accounting of NAs.
164 //
165 // The two specializations, while necessary, are redundant, hence
166 // the macro. The same applies to the 'na.rm = TRUE' variant, and
167 // likewise for colSums, rowMeans, and colMeans.
168 //
169 #define ROW_SUMS_IMPL_KEEPNA(__RTYPE__) \
170  \
171 template <bool NA, typename T, bool NA_RM> \
172 class RowSumsImpl<__RTYPE__, NA, T, NA_RM> : \
173  public Lazy<typename detail::RowSumsReturn<__RTYPE__>::type, RowSumsImpl<__RTYPE__, NA, T, NA_RM> > { \
174 private: \
175  const MatrixBase<__RTYPE__, NA, T>& ref; \
176  \
177  typedef detail::RowSumsReturn<__RTYPE__> return_traits; \
178  typedef typename return_traits::type return_vector; \
179  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
180  \
181  struct bit { \
182  unsigned char x : 1; \
183  }; \
184  \
185 public: \
186  RowSumsImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
187  : ref(ref_) \
188  {} \
189  \
190  return_vector get() const { \
191  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
192  return_vector res(nr); \
193  \
194  std::vector<bit> na_flags(nr); \
195  \
196  for (j = 0; j < nc; j++) { \
197  for (i = 0; i < nr; i++) { \
198  if (detail::check_na(ref(i, j))) { \
199  na_flags[i].x |= 0x1; \
200  } \
201  detail::incr(&res[i], ref(i, j)); \
202  } \
203  } \
204  \
205  for (i = 0; i < nr; i++) { \
206  if (na_flags[i].x) { \
207  res[i] = NA_INTEGER; \
208  } \
209  } \
210  \
211  return res; \
212  } \
213 };
214 
215 ROW_SUMS_IMPL_KEEPNA(LGLSXP)
216 ROW_SUMS_IMPL_KEEPNA(INTSXP)
217 
218 #undef ROW_SUMS_IMPL_KEEPNA
219 
220 // RowSums
221 // na.rm = TRUE
222 // default input
223 // default output
224 //
225 template <int RTYPE, bool NA, typename T>
226 class RowSumsImpl<RTYPE, NA, T, true> :
227  public Lazy<typename detail::RowSumsReturn<RTYPE>::type, RowSumsImpl<RTYPE, NA, T, true> > {
228 private:
229  const MatrixBase<RTYPE, NA, T>& ref;
230 
231  typedef detail::RowSumsReturn<RTYPE> return_traits;
232  typedef typename return_traits::type return_vector;
233  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
234 
235 public:
236  RowSumsImpl(const MatrixBase<RTYPE, NA, T>& ref_)
237  : ref(ref_)
238  {}
239 
240  return_vector get() const {
241  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
242  return_vector res(nr);
243 
244  stored_type current = stored_type();
245 
246  for (j = 0; j < nc; j++) {
247  for (i = 0; i < nr; i++) {
248  current = ref(i, j);
249  if (!detail::check_na(current)) {
250  detail::incr(&res[i], current);
251  }
252  }
253  }
254 
255  return res;
256  }
257 };
258 
259 // RowSums
260 // na.rm = TRUE
261 // LGLSXP / INTSXP input
262 // INTSXP output
263 //
264 #define ROW_SUMS_IMPL_RMNA(__RTYPE__) \
265  \
266 template <bool NA, typename T> \
267 class RowSumsImpl<__RTYPE__, NA, T, true> : \
268  public Lazy<typename detail::RowSumsReturn<__RTYPE__>::type, RowSumsImpl<__RTYPE__, NA, T, true> > { \
269 private: \
270  const MatrixBase<__RTYPE__, NA, T>& ref; \
271  \
272  typedef detail::RowSumsReturn<__RTYPE__> return_traits; \
273  typedef typename return_traits::type return_vector; \
274  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
275  \
276 public: \
277  RowSumsImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
278  : ref(ref_) \
279  {} \
280  \
281  return_vector get() const { \
282  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
283  return_vector res(nr); \
284  \
285  stored_type current = stored_type(); \
286  \
287  for (j = 0; j < nc; j++) { \
288  for (i = 0; i < nr; i++) { \
289  current = ref(i, j); \
290  if (!detail::check_na(current)) { \
291  detail::incr(&res[i], current); \
292  } \
293  } \
294  } \
295  \
296  return res; \
297  } \
298 };
299 
300 ROW_SUMS_IMPL_RMNA(LGLSXP)
301 ROW_SUMS_IMPL_RMNA(INTSXP)
302 
303 #undef ROW_SUMS_IMPL_RMNA
304 
305 // RowSums
306 // Input with template parameter NA = false
307 // RowSumsImpl<..., NA_RM = false>
308 //
309 template <int RTYPE, typename T, bool NA_RM>
310 class RowSumsImpl<RTYPE, false, T, NA_RM>
311  : public RowSumsImpl<RTYPE, false, T, false> {};
312 
313 
314 // ColSums
315 // na.rm = FALSE
316 // default input
317 // default output
318 //
319 template <int RTYPE, bool NA, typename T, bool NA_RM = false>
320 class ColSumsImpl :
321  public Lazy<typename detail::ColSumsReturn<RTYPE>::type, ColSumsImpl<RTYPE, NA, T, NA_RM> > {
322 private:
323  const MatrixBase<RTYPE, NA, T>& ref;
324 
325  typedef detail::ColSumsReturn<RTYPE> return_traits;
326  typedef typename return_traits::type return_vector;
327  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
328 
329 public:
330  ColSumsImpl(const MatrixBase<RTYPE, NA, T>& ref_)
331  : ref(ref_)
332  {}
333 
334  return_vector get() const {
335  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
336  return_vector res(nc);
337 
338  for (j = 0; j < nc; j++) {
339  for (i = 0; i < nr; i++) {
340  detail::incr(&res[j], ref(i, j));
341  }
342  }
343 
344  return res;
345  }
346 };
347 
348 // ColSums
349 // na.rm = FALSE
350 // LGLSXP / INTSXP input
351 // INTSXP output
352 //
353 #define COL_SUMS_IMPL_KEEPNA(__RTYPE__) \
354  \
355 template <bool NA, typename T, bool NA_RM> \
356 class ColSumsImpl<__RTYPE__, NA, T, NA_RM> : \
357  public Lazy<typename detail::ColSumsReturn<__RTYPE__>::type, ColSumsImpl<__RTYPE__, NA, T, NA_RM> > { \
358 private: \
359  const MatrixBase<__RTYPE__, NA, T>& ref; \
360  \
361  typedef detail::ColSumsReturn<__RTYPE__> return_traits; \
362  typedef typename return_traits::type return_vector; \
363  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
364  \
365  struct bit { \
366  unsigned char x : 1; \
367  }; \
368  \
369 public: \
370  ColSumsImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
371  : ref(ref_) \
372  {} \
373  \
374  return_vector get() const { \
375  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
376  return_vector res(nc); \
377  \
378  std::vector<bit> na_flags(nc); \
379  \
380  for (j = 0; j < nc; j++) { \
381  for (i = 0; i < nr; i++) { \
382  if (detail::check_na(ref(i, j))) { \
383  na_flags[j].x |= 0x1; \
384  } \
385  detail::incr(&res[j], ref(i, j)); \
386  } \
387  } \
388  \
389  for (j = 0; j < nc; j++) { \
390  if (na_flags[j].x) { \
391  res[j] = NA_INTEGER; \
392  } \
393  } \
394  \
395  return res; \
396  } \
397 };
398 
399 COL_SUMS_IMPL_KEEPNA(LGLSXP)
400 COL_SUMS_IMPL_KEEPNA(INTSXP)
401 
402 #undef COL_SUMS_IMPL_KEEPNA
403 
404 // ColSums
405 // na.rm = TRUE
406 // default input
407 // default output
408 //
409 template <int RTYPE, bool NA, typename T>
410 class ColSumsImpl<RTYPE, NA, T, true> :
411  public Lazy<typename detail::ColSumsReturn<RTYPE>::type, ColSumsImpl<RTYPE, NA, T, true> > {
412 private:
413  const MatrixBase<RTYPE, NA, T>& ref;
414 
415  typedef detail::ColSumsReturn<RTYPE> return_traits;
416  typedef typename return_traits::type return_vector;
417  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
418 
419 public:
420  ColSumsImpl(const MatrixBase<RTYPE, NA, T>& ref_)
421  : ref(ref_)
422  {}
423 
424  return_vector get() const {
425  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
426  return_vector res(nc);
427 
428  stored_type current = stored_type();
429 
430  for (j = 0; j < nc; j++) {
431  for (i = 0; i < nr; i++) {
432  current = ref(i, j);
433  if (!detail::check_na(current)) {
434  detail::incr(&res[j], current);
435  }
436  }
437  }
438 
439  return res;
440  }
441 };
442 
443 // ColSums
444 // na.rm = TRUE
445 // LGLSXP / INTSXP input
446 // INTSXP output
447 //
448 #define COL_SUMS_IMPL_RMNA(__RTYPE__) \
449  \
450 template <bool NA, typename T> \
451 class ColSumsImpl<__RTYPE__, NA, T, true> : \
452  public Lazy<typename detail::ColSumsReturn<__RTYPE__>::type, ColSumsImpl<__RTYPE__, NA, T, true> > { \
453 private: \
454  const MatrixBase<__RTYPE__, NA, T>& ref; \
455  \
456  typedef detail::ColSumsReturn<__RTYPE__> return_traits; \
457  typedef typename return_traits::type return_vector; \
458  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
459  \
460 public: \
461  ColSumsImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
462  : ref(ref_) \
463  {} \
464  \
465  return_vector get() const { \
466  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
467  return_vector res(nc); \
468  \
469  stored_type current = stored_type(); \
470  \
471  for (j = 0; j < nc; j++) { \
472  for (i = 0; i < nr; i++) { \
473  current = ref(i, j); \
474  if (!detail::check_na(current)) { \
475  detail::incr(&res[j], current); \
476  } \
477  } \
478  } \
479  \
480  return res; \
481  } \
482 };
483 
484 COL_SUMS_IMPL_RMNA(LGLSXP)
485 COL_SUMS_IMPL_RMNA(INTSXP)
486 
487 #undef COL_SUMS_IMPL_RMNA
488 
489 // ColSums
490 // Input with template parameter NA = false
491 // ColSumsImpl<..., NA_RM = false>
492 //
493 template <int RTYPE, typename T, bool NA_RM>
494 class ColSumsImpl<RTYPE, false, T, NA_RM>
495  : public ColSumsImpl<RTYPE, false, T, false> {};
496 
497 
498 // RowMeans
499 // na.rm = FALSE
500 // default input
501 // default output
502 //
503 // All RowMeans and ColMeans variants use a single-pass
504 // mean calculation as in array.c
505 //
506 template <int RTYPE, bool NA, typename T, bool NA_RM = false>
507 class RowMeansImpl :
508  public Lazy<typename detail::RowMeansReturn<RTYPE>::type, RowMeansImpl<RTYPE, NA, T, NA_RM> > {
509 private:
510  const MatrixBase<RTYPE, NA, T>& ref;
511 
512  typedef detail::RowMeansReturn<RTYPE> return_traits;
513  typedef typename return_traits::type return_vector;
514  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
515 
516 public:
517  RowMeansImpl(const MatrixBase<RTYPE, NA, T>& ref_)
518  : ref(ref_)
519  {}
520 
521  return_vector get() const {
522  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
523  return_vector res(nr);
524 
525  for (j = 0; j < nc; j++) {
526  for (i = 0; i < nr; i++) {
527  detail::incr(&res[i], ref(i, j));
528  }
529  }
530 
531  for (i = 0; i < nr; i++) {
532  detail::div(&res[i], nc);
533  }
534 
535  return res;
536  }
537 };
538 
539 // RowMeans
540 // na.rm = FALSE
541 // LGLSXP / INTSXP input
542 // REALSXP output
543 //
544 #define ROW_MEANS_IMPL_KEEPNA(__RTYPE__) \
545  \
546 template <bool NA, typename T, bool NA_RM> \
547 class RowMeansImpl<__RTYPE__, NA, T, NA_RM> : \
548  public Lazy<typename detail::RowMeansReturn<__RTYPE__>::type, RowMeansImpl<__RTYPE__, NA, T, NA_RM> > { \
549 private: \
550  const MatrixBase<__RTYPE__, NA, T>& ref; \
551  \
552  typedef detail::RowMeansReturn<__RTYPE__> return_traits; \
553  typedef typename return_traits::type return_vector; \
554  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
555  \
556  struct bit { \
557  unsigned char x : 1; \
558  }; \
559  \
560 public: \
561  RowMeansImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
562  : ref(ref_) \
563  {} \
564  \
565  return_vector get() const { \
566  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
567  return_vector res(nr); \
568  \
569  std::vector<bit> na_flags(nc); \
570  \
571  for (j = 0; j < nc; j++) { \
572  for (i = 0; i < nr; i++) { \
573  if (detail::check_na(ref(i, j))) { \
574  na_flags[i].x |= 0x1; \
575  } \
576  detail::incr(&res[i], ref(i, j)); \
577  } \
578  } \
579  \
580  for (i = 0; i < nr; i++) { \
581  if (!na_flags[i].x) { \
582  detail::div(&res[i], nc); \
583  } else { \
584  res[i] = NA_REAL; \
585  } \
586  } \
587  \
588  return res; \
589  } \
590 };
591 
592 ROW_MEANS_IMPL_KEEPNA(LGLSXP)
593 ROW_MEANS_IMPL_KEEPNA(INTSXP)
594 
595 #undef ROW_MEANS_IMPL_KEEPNA
596 
597 // RowMeans
598 // na.rm = TRUE
599 // default input
600 // default output
601 //
602 template <int RTYPE, bool NA, typename T>
603 class RowMeansImpl<RTYPE, NA, T, true> :
604  public Lazy<typename detail::RowMeansReturn<RTYPE>::type, RowMeansImpl<RTYPE, NA, T, true> > {
605 private:
606  const MatrixBase<RTYPE, NA, T>& ref;
607 
608  typedef detail::RowMeansReturn<RTYPE> return_traits;
609  typedef typename return_traits::type return_vector;
610  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
611 
612 public:
613  RowMeansImpl(const MatrixBase<RTYPE, NA, T>& ref_)
614  : ref(ref_)
615  {}
616 
617  return_vector get() const {
618  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
619  return_vector res(nr);
620 
621  std::vector<R_xlen_t> n_ok(nr, 0);
622  stored_type current = stored_type();
623 
624 
625  for (j = 0; j < nc; j++) {
626  for (i = 0; i < nr; i++) {
627  current = ref(i, j);
628  if (!detail::check_na(current)) {
629  detail::incr(&res[i], ref(i, j));
630  ++n_ok[i];
631  }
632  }
633  }
634 
635  for (i = 0; i < nr; i++) {
636  if (n_ok[i]) {
637  detail::div(&res[i], n_ok[i]);
638  } else {
639  detail::set_nan(&res[i]);
640  }
641  }
642 
643  return res;
644  }
645 };
646 
647 // RowMeans
648 // na.rm = TRUE
649 // LGLSXP / INTSXP input
650 // REALSXP output
651 //
652 #define ROW_MEANS_IMPL_RMNA(__RTYPE__) \
653  \
654 template <bool NA, typename T> \
655 class RowMeansImpl<__RTYPE__, NA, T, true> : \
656  public Lazy<typename detail::RowMeansReturn<__RTYPE__>::type, RowMeansImpl<__RTYPE__, NA, T, true> > { \
657 private: \
658  const MatrixBase<__RTYPE__, NA, T>& ref; \
659  \
660  typedef detail::RowMeansReturn<__RTYPE__> return_traits; \
661  typedef typename return_traits::type return_vector; \
662  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
663  \
664 public: \
665  RowMeansImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
666  : ref(ref_) \
667  {} \
668  \
669  return_vector get() const { \
670  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
671  return_vector res(nr); \
672  \
673  std::vector<R_xlen_t> n_ok(nr, 0); \
674  \
675  for (j = 0; j < nc; j++) { \
676  for (i = 0; i < nr; i++) { \
677  if (!detail::check_na(ref(i, j))) { \
678  detail::incr(&res[i], ref(i, j)); \
679  ++n_ok[i]; \
680  } \
681  } \
682  } \
683  \
684  for (i = 0; i < nr; i++) { \
685  if (n_ok[i]) { \
686  detail::div(&res[i], n_ok[i]); \
687  } else { \
688  detail::set_nan(&res[i]); \
689  } \
690  } \
691  \
692  return res; \
693  } \
694 };
695 
696 ROW_MEANS_IMPL_RMNA(LGLSXP)
697 ROW_MEANS_IMPL_RMNA(INTSXP)
698 
699 #undef ROW_MEANS_IMPL_RMNA
700 
701 // RowMeans
702 // Input with template parameter NA = false
703 // RowMeansImpl<..., NA_RM = false>
704 //
705 template <int RTYPE, typename T, bool NA_RM>
706 class RowMeansImpl<RTYPE, false, T, NA_RM>
707  : public RowMeansImpl<RTYPE, false, T, false> {};
708 
709 
710 // ColMeans
711 // na.rm = FALSE
712 // default input
713 // default output
714 //
715 template <int RTYPE, bool NA, typename T, bool NA_RM = false>
716 class ColMeansImpl :
717  public Lazy<typename detail::ColMeansReturn<RTYPE>::type, ColMeansImpl<RTYPE, NA, T, NA_RM> > {
718 private:
719  const MatrixBase<RTYPE, NA, T>& ref;
720 
721  typedef detail::ColMeansReturn<RTYPE> return_traits;
722  typedef typename return_traits::type return_vector;
723  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
724 
725 public:
726  ColMeansImpl(const MatrixBase<RTYPE, NA, T>& ref_)
727  : ref(ref_)
728  {}
729 
730  return_vector get() const {
731  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
732  return_vector res(nc);
733 
734  for (j = 0; j < nc; j++) {
735  for (i = 0; i < nr; i++) {
736  detail::incr(&res[j], ref(i, j));
737  }
738  }
739 
740  for (j = 0; j < nc; j++) {
741  detail::div(&res[j], nr);
742  }
743 
744  return res;
745  }
746 };
747 
748 // ColMeans
749 // na.rm = FALSE
750 // LGLSXP / INTSXP input
751 // REALSXP output
752 //
753 #define COL_MEANS_IMPL_KEEPNA(__RTYPE__) \
754  \
755 template <bool NA, typename T, bool NA_RM> \
756 class ColMeansImpl<__RTYPE__, NA, T, NA_RM> : \
757  public Lazy<typename detail::ColMeansReturn<__RTYPE__>::type, ColMeansImpl<__RTYPE__, NA, T, NA_RM> > { \
758 private: \
759  const MatrixBase<__RTYPE__, NA, T>& ref; \
760  \
761  typedef detail::ColMeansReturn<__RTYPE__> return_traits; \
762  typedef typename return_traits::type return_vector; \
763  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
764  \
765  struct bit { \
766  unsigned char x : 1; \
767  }; \
768  \
769 public: \
770  ColMeansImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
771  : ref(ref_) \
772  {} \
773  \
774  return_vector get() const { \
775  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
776  return_vector res(nc); \
777  \
778  std::vector<bit> na_flags(nc); \
779  \
780  for (j = 0; j < nc; j++) { \
781  for (i = 0; i < nr; i++) { \
782  if (detail::check_na(ref(i, j))) { \
783  na_flags[j].x |= 0x1; \
784  } \
785  detail::incr(&res[j], ref(i, j)); \
786  } \
787  } \
788  \
789  for (j = 0; j < nc; j++) { \
790  if (!na_flags[j].x) { \
791  detail::div(&res[j], nr); \
792  } else { \
793  res[j] = NA_REAL; \
794  } \
795  } \
796  \
797  return res; \
798  } \
799 };
800 
801 COL_MEANS_IMPL_KEEPNA(LGLSXP)
802 COL_MEANS_IMPL_KEEPNA(INTSXP)
803 
804 #undef COL_MEANS_IMPL_KEEPNA
805 
806 // ColMeans
807 // na.rm = TRUE
808 // default input
809 // default output
810 //
811 template <int RTYPE, bool NA, typename T>
812 class ColMeansImpl<RTYPE, NA, T, true> :
813  public Lazy<typename detail::ColMeansReturn<RTYPE>::type, ColMeansImpl<RTYPE, NA, T, true> > {
814 private:
815  const MatrixBase<RTYPE, NA, T>& ref;
816 
817  typedef detail::ColMeansReturn<RTYPE> return_traits;
818  typedef typename return_traits::type return_vector;
819  typedef typename traits::storage_type<return_traits::rtype>::type stored_type;
820 
821 public:
822  ColMeansImpl(const MatrixBase<RTYPE, NA, T>& ref_)
823  : ref(ref_)
824  {}
825 
826  return_vector get() const {
827  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol();
828  return_vector res(nc);
829 
830  std::vector<R_xlen_t> n_ok(nc, 0);
831  stored_type current = stored_type();
832 
833 
834  for (j = 0; j < nc; j++) {
835  for (i = 0; i < nr; i++) {
836  current = ref(i, j);
837  if (!detail::check_na(current)) {
838  detail::incr(&res[j], ref(i, j));
839  ++n_ok[j];
840  }
841  }
842  }
843 
844  for (j = 0; j < nc; j++) {
845  if (n_ok[j]) {
846  detail::div(&res[j], n_ok[j]);
847  } else {
848  detail::set_nan(&res[j]);
849  }
850  }
851 
852  return res;
853  }
854 };
855 
856 // ColMeans
857 // na.rm = TRUE
858 // LGLSXP / INTSXP input
859 // REALSXP output
860 //
861 #define COL_MEANS_IMPL_RMNA(__RTYPE__) \
862  \
863 template <bool NA, typename T> \
864 class ColMeansImpl<__RTYPE__, NA, T, true> : \
865  public Lazy<typename detail::ColMeansReturn<__RTYPE__>::type, ColMeansImpl<__RTYPE__, NA, T, true> > { \
866 private: \
867  const MatrixBase<__RTYPE__, NA, T>& ref; \
868  \
869  typedef detail::ColMeansReturn<__RTYPE__> return_traits; \
870  typedef typename return_traits::type return_vector; \
871  typedef typename traits::storage_type<return_traits::rtype>::type stored_type; \
872  \
873 public: \
874  ColMeansImpl(const MatrixBase<__RTYPE__, NA, T>& ref_) \
875  : ref(ref_) \
876  {} \
877  \
878  return_vector get() const { \
879  R_xlen_t i, j, nr = ref.nrow(), nc = ref.ncol(); \
880  return_vector res(nc); \
881  \
882  std::vector<R_xlen_t> n_ok(nc, 0); \
883  \
884  for (j = 0; j < nc; j++) { \
885  for (i = 0; i < nr; i++) { \
886  if (!detail::check_na(ref(i, j))) { \
887  detail::incr(&res[j], ref(i, j)); \
888  ++n_ok[j]; \
889  } \
890  } \
891  } \
892  \
893  for (j = 0; j < nc; j++) { \
894  if (n_ok[j]) { \
895  detail::div(&res[j], n_ok[j]); \
896  } else { \
897  detail::set_nan(&res[j]); \
898  } \
899  } \
900  \
901  return res; \
902  } \
903 };
904 
905 COL_MEANS_IMPL_RMNA(LGLSXP)
906 COL_MEANS_IMPL_RMNA(INTSXP)
907 
908 #undef COL_MEANS_IMPL_RMNA
909 
910 // ColMeans
911 // Input with template parameter NA = false
912 // ColMeansImpl<..., NA_RM = false>
913 //
914 template <int RTYPE, typename T, bool NA_RM>
915 class ColMeansImpl<RTYPE, false, T, NA_RM>
916  : public ColMeansImpl<RTYPE, false, T, false> {};
917 
918 
919 } // sugar
920 
921 
922 template <int RTYPE, bool NA, typename T>
923 inline typename sugar::detail::RowSumsReturn<RTYPE>::type
924 rowSums(const MatrixBase<RTYPE, NA, T>& x, bool na_rm = false) {
925  if (!na_rm) {
926  return sugar::RowSumsImpl<RTYPE, NA, T, false>(x);
927  }
928  return sugar::RowSumsImpl<RTYPE, NA, T, true>(x);
929 }
930 
931 template <int RTYPE, bool NA, typename T>
932 inline typename sugar::detail::ColSumsReturn<RTYPE>::type
933 colSums(const MatrixBase<RTYPE, NA, T>& x, bool na_rm = false) {
934  if (!na_rm) {
935  return sugar::ColSumsImpl<RTYPE, NA, T, false>(x);
936  }
937  return sugar::ColSumsImpl<RTYPE, NA, T, true>(x);
938 }
939 
940 template <int RTYPE, bool NA, typename T>
941 inline typename sugar::detail::RowMeansReturn<RTYPE>::type
942 rowMeans(const MatrixBase<RTYPE, NA, T>& x, bool na_rm = false) {
943  if (!na_rm) {
944  return sugar::RowMeansImpl<RTYPE, NA, T, false>(x);
945  }
946  return sugar::RowMeansImpl<RTYPE, NA, T, true>(x);
947 }
948 
949 template <int RTYPE, bool NA, typename T>
950 inline typename sugar::detail::ColMeansReturn<RTYPE>::type
951 colMeans(const MatrixBase<RTYPE, NA, T>& x, bool na_rm = false) {
952  if (!na_rm) {
953  return sugar::ColMeansImpl<RTYPE, NA, T, false>(x);
954  }
955  return sugar::ColMeansImpl<RTYPE, NA, T, true>(x);
956 }
957 
958 
959 } // Rcpp
960 
961 #endif // Rcpp__sugar__rowSums_h
Rcpp::MatrixBase::nrow
R_xlen_t nrow() const
Definition: MatrixBase.h:48
Rcpp::MatrixBase::ncol
R_xlen_t ncol() const
Definition: MatrixBase.h:49
Rcpp::sugar::ColMeansImpl< RTYPE, NA, T, true >::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:815
Rcpp::sugar::ColMeansImpl< RTYPE, NA, T, true >::ColMeansImpl
ColMeansImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:822
Rcpp::sugar::ColMeansImpl< RTYPE, NA, T, true >::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:819
Rcpp::sugar::ColMeansImpl< RTYPE, NA, T, true >::return_traits
detail::ColMeansReturn< RTYPE > return_traits
Definition: rowSums.h:817
Rcpp::sugar::ColMeansImpl::get
return_vector get() const
Definition: rowSums.h:730
Rcpp::sugar::ColMeansImpl::ColMeansImpl
ColMeansImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:726
Rcpp::sugar::ColMeansImpl::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:719
Rcpp::sugar::ColMeansImpl::return_vector
return_traits::type return_vector
Definition: rowSums.h:722
Rcpp::sugar::ColMeansImpl::return_traits
detail::ColMeansReturn< RTYPE > return_traits
Definition: rowSums.h:721
Rcpp::sugar::ColMeansImpl::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:723
Rcpp::sugar::ColSumsImpl< RTYPE, NA, T, true >::ColSumsImpl
ColSumsImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:420
Rcpp::sugar::ColSumsImpl< RTYPE, NA, T, true >::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:417
Rcpp::sugar::ColSumsImpl< RTYPE, NA, T, true >::return_traits
detail::ColSumsReturn< RTYPE > return_traits
Definition: rowSums.h:415
Rcpp::sugar::ColSumsImpl< RTYPE, NA, T, true >::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:413
Rcpp::sugar::ColSumsImpl::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:323
Rcpp::sugar::ColSumsImpl::ColSumsImpl
ColSumsImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:330
Rcpp::sugar::ColSumsImpl::get
return_vector get() const
Definition: rowSums.h:334
Rcpp::sugar::ColSumsImpl::return_vector
return_traits::type return_vector
Definition: rowSums.h:326
Rcpp::sugar::ColSumsImpl::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:327
Rcpp::sugar::ColSumsImpl::return_traits
detail::ColSumsReturn< RTYPE > return_traits
Definition: rowSums.h:325
Rcpp::sugar::RowMeansImpl< RTYPE, NA, T, true >::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:606
Rcpp::sugar::RowMeansImpl< RTYPE, NA, T, true >::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:610
Rcpp::sugar::RowMeansImpl< RTYPE, NA, T, true >::return_traits
detail::RowMeansReturn< RTYPE > return_traits
Definition: rowSums.h:608
Rcpp::sugar::RowMeansImpl< RTYPE, NA, T, true >::RowMeansImpl
RowMeansImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:613
Rcpp::sugar::RowMeansImpl::return_vector
return_traits::type return_vector
Definition: rowSums.h:513
Rcpp::sugar::RowMeansImpl::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:510
Rcpp::sugar::RowMeansImpl::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:514
Rcpp::sugar::RowMeansImpl::get
return_vector get() const
Definition: rowSums.h:521
Rcpp::sugar::RowMeansImpl::return_traits
detail::RowMeansReturn< RTYPE > return_traits
Definition: rowSums.h:512
Rcpp::sugar::RowMeansImpl::RowMeansImpl
RowMeansImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:517
Rcpp::sugar::RowSumsImpl< RTYPE, NA, T, true >::RowSumsImpl
RowSumsImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:236
Rcpp::sugar::RowSumsImpl< RTYPE, NA, T, true >::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:233
Rcpp::sugar::RowSumsImpl< RTYPE, NA, T, true >::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:229
Rcpp::sugar::RowSumsImpl< RTYPE, NA, T, true >::return_traits
detail::RowSumsReturn< RTYPE > return_traits
Definition: rowSums.h:231
Rcpp::sugar::RowSumsImpl::ref
const MatrixBase< RTYPE, NA, T > & ref
Definition: rowSums.h:131
Rcpp::sugar::RowSumsImpl::get
return_vector get() const
Definition: rowSums.h:142
Rcpp::sugar::RowSumsImpl::stored_type
traits::storage_type< return_traits::rtype >::type stored_type
Definition: rowSums.h:135
Rcpp::sugar::RowSumsImpl::return_traits
detail::RowSumsReturn< RTYPE > return_traits
Definition: rowSums.h:133
Rcpp::sugar::RowSumsImpl::RowSumsImpl
RowSumsImpl(const MatrixBase< RTYPE, NA, T > &ref_)
Definition: rowSums.h:138
Rcpp::sugar::RowSumsImpl::return_vector
return_traits::type return_vector
Definition: rowSums.h:134
Rcpp::sugar::detail::check_na
bool check_na(double x)
Definition: rowSums.h:30
Rcpp::sugar::detail::incr
void incr(double *lhs, double rhs)
Definition: rowSums.h:51
Rcpp::sugar::detail::div
void div(double *lhs, R_xlen_t rhs)
Definition: rowSums.h:65
Rcpp::sugar::detail::set_nan
void set_nan(double *x)
Definition: rowSums.h:75
Rcpp
Rcpp API.
Definition: algo.h:28
Rcpp::colMeans
sugar::detail::ColMeansReturn< RTYPE >::type colMeans(const MatrixBase< RTYPE, NA, T > &x, bool na_rm=false)
Definition: rowSums.h:951
Rcpp::rowSums
sugar::detail::RowSumsReturn< RTYPE >::type rowSums(const MatrixBase< RTYPE, NA, T > &x, bool na_rm=false)
Definition: rowSums.h:924
Rcpp::colSums
sugar::detail::ColSumsReturn< RTYPE >::type colSums(const MatrixBase< RTYPE, NA, T > &x, bool na_rm=false)
Definition: rowSums.h:933
Rcpp::rowMeans
sugar::detail::RowMeansReturn< RTYPE >::type rowMeans(const MatrixBase< RTYPE, NA, T > &x, bool na_rm=false)
Definition: rowSums.h:942
Rcpp::NA
static Na_Proxy NA
Definition: Na_Proxy.h:52
ROW_SUMS_IMPL_RMNA
#define ROW_SUMS_IMPL_RMNA(__RTYPE__)
Definition: rowSums.h:264
ROW_MEANS_IMPL_KEEPNA
#define ROW_MEANS_IMPL_KEEPNA(__RTYPE__)
Definition: rowSums.h:544
ROW_SUMS_IMPL_KEEPNA
#define ROW_SUMS_IMPL_KEEPNA(__RTYPE__)
Definition: rowSums.h:169
COL_SUMS_IMPL_RMNA
#define COL_SUMS_IMPL_RMNA(__RTYPE__)
Definition: rowSums.h:448
ROW_MEANS_IMPL_RMNA
#define ROW_MEANS_IMPL_RMNA(__RTYPE__)
Definition: rowSums.h:652
COL_MEANS_IMPL_KEEPNA
#define COL_MEANS_IMPL_KEEPNA(__RTYPE__)
Definition: rowSums.h:753
COL_SUMS_IMPL_KEEPNA
#define COL_SUMS_IMPL_KEEPNA(__RTYPE__)
Definition: rowSums.h:353
COL_MEANS_IMPL_RMNA
#define COL_MEANS_IMPL_RMNA(__RTYPE__)
Definition: rowSums.h:861
Rcpp::traits::storage_type::type
SEXP type
Definition: storage_type.h:36