dist_test (c++ part)

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
#include <string>
using namespace arma;

void chk_mat(const mat& x, const std::string& varName, const std::string& type){
  if (!is_finite(x))
    Rcpp::stop(varName + " must be numerical.\n");
}

arma::vec quantileCpp(const arma::vec& x, const arma::vec& probs){
  chk_mat(x, "x", "double");
  chk_mat(probs, "probs", "double");
  if (any(probs < 0) || any(probs > 1))
    Rcpp::stop("There is values in probs outside [0, 1].");

  vec x_sort = sort(x);
  vec index = (x.n_elem - 1) * probs;
  uvec lo = conv_to<uvec>::from(floor(index)), hi = conv_to<uvec>::from(ceil(index));
  vec h = (index - conv_to<vec>::from(lo));
  vec qs = (1 - h) % x_sort.elem(lo) + h % x_sort.elem(hi);
  return qs;
}

struct Worker_bs_stat: public RcppParallel::Worker {
  const vec& X;
  const vec& Y;
  const vec& probs;
  mat& bsStat;

  Worker_bs_stat(const vec& X, const vec& Y, const vec& probs, mat& bsStat):
    X(X), Y(Y), probs(probs), bsStat(bsStat) {}

  void operator()(std::size_t begin, std::size_t end)
  {
    for (uword i = begin; i < end; ++i)
    {
      uvec indx_x = randi<uvec>(X.n_elem, distr_param(0, X.n_elem-1));
      uvec indx_y = randi<uvec>(Y.n_elem, distr_param(0, Y.n_elem-1));
      bsStat.col(i) = quantileCpp(X.elem(indx_x), probs) - quantileCpp(Y.elem(indx_y), probs);
    }
  }
};

// [[Rcpp::export]]
arma::vec dist_test(const arma::vec& X, const arma::vec& Y, const arma::vec& probs, const double& bsSize){
  // check data
  chk_mat(X, "X", "double");
  chk_mat(Y, "Y", "double");
  chk_mat(probs, "probs", "double");

  if (bsSize <= 0 || std::abs(bsSize - floor(bsSize)) > 1e-6)
    Rcpp::stop("bsSize must be a positive number.\n");

  mat bsDiff = zeros<mat>(probs.n_elem, bsSize);
  Worker_bs_stat bs_stat(X, Y, probs, bsDiff);
  RcppParallel::parallelFor(0, bsSize, bs_stat);
  vec stat = (quantileCpp(X, probs) - quantileCpp(Y, probs)) / stddev(bsDiff, 0, 1);
  return stat;
}