Create a gist now

Instantly share code, notes, and snippets.

@Puriney /distance.cpp Secret
Created Jan 21, 2018

What would you like to do?
#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));
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment