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