Rcpp Version 0.12.12
median.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 // median.h: Rcpp R/C++ interface class library -- median
4 //
5 // Copyright (C) 2016 Dirk Eddelbuettel, Romain Francois, and 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__median_h
23 #define Rcpp__sugar__median_h
24 
25 namespace Rcpp {
26 namespace sugar {
27 namespace median_detail {
28 
29 // need to return double for integer vectors
30 // (in case of even-length input vector)
31 // and Rcpp::String for STRSXP
32 // also need to return NA_REAL for
33 // integer vector yielding NA result
34 template <int RTYPE>
35 struct result {
37  enum { rtype = RTYPE };
38 };
39 
40 template <>
41 struct result<INTSXP> {
42  typedef double type;
43  enum { rtype = REALSXP };
44 };
45 
46 template <>
47 struct result<STRSXP> {
48  typedef Rcpp::String type;
49  enum { rtype = STRSXP };
50 };
51 
52 // std::nth_element and std::max_element don't
53 // know how to compare Rcomplex values
54 template <typename T>
55 inline bool less(T lhs, T rhs) {
56  return lhs < rhs;
57 }
58 
59 template<>
60 inline bool less<Rcomplex>(Rcomplex lhs, Rcomplex rhs) {
61  if (lhs.r < rhs.r) return true;
62  if (lhs.i < rhs.i) return true;
63  return false;
64 }
65 
66 // compiler does not know how to handle
67 // Rcomplex numerator / double denominator
68 // and need explicit cast for INTSXP case
69 inline double half(double lhs) {
70  return lhs / 2.0;
71 }
72 
73 inline double half(int lhs) {
74  return static_cast<double>(lhs) / 2.0;
75 }
76 
77 inline Rcomplex half(Rcomplex lhs) {
78  lhs.r /= 2.0;
79  lhs.i /= 2.0;
80  return lhs;
81 }
82 
83 } // median_detail
84 
85 // base case
86 template <int RTYPE, bool NA, typename T, bool NA_RM = false>
87 class Median {
88 public:
91  enum { RESULT_RTYPE = median_detail::result<RTYPE>::rtype };
92  typedef T VECTOR;
93 
94 private:
95  VECTOR x;
96 
97 public:
98  Median(const VECTOR& xx)
99  : x(Rcpp::clone(xx)) {}
100 
101  operator result_type() {
102  if (x.size() < 1) {
103  return Rcpp::traits::get_na<RESULT_RTYPE>();
104  }
105 
106  if (Rcpp::any(Rcpp::is_na(x))) {
107  return Rcpp::traits::get_na<RESULT_RTYPE>();
108  }
109 
110  R_xlen_t n = x.size() / 2;
111  std::nth_element(
112  x.begin(), x.begin() + n, x.end(),
113  median_detail::less<stored_type>);
114 
115  if (x.size() % 2) return x[n];
116  return median_detail::half(
117  x[n] + *std::max_element(
118  x.begin(), x.begin() + n,
119  median_detail::less<stored_type>));
120  }
121 };
122 
123 // na.rm = TRUE
124 template <int RTYPE, bool NA, typename T>
125 class Median<RTYPE, NA, T, true> {
126 public:
129  enum { RESULT_RTYPE = median_detail::result<RTYPE>::rtype };
130  typedef T VECTOR;
131 
132 private:
133  VECTOR x;
134 
135 public:
136  Median(const VECTOR& xx)
137  : x(Rcpp::na_omit(Rcpp::clone(xx))) {}
138 
139  operator result_type() {
140  if (!x.size()) {
141  return Rcpp::traits::get_na<RESULT_RTYPE>();
142  }
143 
144  R_xlen_t n = x.size() / 2;
145  std::nth_element(
146  x.begin(), x.begin() + n, x.end(),
147  median_detail::less<stored_type>);
148 
149  if (x.size() % 2) return x[n];
150  return median_detail::half(
151  x[n] + *std::max_element(
152  x.begin(), x.begin() + n,
153  median_detail::less<stored_type>));
154  }
155 };
156 
157 // NA = false
158 template <int RTYPE, typename T, bool NA_RM>
159 class Median<RTYPE, false, T, NA_RM> {
160 public:
163  enum { RESULT_RTYPE = median_detail::result<RTYPE>::rtype };
164  typedef T VECTOR;
165 
166 private:
167  VECTOR x;
168 
169 public:
170  Median(const VECTOR& xx)
171  : x(Rcpp::clone(xx)) {}
172 
173  operator result_type() {
174  if (x.size() < 1) {
175  return Rcpp::traits::get_na<RESULT_RTYPE>();
176  }
177 
178  R_xlen_t n = x.size() / 2;
179  std::nth_element(
180  x.begin(), x.begin() + n, x.end(),
181  median_detail::less<stored_type>);
182 
183  if (x.size() % 2) return x[n];
184  return median_detail::half(
185  x[n] + *std::max_element(
186  x.begin(), x.begin() + n,
187  median_detail::less<stored_type>));
188  }
189 };
190 
191 // specialize for character vector
192 // due to string_proxy's incompatibility
193 // with certain std:: algorithms;
194 // need to return NA for even-length vectors
195 template <bool NA, typename T, bool NA_RM>
196 class Median<STRSXP, NA, T, NA_RM> {
197 public:
200  typedef T VECTOR;
201 
202 private:
203  VECTOR x;
204 
205 public:
206  Median(const VECTOR& xx)
207  : x(Rcpp::clone(xx)) {}
208 
209  operator result_type() {
210  if (!(x.size() % 2)) {
212  }
213 
214  if (Rcpp::any(Rcpp::is_na(x))) {
216  }
217 
218  R_xlen_t n = x.size() / 2;
219  x.sort();
220 
221  return x[n];
222  }
223 };
224 
225 // na.rm = TRUE
226 template <bool NA, typename T>
227 class Median<STRSXP, NA, T, true> {
228 public:
231  typedef T VECTOR;
232 
233 private:
234  VECTOR x;
235 
236 public:
237  Median(const VECTOR& xx)
238  : x(Rcpp::na_omit(Rcpp::clone(xx))) {}
239 
240  operator result_type() {
241  if (!(x.size() % 2)) {
243  }
244 
245  R_xlen_t n = x.size() / 2;
246  x.sort();
247 
248  return x[n];
249  }
250 };
251 
252 // NA = false
253 template <typename T>
254 class Median<STRSXP, false, T, true> {
255 public:
258  typedef T VECTOR;
259 
260 private:
261  VECTOR x;
262 
263 public:
264  Median(const VECTOR& xx)
265  : x(Rcpp::clone(xx)) {}
266 
267  operator result_type() {
268  if (!(x.size() % 2)) {
270  }
271 
272  R_xlen_t n = x.size() / 2;
273  x.sort();
274 
275  return x[n];
276  }
277 };
278 
279 } // sugar
280 
281 template <int RTYPE, bool NA, typename T>
283 median(const Rcpp::VectorBase<RTYPE, NA, T>& x, bool na_rm = false) {
284  if (!na_rm) return sugar::Median<RTYPE, NA, T, false>(x);
286 }
287 
288 } // Rcpp
289 
290 #endif // Rcpp__sugar__median_h
Rcpp::traits::storage_type< RTYPE >::type stored_type
Definition: median.h:128
double half(double lhs)
Definition: median.h:69
Rcpp::traits::storage_type< RTYPE >::type stored_type
Definition: median.h:90
median_detail::result< STRSXP >::type result_type
Definition: median.h:229
Rcpp::traits::storage_type< STRSXP >::type stored_type
Definition: median.h:257
SEXP get_na< STRSXP >()
Definition: get_na.h:50
sugar::median_detail::result< RTYPE >::type median(const Rcpp::VectorBase< RTYPE, NA, T > &x, bool na_rm=false)
Definition: median.h:283
Rcpp::traits::storage_type< RTYPE >::type stored_type
Definition: median.h:162
bool less< Rcomplex >(Rcomplex lhs, Rcomplex rhs)
Definition: median.h:60
median_detail::result< RTYPE >::type result_type
Definition: median.h:127
Rcpp::traits::storage_type< RTYPE >::type type
Definition: median.h:36
bool any(InputIterator first, InputIterator last, const T &value)
Definition: algo.h:89
Rcpp::traits::storage_type< STRSXP >::type stored_type
Definition: median.h:199
median_detail::result< RTYPE >::type result_type
Definition: median.h:161
median_detail::result< RTYPE >::type result_type
Definition: median.h:89
static Na_Proxy NA
Definition: Na_Proxy.h:52
median_detail::result< STRSXP >::type result_type
Definition: median.h:256
sugar::IsNa< RTYPE, NA, T > is_na(const Rcpp::VectorBase< RTYPE, NA, T > &t)
Definition: is_na.h:87
Rcpp::traits::storage_type< STRSXP >::type stored_type
Definition: median.h:230
bool less(T lhs, T rhs)
Definition: median.h:55
Rcpp API.
Definition: algo.h:28
Median(const VECTOR &xx)
Definition: median.h:98
T clone(const T &object)
Definition: clone.h:33
median_detail::result< STRSXP >::type result_type
Definition: median.h:198
Vector< RTYPE > na_omit(const VectorBase< RTYPE, NA, T > &t)
Definition: na_omit.h:75