| 1 | #' computeWerDists |
| 2 | #' |
| 3 | #' Compute the WER distances between the synchrones curves (in columns), which are |
| 4 | #' returned (e.g.) by \code{computeSynchrones()} |
| 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, nvoice, nbytes, endian, |
| 14 | ncores_clust=1, verbose=FALSE, parll=TRUE) |
| 15 | { |
| 16 | n <- length(indices) |
| 17 | L <- length(getSeries(1)) #TODO: not very nice 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 | # Generate n(n-1)/2 pairs for WER distances computations |
| 27 | pairs = list() |
| 28 | V = seq_len(n) |
| 29 | for (i in 1:n) |
| 30 | { |
| 31 | V = V[-1] |
| 32 | pairs = c(pairs, lapply(V, function(v) c(i,v))) |
| 33 | } |
| 34 | |
| 35 | cwt_file = ".cwt.bin" |
| 36 | # Compute the getSeries(indices) CWT, and store the results in the binary file |
| 37 | computeSaveCWT = function(indices) |
| 38 | { |
| 39 | if (parll) |
| 40 | { |
| 41 | require("bigmemory", quietly=TRUE) |
| 42 | require("Rwave", quietly=TRUE) |
| 43 | require("epclust", quietly=TRUE) |
| 44 | } |
| 45 | |
| 46 | # Obtain CWT as big vectors of real part + imaginary part (concatenate) |
| 47 | ts_cwt <- sapply(indices, function(i) { |
| 48 | ts <- scale(ts(getSeries(i)), center=TRUE, scale=FALSE) |
| 49 | ts_cwt <- Rwave::cwt(ts, noctave, nvoice, w0=2*pi, twoD=TRUE, plot=FALSE) |
| 50 | c( as.double(Re(ts_cwt)),as.double(Im(ts_cwt)) ) |
| 51 | }) |
| 52 | |
| 53 | # Serialization |
| 54 | binarize(ts_cwt, cwt_file, nb_cwt_per_chunk, ",", nbytes, endian) |
| 55 | } |
| 56 | |
| 57 | if (parll) |
| 58 | { |
| 59 | cl = parallel::makeCluster(ncores_clust) |
| 60 | Xwer_dist_desc <- bigmemory::describe(Xwer_dist) |
| 61 | parallel::clusterExport(cl, varlist=c("parll","nb_cwt_per_chunk","L", |
| 62 | "Xwer_dist_desc","noctave","nvoice","getCWT"), envir=environment()) |
| 63 | } |
| 64 | |
| 65 | if (verbose) |
| 66 | cat(paste("--- Precompute and serialize synchrones CWT\n", sep="")) |
| 67 | |
| 68 | ignored <- |
| 69 | if (parll) |
| 70 | parallel::parLapply(cl, 1:n, computeSaveCWT) |
| 71 | else |
| 72 | lapply(1:n, computeSaveCWT) |
| 73 | |
| 74 | # Function to retrieve a synchrone CWT from (binary) file |
| 75 | getCWT = function(index, L) |
| 76 | { |
| 77 | flat_cwt <- getDataInFile(index, cwt_file, nbytes, endian) |
| 78 | cwt_length = length(flat_cwt) / 2 |
| 79 | re_part = as.matrix(flat_cwt[1:cwt_length], nrow=L) |
| 80 | im_part = as.matrix(flat_cwt[(cwt_length+1):(2*cwt_length)], nrow=L) |
| 81 | re_part + 1i * im_part |
| 82 | } |
| 83 | |
| 84 | |
| 85 | |
| 86 | |
| 87 | #TODO: better repartition here, |
| 88 | |
| 89 | |
| 90 | |
| 91 | # Compute distance between columns i and j in synchrones |
| 92 | computeDistanceIJ = function(pair) |
| 93 | { |
| 94 | if (parll) |
| 95 | { |
| 96 | # parallel workers start with an empty environment |
| 97 | require("bigmemory", quietly=TRUE) |
| 98 | require("epclust", quietly=TRUE) |
| 99 | Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc) |
| 100 | } |
| 101 | |
| 102 | i = pair[1] ; j = pair[2] |
| 103 | if (verbose && j==i+1 && !parll) |
| 104 | cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep="")) |
| 105 | |
| 106 | # Compute CWT of columns i and j in synchrones |
| 107 | cwt_i <- getCWT(i, L) |
| 108 | cwt_j <- getCWT(j, L) |
| 109 | |
| 110 | # Compute the ratio of integrals formula 5.6 for WER^2 |
| 111 | # in https://arxiv.org/abs/1101.4744v2 ยง5.3 |
| 112 | num <- filterMA(Mod(cwt_i * Conj(cwt_j))) |
| 113 | WX <- filterMA(Mod(cwt_i * Conj(cwt_i))) |
| 114 | WY <- filterMA(Mod(cwt_j * Conj(cwt_j))) |
| 115 | wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY)) |
| 116 | |
| 117 | Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2)) |
| 118 | Xwer_dist[j,i] <- Xwer_dist[i,j] |
| 119 | Xwer_dist[i,i] <- 0. |
| 120 | } |
| 121 | |
| 122 | if (verbose) |
| 123 | cat(paste("--- Compute WER distances\n", sep="")) |
| 124 | |
| 125 | ignored <- |
| 126 | if (parll) |
| 127 | parallel::parLapply(cl, pairs, computeDistanceIJ) |
| 128 | else |
| 129 | lapply(pairs, computeDistanceIJ) |
| 130 | |
| 131 | if (parll) |
| 132 | parallel::stopCluster(cl) |
| 133 | |
| 134 | unlink(cwt_file) |
| 135 | |
| 136 | Xwer_dist[n,n] = 0. |
| 137 | Xwer_dist[,] #~small matrix K1 x K1 |
| 138 | } |