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