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