X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=epclust%2Fsrc%2FcomputeMedoidsIndices.cpp;fp=epclust%2Fsrc%2FcomputeMedoidsIndices.cpp;h=693418118f9f3d38d1007c3923ea534c44703ec9;hb=363ae13430cdee6ba76b42b7316aa4b292b04d93;hp=0000000000000000000000000000000000000000;hpb=777c4b0274c059eeb5d4bd784ef773e819a7f7a2;p=epclust.git

diff --git a/epclust/src/computeMedoidsIndices.cpp b/epclust/src/computeMedoidsIndices.cpp
new file mode 100644
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;
+}