X-Git-Url: https://git.auder.net/?p=epclust.git;a=blobdiff_plain;f=epclust%2FR%2FcomputeWerDists.R;h=061c3609a54bc10928c5cd36702ad326f34b0eb6;hp=a813b8f615d8840762144973d75162d876f34b6f;hb=282342bafdc9ff65c5df98c6e2304d63b33b9fb2;hpb=3c5a4b0880db63367a474a568e1322b3999932fe diff --git a/epclust/R/computeWerDists.R b/epclust/R/computeWerDists.R index a813b8f..061c360 100644 --- a/epclust/R/computeWerDists.R +++ b/epclust/R/computeWerDists.R @@ -10,31 +10,22 @@ #' @return A distances matrix of size K x K where K == length(indices) #' #' @export -computeWerDists = function(indices, getSeries, nb_series_per_chunk, nvoice, nbytes, endian, - ncores_clust=1, verbose=FALSE, parll=TRUE) +computeWerDists <- function(indices, getSeries, nb_series_per_chunk, smooth_lvl, nvoice, + nbytes, endian, ncores_clust=1, verbose=FALSE, parll=TRUE) { n <- length(indices) - L <- length(getSeries(1)) #TODO: not very nice way to get L - noctave = ceiling(log2(L)) #min power of 2 to cover serie range + L <- length(getSeries(1)) #TODO: not very neat way to get L + noctave <- ceiling(log2(L)) #min power of 2 to cover serie range # Since a CWT contains noctave*nvoice complex series, we deduce the number of CWT to # retrieve/put in one chunk. - nb_cwt_per_chunk = max(1, floor(nb_series_per_chunk / (nvoice*noctave*2))) + nb_cwt_per_chunk <- max(1, floor(nb_series_per_chunk / (nvoice*noctave*2))) # Initialize result as a square big.matrix of size 'number of medoids' Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double") - # Generate n(n-1)/2 pairs for WER distances computations - pairs = list() - V = seq_len(n) - for (i in 1:n) - { - V = V[-1] - pairs = c(pairs, lapply(V, function(v) c(i,v))) - } - - cwt_file = ".cwt.bin" + cwt_file <- tempfile(pattern="epclust_cwt.bin_") # Compute the getSeries(indices) CWT, and store the results in the binary file - computeSaveCWT = function(indices) + computeSaveCWT <- function(indices) { if (parll) { @@ -54,42 +45,18 @@ computeWerDists = function(indices, getSeries, nb_series_per_chunk, nvoice, nbyt binarize(ts_cwt, cwt_file, nb_cwt_per_chunk, ",", nbytes, endian) } - if (parll) - { - cl = parallel::makeCluster(ncores_clust) - Xwer_dist_desc <- bigmemory::describe(Xwer_dist) - parallel::clusterExport(cl, varlist=c("parll","nb_cwt_per_chunk","L", - "Xwer_dist_desc","noctave","nvoice","getCWT"), envir=environment()) - } - - if (verbose) - cat(paste("--- Precompute and serialize synchrones CWT\n", sep="")) - - ignored <- - if (parll) - parallel::parLapply(cl, 1:n, computeSaveCWT) - else - lapply(1:n, computeSaveCWT) - # Function to retrieve a synchrone CWT from (binary) file - getCWT = function(index, L) + getCWT <- function(index, L) { 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) + 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 } - - - -#TODO: better repartition here, - - - - # Compute distance between columns i and j in synchrones - computeDistanceIJ = function(pair) + # Compute distance between columns i and j for j>i + computeDistances <- function(i) { if (parll) { @@ -99,40 +66,63 @@ computeWerDists = function(indices, getSeries, nb_series_per_chunk, nvoice, nbyt Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc) } - i = pair[1] ; j = pair[2] - if (verbose && j==i+1 && !parll) - cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep="")) + if (verbose && !parll) + cat(paste(" Distances from ",i," to ",i+1,"...",n,"\n", sep="")) - # Compute CWT of columns i and j in synchrones + # Get CWT of column i, and run computations for columns j>i cwt_i <- getCWT(i, L) - cwt_j <- getCWT(j, L) + WX <- filterMA(Mod(cwt_i * Conj(cwt_i)), smooth_lvl) + + for (j in (i+1):n) + { + cwt_j <- getCWT(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)) + # 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)), smooth_lvl) + WY <- filterMA(Mod(cwt_j * Conj(cwt_j)), smooth_lvl) + wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY)) - Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2)) - Xwer_dist[j,i] <- Xwer_dist[i,j] + Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2)) + Xwer_dist[j,i] <- Xwer_dist[i,j] + } Xwer_dist[i,i] <- 0. } + if (parll) + { + # outfile=="" to see stderr/stdout on terminal + cl <- parallel::makeCluster(ncores_clust, outfile="") + Xwer_dist_desc <- bigmemory::describe(Xwer_dist) + parallel::clusterExport(cl, varlist=c("parll","nb_cwt_per_chunk","n","L", + "Xwer_dist_desc","noctave","nvoice","getCWT"), envir=environment()) + } + + if (verbose) + cat(paste("--- Precompute and serialize synchrones CWT\n", sep="")) + + # Split indices by packets of length at most nb_cwt_per_chunk + indices_cwt <- .splitIndices(seq_len(n), nb_cwt_per_chunk) + ignored <- + if (parll) + parallel::parLapply(cl, indices_cwt, computeSaveCWT) + else + lapply(indices_cwt, computeSaveCWT) + if (verbose) cat(paste("--- Compute WER distances\n", sep="")) ignored <- if (parll) - parallel::parLapply(cl, pairs, computeDistanceIJ) + parallel::parLapply(cl, seq_len(n-1), computeDistances) else - lapply(pairs, computeDistanceIJ) + lapply(seq_len(n-1), computeDistances) + Xwer_dist[n,n] <- 0. if (parll) parallel::stopCluster(cl) - unlink(cwt_file) + unlink(cwt_file) #remove binary file - Xwer_dist[n,n] = 0. Xwer_dist[,] #~small matrix K1 x K1 }