Commit | Line | Data |
---|---|---|
4bcfdbee BA |
1 | #' @name clustering |
2 | #' @rdname clustering | |
eef6f6c9 | 3 | #' @aliases clusteringTask1 clusteringTask2 computeClusters1 computeClusters2 |
4bcfdbee | 4 | #' |
492cd9e7 | 5 | #' @title Two-stage clustering, withing one task (see \code{claws()}) |
4bcfdbee | 6 | #' |
492cd9e7 BA |
7 | #' @description \code{clusteringTask1()} runs one full stage-1 task, which consists in |
8 | #' iterated stage 1 clustering (on nb_curves / ntasks energy contributions, computed | |
bf5c0844 BA |
9 | #' through discrete wavelets coefficients). |
10 | #' \code{clusteringTask2()} runs a full stage-2 task, which consists in synchrones | |
11 | #' and then WER distances computations, before applying the clustering algorithm. | |
12 | #' \code{computeClusters1()} and \code{computeClusters2()} correspond to the atomic | |
13 | #' clustering procedures respectively for stage 1 and 2. The former applies the | |
2b9f5356 BA |
14 | #' first clustering algorithm on a contributions matrix, while the latter clusters |
15 | #' a set of series inside one task (~nb_items_clust) | |
4bcfdbee BA |
16 | #' |
17 | #' @param indices Range of series indices to cluster in parallel (initial data) | |
18 | #' @param getContribs Function to retrieve contributions from initial series indices: | |
19 | #' \code{getContribs(indices)} outpus a contributions matrix | |
20 | #' @param contribs matrix of contributions (e.g. output of \code{curvesToContribs()}) | |
777c4b02 | 21 | #' @param distances matrix of K1 x K1 (WER) distances between synchrones |
4bcfdbee BA |
22 | #' @inheritParams computeSynchrones |
23 | #' @inheritParams claws | |
24 | #' | |
492cd9e7 | 25 | #' @return For \code{clusteringTask1()} and \code{computeClusters1()}, the indices of the |
4bcfdbee | 26 | #' computed (K1) medoids. Indices are irrelevant for stage 2 clustering, thus |
24ed5d83 | 27 | #' \code{computeClusters2()} outputs a big.matrix of medoids |
4bcfdbee BA |
28 | #' (of size limited by nb_series_per_chunk) |
29 | NULL | |
30 | ||
31 | #' @rdname clustering | |
32 | #' @export | |
492cd9e7 | 33 | clusteringTask1 = function( |
2b9f5356 BA |
34 | indices, getContribs, K1, nb_per_chunk, nb_items_clust, ncores_clust=1, |
35 | verbose=FALSE, parll=TRUE) | |
5c652979 | 36 | { |
492cd9e7 | 37 | if (verbose) |
e161499b | 38 | cat(paste("*** Clustering task 1 on ",length(indices)," lines\n", sep="")) |
4bcfdbee | 39 | |
2b9f5356 BA |
40 | |
41 | ||
42 | ||
43 | ||
44 | ||
45 | ##TODO: reviser le spreadIndices, tenant compte de nb_items_clust | |
46 | ||
47 | ##TODO: reviser / harmoniser avec getContribs qui en récupère pt'et + pt'et - !! | |
48 | ||
49 | ||
50 | ||
492cd9e7 | 51 | if (parll) |
7b13d0c2 | 52 | { |
492cd9e7 BA |
53 | cl = parallel::makeCluster(ncores_clust) |
54 | parallel::clusterExport(cl, varlist=c("getContribs","K1","verbose"), envir=environment()) | |
7b13d0c2 | 55 | } |
492cd9e7 BA |
56 | while (length(indices) > K1) |
57 | { | |
58 | indices_workers = .spreadIndices(indices, nb_series_per_chunk) | |
e161499b BA |
59 | indices <- |
60 | if (parll) | |
61 | { | |
62 | unlist( parallel::parLapply(cl, indices_workers, function(inds) { | |
63 | require("epclust", quietly=TRUE) | |
64 | inds[ computeClusters1(getContribs(inds), K1, verbose) ] | |
65 | }) ) | |
66 | } | |
67 | else | |
68 | { | |
69 | unlist( lapply(indices_workers, function(inds) | |
70 | inds[ computeClusters1(getContribs(inds), K1, verbose) ] | |
71 | ) ) | |
72 | } | |
492cd9e7 BA |
73 | } |
74 | if (parll) | |
75 | parallel::stopCluster(cl) | |
76 | ||
56857861 | 77 | indices #medoids |
5c652979 BA |
78 | } |
79 | ||
4bcfdbee BA |
80 | #' @rdname clustering |
81 | #' @export | |
a174b8ea BA |
82 | clusteringTask2 = function(medoids, K2, getRefSeries, nb_ref_curves, |
83 | nb_series_per_chunk, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE) | |
5c652979 | 84 | { |
e161499b BA |
85 | if (verbose) |
86 | cat(paste("*** Clustering task 2 on ",nrow(medoids)," lines\n", sep="")) | |
87 | ||
bf5c0844 BA |
88 | if (nrow(medoids) <= K2) |
89 | return (medoids) | |
492cd9e7 BA |
90 | synchrones = computeSynchrones(medoids, |
91 | getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust, verbose, parll) | |
a174b8ea | 92 | distances = computeWerDists(synchrones, nbytes, endian, ncores_clust, verbose, parll) |
777c4b02 | 93 | medoids[ computeClusters2(distances,K2,verbose), ] |
5c652979 BA |
94 | } |
95 | ||
bf5c0844 BA |
96 | #' @rdname clustering |
97 | #' @export | |
e161499b BA |
98 | computeClusters1 = function(contribs, K1, verbose=FALSE) |
99 | { | |
100 | if (verbose) | |
101 | cat(paste(" computeClusters1() on ",nrow(contribs)," lines\n", sep="")) | |
eef6f6c9 | 102 | cluster::pam( t(contribs) , K1, diss=FALSE)$id.med |
e161499b | 103 | } |
bf5c0844 BA |
104 | |
105 | #' @rdname clustering | |
106 | #' @export | |
e161499b BA |
107 | computeClusters2 = function(distances, K2, verbose=FALSE) |
108 | { | |
109 | if (verbose) | |
110 | cat(paste(" computeClusters2() on ",nrow(distances)," lines\n", sep="")) | |
eef6f6c9 | 111 | cluster::pam( distances , K2, diss=TRUE)$id.med |
e161499b | 112 | } |
bf5c0844 | 113 | |
4bcfdbee BA |
114 | #' computeSynchrones |
115 | #' | |
116 | #' Compute the synchrones curves (sum of clusters elements) from a matrix of medoids, | |
117 | #' using L2 distances. | |
118 | #' | |
24ed5d83 | 119 | #' @param medoids big.matrix of medoids (curves of same length as initial series) |
4bcfdbee BA |
120 | #' @param getRefSeries Function to retrieve initial series (e.g. in stage 2 after series |
121 | #' have been replaced by stage-1 medoids) | |
492cd9e7 | 122 | #' @param nb_ref_curves How many reference series? (This number is known at this stage) |
4bcfdbee BA |
123 | #' @inheritParams claws |
124 | #' | |
eef6f6c9 | 125 | #' @return A big.matrix of size L x K1 where L = length of a serie |
24ed5d83 | 126 | #' |
4bcfdbee | 127 | #' @export |
492cd9e7 BA |
128 | computeSynchrones = function(medoids, getRefSeries, |
129 | nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE) | |
e205f218 | 130 | { |
e161499b BA |
131 | if (verbose) |
132 | cat(paste("--- Compute synchrones\n", sep="")) | |
133 | ||
492cd9e7 | 134 | computeSynchronesChunk = function(indices) |
3eef8d3d | 135 | { |
363ae134 BA |
136 | if (parll) |
137 | { | |
138 | require("bigmemory", quietly=TRUE) | |
6ad3f3fd | 139 | requireNamespace("synchronicity", quietly=TRUE) |
363ae134 BA |
140 | require("epclust", quietly=TRUE) |
141 | synchrones <- bigmemory::attach.big.matrix(synchrones_desc) | |
6ad3f3fd | 142 | counts <- bigmemory::attach.big.matrix(counts_desc) |
363ae134 BA |
143 | medoids <- bigmemory::attach.big.matrix(medoids_desc) |
144 | m <- synchronicity::attach.mutex(m_desc) | |
145 | } | |
146 | ||
6ad3f3fd BA |
147 | ref_series = getRefSeries(indices) |
148 | nb_series = nrow(ref_series) | |
149 | ||
363ae134 | 150 | #get medoids indices for this chunk of series |
2c14dbea | 151 | mi = computeMedoidsIndices(medoids@address, ref_series) |
e161499b BA |
152 | |
153 | for (i in seq_len(nb_series)) | |
56857861 | 154 | { |
492cd9e7 BA |
155 | if (parll) |
156 | synchronicity::lock(m) | |
eef6f6c9 BA |
157 | synchrones[, mi[i] ] = synchrones[, mi[i] ] + ref_series[,i] |
158 | counts[ mi[i] ] = counts[ mi[i] ] + 1 #TODO: remove counts? ...or as arg?! | |
492cd9e7 BA |
159 | if (parll) |
160 | synchronicity::unlock(m) | |
161 | } | |
162 | } | |
163 | ||
e161499b | 164 | K = nrow(medoids) ; L = ncol(medoids) |
492cd9e7 | 165 | # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in // |
24ed5d83 | 166 | # TODO: if size > RAM (not our case), use file-backed big.matrix |
eef6f6c9 | 167 | synchrones = bigmemory::big.matrix(nrow=L, ncol=K, type="double", init=0.) |
e161499b | 168 | counts = bigmemory::big.matrix(nrow=K, ncol=1, type="double", init=0) |
24ed5d83 | 169 | # synchronicity is only for Linux & MacOS; on Windows: run sequentially |
492cd9e7 BA |
170 | parll = (requireNamespace("synchronicity",quietly=TRUE) |
171 | && parll && Sys.info()['sysname'] != "Windows") | |
172 | if (parll) | |
363ae134 | 173 | { |
492cd9e7 | 174 | m <- synchronicity::boost.mutex() |
363ae134 BA |
175 | m_desc <- synchronicity::describe(m) |
176 | synchrones_desc = bigmemory::describe(synchrones) | |
6ad3f3fd | 177 | counts_desc = bigmemory::describe(counts) |
363ae134 | 178 | medoids_desc = bigmemory::describe(medoids) |
24ed5d83 | 179 | cl = parallel::makeCluster(ncores_clust) |
6ad3f3fd BA |
180 | parallel::clusterExport(cl, varlist=c("synchrones_desc","counts_desc","counts", |
181 | "verbose","m_desc","medoids_desc","getRefSeries"), envir=environment()) | |
24ed5d83 BA |
182 | } |
183 | ||
492cd9e7 | 184 | indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk) |
c45fd663 | 185 | ignored <- |
492cd9e7 | 186 | if (parll) |
e161499b | 187 | parallel::parLapply(cl, indices_workers, computeSynchronesChunk) |
492cd9e7 | 188 | else |
c45fd663 | 189 | lapply(indices_workers, computeSynchronesChunk) |
492cd9e7 | 190 | |
24ed5d83 BA |
191 | if (parll) |
192 | parallel::stopCluster(cl) | |
193 | ||
194 | #TODO: can we avoid this loop? ( synchrones = sweep(synchrones, 1, counts, '/') ) | |
492cd9e7 | 195 | for (i in seq_len(K)) |
eef6f6c9 | 196 | synchrones[,i] = synchrones[,i] / counts[i] |
3eef8d3d | 197 | #NOTE: odds for some clusters to be empty? (when series already come from stage 2) |
8702eb86 | 198 | # ...maybe; but let's hope resulting K1' be still quite bigger than K2 |
eef6f6c9 | 199 | noNA_rows = sapply(seq_len(K), function(i) all(!is.nan(synchrones[,i]))) |
24ed5d83 BA |
200 | if (all(noNA_rows)) |
201 | return (synchrones) | |
202 | # Else: some clusters are empty, need to slice synchrones | |
eef6f6c9 | 203 | bigmemory::as.big.matrix(synchrones[,noNA_rows]) |
e205f218 | 204 | } |
1c6f223e | 205 | |
4bcfdbee BA |
206 | #' computeWerDists |
207 | #' | |
208 | #' Compute the WER distances between the synchrones curves (in rows), which are | |
209 | #' returned (e.g.) by \code{computeSynchrones()} | |
210 | #' | |
24ed5d83 BA |
211 | #' @param synchrones A big.matrix of synchrones, in rows. The series have same length |
212 | #' as the series in the initial dataset | |
492cd9e7 | 213 | #' @inheritParams claws |
4bcfdbee | 214 | #' |
777c4b02 | 215 | #' @return A matrix of size K1 x K1 |
24ed5d83 | 216 | #' |
4bcfdbee | 217 | #' @export |
a174b8ea | 218 | computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE) |
d03c0621 | 219 | { |
e161499b BA |
220 | if (verbose) |
221 | cat(paste("--- Compute WER dists\n", sep="")) | |
24ed5d83 | 222 | |
4bcfdbee BA |
223 | n <- nrow(synchrones) |
224 | delta <- ncol(synchrones) | |
db6fc17d | 225 | #TODO: automatic tune of all these parameters ? (for other users) |
d03c0621 | 226 | nvoice <- 4 |
4bcfdbee | 227 | # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones)) |
d7d55bc1 BA |
228 | noctave = 13 |
229 | # 4 here represent 2^5 = 32 half-hours ~ 1 day | |
db6fc17d | 230 | #NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?) |
24ed5d83 | 231 | scalevector <- 2^(4:(noctave * nvoice) / nvoice + 1) |
db6fc17d | 232 | #condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1 |
a174b8ea BA |
233 | s0 = 2 |
234 | w0 = 2*pi | |
db6fc17d BA |
235 | scaled=FALSE |
236 | s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 ) | |
237 | totnoct = noctave + as.integer(s0log/nvoice) + 1 | |
238 | ||
e161499b | 239 | Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double") |
e161499b | 240 | |
4204e877 BA |
241 | cwt_file = ".epclust_bin/cwt" |
242 | #TODO: args, nb_per_chunk, nbytes, endian | |
243 | ||
e161499b BA |
244 | # Generate n(n-1)/2 pairs for WER distances computations |
245 | pairs = list() | |
4204e877 BA |
246 | V = seq_len(n) |
247 | for (i in 1:n) | |
e161499b BA |
248 | { |
249 | V = V[-1] | |
4204e877 BA |
250 | pairs = c(pairs, lapply(V, function(v) c(i,v))) |
251 | } | |
a174b8ea | 252 | |
4204e877 BA |
253 | computeSaveCWT = function(index) |
254 | { | |
255 | ts <- scale(ts(synchrones[index,]), center=TRUE, scale=scaled) | |
256 | totts.cwt = Rwave::cwt(ts, totnoct, nvoice, w0, plot=FALSE) | |
257 | ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)] | |
258 | #Normalization | |
259 | sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0) | |
260 | sqres <- sweep(ts.cwt,2,sqs,'*') | |
261 | res <- sqres / max(Mod(sqres)) | |
262 | #TODO: serializer les CWT, les récupérer via getDataInFile ; | |
263 | #--> OK, faut juste stocker comme séries simples de taille delta*ncol (53*17519) | |
a174b8ea | 264 | binarize(c(as.double(Re(res)),as.double(Im(res))), cwt_file, ncol(res), ",", nbytes, endian) |
4204e877 BA |
265 | } |
266 | ||
267 | if (parll) | |
268 | { | |
269 | cl = parallel::makeCluster(ncores_clust) | |
270 | synchrones_desc <- bigmemory::describe(synchrones) | |
271 | Xwer_dist_desc <- bigmemory::describe(Xwer_dist) | |
272 | parallel::clusterExport(cl, varlist=c("synchrones_desc","Xwer_dist_desc","totnoct", | |
273 | "nvoice","w0","s0log","noctave","s0","verbose","getCWT"), envir=environment()) | |
274 | } | |
275 | ||
276 | #precompute and serialize all CWT | |
277 | ignored <- | |
278 | if (parll) | |
279 | parallel::parLapply(cl, 1:n, computeSaveCWT) | |
280 | else | |
281 | lapply(1:n, computeSaveCWT) | |
282 | ||
283 | getCWT = function(index) | |
284 | { | |
285 | #from cwt_file ... | |
a174b8ea | 286 | res <- getDataInFile(c(2*index-1,2*index), cwt_file, nbytes, endian) |
eef6f6c9 | 287 | ###############TODO: |
4204e877 | 288 | } |
e161499b | 289 | |
777c4b02 | 290 | # Distance between rows i and j |
e161499b BA |
291 | computeDistancesIJ = function(pair) |
292 | { | |
2c14dbea | 293 | if (parll) |
363ae134 | 294 | { |
2c14dbea BA |
295 | require("bigmemory", quietly=TRUE) |
296 | require("epclust", quietly=TRUE) | |
297 | synchrones <- bigmemory::attach.big.matrix(synchrones_desc) | |
298 | Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc) | |
299 | } | |
300 | ||
e161499b BA |
301 | i = pair[1] ; j = pair[2] |
302 | if (verbose && j==i+1) | |
303 | cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep="")) | |
4204e877 BA |
304 | cwt_i <- getCWT(i) |
305 | cwt_j <- getCWT(j) | |
2c14dbea | 306 | |
363ae134 BA |
307 | num <- epclustFilter(Mod(cwt_i * Conj(cwt_j))) |
308 | WX <- epclustFilter(Mod(cwt_i * Conj(cwt_i))) | |
4204e877 | 309 | WY <- epclustFilter(Mod(cwt_j * Conj(cwt_j))) |
e161499b | 310 | wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY)) |
2c14dbea | 311 | Xwer_dist[i,j] <- sqrt(delta * ncol(cwt_i) * max(1 - wer2, 0.)) #FIXME: wer2 should be < 1 |
e161499b BA |
312 | Xwer_dist[j,i] <- Xwer_dist[i,j] |
313 | Xwer_dist[i,i] = 0. | |
314 | } | |
315 | ||
e161499b | 316 | ignored <- |
492cd9e7 | 317 | if (parll) |
e161499b | 318 | parallel::parLapply(cl, pairs, computeDistancesIJ) |
492cd9e7 | 319 | else |
e161499b | 320 | lapply(pairs, computeDistancesIJ) |
492cd9e7 BA |
321 | |
322 | if (parll) | |
323 | parallel::stopCluster(cl) | |
6ad3f3fd | 324 | |
492cd9e7 | 325 | Xwer_dist[n,n] = 0. |
777c4b02 BA |
326 | distances <- Xwer_dist[,] |
327 | rm(Xwer_dist) ; gc() | |
328 | distances #~small matrix K1 x K1 | |
492cd9e7 BA |
329 | } |
330 | ||
331 | # Helper function to divide indices into balanced sets | |
332 | .spreadIndices = function(indices, nb_per_chunk) | |
333 | { | |
334 | L = length(indices) | |
335 | nb_workers = floor( L / nb_per_chunk ) | |
336 | if (nb_workers == 0) | |
337 | { | |
338 | # L < nb_series_per_chunk, simple case | |
339 | indices_workers = list(indices) | |
340 | } | |
341 | else | |
342 | { | |
343 | indices_workers = lapply( seq_len(nb_workers), function(i) | |
344 | indices[(nb_per_chunk*(i-1)+1):(nb_per_chunk*i)] ) | |
345 | # Spread the remaining load among the workers | |
346 | rem = L %% nb_per_chunk | |
347 | while (rem > 0) | |
348 | { | |
349 | index = rem%%nb_workers + 1 | |
350 | indices_workers[[index]] = c(indices_workers[[index]], indices[L-rem+1]) | |
351 | rem = rem - 1 | |
d03c0621 | 352 | } |
1c6f223e | 353 | } |
492cd9e7 | 354 | indices_workers |
1c6f223e | 355 | } |