From e161499b97c782aadfc287c22b55f85724f86fae Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Wed, 8 Mar 2017 21:50:53 +0100
Subject: [PATCH] improvements

---
 .gitignore                               |   4 +
 epclust/DESCRIPTION                      |   3 +-
 epclust/R/a_NAMESPACE.R                  |   7 +-
 epclust/R/clustering.R                   | 168 ++++++++++++-----------
 epclust/R/main.R                         |   9 +-
 epclust/src/WER.c                        | 117 ----------------
 epclust/src/computeMedoidsIndices.c      |  51 +++++++
 epclust/src/filter.c                     |  37 +++++
 epclust/tests/testthat/test.clustering.R | 109 ++++++++++-----
 9 files changed, 266 insertions(+), 239 deletions(-)
 delete mode 100644 epclust/src/WER.c
 create mode 100644 epclust/src/computeMedoidsIndices.c
 create mode 100644 epclust/src/filter.c

diff --git a/.gitignore b/.gitignore
index af4c22c..3a326ea 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,3 +32,7 @@
 
 #ignore jupyter generated file (HTML vignette, and reports)
 *.ipynb.html
+
+#ignore object files
+*.o
+*.so
diff --git a/epclust/DESCRIPTION b/epclust/DESCRIPTION
index 304fdff..1f2b5ea 100644
--- a/epclust/DESCRIPTION
+++ b/epclust/DESCRIPTION
@@ -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'
diff --git a/epclust/R/a_NAMESPACE.R b/epclust/R/a_NAMESPACE.R
index 89aa453..e9aa830 100644
--- a/epclust/R/a_NAMESPACE.R
+++ b/epclust/R/a_NAMESPACE.R
@@ -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
diff --git a/epclust/R/clustering.R b/epclust/R/clustering.R
index cda7fbe..92adda2 100644
--- a/epclust/R/clustering.R
+++ b/epclust/R/clustering.R
@@ -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
 }
diff --git a/epclust/R/main.R b/epclust/R/main.R
index 2037dbe..977e61b 100644
--- a/epclust/R/main.R
+++ b/epclust/R/main.R
@@ -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
index 36bfba7..0000000
--- a/epclust/src/WER.c
+++ /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
index 0000000..98a0111
--- /dev/null
+++ b/epclust/src/computeMedoidsIndices.c
@@ -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
index 0000000..382910a
--- /dev/null
+++ b/epclust/src/filter.c
@@ -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;
+}
diff --git a/epclust/tests/testthat/test.clustering.R b/epclust/tests/testthat/test.clustering.R
index 49afe60..6c94f92 100644
--- a/epclust/tests/testthat/test.clustering.R
+++ b/epclust/tests/testthat/test.clustering.R
@@ -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),]) )
+#})
-- 
2.44.0