| 1 | #' computeWerDists |
| 2 | #' |
| 3 | #' Compute the WER distances between the series at specified indices, which are |
| 4 | #' obtaind by \code{getSeries(indices)} |
| 5 | #' |
| 6 | #' @param indices Range of series indices to cluster |
| 7 | #' @inheritParams claws |
| 8 | #' @inheritParams computeSynchrones |
| 9 | #' |
| 10 | #' @return A distances matrix of size K x K where K == length(indices) |
| 11 | #' |
| 12 | #' @export |
| 13 | computeWerDists <- function(indices, getSeries, nb_series_per_chunk, smooth_lvl, nvoice, |
| 14 | nbytes, endian, ncores_clust=3, verbose=FALSE, parll=TRUE) |
| 15 | { |
| 16 | n <- length(indices) |
| 17 | L <- length(getSeries(1)) #TODO: not very neat way to get L |
| 18 | noctave <- ceiling(log2(L)) #min power of 2 to cover serie range |
| 19 | # Since a CWT contains noctave*nvoice complex series, we deduce the number of CWT to |
| 20 | # retrieve/put in one chunk. |
| 21 | nb_cwt_per_chunk <- max(1, floor(nb_series_per_chunk / (nvoice*noctave*2))) |
| 22 | |
| 23 | # Initialize result as a square big.matrix of size 'number of medoids' |
| 24 | Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double") |
| 25 | |
| 26 | cwt_file <- tempfile(pattern="epclust_cwt.bin_") |
| 27 | # Compute the getSeries(indices) CWT, and store the results in the binary file |
| 28 | computeSaveCWT <- function(indices) |
| 29 | { |
| 30 | if (parll) |
| 31 | { |
| 32 | # parallel workers start with an empty environment |
| 33 | require("epclust", quietly=TRUE) |
| 34 | } |
| 35 | |
| 36 | # Obtain CWT as big vectors of real part + imaginary part (concatenate) |
| 37 | ts_cwt <- sapply(indices, function(i) { |
| 38 | ts <- scale(ts(getSeries(i)), center=TRUE, scale=FALSE) |
| 39 | ts_cwt <- Rwave::cwt(ts, noctave, nvoice, w0=2*pi, twoD=TRUE, plot=FALSE) |
| 40 | c( as.double(Re(ts_cwt)),as.double(Im(ts_cwt)) ) |
| 41 | }) |
| 42 | |
| 43 | # Serialization |
| 44 | binarize(ts_cwt, cwt_file, nb_cwt_per_chunk, ",", nbytes, endian) |
| 45 | } |
| 46 | |
| 47 | # Function to retrieve a synchrone CWT from (binary) file |
| 48 | getCWT <- function(index, L) |
| 49 | { |
| 50 | flat_cwt <- getDataInFile(index, cwt_file, nbytes, endian) |
| 51 | cwt_length <- length(flat_cwt) / 2 |
| 52 | re_part <- as.matrix(flat_cwt[1:cwt_length], nrow=L) |
| 53 | im_part <- as.matrix(flat_cwt[(cwt_length+1):(2*cwt_length)], nrow=L) |
| 54 | re_part + 1i * im_part |
| 55 | } |
| 56 | |
| 57 | # Compute distance between columns i and j for j>i |
| 58 | computeDistances <- function(i) |
| 59 | { |
| 60 | if (parll) |
| 61 | { |
| 62 | # parallel workers start with an empty environment |
| 63 | require("epclust", quietly=TRUE) |
| 64 | Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc) |
| 65 | } |
| 66 | |
| 67 | if (verbose && !parll) |
| 68 | cat(paste(" Distances from ",i," to ",i+1,"...",n,"\n", sep="")) |
| 69 | |
| 70 | # Get CWT of column i, and run computations for columns j>i |
| 71 | cwt_i <- getCWT(i, L) |
| 72 | WX <- filterMA(Mod(cwt_i * Conj(cwt_i)), smooth_lvl) |
| 73 | |
| 74 | for (j in (i+1):n) |
| 75 | { |
| 76 | cwt_j <- getCWT(j, L) |
| 77 | |
| 78 | # Compute the ratio of integrals formula 5.6 for WER^2 |
| 79 | # in https://arxiv.org/abs/1101.4744v2 paragraph 5.3 |
| 80 | num <- filterMA(Mod(cwt_i * Conj(cwt_j)), smooth_lvl) |
| 81 | WY <- filterMA(Mod(cwt_j * Conj(cwt_j)), smooth_lvl) |
| 82 | wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY)) |
| 83 | |
| 84 | Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2)) |
| 85 | Xwer_dist[j,i] <- Xwer_dist[i,j] |
| 86 | } |
| 87 | Xwer_dist[i,i] <- 0. |
| 88 | } |
| 89 | |
| 90 | if (parll) |
| 91 | { |
| 92 | # outfile=="" to see stderr/stdout on terminal |
| 93 | cl <- |
| 94 | if (verbose) |
| 95 | parallel::makeCluster(ncores_clust, outfile="") |
| 96 | else |
| 97 | parallel::makeCluster(ncores_clust) |
| 98 | Xwer_dist_desc <- bigmemory::describe(Xwer_dist) |
| 99 | parallel::clusterExport(cl, varlist=c("parll","nb_cwt_per_chunk","n","L", |
| 100 | "Xwer_dist_desc","noctave","nvoice","getCWT"), envir=environment()) |
| 101 | } |
| 102 | |
| 103 | if (verbose) |
| 104 | cat(paste("--- Precompute and serialize synchrones CWT\n", sep="")) |
| 105 | |
| 106 | # Split indices by packets of length at most nb_cwt_per_chunk |
| 107 | indices_cwt <- .splitIndices(seq_len(n), nb_cwt_per_chunk) |
| 108 | ignored <- |
| 109 | if (parll) |
| 110 | parallel::parLapply(cl, indices_cwt, computeSaveCWT) |
| 111 | else |
| 112 | lapply(indices_cwt, computeSaveCWT) |
| 113 | |
| 114 | if (verbose) |
| 115 | cat(paste("--- Compute WER distances\n", sep="")) |
| 116 | |
| 117 | ignored <- |
| 118 | if (parll) |
| 119 | parallel::parLapply(cl, seq_len(n-1), computeDistances) |
| 120 | else |
| 121 | lapply(seq_len(n-1), computeDistances) |
| 122 | Xwer_dist[n,n] <- 0. |
| 123 | |
| 124 | if (parll) |
| 125 | parallel::stopCluster(cl) |
| 126 | |
| 127 | unlink(cwt_file) #remove binary file |
| 128 | |
| 129 | Xwer_dist[,] #~small matrix K1 x K1 |
| 130 | } |