WER distances is a regular matrix, fix doc (weird latex error?)
[epclust.git] / epclust / src / computeMedoidsIndices.c
CommitLineData
e161499b
BA
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
11SEXP 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}