epclust/src/computeMedoidsIndices.cpp

   1 #include <Rcpp.h>
   2
   3 // [[Rcpp::depends(BH, bigmemory)]]
   4 #include <bigmemory/MatrixAccessor.hpp>
   5
   6 #include <numeric>
   7 #include <cfloat>
   8
   9 using namespace Rcpp;
  10
  11 //' computeMedoidsIndices
  12 //'
  13 //' For each column of the 'series' matrix input, search for the closest medoid
  14 //' (euclidian distance) and store its index
  15 //'
  16 //' @param pMedoids External pointer (a big.matrix 'address' slot in R)
  17 //' @param series (reference) series, a matrix of size Lxn
  18 //'
  19 //' @return An integer vector of the closest medoids indices, for each (column) serie
  20 // [[Rcpp::export]]
  21 IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix series)
  22 {
  23     // Turn SEXP external pointer into BigMatrix (description) object
  24     XPtr<BigMatrix> pMed(pMedoids);
  25     // medoids: access to the content of the BigMatrix object
  26     MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);
  27
  28     int nb_series = series.ncol(),
  29         K = pMed->ncol(),
  30         L = pMed->nrow();
  31     IntegerVector mi(nb_series);
  32
  33     for (int i=0; i<nb_series ; i++) //column index in series
  34     {
  35         // In R: mi[i] <- which.min( rowSums( sweep(medoids, 2, series[i,], '-')^2 ) )
  36         // In C(++), computations must be unrolled
  37         mi[i] = 0;
  38         double best_dist = DBL_MAX;
  39         for (int j=0; j<K; j++) //column index in medoids
  40         {
  41             double dist_ij = 0.;
  42             for (int k=0; k<L; k++) //common row index (medoids, series)
  43             {
  44                 // Accessing values for a big matrix is a bit weird; see Rcpp gallery/bigmemory
  45                 double delta = series(k,i) - *(medoids[j]+k);
  46                 dist_ij += delta * delta;
  47             }
  48             if (dist_ij < best_dist)
  49             {
  50                 mi[i] = j+1; //R indices start at 1
  51                 best_dist = dist_ij;
  52             }
  53         }
  54     }
  55
  56     return mi;
  57 }