From 0fe757f750f51e580d2c5a7b7f7df87cc405d12d Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Fri, 10 Mar 2017 22:55:26 +0100
Subject: [PATCH] save state; test clustering not OK, all others OK

---
 epclust/R/clustering.R                        |  6 +--
 epclust/R/main.R                              |  3 +-
 epclust/src/computeMedoidsIndices.cpp         |  8 ++--
 epclust/src/filter.cpp                        | 45 ++++++++++++-------
 epclust/tests/testthat/test.clustering.R      |  9 ++--
 .../testthat/test.computeMedoidsIndices.R     | 11 ++---
 epclust/tests/testthat/test.filter.R          | 12 ++---
 epclust/tests/testthat/test.wavelets.R        |  3 ++
 8 files changed, 58 insertions(+), 39 deletions(-)
 create mode 100644 epclust/tests/testthat/test.wavelets.R

diff --git a/epclust/R/clustering.R b/epclust/R/clustering.R
index b09c1bc..8662f89 100644
--- a/epclust/R/clustering.R
+++ b/epclust/R/clustering.R
@@ -318,16 +318,16 @@ computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALS
 {
 	L = length(indices)
 	nb_workers = floor( L / nb_per_set )
-	rem = L %% max_nb_per_set
+	rem = L %% nb_per_set
 	if (nb_workers == 0 || (nb_workers==1 && rem==0))
 	{
-		# L <= max_nb_per_set, simple case
+		# L <= nb_per_set, simple case
 		indices_workers = list(indices)
 	}
 	else
 	{
 		indices_workers = lapply( seq_len(nb_workers), function(i)
-			indices[(nb_per_chunk*(i-1)+1):(nb_per_set*i)] )
+			indices[(nb_per_set*(i-1)+1):(nb_per_set*i)] )
 		# Spread the remaining load among the workers
 		rem = L %% nb_per_set
 		while (rem > 0)
diff --git a/epclust/R/main.R b/epclust/R/main.R
index 9e9b641..e003933 100644
--- a/epclust/R/main.R
+++ b/epclust/R/main.R
@@ -348,9 +348,10 @@ claws <- function(getSeries, K1, K2, nb_series_per_chunk,
 #' @return A [big.]matrix of size log(L) x n containing contributions in columns
 #'
 #' @export
-curvesToContribs = function(series, wav_filt, contrib_type)
+curvesToContribs = function(series, wav_filt, contrib_type, coin=FALSE)
 {
 	L = nrow(series)
+	if (coin) browser()
 	D = ceiling( log2(L) )
 	nb_sample_points = 2^D
 	apply(series, 2, function(x) {
diff --git a/epclust/src/computeMedoidsIndices.cpp b/epclust/src/computeMedoidsIndices.cpp
index 6934181..f247584 100644
--- a/epclust/src/computeMedoidsIndices.cpp
+++ b/epclust/src/computeMedoidsIndices.cpp
@@ -21,9 +21,9 @@ IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix ref_series)
 {
 	XPtr<BigMatrix> pMed(pMedoids);
 	MatrixAccessor<double> medoids = MatrixAccessor<double>(*pMed);
-	int nb_series = ref_series.nrow(),
-		K = pMed->nrow(),
-		L = pMed->ncol();
+	int nb_series = ref_series.ncol(),
+		K = pMed->ncol(),
+		L = pMed->nrow();
 	IntegerVector mi(nb_series);
 
 	for (int i=0; i<nb_series ; i++)
@@ -36,7 +36,7 @@ IntegerVector computeMedoidsIndices(SEXP pMedoids, NumericMatrix ref_series)
 			double dist_ij = 0.;
 			for (int k=0; k<L; k++)
 			{
-				double delta = ref_series(i,k) - *(medoids[k]+j);
+				double delta = ref_series(k,i) - *(medoids[j]+k);
 				dist_ij += delta * delta;
 			}
 			if (dist_ij < best_dist)
diff --git a/epclust/src/filter.cpp b/epclust/src/filter.cpp
index ee24af6..cb5a8a0 100644
--- a/epclust/src/filter.cpp
+++ b/epclust/src/filter.cpp
@@ -1,5 +1,12 @@
 #include <Rcpp.h>
 
+#include <R.h>           //  Rprintf()
+//#include <R_ext/Utils.h> //  user interrupts
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+#include <stdio.h>
+
 using namespace Rcpp;
 
 //' filter
@@ -10,31 +17,35 @@ using namespace Rcpp;
 //'
 //' @return The filtered CWT
 // [[Rcpp::export]]
-NumericMatrix epclustFilter(NumericMatrix cwt)
+RcppExport SEXP epclustFilter(SEXP cwt_)
 {
-	int L = cwt.nrow(),
-		D = cwt.ncol();
-	NumericMatrix fcwt(L, D); //fill with 0... TODO: back to SEXP C-style?
-	double *cwt_c = cwt.begin(),
-		*fcwt_c = fcwt.begin();
+	int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0],
+		D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1];
+	double *cwt = REAL(cwt_);
+	SEXP fcwt_;
+	PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
+		double* fcwt = REAL(fcwt_); //(double*)malloc(L*D*sizeof(double));
 
 	//TODO: coding style is terrible... no time for now.
 	for (int col=0; col<D; col++)
 	{
-		double v1 = cwt_c[0];
-		double ma = v1 + cwt[1] + cwt_c[2];
+		double v1 = cwt[0];
+		double ma = v1 + cwt[1] + cwt[2];
 		for (int i=1; i<L-2; i++)
 		{
-			fcwt_c[i] = ma / 3.;
-			ma = ma - v1 + cwt_c[i+2];
-			v1 = cwt_c[i];
+			fcwt[i] = ma / 3.;
+			ma = ma - v1 + cwt[i+2];
+			v1 = cwt[i];
 		}
-		fcwt_c[0] = cwt_c[0];
-		fcwt_c[L-2] = ma / 3.;
-		fcwt_c[L-1] = cwt_c[L-1];
-		cwt_c = cwt_c + L;
-		fcwt_c = fcwt_c + L;
+		fcwt[0] = cwt[0];
+		fcwt[L-2] = ma / 3.;
+		fcwt[L-1] = cwt[L-1];
+		cwt = cwt + L;
+		fcwt = fcwt + L;
 	}
 
-	return fcwt;
+//	REAL(fcwt_) = fcwt;
+	UNPROTECT(1);
+
+	return fcwt_;
 }
diff --git a/epclust/tests/testthat/test.clustering.R b/epclust/tests/testthat/test.clustering.R
index a5dc3bd..77faeb9 100644
--- a/epclust/tests/testthat/test.clustering.R
+++ b/epclust/tests/testthat/test.clustering.R
@@ -62,7 +62,7 @@ test_that("clusteringTask1 behave as expected",
 	L = length(x) #96 1/4h
 	K1 = 60
 	s = lapply( seq_len(K1), function(i) x^(1+i/30)*cos(x+i) )
-	series = matrix(nrow=n, ncol=L)
+	series = matrix(nrow=L, ncol=n)
 	for (i in seq_len(n))
 		series[,i] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
 	getSeries = function(indices) {
@@ -72,7 +72,10 @@ test_that("clusteringTask1 behave as expected",
 	wf = "haar"
 	ctype = "absolute"
 	getContribs = function(indices) curvesToContribs(series[,indices],wf,ctype)
-	indices1 = clusteringTask1(1:n, getContribs, K1, 75, verbose=TRUE, parll=FALSE)
+	require("cluster", quietly=TRUE)
+	browser()
+	algoClust1 = function(contribs,K) cluster::pam(contribs,K,diss=FALSE)$id.med
+	indices1 = clusteringTask1(1:n, getContribs, K1, algoClust1, 75, verbose=TRUE, parll=FALSE)
 	medoids_K1 = getSeries(indices1)
 
 	expect_equal(dim(medoids_K1), c(L,K1))
@@ -92,7 +95,7 @@ test_that("clusteringTask2 behave as expected",
 	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)
+	series = matrix(nrow=L, ncol=n)
 	for (i in seq_len(n))
 		series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01)
 	getRefSeries = function(indices) {
diff --git a/epclust/tests/testthat/test.computeMedoidsIndices.R b/epclust/tests/testthat/test.computeMedoidsIndices.R
index efd6af9..8347fb6 100644
--- a/epclust/tests/testthat/test.computeMedoidsIndices.R
+++ b/epclust/tests/testthat/test.computeMedoidsIndices.R
@@ -8,17 +8,18 @@ test_that("computeMedoidsIndices behave as expected",
 	medoids = cbind( rep(0,L), rep(-5,L), rep(5,L) )
 	# short series...
 	series = t( rbind( MASS::mvrnorm(n/3, medoids[,1], diag(L)),
-		MASS::mvrnorm(n/3, medoids[,2], diag(L),
-		MASS::mvrnorm(n/3, medoids[,3], diag(L))) ) )
+		MASS::mvrnorm(n/3, medoids[,2], diag(L)),
+		MASS::mvrnorm(n/3, medoids[,3], diag(L)) ) )
 
 	# With high probability, medoids indices should resemble 1,1,1,...,2,2,2,...,3,3,3,...
-	mi = epclust:::.computeMedoidsIndices(medoids, series)
+	require("bigmemory", quietly=TRUE)
+	mi = epclust:::computeMedoidsIndices(bigmemory::as.big.matrix(medoids)@address, series)
 	mi_ref = rep(1:3, each=n/3)
-	expect_that( mean(mi != mi_ref) < 0.01 )
+	expect_lt( mean(mi != mi_ref), 0.01 )
 
 	# Now with a random matrix, compare with (trusted) R version
 	series = matrix(runif(n*L, min=-7, max=7), nrow=L)
-	mi = epclust:::.computeMedoidsIndices(medoids, series)
+	mi = epclust:::computeMedoidsIndices(bigmemory::as.big.matrix(medoids)@address, series)
 	mi_ref = R_computeMedoidsIndices(medoids, series)
 	expect_equal(mi, mi_ref)
 })
diff --git a/epclust/tests/testthat/test.filter.R b/epclust/tests/testthat/test.filter.R
index 9d1916d..a0263a4 100644
--- a/epclust/tests/testthat/test.filter.R
+++ b/epclust/tests/testthat/test.filter.R
@@ -6,13 +6,13 @@ test_that("[time-]serie filtering behave as expected",
 	# Currently just a mean of 3 values
 	M = matrix(runif(1000,min=-7,max=7), ncol=10)
 	ref_fM = stats::filter(M, c(1/3,1/3,1/3), circular=FALSE)
-	fM = epclust:::.epclustFilter(M)
+	fM = epclust:::epclustFilter(M)
 
 	#Expect an agreement on all inner values
 	expect_equal(dim(fM), c(100,10))
-	expect_equal(fM[2:99,], ref_fM[,2:99])
+	expect_equal(fM[2:99,], ref_fM[2:99,])
 
-	#For border values, just apply a "by-hand" mean
-	expect_equal(fM[1,], colMeans(rbind(M[1,],M[2,],M[100,])))
-	expect_equal(fM[100,], colMeans(rbind(M[100,],M[99,],M[1,])))
-}
+	#Border values should be unchanged
+	expect_equal(fM[1,], M[1,])
+	expect_equal(fM[100,], M[100,])
+})
diff --git a/epclust/tests/testthat/test.wavelets.R b/epclust/tests/testthat/test.wavelets.R
new file mode 100644
index 0000000..5eb19c5
--- /dev/null
+++ b/epclust/tests/testthat/test.wavelets.R
@@ -0,0 +1,3 @@
+#TODO: test computeWerDists (reference result? Jairo?)
+
+#TODO: test sur curvesToCoefs! ref results ?
-- 
2.44.0