code seems OK; still wavelets test to write

[epclust.git] / epclust / src / computeMedoidsIndices.cpp

#include <Rcpp.h>

// [[Rcpp::depends(BH, bigmemory)]]
#include <bigmemory/MatrixAccessor.hpp>

#include <numeric>
#include <cfloat>

using namespace Rcpp;

//' computeMedoidsIndices
//'
//' For each column of the 'series' matrix input, search for the closest medoid
//' (euclidian distance) and store its index
//'
//' @param pMedoids External pointer (a big.matrix 'address' slot in R)
//' @param series (reference) series, a matrix of size Lxn
//'
//' @return An integer vector of the closest medoids indices, for each (column) serie
// [[Rcpp::export]]
IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix series)
{
	// Turn SEXP external pointer into BigMatrix (description) object
	XPtr<BigMatrix> pMed(pMedoids);
	// medoids: access to the content of the BigMatrix object
	MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);

	int nb_series = series.ncol(),
		K = pMed->ncol(),
		L = pMed->nrow();
	IntegerVector mi(nb_series);

	for (int i=0; i<nb_series ; i++) //column index in series
	{
		// In R: mi[i] <- which.min( rowSums( sweep(medoids, 2, series[i,], '-')^2 ) )
		// In C(++), computations must be unrolled
		mi[i] = 0;
		double best_dist = DBL_MAX;
		for (int j=0; j<K; j++) //column index in medoids
		{
			double dist_ij = 0.;
			for (int k=0; k<L; k++) //common row index (medoids, series)
			{
				// Accessing values for a big matrix is a bit weird; see Rcpp gallery/bigmemory
				double delta = series(k,i) - *(medoids[j]+k);
				dist_ij += delta * delta;
			}
			if (dist_ij < best_dist)
			{
				mi[i] = j+1; //R indices start at 1
				best_dist = dist_ij;
			}
		}
	}

	return mi;
}