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 //' Compute medoids indices
  14 //'
  15 //' @param pMedoids External pointer
  16 //' @param ref_series reference series
  17 //'
  18 //' @return A map serie number -> medoid index
  19 // [[Rcpp::export]]
  20 IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix ref_series)
  21 {
  22     XPtr<BigMatrix> pMed(pMedoids);
  23     MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);
  24     int nb_series = ref_series.nrow(),
  25         K = pMed->nrow(),
  26         L = pMed->ncol();
  27     IntegerVector mi(nb_series);
  28
  29     for (int i=0; i<nb_series ; i++)
  30     {
  31         // mi[i] <- which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
  32         mi[i] = 0;
  33         double best_dist = DBL_MAX;
  34         for (int j=0; j<K; j++)
  35         {
  36             double dist_ij = 0.;
  37             for (int k=0; k<L; k++)
  38             {
  39                 double delta = ref_series(i,k) - *(medoids[k]+j);
  40                 dist_ij += delta * delta;
  41             }
  42             if (dist_ij < best_dist)
  43             {
  44                 mi[i] = j+1; //R indices start at 1
  45                 best_dist = dist_ij;
  46             }
  47         }
  48     }
  49
  50     return mi;
  51 }