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