-#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;
-}