From: Benjamin Auder Date: Wed, 8 Mar 2017 00:52:24 +0000 (+0100) Subject: option to run sequentially. various fixes. R CMD check OK X-Git-Url: https://git.auder.net/doc/html/img/cross.svg?a=commitdiff_plain;h=492cd9e74a79cbcc0ecde55fa3071a44b7e463dc;p=epclust.git option to run sequentially. various fixes. R CMD check OK --- diff --git a/TODO b/TODO index 275c10d..1788ce6 100644 --- a/TODO +++ b/TODO @@ -38,3 +38,10 @@ utiliser Rcpp ? # x <- x / sqrt(rowSums(x^2)) # x %*% t(x) #} + +#TODO: soften condition clustering.R line 37 ? +#regarder mapply et mcmapply pour le // (pas OK pour Windows ou GUI... mais ?) +#TODO: map-reduce more appropriate R/clustering.R ligne 88 +#Alternative: use bigmemory to share series when CSV or matrix(...) + +#' @importFrom synchronicity boost.mutex lock unlock diff --git a/epclust/DESCRIPTION b/epclust/DESCRIPTION index 66d1712..304fdff 100644 --- a/epclust/DESCRIPTION +++ b/epclust/DESCRIPTION @@ -18,8 +18,10 @@ Imports: parallel, cluster, wavelets, + bigmemory, Rwave Suggests: + synchronicity, devtools, testthat, MASS, diff --git a/epclust/R/a_NAMESPACE.R b/epclust/R/a_NAMESPACE.R index 3453108..89aa453 100644 --- a/epclust/R/a_NAMESPACE.R +++ b/epclust/R/a_NAMESPACE.R @@ -9,5 +9,5 @@ #' @importFrom stats filter spline #' @importFrom utils tail #' @importFrom methods is +#' @importFrom bigmemory as.big.matrix is.big.matrix NULL - diff --git a/epclust/R/clustering.R b/epclust/R/clustering.R index 6408370..74d009e 100644 --- a/epclust/R/clustering.R +++ b/epclust/R/clustering.R @@ -1,15 +1,16 @@ #' @name clustering #' @rdname clustering -#' @aliases clusteringTask computeClusters1 computeClusters2 +#' @aliases clusteringTask1 computeClusters1 computeClusters2 #' -#' @title Two-stages clustering, withing one task (see \code{claws()}) +#' @title Two-stage clustering, withing one task (see \code{claws()}) #' -#' @description \code{clusteringTask()} runs one full task, which consists in iterated stage 1 -#' clustering (on nb_curves / ntasks energy contributions, computed through discrete -#' wavelets coefficients). \code{computeClusters1()} and \code{computeClusters2()} -#' correspond to the atomic clustering procedures respectively for stage 1 and 2. -#' The former applies the clustering algorithm (PAM) on a contributions matrix, while -#' the latter clusters a chunk of series inside one task (~max nb_series_per_chunk) +#' @description \code{clusteringTask1()} runs one full stage-1 task, which consists in +#' iterated stage 1 clustering (on nb_curves / ntasks energy contributions, computed +#' through discrete wavelets coefficients). \code{computeClusters1()} and +#' \code{computeClusters2()} correspond to the atomic clustering procedures respectively +#' for stage 1 and 2. The former applies the clustering algorithm (PAM) on a +#' contributions matrix, while the latter clusters a chunk of series inside one task +#' (~max nb_series_per_chunk) #' #' @param indices Range of series indices to cluster in parallel (initial data) #' @param getContribs Function to retrieve contributions from initial series indices: @@ -18,7 +19,7 @@ #' @inheritParams computeSynchrones #' @inheritParams claws #' -#' @return For \code{clusteringTask()} and \code{computeClusters1()}, the indices of the +#' @return For \code{clusteringTask1()} and \code{computeClusters1()}, the indices of the #' computed (K1) medoids. Indices are irrelevant for stage 2 clustering, thus #' \code{computeClusters2()} outputs a matrix of medoids #' (of size limited by nb_series_per_chunk) @@ -26,42 +27,36 @@ NULL #' @rdname clustering #' @export -clusteringTask = function(indices, getContribs, K1, nb_series_per_chunk, ncores_clust) +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="")) -#NOTE: comment out parallel sections for debugging -#propagate verbose arg ?! + 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) ] + } -# cl = parallel::makeCluster(ncores_clust) -# parallel::clusterExport(cl, varlist=c("getContribs","K1"), envir=environment()) - repeat + if (parll) { - -print(length(indices)) - - nb_workers = max( 1, floor( length(indices) / nb_series_per_chunk ) ) - indices_workers = lapply( seq_len(nb_workers), function(i) - indices[(nb_series_per_chunk*(i-1)+1):(nb_series_per_chunk*i)] ) - # Spread the remaining load among the workers - rem = length(indices) %% nb_series_per_chunk - while (rem > 0) - { - index = rem%%nb_workers + 1 - indices_workers[[index]] = c(indices_workers[[index]], tail(indices,rem)) - rem = rem - 1 - } -# indices = unlist( parallel::parLapply( cl, indices_workers, function(inds) { - indices = unlist( lapply( indices_workers, function(inds) { -# require("epclust", quietly=TRUE) - -print(paste(" ",length(inds))) ## PROBLEME ICI : 21104 ??! - - inds[ computeClusters1(getContribs(inds), K1) ] - } ) ) - if (length(indices) == K1) - break + cl = parallel::makeCluster(ncores_clust) + parallel::clusterExport(cl, varlist=c("getContribs","K1","verbose"), envir=environment()) } -# parallel::stopCluster(cl) + 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) ) + } + if (parll) + parallel::stopCluster(cl) + indices #medoids } @@ -72,10 +67,13 @@ computeClusters1 = function(contribs, K1) #' @rdname clustering #' @export -computeClusters2 = function(medoids, K2, getRefSeries, nb_series_per_chunk) +computeClusters2 = function(medoids, K2, + getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE) { - synchrones = computeSynchrones(medoids, getRefSeries, nb_series_per_chunk) - medoids[ cluster::pam(computeWerDists(synchrones), K2, diss=TRUE)$medoids , ] + synchrones = computeSynchrones(medoids, + getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust, verbose, parll) + distances = computeWerDists(synchrones, ncores_clust, verbose, parll) + medoids[ cluster::pam(distances, K2, diss=TRUE)$medoids , ] } #' computeSynchrones @@ -86,34 +84,67 @@ computeClusters2 = function(medoids, K2, getRefSeries, nb_series_per_chunk) #' @param medoids Matrix of medoids (curves of same legnth as initial series) #' @param getRefSeries Function to retrieve initial series (e.g. in stage 2 after series #' have been replaced by stage-1 medoids) +#' @param nb_ref_curves How many reference series? (This number is known at this stage) #' @inheritParams claws #' #' @export -computeSynchrones = function(medoids, getRefSeries, nb_series_per_chunk) +computeSynchrones = function(medoids, getRefSeries, + nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE) { - K = nrow(medoids) - synchrones = matrix(0, nrow=K, ncol=ncol(medoids)) - counts = rep(0,K) - index = 1 - repeat + computeSynchronesChunk = function(indices) { - range = (index-1) + seq_len(nb_series_per_chunk) - ref_series = getRefSeries(range) - if (is.null(ref_series)) - break + ref_series = getRefSeries(indices) #get medoids indices for this chunk of series for (i in seq_len(nrow(ref_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] = counts[j] + 1 + counts[j,1] = counts[j,1] + 1 + if (parll) + synchronicity::unlock(m) + } + } + + K = nrow(medoids) + # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in // + 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) + # Fork (// run) only on Linux & MacOS; on Windows: run sequentially + parll = (requireNamespace("synchronicity",quietly=TRUE) + && parll && Sys.info()['sysname'] != "Windows") + if (parll) + m <- synchronicity::boost.mutex() + + indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk) + for (inds in indices_workers) + { + if (verbose) + { + cat(paste("--- Compute synchrones for indices range ", + min(inds)," -> ",max(inds),"\n", sep="")) } - index = index + nb_series_per_chunk + if (parll) + ignored <- parallel::mcparallel(computeSynchronesChunk(inds)) + else + computeSynchronesChunk(inds) + } + if (parll) + parallel::mccollect() + + mat_syncs = matrix(nrow=K, ncol=ncol(medoids)) + vec_count = rep(NA, K) + #TODO: can we avoid this loop? + for (i in seq_len(K)) + { + mat_syncs[i,] = synchrones[i,] + vec_count[i] = counts[i,1] } #NOTE: odds for some clusters to be empty? (when series already come from stage 2) # ...maybe; but let's hope resulting K1' be still quite bigger than K2 - synchrones = sweep(synchrones, 1, counts, '/') - synchrones[ sapply(seq_len(K), function(i) all(!is.nan(synchrones[i,]))) , ] + mat_syncs = sweep(mat_syncs, 1, vec_count, '/') + mat_syncs[ sapply(seq_len(K), function(i) all(!is.nan(mat_syncs[i,]))) , ] } #' computeWerDists @@ -122,10 +153,11 @@ computeSynchrones = function(medoids, getRefSeries, nb_series_per_chunk) #' returned (e.g.) by \code{computeSynchrones()} #' #' @param synchrones A matrix of synchrones, in rows. The series have same length as the -#' series in the initial dataset +#' series in the initial dataset +#' @inheritParams claws #' #' @export -computeWerDists = function(synchrones) +computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE) { n <- nrow(synchrones) delta <- ncol(synchrones) @@ -143,8 +175,10 @@ computeWerDists = function(synchrones) s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 ) totnoct = noctave + as.integer(s0log/nvoice) + 1 - # (normalized) observations node with CWT - Xcwt4 <- lapply(seq_len(n), function(i) { + 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)] @@ -152,24 +186,96 @@ computeWerDists = function(synchrones) sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) sqres <- sweep(ts.cwt,MARGIN=2,sqs,'*') sqres / max(Mod(sqres)) - }) + } - Xwer_dist <- matrix(0., n, n) + if (parll) + { + cl = parallel::makeCluster(ncores_clust) + parallel::clusterExport(cl, varlist=c("getContribs","K1","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?!) - for (i in 1:(n-1)) + if (verbose) + cat("*** Compute WER distances from CWT\n") + + computeDistancesLineI = function(i) { + if (verbose) + cat(paste(" Line ",i,"\n", sep="")) for (j in (i+1):n) { - #TODO: later, compute CWT here (because not enough storage space for 200k series) - # 'circular=TRUE' is wrong, should just take values on the sides; to rewrite in C + #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() + + for (i in 1:(n-1)) + { + if (parll) + ignored <- parallel::mcparallel(computeDistancesLineI(i)) + else + computeDistancesLineI(i) + } + Xwer_dist[n,n] = 0. + + if (parll) + parallel::mccollect() + + mat_dists = matrix(nrow=n, ncol=n) + #TODO: avoid this loop? + for (i in 1:n) + mat_dists[i,] = Xwer_dist[i,] + mat_dists +} + +# Helper function to divide indices into balanced sets +.spreadIndices = function(indices, nb_per_chunk) +{ + L = length(indices) + nb_workers = floor( L / nb_per_chunk ) + if (nb_workers == 0) + { + # L < nb_series_per_chunk, 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_chunk*i)] ) + # Spread the remaining load among the workers + rem = L %% nb_per_chunk + while (rem > 0) + { + index = rem%%nb_workers + 1 + indices_workers[[index]] = c(indices_workers[[index]], indices[L-rem+1]) + rem = rem - 1 } } - diag(Xwer_dist) <- numeric(n) - Xwer_dist + indices_workers } diff --git a/epclust/R/de_serialize.R b/epclust/R/de_serialize.R index 8dde258..b6684d2 100644 --- a/epclust/R/de_serialize.R +++ b/epclust/R/de_serialize.R @@ -1,12 +1,15 @@ #' @name de_serialize #' @rdname de_serialize -#' @aliases binarize getDataInFile +#' @aliases binarize binarizeTransform getDataInFile #' -#' @title (De)Serialization of a matrix +#' @title (De)Serialization of a [big]matrix or data stream #' #' @description \code{binarize()} serializes a matrix or CSV file with minimal overhead, #' into a binary file. \code{getDataInFile()} achieves the inverse task: she retrieves -#' (ASCII) data rows from indices in the binary file +#' (ASCII) data rows from indices in the binary file. Finally, +#' \code{binarizeTransform()} serialize transformations of all data chunks; to use it, +#' a data-retrieval function must be provided, thus \code{binarize} will most likely be +#' used first (and then a function defined to seek in generated binary file) #' #' @param data_ascii Either a matrix or CSV file, with items in rows #' @param indices Indices of the lines to retrieve @@ -14,9 +17,12 @@ #' or intput (\code{getDataInFile}) #' @param nb_per_chunk Number of lines to process in one batch #' @inheritParams claws +#' @param getData Function to retrieve data chunks +#' @param transform Transformation function to apply on data chunks #' #' @return For \code{getDataInFile()}, the matrix with rows corresponding to the -#' requested indices +#' requested indices. \code{binarizeTransform} returns the number of processed lines. +#' \code{binarize} is designed to serialize in several calls, thus returns nothing. NULL #' @rdname de_serialize @@ -28,6 +34,7 @@ binarize = function(data_ascii, data_bin_file, nb_per_chunk, data_ascii = file(data_ascii, open="r") else if (methods::is(data_ascii,"connection") && !isOpen(data_ascii)) open(data_ascii) + is_matrix = !methods::is(data_ascii,"connection") first_write = (!file.exists(data_bin_file) || file.info(data_bin_file)$size == 0) data_bin = file(data_bin_file, open=ifelse(first_write,"wb","ab")) @@ -37,9 +44,9 @@ binarize = function(data_ascii, data_bin_file, nb_per_chunk, { #number of items always on 8 bytes writeBin(0L, data_bin, size=8, endian=endian) - if (is.matrix(data_ascii)) + if ( is_matrix ) data_length = ncol(data_ascii) - else #if (is(data, "connection")) + else #connection { data_line = scan(data_ascii, double(), sep=sep, nlines=1, quiet=TRUE) writeBin(data_line, data_bin, size=nbytes, endian=endian) @@ -47,18 +54,17 @@ binarize = function(data_ascii, data_bin_file, nb_per_chunk, } } - if (is.matrix(data_ascii)) + if (is_matrix) index = 1 repeat { - if (is.matrix(data_ascii)) + if ( is_matrix ) { - range = index:min(nrow(data_ascii),index+nb_per_chunk) data_chunk = - if (range[1] <= nrow(data_ascii)) - as.double(t(data_ascii[range,])) + if (index <= nrow(data_ascii)) + as.double(t(data_ascii[index:min(nrow(data_ascii),index+nb_per_chunk-1),])) else - integer(0) + double(0) index = index + nb_per_chunk } else @@ -70,16 +76,36 @@ binarize = function(data_ascii, data_bin_file, nb_per_chunk, if (first_write) { - #ecrire file_size-1 / (nbytes*nbWritten) en 0 dans bin_data ! ignored == file_size + # Write data_length, = (file_size-1) / (nbytes*nbWritten) at offset 0 in data_bin ignored = seek(data_bin, 0) writeBin(data_length, data_bin, size=8, endian=endian) } close(data_bin) - if (methods::is(data_ascii,"connection")) + if ( ! is_matrix ) close(data_ascii) } +#' @rdname de_serialize +#' @export +binarizeTransform = function(getData, transform, data_bin_file, nb_per_chunk, + nbytes=4, endian=.Platform$endian) +{ + nb_items = 0 + index = 1 + repeat + { + data_chunk = getData((index-1)+seq_len(nb_per_chunk)) + if (is.null(data_chunk)) + break + transformed_chunk = transform(data_chunk) + binarize(transformed_chunk, data_bin_file, nb_per_chunk, ",", nbytes, endian) + index = index + nb_per_chunk + nb_items = nb_items + nrow(data_chunk) + } + nb_items +} + #' @rdname de_serialize #' @export getDataInFile = function(indices, data_bin_file, nbytes=4, endian=.Platform$endian) diff --git a/epclust/R/main.R b/epclust/R/main.R index a982f4c..9064dfa 100644 --- a/epclust/R/main.R +++ b/epclust/R/main.R @@ -30,6 +30,7 @@ #' @param nbytes Number of bytes to serialize a floating-point number; 4 or 8 #' @param endian Endianness to use for (de)serialization. Use "little" or "big" for portability #' @param verbose Level of verbosity (0/FALSE for nothing or 1/TRUE for all; devel stage) +#' @param parll TRUE to fully parallelize; otherwise run sequentially (debug, comparison) #' #' @return A matrix of the final medoids curves (K2) in rows #' @@ -104,13 +105,14 @@ claws = function(getSeries, K1, K2, nb_series_per_chunk=50*K1, min_series_per_chunk=5*K1, #chunk size sep=",", #ASCII input separator nbytes=4, endian=.Platform$endian, #serialization (write,read) - verbose=FALSE) + verbose=FALSE, parll=TRUE) { # Check/transform arguments - if (!is.matrix(getSeries) && !is.function(getSeries) && - !methods::is(getSeries, "connection" && !is.character(getSeries))) + if (!is.matrix(getSeries) && !bigmemory::is.big.matrix(getSeries) + && !is.function(getSeries) + && !methods::is(getSeries,"connection") && !is.character(getSeries)) { - stop("'getSeries': matrix, function, file or valid connection (no NA)") + stop("'getSeries': [big]matrix, function, file or valid connection (no NA)") } K1 = .toInteger(K1, function(x) x>=2) K2 = .toInteger(K2, function(x) x>=2) @@ -148,16 +150,9 @@ claws = function(getSeries, K1, K2, nb_curves = 0 if (verbose) cat("...Compute contributions and serialize them\n") - repeat - { - series = getSeries((index-1)+seq_len(nb_series_per_chunk)) - if (is.null(series)) - break - contribs_chunk = curvesToContribs(series, wf, ctype) - binarize(contribs_chunk, contribs_file, nb_series_per_chunk, sep, nbytes, endian) - index = index + nb_series_per_chunk - nb_curves = nb_curves + nrow(contribs_chunk) - } + nb_curves = binarizeTransform(getSeries, + function(series) curvesToContribs(series, wf, ctype), + contribs_file, nb_series_per_chunk, nbytes, endian) getContribs = function(indices) getDataInFile(indices, contribs_file, nbytes, endian) if (nb_curves < min_series_per_chunk) @@ -174,28 +169,37 @@ claws = function(getSeries, K1, K2, }) if (verbose) cat(paste("...Run ",ntasks," x stage 1 in parallel\n",sep="")) -# cl = parallel::makeCluster(ncores_tasks) -# parallel::clusterExport(cl, varlist=c("getSeries","getContribs","K1","K2", -# "nb_series_per_chunk","ncores_clust","synchrones_file","sep","nbytes","endian"), -# envir = environment()) - # 1000*K1 indices [if WER=="end"], or empty vector [if WER=="mix"] --> series on file -# indices = unlist( parallel::parLapply(cl, indices_tasks, function(inds) { - indices = unlist( lapply(indices_tasks, function(inds) { -# require("epclust", quietly=TRUE) - - browser() #TODO: CONTINUE DEBUG HERE + if (parll) + { + cl = parallel::makeCluster(ncores_tasks) + parallel::clusterExport(cl, varlist=c("getSeries","getContribs","K1","K2","verbose","parll", + "nb_series_per_chunk","ncores_clust","synchrones_file","sep","nbytes","endian"), + envir = environment()) + } - indices_medoids = clusteringTask(inds,getContribs,K1,nb_series_per_chunk,ncores_clust) + runTwoStepClustering = function(inds) + { + if (parll) + require("epclust", quietly=TRUE) + indices_medoids = clusteringTask1( + inds, getContribs, K1, nb_series_per_chunk, ncores_clust, verbose, parll) if (WER=="mix") { - medoids2 = computeClusters2( - getSeries(indices_medoids), K2, getSeries, nb_series_per_chunk) + medoids2 = computeClusters2(getSeries(indices_medoids), + K2, getSeries, nb_curves, nb_series_per_chunk, ncores_clust, verbose, parll) binarize(medoids2, synchrones_file, nb_series_per_chunk, sep, nbytes, endian) return (vector("integer",0)) } indices_medoids - }) ) -# parallel::stopCluster(cl) + } + + # 1000*K1 indices [if WER=="end"], or empty vector [if WER=="mix"] --> series on file + if (parll) + indices = unlist( parallel::parLapply(cl, indices_tasks, runTwoStepClustering) ) + else + indices = unlist( lapply(indices_tasks, runTwoStepClustering) ) + if (parll) + parallel::stopCluster(cl) getRefSeries = getSeries synchrones_file = paste(bin_dir,"synchrones",sep="") ; unlink(synchrones_file) @@ -209,23 +213,18 @@ claws = function(getSeries, K1, K2, index = 1 if (verbose) cat("...Serialize contributions computed on synchrones\n") - repeat - { - series = getSeries((index-1)+seq_len(nb_series_per_chunk)) - if (is.null(series)) - break - contribs_chunk = curvesToContribs(series, wf, ctype) - binarize(contribs_chunk, contribs_file, nb_series_per_chunk, sep, nbytes, endian) - index = index + nb_series_per_chunk - } + ignored = binarizeTransform(getSeries, + function(series) curvesToContribs(series, wf, ctype), + contribs_file, nb_series_per_chunk, nbytes, endian) } # Run step2 on resulting indices or series (from file) if (verbose) cat("...Run final // stage 1 + stage 2\n") - indices_medoids = clusteringTask( - indices, getContribs, K1, nb_series_per_chunk, ncores_tasks*ncores_clust) - medoids = computeClusters2(getSeries(indices_medoids),K2,getRefSeries,nb_series_per_chunk) + indices_medoids = clusteringTask1( + indices, getContribs, K1, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose) + medoids = computeClusters2(getSeries(indices_medoids), + K2, getRefSeries, nb_curves, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose) # Cleanup unlink(bin_dir, recursive=TRUE) @@ -259,7 +258,7 @@ curvesToContribs = function(series, wf, ctype) }) ) } -# Helper for main function: check integer arguments with functiional conditions +# Check integer arguments with functional conditions .toInteger <- function(x, condition) { if (!is.integer(x)) diff --git a/epclust/tests/testthat/test.clustering.R b/epclust/tests/testthat/test.clustering.R index 9333876..49afe60 100644 --- a/epclust/tests/testthat/test.clustering.R +++ b/epclust/tests/testthat/test.clustering.R @@ -63,10 +63,11 @@ test_that("computeSynchrones behave as expected", for (i in seq_len(n)) series[i,] = s[[I(i,K)]] + rnorm(L,sd=0.01) getRefSeries = function(indices) { - indices = indices[indices < n] + indices = indices[indices <= n] if (length(indices)>0) series[indices,] else NULL } - synchrones = computeSynchrones(rbind(s1,s2,s3), getRefSeries, 100) + synchrones = computeSynchrones(rbind(s1,s2,s3), getRefSeries, n, 100, + verbose=TRUE, parll=FALSE) expect_equal(dim(synchrones), c(K,L)) for (i in 1:K) @@ -95,23 +96,23 @@ test_that("computeClusters2 behave as expected", for (i in seq_len(n)) series[i,] = s[[I(i,K1)]] + rnorm(L,sd=0.01) getRefSeries = function(indices) { - indices = indices[indices < n] + 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, 75) + medoids_K2 = computeClusters2(medoids_K1, K2, getRefSeries, n, 75, + verbose=TRUE, parll=FALSE) 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),]) ) }) -test_that("clusteringTask + computeClusters2 behave as expected", +test_that("clusteringTask1 + computeClusters2 behave as expected", { n = 900 x = seq(0,9.5,0.1) @@ -129,8 +130,10 @@ test_that("clusteringTask + computeClusters2 behave as expected", wf = "haar" ctype = "absolute" getContribs = function(indices) curvesToContribs(series[indices,],wf,ctype) - medoids_K1 = getSeries( clusteringTask(1:n, getContribs, K1, 75, 4) ) - medoids_K2 = computeClusters2(medoids_K1, K2, getSeries, 120) + 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) expect_equal(dim(medoids_K1), c(K1,L)) expect_equal(dim(medoids_K2), c(K2,L)) diff --git a/epclust/tests/testthat/test.de_serialize.R b/epclust/tests/testthat/test.de_serialize.R index a2fae5e..8403e6d 100644 --- a/epclust/tests/testthat/test.de_serialize.R +++ b/epclust/tests/testthat/test.de_serialize.R @@ -32,6 +32,33 @@ test_that("serialization + getDataInFile retrieve original data / from matrix", unlink(data_bin_file) }) +test_that("serialization + transform + getDataInFile retrieve original transforms", +{ + data_bin_file = "/tmp/epclust_test_t.bin" + unlink(data_bin_file) + + #dataset 200 lignes / 30 columns + data_ascii = matrix(runif(200*30,-10,10),ncol=30) + nbytes = 8 + endian = "little" + + binarize(data_ascii, data_bin_file, 500, ",", nbytes, endian) + # Serialize transformation (just compute range) into a new binary file + trans_bin_file = "/tmp/epclust_test_t_trans.bin" + unlink(trans_bin_file) + getSeries = function(inds) getDataInFile(inds, data_bin_file, nbytes, endian) + binarizeTransform(getSeries, function(series) t(apply(series, 1, range)), + trans_bin_file, 250, nbytes, endian) + unlink(data_bin_file) + expect_equal(file.info(trans_bin_file)$size, 2*nrow(data_ascii)*nbytes+8) + for (indices in list(c(1,3,5), 3:13, c(5,20,50), c(75,130:135), 196:200)) + { + trans_lines = getDataInFile(indices, trans_bin_file, nbytes, endian) + expect_equal(trans_lines, t(apply(data_ascii[indices,],1,range)), tolerance=1e-6) + } + unlink(trans_bin_file) +}) + test_that("serialization + getDataInFile retrieve original data / from connection", { data_bin_file = "/tmp/epclust_test_c.bin"