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