save state; test clustering not OK, all others OK

diff --git a/epclust/R/clustering.R b/epclust/R/clustering.R
index b09c1bc..8662f89 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
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 ?