improvements

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

diff --git a/epclust/src/computeMedoidsIndices.c b/epclust/src/computeMedoidsIndices.c
new file mode 100644 (file)
index 0000000..98a0111
--- /dev/null
+++ b/epclust/src/computeMedoidsIndices.c
@@ -0,0 +1,51 @@
+#include <stdlib.h>
+#include <float.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+
+#include <stdio.h>
+
+// (K,L): dim(medoids)
+// mi: medoids indices
+SEXP computeMedoidsIndices(SEXP medoids_, SEXP ref_series_)
+{
+       double *medoids = (double*) R_ExternalPtrAddr(medoids_),
+               *ref_series = REAL(ref_series_);
+       int nb_series = INTEGER(getAttrib(ref_series_, R_DimSymbol))[0],
+               K = INTEGER(getAttrib(medoids_, R_DimSymbol))[0],
+               L = INTEGER(getAttrib(ref_series_, R_DimSymbol))[1],
+               *mi = (int*)malloc(nb_series*sizeof(int));
+
+       //TODO: debug this: medoids have same addresses on both sides, but this side fails
+       printf("MED: %x\n",medoids);
+
+       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[k*K+i] - medoids[k*K+j];
+                               dist_ij += delta * delta;
+                       }
+                       if (dist_ij < best_dist)
+                       {
+                               mi[i] = j+1; //R indices start at 1
+                               best_dist = dist_ij;
+                       }
+               }
+       }
+
+       SEXP R_mi;
+       PROTECT(R_mi = allocVector(INTSXP, nb_series));
+       for (int i = 0; i < nb_series; i++)
+               INTEGER(R_mi)[i] = mi[i];
+       UNPROTECT(1);
+       free(mi);
+       return R_mi;
+}