improvements
authorBenjamin Auder <benjamin.auder@somewhere>
Wed, 8 Mar 2017 20:50:53 +0000 (21:50 +0100)
committerBenjamin Auder <benjamin.auder@somewhere>
Wed, 8 Mar 2017 20:50:53 +0000 (21:50 +0100)
.gitignore
epclust/DESCRIPTION
epclust/R/a_NAMESPACE.R
epclust/R/clustering.R
epclust/R/main.R
epclust/src/WER.c [deleted file]
epclust/src/computeMedoidsIndices.c [new file with mode: 0644]
epclust/src/filter.c [new file with mode: 0644]
epclust/tests/testthat/test.clustering.R

index af4c22c..3a326ea 100644 (file)
@@ -32,3 +32,7 @@
 
 #ignore jupyter generated file (HTML vignette, and reports)
 *.ipynb.html
+
+#ignore object files
+*.o
+*.so
index 304fdff..1f2b5ea 100644 (file)
@@ -30,8 +30,9 @@ Suggests:
     DBI
 License: MIT + file LICENSE
 RoxygenNote: 6.0.1
-Collate:
+Collate: 
     'main.R'
     'clustering.R'
     'de_serialize.R'
     'a_NAMESPACE.R'
+    'plot.R'
index 89aa453..e9aa830 100644 (file)
@@ -2,12 +2,13 @@
 #' @include clustering.R
 #' @include main.R
 #'
+#' @useDynLib epclust
+#'
 #' @importFrom Rwave cwt
 #' @importFrom cluster pam
 #' @importFrom parallel makeCluster clusterExport parLapply stopCluster
 #' @importFrom wavelets dwt wt.filter
-#' @importFrom stats filter spline
-#' @importFrom utils tail
+#' @importFrom stats spline
 #' @importFrom methods is
-#' @importFrom bigmemory as.big.matrix is.big.matrix
+#' @importFrom bigmemory big.matrix as.big.matrix is.big.matrix
 NULL
