| 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 | } |