distance.cpp

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(openmp)]]
// [[Rcpp::plugins(cpp11)]]
using namespace Rcpp;


// [[Rcpp::export]]
// Velocity.R colEuclid function
arma::mat  colEuclid(const arma::mat& e, const arma::mat& p, int nthreads=1) {
  arma::mat rm(e.n_cols,e.n_cols);

#pragma omp parallel for shared(rm) num_threads(nthreads)
  for(int i=0;i<e.n_cols;i++) {
    arma::mat t(e); t.each_col() -= p.col(i);
    t%=t;
    arma::rowvec v=sqrt(sum(t,0));
    rm.col(i)=trans(v);
  }
  return(rm);
}

//' Pair-wise euclidean distance among two matrices.
//'
//' @param Ar A matrix with samples on rows.
//' @param Br A matrix with samples on rows.
//' @return A distance matrix in a shape of nrow(\code{Ar}) x nrow(\code{Br}).
//' @references \url{http://blog.felixriedel.com/2013/05/pairwise-distances-in-r/}
//' @examples
//' x <- matrix(1:12, 3)
//' all.equal(c(as.matrix(dist(x))), c(dist_euclidean(x, x)))
// [[Rcpp::export]]
NumericMatrix dist_euclidean(NumericMatrix Ar, NumericMatrix Br) {

  int m = Ar.nrow(),
    n = Br.nrow(),
    k = Ar.ncol();
  arma::mat A = arma::mat(Ar.begin(), m, k, false);
  arma::mat B = arma::mat(Br.begin(), n, k, false);

  arma::colvec An =  sum(square(A), 1);  // rowSum
  arma::colvec Bn =  sum(square(B), 1);

  arma::mat C = -2 * (A * B.t());
  C.each_col() += An; checkUserInterrupt();
  C.each_row() += Bn.t(); checkUserInterrupt();

  return wrap(sqrt(C));
}