index cda7fbe..92adda2 100644 (file)
@@ -33,15 +33,7 @@ clusteringTask1 = function(
        indices, getContribs, K1, nb_series_per_chunk, ncores_clust=1, verbose=FALSE, parll=TRUE)
 {
        if (verbose)
-               cat(paste("*** Clustering task on ",length(indices)," lines\n", sep=""))
-
-       wrapComputeClusters1 = function(inds) {
-               if (parll)
-                       require("epclust", quietly=TRUE)
-               if (verbose)
-                       cat(paste("   computeClusters1() on ",length(inds)," lines\n", sep=""))
-               inds[ computeClusters1(getContribs(inds), K1) ]
-       }
+               cat(paste("*** Clustering task 1 on ",length(indices)," lines\n", sep=""))
 
        if (parll)
        {
@@ -51,10 +43,20 @@ clusteringTask1 = function(
        while (length(indices) > K1)
        {
                indices_workers = .spreadIndices(indices, nb_series_per_chunk)
-               if (parll)
-                       indices = unlist( parallel::parLapply(cl, indices_workers, wrapComputeClusters1) )
-               else
-                       indices = unlist( lapply(indices_workers, wrapComputeClusters1) )
+               indices <-
+                       if (parll)
+                       {
+                               unlist( parallel::parLapply(cl, indices_workers, function(inds) {
+                                       require("epclust", quietly=TRUE)
+                                       inds[ computeClusters1(getContribs(inds), K1, verbose) ]
+                               }) )
+                       }
+                       else
+                       {
+                               unlist( lapply(indices_workers, function(inds)
+                                       inds[ computeClusters1(getContribs(inds), K1, verbose) ]
+                               ) )
+                       }
        }
        if (parll)
                parallel::stopCluster(cl)
@@ -67,27 +69,35 @@ clusteringTask1 = function(
 clusteringTask2 = function(medoids, K2,
        getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE)
 {
+       if (verbose)
+               cat(paste("*** Clustering task 2 on ",nrow(medoids)," lines\n", sep=""))
+
        if (nrow(medoids) <= K2)
                return (medoids)
        synchrones = computeSynchrones(medoids,
                getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust, verbose, parll)
        distances = computeWerDists(synchrones, ncores_clust, verbose, parll)
        # PAM in package 'cluster' cannot take big.matrix in input: need to cast it
-       mat_dists = matrix(nrow=K1, ncol=K1)
-       for (i in seq_len(K1))
-               mat_dists[i,] = distances[i,]
-       medoids[ computeClusters2(mat_dists,K2), ]
+       medoids[ computeClusters2(distances[,],K2,verbose), ]
 }
 
 #' @rdname clustering
 #' @export
-computeClusters1 = function(contribs, K1)
+computeClusters1 = function(contribs, K1, verbose=FALSE)
+{
+       if (verbose)
+               cat(paste("   computeClusters1() on ",nrow(contribs)," lines\n", sep=""))
        cluster::pam(contribs, K1, diss=FALSE)$id.med
+}
 
 #' @rdname clustering
 #' @export
-computeClusters2 = function(distances, K2)
+computeClusters2 = function(distances, K2, verbose=FALSE)
+{
+       if (verbose)
+               cat(paste("   computeClusters2() on ",nrow(distances)," lines\n", sep=""))
        cluster::pam(distances, K2, diss=TRUE)$id.med
+}
 
 #' computeSynchrones
 #'
@@ -106,29 +116,41 @@ computeClusters2 = function(distances, K2)
 computeSynchrones = function(medoids, getRefSeries,
        nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE)
 {
+       if (verbose)
+               cat(paste("--- Compute synchrones\n", sep=""))
+
        computeSynchronesChunk = function(indices)
        {
-               if (verbose)
-                       cat(paste("--- Compute synchrones for ",length(indices)," lines\n", sep=""))
                ref_series = getRefSeries(indices)
+               nb_series = nrow(ref_series)
                #get medoids indices for this chunk of series
-               for (i in seq_len(nrow(ref_series)))
+
+               #TODO: debug this (address is OK but values are garbage: why?)
+#               mi = .Call("computeMedoidsIndices", medoids@address, ref_series, PACKAGE="epclust")
+
+               #R-equivalent, requiring a matrix (thus potentially breaking "fit-in-memory" hope)
+               mat_meds = medoids[,]
+               mi = rep(NA,nb_series)
+               for (i in 1:nb_series)
+                       mi[i] <- which.min( rowSums( sweep(mat_meds, 2, ref_series[i,], '-')^2 ) )
+               rm(mat_meds); gc()
+
+               for (i in seq_len(nb_series))
                {
-                       j = which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
                        if (parll)
                                synchronicity::lock(m)
-                       synchrones[j,] = synchrones[j,] + ref_series[i,]
-                       counts[j,1] = counts[j,1] + 1
+                       synchrones[mi[i],] = synchrones[mi[i],] + ref_series[i,]
+                       counts[mi[i],1] = counts[mi[i],1] + 1
                        if (parll)
                                synchronicity::unlock(m)
                }
        }
 
-       K = nrow(medoids)
+       K = nrow(medoids) ; L = ncol(medoids)
        # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in //
        # TODO: if size > RAM (not our case), use file-backed big.matrix
-       synchrones = bigmemory::big.matrix(nrow=K,ncol=ncol(medoids),type="double",init=0.)
-       counts = bigmemory::big.matrix(nrow=K,ncol=1,type="double",init=0)
+       synchrones = bigmemory::big.matrix(nrow=K, ncol=L, type="double", init=0.)
+       counts = bigmemory::big.matrix(nrow=K, ncol=1, type="double", init=0)
        # synchronicity is only for Linux & MacOS; on Windows: run sequentially
        parll = (requireNamespace("synchronicity",quietly=TRUE)
                && parll && Sys.info()['sysname'] != "Windows")
@@ -144,9 +166,10 @@ computeSynchrones = function(medoids, getRefSeries,
        }
 
        indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk)
+               browser()
        ignored <-
                if (parll)
-                       parallel::parLapply(indices_workers, computeSynchronesChunk)
+                       parallel::parLapply(cl, indices_workers, computeSynchronesChunk)
                else
                        lapply(indices_workers, computeSynchronesChunk)
 
@@ -179,11 +202,8 @@ computeSynchrones = function(medoids, getRefSeries,
 #' @export
 computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE)
 {
-
-
-
-#TODO: re-organize to call computeWerDist(x,y) [C] (in //?) from two indices + big.matrix
-
+       if (verbose)
+               cat(paste("--- Compute WER dists\n", sep=""))
 
        n <- nrow(synchrones)
        delta <- ncol(synchrones)
@@ -201,10 +221,20 @@ computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE)
        s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 )
        totnoct = noctave + as.integer(s0log/nvoice) + 1
 
+       Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
+       fcoefs = rep(1/3, 3) #moving average on 3 values
+
+       # Generate n(n-1)/2 pairs for WER distances computations
+       pairs = list()
+       V = seq_len(n)
+       for (i in 1:n)
+       {
+               V = V[-1]
+               pairs = c(pairs, lapply(V, function(v) c(i,v)))
+       }
+
        computeCWT = function(i)
        {
-               if (verbose)
-                       cat(paste("+++ Compute Rwave::cwt() on serie ",i,"\n", sep=""))
                ts <- scale(ts(synchrones[i,]), center=TRUE, scale=scaled)
                totts.cwt = Rwave::cwt(ts, totnoct, nvoice, w0, plot=FALSE)
                ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
@@ -214,6 +244,22 @@ computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE)
                sqres / max(Mod(sqres))
        }
 
+       computeDistancesIJ = function(pair)
+       {
+               i = pair[1] ; j = pair[2]
+               if (verbose && j==i+1)
+                       cat(paste("   Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep=""))
+               cwt_i = computeCWT(i)
+               cwt_j = computeCWT(j)
+               num <- .Call("filter", Mod(cwt_i * Conj(cwt_j)), PACKAGE="epclust")
+               WX  <- .Call("filter", Mod(cwt_i * Conj(cwt_i)), PACKAGE="epclust")
+               WY  <- .Call("filter", Mod(cwt_j * Conj(cwt_j)), PACKAGE="epclust")
+               wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY))
+               Xwer_dist[i,j] <- sqrt(delta * ncol(cwt_i) * (1 - wer2))
+               Xwer_dist[j,i] <- Xwer_dist[i,j]
+               Xwer_dist[i,i] = 0.
+       }
+
        if (parll)
        {
                cl = parallel::makeCluster(ncores_clust)
@@ -222,59 +268,15 @@ computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE)
                        envir=environment())
        }
 
