drop enercast submodule; drop Rcpp requirement; fix doc, complete code, fix fix fix
[epclust.git] / epclust / R / computeWerDists.R
index a813b8f..061c360 100644 (file)
 #' @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
 }