epclust/src/computeMedoidsIndices.c

   1 #include <stdlib.h>
   2 #include <float.h>
   3 #include <R.h>
   4 #include <Rinternals.h>
   5 #include <Rmath.h>
   6
   7 #include <stdio.h>
   8
   9 // (K,L): dim(medoids)
  10 // mi: medoids indices
  11 SEXP computeMedoidsIndices(SEXP medoids_, SEXP ref_series_)
  12 {
  13     double *medoids = (double*) R_ExternalPtrAddr(medoids_),
  14         *ref_series = REAL(ref_series_);
  15     int nb_series = INTEGER(getAttrib(ref_series_, R_DimSymbol))[0],
  16         K = INTEGER(getAttrib(medoids_, R_DimSymbol))[0],
  17         L = INTEGER(getAttrib(ref_series_, R_DimSymbol))[1],
  18         *mi = (int*)malloc(nb_series*sizeof(int));
  19
  20     //TODO: debug this: medoids have same addresses on both sides, but this side fails
  21     printf("MED: %x\n",medoids);
  22
  23     for (int i=0; i<nb_series ; i++)
  24     {
  25         // mi[i] <- which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
  26         mi[i] = 0;
  27         double best_dist = DBL_MAX;
  28         for (int j=0; j<K; j++)
  29         {
  30             double dist_ij = 0.;
  31             for (int k=0; k<L; k++)
  32             {
  33                 double delta = ref_series[k*K+i] - medoids[k*K+j];
  34                 dist_ij += delta * delta;
  35             }
  36             if (dist_ij < best_dist)
  37             {
  38                 mi[i] = j+1; //R indices start at 1
  39                 best_dist = dist_ij;
  40             }
  41         }
  42     }
  43
  44     SEXP R_mi;
  45     PROTECT(R_mi = allocVector(INTSXP, nb_series));
  46     for (int i = 0; i < nb_series; i++)
  47         INTEGER(R_mi)[i] = mi[i];
  48     UNPROTECT(1);
  49     free(mi);
  50     return R_mi;
  51 }