-       # list of CWT from synchrones
-       # TODO: fit in RAM, OK? If not, 2 options: serialize, compute individual distances
-       Xcwt4 <-
+       ignored <-
                if (parll)
-                       parallel::parLapply(cl, seq_len(n), computeCWT)
+                       parallel::parLapply(cl, pairs, computeDistancesIJ)
                else
-                       lapply(seq_len(n), computeCWT)
+                       lapply(pairs, computeDistancesIJ)
 
        if (parll)
                parallel::stopCluster(cl)
 
-       Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
-       fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!)
-       if (verbose)
-               cat("*** Compute WER distances from CWT\n")
-
-       #TODO: computeDistances(i,j), et répartir les n(n-1)/2 couples d'indices
-       #là c'est trop déséquilibré
-
-       computeDistancesLineI = function(i)
-       {
-               if (verbose)
-                       cat(paste("   Line ",i,"\n", sep=""))
-               for (j in (i+1):n)
-               {
-                       #TODO: 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
-                       num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
-                       WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)
-                       WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
-                       wer2 <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) )
-                       if (parll)
-                               synchronicity::lock(m)
-                       Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2))
-                       Xwer_dist[j,i] <- Xwer_dist[i,j]
-                       if (parll)
-                               synchronicity::unlock(m)
-               }
-               Xwer_dist[i,i] = 0.
-       }
-
-       parll = (requireNamespace("synchronicity",quietly=TRUE)
-               && parll && Sys.info()['sysname'] != "Windows")
-       if (parll)
-               m <- synchronicity::boost.mutex()
-
-       ignored <-
-               if (parll)
-               {
-                       parallel::mclapply(seq_len(n-1), computeDistancesLineI,
-                               mc.cores=ncores_clust, mc.allow.recursive=FALSE)
-               }
-               else
-                       lapply(seq_len(n-1), computeDistancesLineI)
        Xwer_dist[n,n] = 0.
        Xwer_dist
 }
index 2037dbe..977e61b 100644 (file)
@@ -186,7 +186,12 @@ claws = function(getSeries, K1, K2,
                indices_all[((i-1)*nb_series_per_task+1):upper_bound]
        })
        if (verbose)
-               cat(paste("...Run ",ntasks," x stage 1 in parallel\n",sep=""))
+       {
+               message = paste("...Run ",ntasks," x stage 1", sep="")
+               if (WER=="mix")
+                       message = paste(message," + stage 2", sep="")
+               cat(paste(message,"\n", sep=""))
+       }
        if (WER=="mix")
                {synchrones_file = paste(bin_dir,"synchrones",sep="") ; unlink(synchrones_file)}
        if (parll && ntasks>1)
@@ -229,7 +234,7 @@ claws = function(getSeries, K1, K2,
        indices_medoids = clusteringTask1(
                indices, getContribs, K1, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose, parll)
        medoids1 = bigmemory::as.big.matrix( getSeries(indices_medoids) )
-       medoids2 = computeClusters2(medoids1, K2,
+       medoids2 = clusteringTask2(medoids1, K2,
                getRefSeries, nb_curves, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose, parll)
 
        # Cleanup
diff --git a/epclust/src/WER.c b/epclust/src/WER.c
deleted file mode 100644 (file)
index 36bfba7..0000000
+++ /dev/null
@@ -1,117 +0,0 @@
-#include <stdlib.h>
-#include <math.h>
-#include <stdbool.h>
-
-#ifndef M_PI
-#define M_PI 3.14159265358979323846
-#endif
-
-// n: number of synchrones, m: length of a synchrone
-float computeWerDist(float* s1, float* s2, int n, int m)
-{
-       //TODO: automatic tune of all these parameters ? (for other users)
-       int nvoice = 4;
-       //noctave 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones))
-       int noctave = 13
-       // 4 here represent 2^5 = 32 half-hours ~ 1 day
-       //NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?)
-       //R: scalevector <- 2^(4:(noctave * nvoice) / nvoice + 1)
-       int* scalevector = (int*)malloc( (noctave*nvoice-4 + 1) * sizeof(int))
-       for (int i=4; i<=noctave*nvoice; i++)
-               scalevector[i-4] = pow(2., (float)i/nvoice + 1.);
-       //condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1
-       int s0 = 2;
-       double w0 = 2*M_PI;
-       bool scaled = false;
-       int s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 )
-       int totnoct = noctave + as.integer(s0log/nvoice) + 1
-
-
-
-
-
-///TODO: continue
-
-
-
-       computeCWT = function(i)
-       {
-               if (verbose)
-                       cat(paste("+++ Compute Rwave::cwt() on serie ",i,"\n", sep=""))
-               ts <- scale(ts(synchrones[i,]), center=TRUE, scale=scaled)
-               totts.cwt = Rwave::cwt(ts,totnoct,nvoice,w0,plot=0)
-               ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
-               #Normalization
-               sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
-               sqres <- sweep(ts.cwt,2,sqs,'*')
-               sqres / max(Mod(sqres))
-       }
-
-       if (parll)
-       {
-               cl = parallel::makeCluster(ncores_clust)
-               parallel::clusterExport(cl,
-                       varlist=c("synchrones","totnoct","nvoice","w0","s0log","noctave","s0","verbose"),
-                       envir=environment())
-       }
-
-       # (normalized) observations node with CWT
-       Xcwt4 <-
-               if (parll)
-                       parallel::parLapply(cl, seq_len(n), computeCWT)
-               else
-                       lapply(seq_len(n), computeCWT)
-
-       if (parll)
-               parallel::stopCluster(cl)
-
-       Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
-       fcoefs = rep(1/3, 3) #moving average on 3 values (TODO: very slow! correct?!)
-       if (verbose)
-               cat("*** Compute WER distances from CWT\n")
-
-       #TODO: computeDistances(i,j), et répartir les n(n-1)/2 couples d'indices
-       #là c'est trop déséquilibré
-
-       computeDistancesLineI = function(i)
-       {
-               if (verbose)
-                       cat(paste("   Line ",i,"\n", sep=""))
-               for (j in (i+1):n)
-               {
-                       #TODO: 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C
-                       num <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
-                       WX <- filter(Mod(Xcwt4[[i]] * Conj(Xcwt4[[i]])), fcoefs, circular=TRUE)
-                       WY <- filter(Mod(Xcwt4[[j]] * Conj(Xcwt4[[j]])), fcoefs, circular=TRUE)
-                       wer2    <- sum(colSums(num)^2) / sum( sum(colSums(WX) * colSums(WY)) )
-                       if (parll)
-                               synchronicity::lock(m)
-                       Xwer_dist[i,j] <- sqrt(delta * ncol(Xcwt4[[1]]) * (1 - wer2))
-                       Xwer_dist[j,i] <- Xwer_dist[i,j]
-                       if (parll)
-                               synchronicity::unlock(m)
-               }
-               Xwer_dist[i,i] = 0.
-       }
-
-       parll = (requireNamespace("synchronicity",quietly=TRUE)
-               && parll && Sys.info()['sysname'] != "Windows")
-       if (parll)
-               m <- synchronicity::boost.mutex()
-
-       ignored <-
-               if (parll)
-               {
-                       parallel::mclapply(seq_len(n-1), computeDistancesLineI,
-                               mc.cores=ncores_clust, mc.allow.recursive=FALSE)
-               }
-               else
-                       lapply(seq_len(n-1), computeDistancesLineI)
-       Xwer_dist[n,n] = 0.
-
-       mat_dists = matrix(nrow=n, ncol=n)
-       #TODO: avoid this loop?
-       for (i in 1:n)
-               mat_dists[i,] = Xwer_dist[i,]
-       mat_dists
-
diff --git a/epclust/src/computeMedoidsIndices.c b/epclust/src/computeMedoidsIndices.c
new file mode 100644 (file)
index 0000000..98a0111
--- /dev/null
@@ -0,0 +1,51 @@
+#include <stdlib.h>
+#include <float.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+
+#include <stdio.h>
+
+// (K,L): dim(medoids)
+// mi: medoids indices
+SEXP computeMedoidsIndices(SEXP medoids_, SEXP ref_series_)
+{
+       double *medoids = (double*) R_ExternalPtrAddr(medoids_),
+               *ref_series = REAL(ref_series_);
+       int nb_series = INTEGER(getAttrib(ref_series_, R_DimSymbol))[0],
+               K = INTEGER(getAttrib(medoids_, R_DimSymbol))[0],
+               L = INTEGER(getAttrib(ref_series_, R_DimSymbol))[1],
+               *mi = (int*)malloc(nb_series*sizeof(int));
+
+       //TODO: debug this: medoids have same addresses on both sides, but this side fails
+       printf("MED: %x\n",medoids);
+
+       for (int i=0; i<nb_series ; i++)
+       {
+               // mi[i] <- which.min( rowSums( sweep(medoids, 2, ref_series[i,], '-')^2 ) )
+               mi[i] = 0;
+               double best_dist = DBL_MAX;
+               for (int j=0; j<K; j++)
+               {
+                       double dist_ij = 0.;
+                       for (int k=0; k<L; k++)
+                       {
+                               double delta = ref_series[k*K+i] - medoids[k*K+j];
+                               dist_ij += delta * delta;
+                       }
+                       if (dist_ij < best_dist)
+                       {
+                               mi[i] = j+1; //R indices start at 1
+                               best_dist = dist_ij;
+                       }
+               }
+       }
+
+       SEXP R_mi;
+       PROTECT(R_mi = allocVector(INTSXP, nb_series));
+       for (int i = 0; i < nb_series; i++)
+               INTEGER(R_mi)[i] = mi[i];
+       UNPROTECT(1);
+       free(mi);
+       return R_mi;
+}
diff --git a/epclust/src/filter.c b/epclust/src/filter.c
new file mode 100644 (file)
index 0000000..382910a
--- /dev/null
@@ -0,0 +1,37 @@
+#include <stdlib.h>
+#include <R.h>
+#include <Rinternals.h>
+#include <Rmath.h>
+
+#include <stdio.h>
+
+SEXP filter(SEXP cwt_)
+{
+       int L = INTEGER(getAttrib(cwt_, R_DimSymbol))[0],
+               D = INTEGER(getAttrib(cwt_, R_DimSymbol))[1];
+       double *cwt = REAL(cwt_);
+       SEXP R_fcwt;
+       PROTECT(R_fcwt = allocMatrix(REALSXP, L, D));
+       double* fcwt = REAL(R_fcwt);
+
+       //TODO: coding style is terrible... no time for now.
+       for (int col=0; col<D; col++)
+       {
+               double v1 = cwt[0];
+               double ma = v1 + cwt[1] + cwt[2];
+               for (int i=1; i<L-2; i++)
+               {
+                       fcwt[i] = ma / 3.;
+                       ma = ma - v1 + cwt[i+2];
+                       v1 = cwt[i];
+               }
+               fcwt[0] = cwt[0];
+               fcwt[L-2] = ma / 3.;
+               fcwt[L-1] = cwt[L-1];
+               cwt = cwt + L;
+               fcwt = fcwt + L;
+       }
+
+       UNPROTECT(1);
+       return R_fcwt;
+}
index 49afe60..6c94f92 100644 (file)
@@ -4,7 +4,7 @@ context("clustering")
 I = function(i, base)
        (i-1) %% base + 1
 
-test_that("computeClusters1 behave as expected",
+test_that("computeClusters1&2 behave as expected",
 {
        require("MASS", quietly=TRUE)
        if (!require("clue", quietly=TRUE))
@@ -21,27 +21,38 @@ test_that("computeClusters1 behave as expected",
                Id = diag(d)
                coefs = do.call(rbind,
                        lapply(1:K, function(i) MASS::mvrnorm(cs, c(rep(0,(i-1)),5,rep(0,d-i)), Id)))
-               indices_medoids = computeClusters1(coefs, K)
+               indices_medoids1 = computeClusters1(coefs, K, verbose=TRUE)
+               indices_medoids2 = computeClusters2(dist(coefs), K, verbose=TRUE)
                # Get coefs assignments (to medoids)
-               assignment = sapply(seq_len(n), function(i)
-                       which.min( rowSums( sweep(coefs[indices_medoids,],2,coefs[i,],'-')^2 ) ) )
+               assignment1 = sapply(seq_len(n), function(i)
+                       which.min( rowSums( sweep(coefs[indices_medoids1,],2,coefs[i,],'-')^2 ) ) )
+               assignment2 = sapply(seq_len(n), function(i)
+                       which.min( rowSums( sweep(coefs[indices_medoids2,],2,coefs[i,],'-')^2 ) ) )
                for (i in 1:K)
-                       expect_equal(sum(assignment==i), cs, tolerance=5)
+               {
+                       expect_equal(sum(assignment1==i), cs, tolerance=5)
+                       expect_equal(sum(assignment2==i), cs, tolerance=5)
+               }
 
-               costs_matrix = matrix(nrow=K,ncol=K)
+               costs_matrix1 = matrix(nrow=K,ncol=K)
+               costs_matrix2 = matrix(nrow=K,ncol=K)
                for (i in 1:K)
                {
                        for (j in 1:K)
                        {
                                # assign i (in result) to j (order 1,2,3)
-                               costs_matrix[i,j] = abs( mean(assignment[((i-1)*cs+1):(i*cs)]) - j )
+                               costs_matrix1[i,j] = abs( mean(assignment1[((i-1)*cs+1):(i*cs)]) - j )
+                               costs_matrix2[i,j] = abs( mean(assignment2[((i-1)*cs+1):(i*cs)]) - j )
                        }
                }
-               permutation = as.integer( clue::solve_LSAP(costs_matrix) )
+               permutation1 = as.integer( clue::solve_LSAP(costs_matrix1) )
+               permutation2 = as.integer( clue::solve_LSAP(costs_matrix2) )
                for (i in 1:K)
                {
                        expect_equal(
-                               mean(assignment[((i-1)*cs+1):(i*cs)]), permutation[i], tolerance=0.05)
+                               mean(assignment1[((i-1)*cs+1):(i*cs)]), permutation1[i], tolerance=0.05)
+                       expect_equal(
+                               mean(assignment2[((i-1)*cs+1):(i*cs)]), permutation2[i], tolerance=0.05)
                }
        }
 })
@@ -66,76 +77,75 @@ test_that("computeSynchrones behave as expected",
                indices = indices[indices <= n]
                if (length(indices)>0) series[indices,] else NULL
        }
-       synchrones = computeSynchrones(rbind(s1,s2,s3), getRefSeries, n, 100,
-               verbose=TRUE, parll=FALSE)
+       synchrones = computeSynchrones(bigmemory::as.big.matrix(rbind(s1,s2,s3)), getRefSeries,
+               n, 100, verbose=TRUE, parll=FALSE)
 
        expect_equal(dim(synchrones), c(K,L))
        for (i in 1:K)
                expect_equal(synchrones[i,], s[[i]], tolerance=0.01)
 })
 
