'update'
[epclust.git] / epclust / R / computeWerDists.R
CommitLineData
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
3c5a4b08
BA
13computeWerDists = function(indices, getSeries, nb_series_per_chunk, nvoice, nbytes, endian,
14 ncores_clust=1, verbose=FALSE, parll=TRUE)
40f12a2f 15{
3c5a4b08
BA
16 n <- length(indices)
17 L <- length(getSeries(1)) #TODO: not very nice way to get L
40f12a2f 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.
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
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"
3c5a4b08 36 # Compute the getSeries(indices) CWT, and store the results in the binary file
40f12a2f
BA
37 computeSaveCWT = function(indices)
38 {
39 if (parll)
40 {
41 require("bigmemory", quietly=TRUE)
42 require("Rwave", quietly=TRUE)
43 require("epclust", quietly=TRUE)
40f12a2f
BA
44 }
45
46 # Obtain CWT as big vectors of real part + imaginary part (concatenate)
47 ts_cwt <- sapply(indices, function(i) {
3c5a4b08 48 ts <- scale(ts(getSeries(i)), center=TRUE, scale=FALSE)
40f12a2f
BA
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
3c5a4b08 54 binarize(ts_cwt, cwt_file, nb_cwt_per_chunk, ",", nbytes, endian)
40f12a2f
BA
55 }
56
57 if (parll)
58 {
59 cl = parallel::makeCluster(ncores_clust)
40f12a2f 60 Xwer_dist_desc <- bigmemory::describe(Xwer_dist)
3c5a4b08
BA
61 parallel::clusterExport(cl, varlist=c("parll","nb_cwt_per_chunk","L",
62 "Xwer_dist_desc","noctave","nvoice","getCWT"), envir=environment())
40f12a2f
BA
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
3c5a4b08 75 getCWT = function(index, L)
40f12a2f
BA
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,
40f12a2f
BA
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)
40f12a2f
BA
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
3c5a4b08
BA
107 cwt_i <- getCWT(i, L)
108 cwt_j <- getCWT(j, L)
40f12a2f
BA
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}