Sugar median with unit tests

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

@@ -80,5 +80,7 @@
 #include <Rcpp/sugar/functions/cummin.h>
 #include <Rcpp/sugar/functions/cummax.h>
 
+#include <Rcpp/sugar/functions/median.h>
+
 #endif
 

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

@@ -0,0 +1,290 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+//
+// median.h: Rcpp R/C++ interface class library -- median
+//
+// Copyright (C) 2016 Dirk Eddelbuettel, Romain Francois, and Nathan Russell
+//
+// This file is part of Rcpp.
+//
+// Rcpp is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// Rcpp is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Rcpp.  If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef Rcpp__sugar__median_h
+#define Rcpp__sugar__median_h
+
+namespace Rcpp {
+namespace sugar {
+namespace median_detail {
+
+// need to return double for integer vectors 
+// (in case of even-length input vector)
+// and Rcpp::String for STRSXP
+// also need to return NA_REAL for 
+// integer vector yielding NA result
+template <int RTYPE>
+struct result {
+    typedef typename Rcpp::traits::storage_type<RTYPE>::type type;
+    enum { rtype = RTYPE };
+};
+
+template <>
+struct result<INTSXP> {
+    typedef double type;
+    enum { rtype = REALSXP };
+};
+
+template <>
+struct result<STRSXP> {
+    typedef Rcpp::String type;
+    enum { rtype = STRSXP };
+};
+
+// std::nth_element and std::max_element don't 
+// know how to compare Rcomplex values
+template <typename T>
+inline bool less(T lhs, T rhs) {
+    return lhs < rhs;
+}
+
+template<>
+inline bool less<Rcomplex>(Rcomplex lhs, Rcomplex rhs) {
+    if (lhs.r < rhs.r) return true;
+    if (lhs.i < rhs.i) return true;
+    return false; 
+}
+
+// compiler does not know how to handle 
+// Rcomplex numerator / double denominator
+// and need explicit cast for INTSXP case
+inline double half(double lhs) {
+    return lhs / 2.0;
+}
+
+inline double half(int lhs) {
+    return static_cast<double>(lhs) / 2.0;
+}
+
+inline Rcomplex half(Rcomplex lhs) {
+    lhs.r /= 2.0;
+    lhs.i /= 2.0;
+    return lhs;
+}
+
+} // median_detail
+
+// base case
+template <int RTYPE, bool NA, typename T, bool NA_RM = false>
+class Median {
+public:
+    typedef typename median_detail::result<RTYPE>::type result_type;
+    typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
+    enum { RESULT_RTYPE = median_detail::result<RTYPE>::rtype };
+    typedef T VECTOR;
+
+private:
+    VECTOR x;
+
+public:
+    Median(const VECTOR& xx)
+        : x(Rcpp::clone(xx)) {}
+
+    operator result_type() {
+        if (x.size() < 1) {
+            return Rcpp::traits::get_na<RESULT_RTYPE>();
+        }
+
+        if (Rcpp::any(Rcpp::is_na(x))) {
+            return Rcpp::traits::get_na<RESULT_RTYPE>();
+        }
+
+        R_xlen_t n = x.size() / 2;
+        std::nth_element(
+            x.begin(), x.begin() + n, x.end(), 
+            median_detail::less<stored_type>);
+
+        if (x.size() % 2) return x[n];
+        return median_detail::half(
+            x[n] + *std::max_element(
+                    x.begin(), x.begin() + n, 
+                    median_detail::less<stored_type>));
+    }
+};
+
+// na.rm = TRUE
+template <int RTYPE, bool NA, typename T>
+class Median<RTYPE, NA, T, true> {
+public:
+    typedef typename median_detail::result<RTYPE>::type result_type;
+    typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
+    enum { RESULT_RTYPE = median_detail::result<RTYPE>::rtype };
+    typedef T VECTOR;
+
+private:
+    VECTOR x;
+
+public:
+    Median(const VECTOR& xx)
+        : x(Rcpp::na_omit(Rcpp::clone(xx))) {}
+
+    operator result_type() {
+        if (!x.size()) {
+            return Rcpp::traits::get_na<RESULT_RTYPE>();
+        }
+
+        R_xlen_t n = x.size() / 2;
+        std::nth_element(
+            x.begin(), x.begin() + n, x.end(), 
+            median_detail::less<stored_type>);
+
+        if (x.size() % 2) return x[n];
+        return median_detail::half(
+            x[n] + *std::max_element(
+                    x.begin(), x.begin() + n, 
+                    median_detail::less<stored_type>));
+    }
+};
+
+// NA = false
+template <int RTYPE, typename T, bool NA_RM>
+class Median<RTYPE, false, T, NA_RM> {
+public:
+    typedef typename median_detail::result<RTYPE>::type result_type;
+    typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
+    enum { RESULT_RTYPE = median_detail::result<RTYPE>::rtype };
+    typedef T VECTOR;
+
+private:
+    VECTOR x;
+
+public:
+    Median(const VECTOR& xx)
+        : x(Rcpp::clone(xx)) {}
+
+    operator result_type() {
+        if (x.size() < 1) {
+            return Rcpp::traits::get_na<RESULT_RTYPE>();
+        }
+
+        R_xlen_t n = x.size() / 2;
+        std::nth_element(
+            x.begin(), x.begin() + n, x.end(), 
+            median_detail::less<stored_type>);
+
+        if (x.size() % 2) return x[n];
+        return median_detail::half(
+            x[n] + *std::max_element(
+                    x.begin(), x.begin() + n, 
+                    median_detail::less<stored_type>));
+    }
+};
+
+// specialize for character vector
+// due to string_proxy's incompatibility
+// with certain std:: algorithms;
+// need to return NA for even-length vectors
+template <bool NA, typename T, bool NA_RM>
+class Median<STRSXP, NA, T, NA_RM> {
+public:
+    typedef typename median_detail::result<STRSXP>::type result_type;
+    typedef typename Rcpp::traits::storage_type<STRSXP>::type stored_type;
+    typedef T VECTOR;
+
+private:
+    VECTOR x;
+
+public:
+    Median(const VECTOR& xx)
+        : x(Rcpp::clone(xx)) {}
+
+    operator result_type() {
+        if (!(x.size() % 2)) {
+            return Rcpp::traits::get_na<STRSXP>();
+        }
+
+        if (Rcpp::any(Rcpp::is_na(x))) {
+            return Rcpp::traits::get_na<STRSXP>();
+        }
+
+        R_xlen_t n = x.size() / 2;
+        x.sort();
+
+        return x[n];
+    }
+};
+
+// na.rm = TRUE
+template <bool NA, typename T>
+class Median<STRSXP, NA, T, true> {
+public:
+    typedef typename median_detail::result<STRSXP>::type result_type;
+    typedef typename Rcpp::traits::storage_type<STRSXP>::type stored_type;
+    typedef T VECTOR;
+
+private:
+    VECTOR x;
+
+public:
+    Median(const VECTOR& xx)
+        : x(Rcpp::na_omit(Rcpp::clone(xx))) {}
+
+    operator result_type() {
+        if (!(x.size() % 2)) {
+            return Rcpp::traits::get_na<STRSXP>();
+        }
+
+        R_xlen_t n = x.size() / 2;
+        x.sort();
+
+        return x[n];
+    }
+};
+
+// NA = false
+template <typename T>
+class Median<STRSXP, false, T, true> {
+public:
+    typedef typename median_detail::result<STRSXP>::type result_type;
+    typedef typename Rcpp::traits::storage_type<STRSXP>::type stored_type;
+    typedef T VECTOR;
+
+private:
+    VECTOR x;
+
+public:
+    Median(const VECTOR& xx)
+        : x(Rcpp::clone(xx)) {}
+
+    operator result_type() {
+        if (!(x.size() % 2)) {
+            return Rcpp::traits::get_na<STRSXP>();
+        }
+
+        R_xlen_t n = x.size() / 2;
+        x.sort();
+
+        return x[n];
+    }
+};
+
+} // sugar
+
+template <int RTYPE, bool NA, typename T>
+inline typename sugar::median_detail::result<RTYPE>::type
+median(const Rcpp::VectorBase<RTYPE, NA, T>& x, bool na_rm = false) {
+    if (!na_rm) return sugar::Median<RTYPE, NA, T, false>(x);
+    return sugar::Median<RTYPE, NA, T, true>(x);
+}
+
+} // Rcpp
+
+#endif // Rcpp__sugar__median_h

inst/unitTests/cpp/sugar.cpp

@@ -716,3 +716,25 @@ IntegerVector runit_cummax_iv(IntegerVector xx){
 }
 
 
+// 18 January 2016: median 
+
+// [[Rcpp::export]]
+double median_dbl(Rcpp::NumericVector x, bool na_rm = false) {
+    return Rcpp::median(x, na_rm);
+}
+
+// [[Rcpp::export]]
+double median_int(Rcpp::IntegerVector x, bool na_rm = false) {
+    return Rcpp::median(x, na_rm);
+}
+
+// [[Rcpp::export]]
+Rcomplex median_cx(Rcpp::ComplexVector x, bool na_rm = false) {
+    return Rcpp::median(x, na_rm);
+}
+
+// [[Rcpp::export]]
+Rcpp::String median_ch(Rcpp::CharacterVector x, bool na_rm = false) {
+    return Rcpp::median(x, na_rm);
+}
+