code seems OK; still wavelets test to write
[epclust.git] / epclust / src / computeMedoidsIndices.cpp
... / ...
CommitLineData
1#include <Rcpp.h>
2
3// [[Rcpp::depends(BH, bigmemory)]]
4#include <bigmemory/MatrixAccessor.hpp>
5
6#include <numeric>
7#include <cfloat>
8
9using namespace Rcpp;
10
11//' computeMedoidsIndices
12//'
13//' For each column of the 'series' matrix input, search for the closest medoid
14//' (euclidian distance) and store its index
15//'
16//' @param pMedoids External pointer (a big.matrix 'address' slot in R)
17//' @param series (reference) series, a matrix of size Lxn
18//'
19//' @return An integer vector of the closest medoids indices, for each (column) serie
20// [[Rcpp::export]]
21IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix series)
22{
23 // Turn SEXP external pointer into BigMatrix (description) object
24 XPtr<BigMatrix> pMed(pMedoids);
25 // medoids: access to the content of the BigMatrix object
26 MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);
27
28 int nb_series = series.ncol(),
29 K = pMed->ncol(),
30 L = pMed->nrow();
31 IntegerVector mi(nb_series);
32
33 for (int i=0; i<nb_series ; i++) //column index in series
34 {
35 // In R: mi[i] <- which.min( rowSums( sweep(medoids, 2, series[i,], '-')^2 ) )
36 // In C(++), computations must be unrolled
37 mi[i] = 0;
38 double best_dist = DBL_MAX;
39 for (int j=0; j<K; j++) //column index in medoids
40 {
41 double dist_ij = 0.;
42 for (int k=0; k<L; k++) //common row index (medoids, series)
43 {
44 // Accessing values for a big matrix is a bit weird; see Rcpp gallery/bigmemory
45 double delta = series(k,i) - *(medoids[j]+k);
46 dist_ij += delta * delta;
47 }
48 if (dist_ij < best_dist)
49 {
50 mi[i] = j+1; //R indices start at 1
51 best_dist = dist_ij;
52 }
53 }
54 }
55
56 return mi;
57}