X-Git-Url: https://git.auder.net/?p=epclust.git;a=blobdiff_plain;f=epclust%2FR%2Fclustering.R;h=8be871531f22e7a7dfe7b8828744b59f4c058c79;hp=2ce4267ef35e717a7bea0a8667f648367d38fbcf;hb=a52836b23adb4bfa6722642ec6426fb7b5f39650;hpb=d9bb53c5e1392018bf67f92140edb10137f3423c diff --git a/epclust/R/clustering.R b/epclust/R/clustering.R index 2ce4267..8be8715 100644 --- a/epclust/R/clustering.R +++ b/epclust/R/clustering.R @@ -66,7 +66,7 @@ clusteringTask1 = function(indices, getContribs, K1, algoClust1, nb_items_clust1 #' @rdname clustering #' @export clusteringTask2 = function(medoids, K2, algoClust2, getRefSeries, nb_ref_curves, - nb_series_per_chunk, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE) + nb_series_per_chunk, nvoice, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE) { if (verbose) cat(paste("*** Clustering task 2 on ",ncol(medoids)," synchrones\n", sep="")) @@ -80,11 +80,12 @@ clusteringTask2 = function(medoids, K2, algoClust2, getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust, verbose, parll) # B) Compute the WER distances (Wavelets Extended coefficient of deteRmination) - distances = computeWerDists(synchrones, nbytes, endian, ncores_clust, verbose, parll) + distances = computeWerDists( + synchrones, nvoice, nbytes, endian, ncores_clust, verbose, parll) # C) Apply clustering algorithm 2 on the WER distances matrix if (verbose) - cat(paste(" algoClust2() on ",nrow(distances)," items\n", sep="")) + cat(paste("*** algoClust2() on ",nrow(distances)," items\n", sep="")) medoids[ ,algoClust2(distances,K2) ] } @@ -135,14 +136,15 @@ computeSynchrones = function(medoids, getRefSeries, nb_ref_curves, if (parll) synchronicity::unlock(m) } + NULL } K = ncol(medoids) ; L = nrow(medoids) # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in // synchrones = bigmemory::big.matrix(nrow=L, ncol=K, type="double", init=0.) # NOTE: synchronicity is only for Linux & MacOS; on Windows: run sequentially - parll = (requireNamespace("synchronicity",quietly=TRUE) - && parll && Sys.info()['sysname'] != "Windows") + parll = (parll && requireNamespace("synchronicity",quietly=TRUE) + && Sys.info()['sysname'] != "Windows") if (parll) { m <- synchronicity::boost.mutex() #for lock/unlock, see computeSynchronesChunk @@ -176,31 +178,26 @@ computeSynchrones = function(medoids, getRefSeries, nb_ref_curves, #' computeWerDists #' -#' Compute the WER distances between the synchrones curves (in rows), which are +#' Compute the WER distances between the synchrones curves (in columns), which are #' returned (e.g.) by \code{computeSynchrones()} #' -#' @param synchrones A big.matrix of synchrones, in rows. The series have same length -#' as the series in the initial dataset +#' @param synchrones A big.matrix of synchrones, in columns. The series have same +#' length as the series in the initial dataset #' @inheritParams claws #' -#' @return A matrix of size K1 x K1 +#' @return A distances matrix of size K1 x K1 #' #' @export -computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE) +computeWerDists = function(synchrones, nvoice, nbytes,endian,ncores_clust=1, + verbose=FALSE,parll=TRUE) { n <- ncol(synchrones) L <- nrow(synchrones) - #TODO: automatic tune of all these parameters ? (for other users) - # 4 here represent 2^5 = 32 half-hours ~ 1 day - nvoice <- 4 - # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones)) - noctave = 13 + noctave = ceiling(log2(L)) #min power of 2 to cover serie range + # Initialize result as a square big.matrix of size 'number of synchrones' Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double") - cwt_file = ".epclust_bin/cwt" - #TODO: args, nb_per_chunk, nbytes, endian - # Generate n(n-1)/2 pairs for WER distances computations pairs = list() V = seq_len(n) @@ -210,18 +207,23 @@ computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALS pairs = c(pairs, lapply(V, function(v) c(i,v))) } + cwt_file = ".cwt.bin" + # Compute the synchrones[,index] CWT, and store it in the binary file above computeSaveCWT = function(index) { + if (parll && !exists(synchrones)) #avoid going here after first call on a worker + { + require("bigmemory", quietly=TRUE) + require("Rwave", quietly=TRUE) + require("epclust", quietly=TRUE) + synchrones <- bigmemory::attach.big.matrix(synchrones_desc) + } ts <- scale(ts(synchrones[,index]), center=TRUE, scale=FALSE) - totts.cwt = Rwave::cwt(ts, totnoct, nvoice, w0=2*pi, twoD=TRUE, plot=FALSE) - ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)] - #Normalization - sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) - sqres <- sweep(ts.cwt,2,sqs,'*') - res <- sqres / max(Mod(sqres)) - #TODO: serializer les CWT, les récupérer via getDataInFile ; - #--> OK, faut juste stocker comme séries simples de taille L*n' (53*17519) - binarize(c(as.double(Re(res)),as.double(Im(res))), cwt_file, ncol(res), ",", nbytes, endian) + ts_cwt = Rwave::cwt(ts, noctave, nvoice, w0=2*pi, twoD=TRUE, plot=FALSE) + + # Serialization + binarize(as.matrix(c(as.double(Re(ts_cwt)),as.double(Im(ts_cwt)))), cwt_file, 1, + ",", nbytes, endian) } if (parll) @@ -229,35 +231,35 @@ computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALS cl = parallel::makeCluster(ncores_clust) synchrones_desc <- bigmemory::describe(synchrones) Xwer_dist_desc <- bigmemory::describe(Xwer_dist) - parallel::clusterExport(cl, envir=environment(), - varlist=c("synchrones_desc","Xwer_dist_desc","totnoct","nvoice","w0","s0log", - "noctave","s0","verbose","getCWT")) + parallel::clusterExport(cl, varlist=c("parll","synchrones_desc","Xwer_dist_desc", + "noctave","nvoice","verbose","getCWT"), envir=environment()) } - + if (verbose) - { - cat(paste("--- Compute WER dists\n", sep="")) - # precompute save all CWT........ - } - #precompute and serialize all CWT + cat(paste("--- Precompute and serialize synchrones CWT\n", sep="")) + ignored <- if (parll) parallel::parLapply(cl, 1:n, computeSaveCWT) else lapply(1:n, computeSaveCWT) - getCWT = function(index) + # Function to retrieve a synchrone CWT from (binary) file + getSynchroneCWT = function(index, L) { - #from cwt_file ... - res <- getDataInFile(c(2*index-1,2*index), cwt_file, nbytes, endian) - ###############TODO: + flat_cwt <- getDataInFile(index, cwt_file, nbytes, endian) + cwt_length = length(flat_cwt) / 2 + re_part = as.matrix(flat_cwt[1:cwt_length], nrow=L) + im_part = as.matrix(flat_cwt[(cwt_length+1):(2*cwt_length)], nrow=L) + re_part + 1i * im_part } - # Distance between rows i and j - computeDistancesIJ = function(pair) + # Compute distance between columns i and j in synchrones + computeDistanceIJ = function(pair) { if (parll) { + # parallel workers start with an empty environment require("bigmemory", quietly=TRUE) require("epclust", quietly=TRUE) synchrones <- bigmemory::attach.big.matrix(synchrones_desc) @@ -265,37 +267,42 @@ computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALS } i = pair[1] ; j = pair[2] - if (verbose && j==i+1) + if (verbose && j==i+1 && !parll) cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep="")) - cwt_i <- getCWT(i) - cwt_j <- getCWT(j) - num <- epclustFilter(Mod(cwt_i * Conj(cwt_j))) - WX <- epclustFilter(Mod(cwt_i * Conj(cwt_i))) - WY <- epclustFilter(Mod(cwt_j * Conj(cwt_j))) + # Compute CWT of columns i and j in synchrones + L = nrow(synchrones) + cwt_i <- getSynchroneCWT(i, L) + cwt_j <- getSynchroneCWT(j, L) + + # Compute the ratio of integrals formula 5.6 for WER^2 + # in https://arxiv.org/abs/1101.4744v2 §5.3 + num <- filterMA(Mod(cwt_i * Conj(cwt_j))) + WX <- filterMA(Mod(cwt_i * Conj(cwt_i))) + WY <- filterMA(Mod(cwt_j * Conj(cwt_j))) wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY)) - Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * max(1 - wer2, 0.)) + + Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2)) Xwer_dist[j,i] <- Xwer_dist[i,j] - Xwer_dist[i,i] = 0. + Xwer_dist[i,i] <- 0. } if (verbose) - { - cat(paste("--- Compute WER dists\n", sep="")) - } + cat(paste("--- Compute WER distances\n", sep="")) + ignored <- if (parll) - parallel::parLapply(cl, pairs, computeDistancesIJ) + parallel::parLapply(cl, pairs, computeDistanceIJ) else - lapply(pairs, computeDistancesIJ) + lapply(pairs, computeDistanceIJ) if (parll) parallel::stopCluster(cl) + unlink(cwt_file) + Xwer_dist[n,n] = 0. - distances <- Xwer_dist[,] - rm(Xwer_dist) ; gc() - distances #~small matrix K1 x K1 + Xwer_dist[,] #~small matrix K1 x K1 } # Helper function to divide indices into balanced sets @@ -318,7 +325,7 @@ computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALS if (max) { # Sets are not so well balanced, but size is supposed to be critical - return ( c( indices_workers, (L-rem+1):L ) ) + return ( c( indices_workers, if (rem>0) list((L-rem+1):L) else NULL ) ) } # Spread the remaining load among the workers