Skip to content

Commit 7bad63e

Browse files
nathan-russellnathan-russell
committed
Sugar function median with unit tests
1 parent a8bbece commit 7bad63e

File tree

4 files changed

+436
-0
lines changed

4 files changed

+436
-0
lines changed

‎inst/include/Rcpp/sugar/functions/functions.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,5 +80,7 @@
8080
#include <Rcpp/sugar/functions/cummin.h>
8181
#include <Rcpp/sugar/functions/cummax.h>
8282

83+
#include <Rcpp/sugar/functions/median.h>
84+
8385
#endif
8486

‎inst/include/Rcpp/sugar/functions/median.h

Lines changed: 297 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,297 @@
1+
// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
2+
//
3+
// median.h: Rcpp R/C++ interface class library -- median
4+
//
5+
// Copyright (C) 2012 - 2013 Dirk Eddelbuettel and Romain Francois
6+
// Copyright (C) 2016 - Nathan Russell
7+
//
8+
// This file is part of Rcpp.
9+
//
10+
// Rcpp is free software: you can redistribute it and/or modify it
11+
// under the terms of the GNU General Public License as published by
12+
// the Free Software Foundation, either version 2 of the License, or
13+
// (at your option) any later version.
14+
//
15+
// Rcpp is distributed in the hope that it will be useful, but
16+
// WITHOUT ANY WARRANTY; without even the implied warranty of
17+
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18+
// GNU General Public License for more details.
19+
//
20+
// You should have received a copy of the GNU General Public License
21+
// along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
22+
23+
#ifndef Rcpp__sugar__median_h
24+
#define Rcpp__sugar__median_h
25+
26+
namespace Rcpp {
27+
namespace sugar {
28+
namespace median_detail {
29+
30+
// need to return double for integer vectors
31+
// (in case of even-length input vector)
32+
// and Rcpp::String for STRSXP
33+
// also need to return NA_REAL for
34+
// integer vector yielding NA result
35+
template <int RTYPE>
36+
struct result {
37+
typedef typename Rcpp::traits::storage_type<RTYPE>::type type;
38+
enum { rtype = RTYPE };
39+
};
40+
41+
template <>
42+
struct result<INTSXP> {
43+
typedef double type;
44+
enum { rtype = REALSXP };
45+
};
46+
47+
template <>
48+
struct result<STRSXP> {
49+
typedef Rcpp::String type;
50+
enum { rtype = STRSXP };
51+
};
52+
53+
// std::nth_element and std::max_element don't
54+
// know how to compare Rcomplex values
55+
template <typename T>
56+
inline bool less(T lhs, T rhs) {
57+
return lhs < rhs;
58+
}
59+
60+
template<>
61+
inline bool less<Rcomplex>(Rcomplex lhs, Rcomplex rhs) {
62+
if (lhs.r < rhs.r) return true;
63+
if (lhs.i < rhs.i) return true;
64+
return false;
65+
}
66+
67+
// compiler does not know how to handle
68+
// Rcomplex numerator / double denominator
69+
// and need explicit cast for INTSXP case
70+
inline double half(double lhs) {
71+
return lhs / 2.0;
72+
}
73+
74+
inline double half(int lhs) {
75+
return static_cast<double>(lhs) / 2.0;
76+
}
77+
78+
inline Rcomplex half(Rcomplex lhs) {
79+
lhs.r /= 2.0;
80+
lhs.i /= 2.0;
81+
return lhs;
82+
}
83+
84+
} // median_detail
85+
86+
// base case
87+
template <int RTYPE, bool NA, typename T, bool NA_RM = false>
88+
class Median {
89+
public:
90+
typedef typename median_detail::result<RTYPE>::type result_type;
91+
typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
92+
enum { RTYPE2 = median_detail::result<RTYPE>::rtype };
93+
typedef T VECTOR;
94+
95+
private:
96+
VECTOR x;
97+
98+
public:
99+
Median(const VECTOR& xx)
100+
: x(Rcpp::clone(xx)) {}
101+
102+
operator result_type() {
103+
if (x.size() < 1) {
104+
return Rcpp::traits::get_na<RTYPE2>();
105+
}
106+
107+
if (Rcpp::any(Rcpp::is_na(x))) {
108+
return Rcpp::traits::get_na<RTYPE2>();
109+
}
110+
111+
R_xlen_t n = x.size() / 2;
112+
std::nth_element(
113+
x.begin(), x.begin() + n, x.end(),
114+
median_detail::less<stored_type>);
115+
116+
if (x.size() % 2) return x[n];
117+
return median_detail::half(
118+
x[n] + *std::max_element(
119+
x.begin(), x.begin() + n,
120+
median_detail::less<stored_type>));
121+
}
122+
};
123+
124+
// na.rm = TRUE
125+
template <int RTYPE, bool NA, typename T>
126+
class Median<RTYPE, NA, T, true> {
127+
public:
128+
typedef typename median_detail::result<RTYPE>::type result_type;
129+
typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
130+
enum { RTYPE2 = median_detail::result<RTYPE>::rtype };
131+
typedef T VECTOR;
132+
133+
private:
134+
VECTOR x;
135+
136+
public:
137+
Median(const VECTOR& xx)
138+
: x(Rcpp::na_omit(Rcpp::clone(xx))) {}
139+
140+
operator result_type() {
141+
if (!x.size()) {
142+
return Rcpp::traits::get_na<RTYPE2>();
143+
}
144+
145+
R_xlen_t n = x.size() / 2;
146+
std::nth_element(
147+
x.begin(), x.begin() + n, x.end(),
148+
median_detail::less<stored_type>);
149+
150+
if (x.size() % 2) return x[n];
151+
return median_detail::half(
152+
x[n] + *std::max_element(
153+
x.begin(), x.begin() + n,
154+
median_detail::less<stored_type>));
155+
}
156+
};
157+
158+
// NA = false
159+
template <int RTYPE, typename T, bool NA_RM>
160+
class Median<RTYPE, false, T, NA_RM> {
161+
public:
162+
typedef typename median_detail::result<RTYPE>::type result_type;
163+
typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
164+
enum { RTYPE2 = median_detail::result<RTYPE>::rtype };
165+
typedef T VECTOR;
166+
167+
private:
168+
VECTOR x;
169+
170+
public:
171+
Median(const VECTOR& xx)
172+
: x(Rcpp::clone(xx)) {}
173+
174+
operator result_type() {
175+
if (x.size() < 1) {
176+
return Rcpp::traits::get_na<RTYPE2>();
177+
}
178+
179+
R_xlen_t n = x.size() / 2;
180+
std::nth_element(
181+
x.begin(), x.begin() + n, x.end(),
182+
median_detail::less<stored_type>);
183+
184+
if (x.size() % 2) return x[n];
185+
return median_detail::half(
186+
x[n] + *std::max_element(
187+
x.begin(), x.begin() + n,
188+
median_detail::less<stored_type>));
189+
}
190+
};
191+
192+
// specialize for character vector
193+
// due to string_proxy's incompatibility
194+
// with certain std:: algorithms;
195+
// need to return NA for even-length vectors
196+
template <bool NA, typename T, bool NA_RM>
197+
class Median<STRSXP, NA, T, NA_RM> {
198+
public:
199+
typedef typename median_detail::result<STRSXP>::type result_type;
200+
typedef typename Rcpp::traits::storage_type<STRSXP>::type stored_type;
201+
typedef T VECTOR;
202+
203+
private:
204+
VECTOR x;
205+
206+
public:
207+
Median(const VECTOR& xx)
208+
: x(Rcpp::clone(xx)) {}
209+
210+
operator result_type() {
211+
if (!(x.size() % 2)) {
212+
return Rcpp::traits::get_na<STRSXP>();
213+
}
214+
215+
if (Rcpp::any(Rcpp::is_na(x))) {
216+
return Rcpp::traits::get_na<STRSXP>();
217+
}
218+
219+
R_xlen_t n = x.size() / 2;
220+
x.sort();
221+
222+
return x[n];
223+
}
224+
};
225+
226+
// na.rm = TRUE
227+
template <bool NA, typename T>
228+
class Median<STRSXP, NA, T, true> {
229+
public:
230+
typedef typename median_detail::result<STRSXP>::type result_type;
231+
typedef typename Rcpp::traits::storage_type<STRSXP>::type stored_type;
232+
typedef T VECTOR;
233+
234+
private:
235+
VECTOR x;
236+
237+
public:
238+
Median(const VECTOR& xx)
239+
: x(Rcpp::na_omit(Rcpp::clone(xx))) {}
240+
241+
operator result_type() {
242+
if (!(x.size() % 2)) {
243+
return Rcpp::traits::get_na<STRSXP>();
244+
}
245+
246+
R_xlen_t n = x.size() / 2;
247+
x.sort();
248+
249+
return x[n];
250+
}
251+
};
252+
253+
// NA = false
254+
template <typename T>
255+
class Median<STRSXP, false, T, true> {
256+
public:
257+
typedef typename median_detail::result<STRSXP>::type result_type;
258+
typedef typename Rcpp::traits::storage_type<STRSXP>::type stored_type;
259+
typedef T VECTOR;
260+
261+
private:
262+
VECTOR x;
263+
264+
public:
265+
Median(const VECTOR& xx)
266+
: x(Rcpp::clone(xx)) {}
267+
268+
operator result_type() {
269+
if (!(x.size() % 2)) {
270+
return Rcpp::traits::get_na<STRSXP>();
271+
}
272+
273+
R_xlen_t n = x.size() / 2;
274+
x.sort();
275+
276+
return x[n];
277+
}
278+
};
279+
280+
} // sugar
281+
282+
template <int RTYPE, bool NA, typename T>
283+
inline typename sugar::median_detail::result<RTYPE>::type
284+
median(const Rcpp::VectorBase<RTYPE, NA, T>& x, bool na_rm = false) {
285+
switch (static_cast<int>(na_rm)) {
286+
case 0:
287+
return sugar::Median<RTYPE, NA, T, false>(x);
288+
case 1:
289+
return sugar::Median<RTYPE, NA, T, true>(x);
290+
default:
291+
return Rcpp::traits::get_na<RTYPE>();
292+
}
293+
}
294+
295+
} // Rcpp
296+
297+
#endif // Rcpp__sugar__median_h

‎inst/unitTests/cpp/sugar.cpp

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -716,3 +716,25 @@ IntegerVector runit_cummax_iv(IntegerVector xx){
716716
}
717717

718718

719+
// 18 January 2016: median
720+
721+
// [[Rcpp::export]]
722+
double median_dbl(Rcpp::NumericVector x, bool na_rm = false) {
723+
return Rcpp::median(x, na_rm);
724+
}
725+
726+
// [[Rcpp::export]]
727+
double median_int(Rcpp::IntegerVector x, bool na_rm = false) {
728+
return Rcpp::median(x, na_rm);
729+
}
730+
731+
// [[Rcpp::export]]
732+
Rcomplex median_cx(Rcpp::ComplexVector x, bool na_rm = false) {
733+
return Rcpp::median(x, na_rm);
734+
}
735+
736+
// [[Rcpp::export]]
737+
Rcpp::String median_ch(Rcpp::CharacterVector x, bool na_rm = false) {
738+
return Rcpp::median(x, na_rm);
739+
}
740+

0 commit comments

Comments
 (0)