//' computeMedoidsIndices
//'
-//' Compute medoids indices
+//' For each column of the 'series' matrix input, search for the closest medoid
+//' (euclidian distance) and store its index
//'
-//' @param pMedoids External pointer
-//' @param ref_series reference series
+//' @param pMedoids External pointer (a big.matrix 'address' slot in R)
+//' @param series (reference) series, a matrix of size Lxn
//'
-//' @return A map serie number -> medoid index
+//' @return An integer vector of the closest medoids indices, for each (column) serie
// [[Rcpp::export]]
-IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix ref_series)
+IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix series)
{
+ // Turn SEXP external pointer into BigMatrix (description) object
XPtr<BigMatrix> pMed(pMedoids);
+ // medoids: access to the content of the BigMatrix object
MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);
- int nb_series = ref_series.nrow(),
- K = pMed->nrow(),
- L = pMed->ncol();
+
+ int nb_series = series.ncol(),
+ K = pMed->ncol(),
+ L = pMed->nrow();
IntegerVector mi(nb_series);
- for (int i=0; i<nb_series ; i++)
+ for (int i=0; i<nb_series ; i++) //column index in series
{
- // mi[i] <- which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
+ // In R: mi[i] <- which.min( rowSums( sweep(medoids, 2, series[i,], '-')^2 ) )
+ // In C(++), computations must be unrolled
mi[i] = 0;
double best_dist = DBL_MAX;
- for (int j=0; j<K; j++)
+ for (int j=0; j<K; j++) //column index in medoids
{
double dist_ij = 0.;
- for (int k=0; k<L; k++)
+ for (int k=0; k<L; k++) //common row index (medoids, series)
{
- double delta = ref_series(i,k) - *(medoids[k]+j);
+ // Accessing values for a big matrix is a bit weird; see Rcpp gallery/bigmemory
+ double delta = series(k,i) - *(medoids[j]+k);
dist_ij += delta * delta;
}
if (dist_ij < best_dist)