save intermediate
[epclust.git] / epclust / R / clustering.R
... / ...
CommitLineData
1#' @name clustering
2#' @rdname clustering
3#' @aliases clusteringTask1 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)
29NULL
30
31#' @rdname clustering
32#' @export
33clusteringTask1 = function(
34 indices, getContribs, K1, nb_series_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
70clusteringTask2 = function(medoids, K2,
71 getRefSeries, nb_ref_curves, nb_series_per_chunk, 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, ncores_clust, verbose, parll)
81 medoids[ computeClusters2(distances,K2,verbose), ]
82}
83
84#' @rdname clustering
85#' @export
86computeClusters1 = function(contribs, K1, verbose=FALSE)
87{
88 if (verbose)
89 cat(paste(" computeClusters1() on ",nrow(contribs)," lines\n", sep=""))
90 cluster::pam(contribs, K1, diss=FALSE)$id.med
91}
92
93#' @rdname clustering
94#' @export
95computeClusters2 = 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 K1 x L where L = data_length
114#'
115#' @export
116computeSynchrones = 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?
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=K, ncol=L, 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,1]
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 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
206computeWerDists = function(synchrones, ncores_clust=1,verbose=FALSE,parll=TRUE)
207{
208 if (verbose)
209 cat(paste("--- Compute WER dists\n", sep=""))
210
211
212
213
214#TODO: serializer les CWT, les récupérer via getDataInFile
215#--> OK, faut juste stocker comme séries simples de taille delta*ncol (53*17519)
216
217
218
219
220 n <- nrow(synchrones)
221 delta <- ncol(synchrones)
222 #TODO: automatic tune of all these parameters ? (for other users)
223 nvoice <- 4
224 # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones))
225 noctave = 13
226 # 4 here represent 2^5 = 32 half-hours ~ 1 day
227 #NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?)
228 scalevector <- 2^(4:(noctave * nvoice) / nvoice + 1)
229 #condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1
230 s0=2
231 w0=2*pi
232 scaled=FALSE
233 s0log = as.integer( (log2( s0*w0/(2*pi) ) - 1) * nvoice + 1.5 )
234 totnoct = noctave + as.integer(s0log/nvoice) + 1
235
236 Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
237
238 # Generate n(n-1)/2 pairs for WER distances computations
239# pairs = list()
240# V = seq_len(n)
241# for (i in 1:n)
242# {
243# V = V[-1]
244# pairs = c(pairs, lapply(V, function(v) c(i,v)))
245# }
246 # Generate "smart" pairs for WER distances computations
247 pairs = list()
248 F = floor(2*n/3)
249 for (i in 1:F)
250 pairs = c(pairs, lapply((i+1):n, function(v) c(i,v)))
251 V = (F+1):n
252 for (i in (F+1):(n-1))
253 {
254 V = V[-1]
255 pairs = c(pairs,
256
257 # Distance between rows i and j
258 computeDistancesIJ = function(pair)
259 {
260 if (parll)
261 {
262 require("bigmemory", quietly=TRUE)
263 require("epclust", quietly=TRUE)
264 synchrones <- bigmemory::attach.big.matrix(synchrones_desc)
265 Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc)
266 }
267
268 computeCWT = function(index)
269 {
270 ts <- scale(ts(synchrones[index,]), center=TRUE, scale=scaled)
271 totts.cwt = Rwave::cwt(ts, totnoct, nvoice, w0, plot=FALSE)
272 ts.cwt = totts.cwt[,s0log:(s0log+noctave*nvoice)]
273 #Normalization
274 sqs <- sqrt(2^(0:(noctave*nvoice)/nvoice)*s0)
275 sqres <- sweep(ts.cwt,2,sqs,'*')
276 sqres / max(Mod(sqres))
277 }
278
279 i = pair[1] ; j = pair[2]
280 if (verbose && j==i+1)
281 cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep=""))
282 cwt_i <- computeCWT(i)
283 cwt_j <- computeCWT(j)
284
285#print(system.time( {
286 num <- epclustFilter(Mod(cwt_i * Conj(cwt_j)))
287 WX <- epclustFilter(Mod(cwt_i * Conj(cwt_i)))
288 WY <- epclustFilter(Mod(cwt_j * Conj(cwt_j)))
289 wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY))
290 Xwer_dist[i,j] <- sqrt(delta * ncol(cwt_i) * max(1 - wer2, 0.)) #FIXME: wer2 should be < 1
291 Xwer_dist[j,i] <- Xwer_dist[i,j]
292#} ) )
293 Xwer_dist[i,i] = 0.
294 }
295
296 if (parll)
297 {
298 cl = parallel::makeCluster(ncores_clust)
299 synchrones_desc <- bigmemory::describe(synchrones)
300 Xwer_dist_desc <- bigmemory::describe(Xwer_dist)
301
302 parallel::clusterExport(cl, varlist=c("synchrones_desc","Xwer_dist_desc","totnoct",
303 "nvoice","w0","s0log","noctave","s0","verbose"), envir=environment())
304 }
305
306 ignored <-
307 if (parll)
308 parallel::parLapply(cl, pairs, computeDistancesIJ)
309 else
310 lapply(pairs, computeDistancesIJ)
311
312 if (parll)
313 parallel::stopCluster(cl)
314
315 Xwer_dist[n,n] = 0.
316 distances <- Xwer_dist[,]
317 rm(Xwer_dist) ; gc()
318 distances #~small matrix K1 x K1
319}
320
321# Helper function to divide indices into balanced sets
322.spreadIndices = function(indices, nb_per_chunk)
323{
324 L = length(indices)
325 nb_workers = floor( L / nb_per_chunk )
326 if (nb_workers == 0)
327 {
328 # L < nb_series_per_chunk, simple case
329 indices_workers = list(indices)
330 }
331 else
332 {
333 indices_workers = lapply( seq_len(nb_workers), function(i)
334 indices[(nb_per_chunk*(i-1)+1):(nb_per_chunk*i)] )
335 # Spread the remaining load among the workers
336 rem = L %% nb_per_chunk
337 while (rem > 0)
338 {
339 index = rem%%nb_workers + 1
340 indices_workers[[index]] = c(indices_workers[[index]], indices[L-rem+1])
341 rem = rem - 1
342 }
343 }
344 indices_workers
345}