Sugar Median

Median-With-NA-Option.cpp

#include <Rcpp.h>

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 { RTYPE2 = 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<RTYPE2>();
        }
        
        if (Rcpp::any(Rcpp::is_na(x))) {
            return Rcpp::traits::get_na<RTYPE2>();
        }
        
        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 { RTYPE2 = 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<RTYPE2>();
        }
        
        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 { RTYPE2 = 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<RTYPE2>();
        }
        
        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];
    }
};

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];
    }
};

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) {
    switch (static_cast<int>(na_rm)) {
        case 0:
            return sugar::Median<RTYPE, NA, T, false>(x);
        case 1:
            return sugar::Median<RTYPE, NA, T, true>(x);
        default:
            return Rcpp::traits::get_na<RTYPE>();
    }
}

// [[Rcpp::export]]
double median_dbl(Rcpp::NumericVector x, bool na_rm = false) {
    return median(x, na_rm);
}

// [[Rcpp::export]]
double median_int(Rcpp::IntegerVector x, bool na_rm = false) {
    return median(x, na_rm);
}

// [[Rcpp::export]]
Rcomplex median_cx(Rcpp::ComplexVector x, bool na_rm = false) {
    return median(x, na_rm);
}

// [[Rcpp::export]]
Rcpp::String median_ch(Rcpp::CharacterVector x, bool na_rm = false) {
    return median(x, na_rm);
}

Median-Without-NA-Option.cpp

#include <Rcpp.h>

namespace sugar {
namespace median_detail {

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>
class Median {
public:
    typedef typename median_detail::result<RTYPE>::type result_type;
    typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
    enum { RTYPE2 = 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<RTYPE2>();
        }
        
        if (Rcpp::any(Rcpp::is_na(x))) {
            return Rcpp::traits::get_na<RTYPE2>();
        }
        
        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>
class Median<RTYPE, false, T> {
public:
    typedef typename median_detail::result<RTYPE>::type result_type;
    typedef typename Rcpp::traits::storage_type<RTYPE>::type stored_type;
    enum { RTYPE2 = 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<RTYPE2>();
        }
        
        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>
class Median<STRSXP, NA, T> {
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];
    }
};

template <typename T>
class Median<STRSXP, false, T> {
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 sugar::Median<RTYPE, NA, T>
median(const Rcpp::VectorBase<RTYPE, NA, T>& x) {
    return sugar::Median<RTYPE, NA, T>(x);
}


// [[Rcpp::export]]
double median_dbl(Rcpp::NumericVector x) {
    return median(x);
}

// [[Rcpp::export]]
double median_int(Rcpp::IntegerVector x) {
    return median(x);
}

// [[Rcpp::export]]
Rcomplex median_cx(Rcpp::ComplexVector x) {
    return median(x);
}

// [[Rcpp::export]]
Rcpp::String median_ch(Rcpp::CharacterVector x) {
    return median(x);
}