use Rcpp; ongoing debug for parallel synchrones computation

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

diff --git a/epclust/src/computeMedoidsIndices.cpp b/epclust/src/computeMedoidsIndices.cpp
new file mode 100644 (file)
index 0000000..6934181
--- /dev/null
+++ b/epclust/src/computeMedoidsIndices.cpp
@@ -0,0 +1,51 @@
+#include <Rcpp.h>
+
+// [[Rcpp::depends(BH, bigmemory)]]
+#include <bigmemory/MatrixAccessor.hpp>
+
+#include <numeric>
+#include <cfloat>
+
+using namespace Rcpp;
+
+//' computeMedoidsIndices
+//'
+//' Compute medoids indices
+//'
+//' @param pMedoids External pointer
+//' @param ref_series reference series
+//'
+//' @return A map serie number -> medoid index
+// [[Rcpp::export]]
+IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix ref_series)
+{
+       XPtr<BigMatrix> pMed(pMedoids);
+       MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);
+       int nb_series = ref_series.nrow(),
+               K = pMed->nrow(),
+               L = pMed->ncol();
+       IntegerVector mi(nb_series);
+
+       for (int i=0; i<nb_series ; i++)
+       {
+               // mi[i] <- which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
+               mi[i] = 0;
+               double best_dist = DBL_MAX;
+               for (int j=0; j<K; j++)
+               {
+                       double dist_ij = 0.;
+                       for (int k=0; k<L; k++)
+                       {
+                               double delta = ref_series(i,k) - *(medoids[k]+j);
+                               dist_ij += delta * delta;
+                       }
+                       if (dist_ij < best_dist)
+                       {
+                               mi[i] = j+1; //R indices start at 1
+                               best_dist = dist_ij;
+                       }
+               }
+       }
+
+       return mi;
+}