+# NOTE: medoids can be a big.matrix
 computeDistortion = function(series, medoids)
 {
        n = nrow(series) ; L = ncol(series)
        distortion = 0.
+       if (bigmemory::is.big.matrix(medoids))
+               medoids = medoids[,]
        for (i in seq_len(n))
                distortion = distortion + min( rowSums( sweep(medoids,2,series[i,],'-')^2 ) / L )
        distortion / n
 }
 
-test_that("computeClusters2 behave as expected",
+test_that("clusteringTask1 behave as expected",
 {
        n = 900
        x = seq(0,9.5,0.1)
        L = length(x) #96 1/4h
        K1 = 60
-       K2 = 3
-       #for (i in 1:60) {plot(x^(1+i/30)*cos(x+i),type="l",col=i,ylim=c(-50,50)); par(new=TRUE)}
        s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
        series = matrix(nrow=n, ncol=L)
        for (i in seq_len(n))
                series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
-       getRefSeries = function(indices) {
+       getSeries = function(indices) {
                indices = indices[indices <= n]
                if (length(indices)>0) series[indices,] else NULL
        }
-       # Artificially simulate 60 medoids - perfect situation, all equal to one of the refs
-       medoids_K1 = do.call(rbind, lapply( 1:K1, function(i) s[[I(i,K1)]] ) )
-       medoids_K2 = computeClusters2(medoids_K1, K2, getRefSeries, n, 75,
-               verbose=TRUE, parll=FALSE)
+       wf = "haar"
+       ctype = "absolute"
+       getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
+       indices1 = clusteringTask1(1:n, getContribs, K1, 75, verbose=TRUE, parll=FALSE)
+       medoids_K1 = getSeries(indices1)
 
-       expect_equal(dim(medoids_K2), c(K2,L))
+       expect_equal(dim(medoids_K1), c(K1,L))
        # Not easy to evaluate result: at least we expect it to be better than random selection of
-       # medoids within 1...K1 (among references)
-       distorGood = computeDistortion(series, medoids_K2)
+       # medoids within initial series
+       distorGood = computeDistortion(series, medoids_K1)
        for (i in 1:3)
-               expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
+               expect_lte( distorGood, computeDistortion(series,series[sample(1:n, K1),]) )
 })
 
-test_that("clusteringTask1 + computeClusters2 behave as expected",
+test_that("clusteringTask2 behave as expected",
 {
        n = 900
        x = seq(0,9.5,0.1)
        L = length(x) #96 1/4h
        K1 = 60
        K2 = 3
+       #for (i in 1:60) {plot(x^(1+i/30)*cos(x+i),type="l",col=i,ylim=c(-50,50)); par(new=TRUE)}
        s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
        series = matrix(nrow=n, ncol=L)
        for (i in seq_len(n))
                series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
-       getSeries = function(indices) {
+       getRefSeries = function(indices) {
                indices = indices[indices <= n]
                if (length(indices)>0) series[indices,] else NULL
        }
-       wf = "haar"
-       ctype = "absolute"
-       getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
-       medoids_K1 = getSeries( clusteringTask1(1:n, getContribs, K1, 75,
-               verbose=TRUE, parll=FALSE) )
-       medoids_K2 = computeClusters2(medoids_K1, K2, getSeries, n, 120,
-               verbose=TRUE, parll=FALSE)
+       # Artificially simulate 60 medoids - perfect situation, all equal to one of the refs
+       medoids_K1 = bigmemory::as.big.matrix(
+               do.call(rbind, lapply( 1:K1, function(i) s[[I(i,K1)]] ) ) )
+       medoids_K2 = clusteringTask2(medoids_K1, K2, getRefSeries, n, 75, verbose=TRUE, parll=FALSE)
 
-       expect_equal(dim(medoids_K1), c(K1,L))
        expect_equal(dim(medoids_K2), c(K2,L))
        # Not easy to evaluate result: at least we expect it to be better than random selection of
        # medoids within 1...K1 (among references)
@@ -143,3 +153,36 @@ test_that("clusteringTask1 + computeClusters2 behave as expected",
        for (i in 1:3)
                expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
 })
+
+#NOTE: rather redundant test
+#test_that("clusteringTask1 + clusteringTask2 behave as expected",
+#{
+#      n = 900
+#      x = seq(0,9.5,0.1)
+#      L = length(x) #96 1/4h
+#      K1 = 60
+#      K2 = 3
+#      s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
+#      series = matrix(nrow=n, ncol=L)
+#      for (i in seq_len(n))
+#              series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
+#      getSeries = function(indices) {
+#              indices = indices[indices <= n]
+#              if (length(indices)>0) series[indices,] else NULL
+#      }
+#      wf = "haar"
+#      ctype = "absolute"
+#      getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype)
+#      require("bigmemory", quietly=TRUE)
+#      indices1 = clusteringTask1(1:n, getContribs, K1, 75, verbose=TRUE, parll=FALSE)
+#      medoids_K1 = bigmemory::as.big.matrix( getSeries(indices1) )
+#      medoids_K2 = clusteringTask2(medoids_K1, K2, getSeries, n, 120, verbose=TRUE, parll=FALSE)
+#
+#      expect_equal(dim(medoids_K1), c(K1,L))
+#      expect_equal(dim(medoids_K2), c(K2,L))
+#      # Not easy to evaluate result: at least we expect it to be better than random selection of
+#      # medoids within 1...K1 (among references)
+#      distorGood = computeDistortion(series, medoids_K2)
+#      for (i in 1:3)
+#              expect_lte( distorGood, computeDistortion(series,medoids_K1[sample(1:K1, K2),]) )
+#})