Rcpp sd test

// based on http://stackoverflow.com/questions/24253228/performance-of-r-statssd-vs-armastddev-vs-rcpp-implementation

library(Rcpp)
library(RcppArmadillo)
library(RcppParallel)

sourceCpp(code = '
#include <vector>
#include <cmath>
#include <numeric>
// [[Rcpp::plugins(cpp11)]]
//[[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace Rcpp;

struct Sum : public RcppParallel::Worker {
  const RcppParallel::RVector<double> input;
  double value;

  Sum(const Rcpp::NumericVector input)
  : input(input), value(0) {}
  Sum(const Sum& sum, RcppParallel::Split)
  : input(sum.input), value(0) {}

  void operator()(std::size_t begin, std::size_t end) {
    value += std::accumulate(input.begin() + begin,
                             input.begin() + end,
                             0.0);
  }

  void join(const Sum& rhs) {
    value += rhs.value;
  }
};


struct InnerProduct : public RcppParallel::Worker {
  const RcppParallel::RVector<double> x;
  const RcppParallel::RVector<double> y;
  double mean;
  double product;

  InnerProduct(const Rcpp::NumericVector x,
               const Rcpp::NumericVector y,
               const double mean)
  : x(x), y(y), mean(mean), product(0) {}

  InnerProduct(const InnerProduct& innerProduct,
               RcppParallel::Split)
  : x(innerProduct.x), y(innerProduct.y),
    mean(innerProduct.mean), product(0) {}

  void operator()(std::size_t begin, std::size_t end) {
    product += std::inner_product(x.begin() + begin,
                                  x.begin() + end,
                                  y.begin() + begin,
                                  0.0, std::plus<double>(),
                [&](double lhs, double rhs)->double {
                  return ( (lhs-mean)*(rhs-mean) );
                });
  }

  void join(const InnerProduct& rhs) {
    product += rhs.product;
  }
};

// [[Rcpp::export]]
double sd_RP(const Rcpp::NumericVector& x_)
{
  int N = x_.size();
  Rcpp::NumericVector y_(x_);
  Sum sum(x_);
  RcppParallel::parallelReduce(0, x_.length(), sum);

  double mean = sum.value / N;
  InnerProduct innerProduct(x_, y_, mean);
  RcppParallel::parallelReduce(0, x_.length(), innerProduct);

  return std::sqrt( innerProduct.product / (N-1) );
}

// [[Rcpp::export]]
double sd_cpp(Rcpp::NumericVector rinVec) {
  std::vector<double> inVec(rinVec.begin(),rinVec.end());
  int n = inVec.size();
  double sum = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  double mean = sum / inVec.size();

  for(std::vector<double>::iterator iter = inVec.begin();
      iter != inVec.end();
      ++iter){
    double temp = (*iter - mean)*(*iter - mean);
    *iter = temp;
  }

  double sd = std::accumulate(inVec.begin(), inVec.end(), 0.0);
  return std::sqrt( sd / (n-1) );
}

//[[Rcpp::export]]
double sd_arma(arma::vec& xin){
  return arma::stddev(xin);
}
//[[Rcpp::export]]
double sd_arma2(arma::vec xin){
  return arma::stddev(xin);
}
//[[Rcpp::export]]
double sd_arma3(NumericVector xin){
  arma::vec x = arma::vec(xin.begin(), xin.size(), false);
  return arma::stddev(x);
}
//[[Rcpp::export]]
double sd_sugar(NumericVector& xin){
  return sd(xin);
}')

library(microbenchmark)
set.seed(123)
x <- rnorm(6e6)
microbenchmark(
  sd(x),
  sd_cpp(x),
  sd_arma(x),
  sd_arma2(x),
  sd_arma3(x),
  sd_sugar(x),
  sd_RP(x),
  times=300L,
  control=list(warmup=25))

# Unit: milliseconds
#         expr       min        lq      mean   median        uq       max neval
#        sd(x) 22.511043 23.450942 23.765945 23.76388 23.994588 26.869210   300
#    sd_cpp(x) 28.855132 29.302430 29.601853 29.45789 29.753917 31.651868   300
#   sd_arma(x) 22.628647 23.837267 24.102153 24.04160 24.224056 28.056108   300
#  sd_arma2(x) 36.968243 38.796221 39.328432 39.06215 39.688645 42.672182   300
#  sd_arma3(x) 23.078272 23.830440 24.117261 24.03012 24.280841 26.794738   300
#  sd_sugar(x) 14.803495 16.029491 16.218793 16.21459 16.404333 20.748227   300
#     sd_RP(x)  2.164655  2.422824  2.848265  2.83661  3.069336  5.128333   300