3 #' @aliases clusteringTask1 clusteringTask2 computeClusters1 computeClusters2
5 #' @title Two-stage clustering, withing one task (see \code{claws()})
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 #' first clustering algorithm on a contributions matrix, while the latter clusters
15 #' a set of series inside one task (~nb_items_clust1)
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 #' @inheritParams computeSynchrones
21 #' @inheritParams claws
23 #' @return For \code{clusteringTask1()}, the indices of the computed (K1) medoids.
24 #' Indices are irrelevant for stage 2 clustering, thus \code{clusteringTask2()}
25 #' outputs a big.matrix of medoids (of size LxK2, K2 = final number of clusters)
30 clusteringTask1 = function(indices, getContribs, K1, algoClust1, nb_items_clust1,
31 ncores_clust=1, verbose=FALSE, parll=TRUE)
35 cl = parallel::makeCluster(ncores_clust, outfile = "")
36 parallel::clusterExport(cl, c("getContribs","K1","verbose"), envir=environment())
38 # Iterate clustering algorithm 1 until K1 medoids are found
39 while (length(indices) > K1)
41 # Balance tasks by splitting the indices set - as evenly as possible
42 indices_workers = .spreadIndices(indices, nb_items_clust1)
44 cat(paste("*** [iterated] Clustering task 1 on ",length(indices)," series\n", sep=""))
48 unlist( parallel::parLapply(cl, indices_workers, function(inds) {
49 require("epclust", quietly=TRUE)
50 inds[ algoClust1(getContribs(inds), K1) ]
55 unlist( lapply(indices_workers, function(inds)
56 inds[ algoClust1(getContribs(inds), K1) ]
61 parallel::stopCluster(cl)
68 clusteringTask2 = function(medoids, K2, algoClust2, getRefSeries, nb_ref_curves,
69 nb_series_per_chunk, nvoice, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE)
72 cat(paste("*** Clustering task 2 on ",ncol(medoids)," synchrones\n", sep=""))
74 if (ncol(medoids) <= K2)
77 # A) Obtain synchrones, that is to say the cumulated power consumptions
78 # for each of the K1 initial groups
79 synchrones = computeSynchrones(medoids, getRefSeries, nb_ref_curves,
80 nb_series_per_chunk, ncores_clust, verbose, parll)
82 # B) Compute the WER distances (Wavelets Extended coefficient of deteRmination)
83 distances = computeWerDists(
84 synchrones, nvoice, nbytes, endian, ncores_clust, verbose, parll)
86 # C) Apply clustering algorithm 2 on the WER distances matrix
88 cat(paste("*** algoClust2() on ",nrow(distances)," items\n", sep=""))
89 medoids[ ,algoClust2(distances,K2) ]
94 #' Compute the synchrones curves (sum of clusters elements) from a matrix of medoids,
95 #' using euclidian distance.
97 #' @param medoids big.matrix of medoids (curves of same length as initial series)
98 #' @param getRefSeries Function to retrieve initial series (e.g. in stage 2 after series
99 #' have been replaced by stage-1 medoids)
100 #' @param nb_ref_curves How many reference series? (This number is known at this stage)
101 #' @inheritParams claws
103 #' @return A big.matrix of size L x K1 where L = length of a serie
106 computeSynchrones = function(medoids, getRefSeries, nb_ref_curves,
107 nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE)
109 # Synchrones computation is embarassingly parallel: compute it by chunks of series
110 computeSynchronesChunk = function(indices)
114 require("bigmemory", quietly=TRUE)
115 requireNamespace("synchronicity", quietly=TRUE)
116 require("epclust", quietly=TRUE)
117 # The big.matrix objects need to be attached to be usable on the workers
118 synchrones <- bigmemory::attach.big.matrix(synchrones_desc)
119 medoids <- bigmemory::attach.big.matrix(medoids_desc)
120 m <- synchronicity::attach.mutex(m_desc)
123 # Obtain a chunk of reference series
124 ref_series = getRefSeries(indices)
125 nb_series = ncol(ref_series)
127 # Get medoids indices for this chunk of series
128 mi = computeMedoidsIndices(medoids@address, ref_series)
130 # Update synchrones using mi above
131 for (i in seq_len(nb_series))
134 synchronicity::lock(m) #locking required because several writes at the same time
135 synchrones[, mi[i] ] = synchrones[, mi[i] ] + ref_series[,i]
137 synchronicity::unlock(m)
142 K = ncol(medoids) ; L = nrow(medoids)
143 # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in //
144 synchrones = bigmemory::big.matrix(nrow=L, ncol=K, type="double", init=0.)
145 # NOTE: synchronicity is only for Linux & MacOS; on Windows: run sequentially
146 parll = (parll && requireNamespace("synchronicity",quietly=TRUE)
147 && Sys.info()['sysname'] != "Windows")
150 m <- synchronicity::boost.mutex() #for lock/unlock, see computeSynchronesChunk
151 # mutex and big.matrix objects cannot be passed directly:
152 # they will be accessed from their description
153 m_desc <- synchronicity::describe(m)
154 synchrones_desc = bigmemory::describe(synchrones)
155 medoids_desc = bigmemory::describe(medoids)
156 cl = parallel::makeCluster(ncores_clust)
157 parallel::clusterExport(cl, envir=environment(),
158 varlist=c("synchrones_desc","m_desc","medoids_desc","getRefSeries"))
162 cat(paste("--- Compute ",K," synchrones with ",nb_ref_curves," series\n", sep=""))
164 # Balance tasks by splitting the indices set - maybe not so evenly, but
165 # max==TRUE in next call ensures that no set has more than nb_series_per_chunk items.
166 indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk, max=TRUE)
169 parallel::parLapply(cl, indices_workers, computeSynchronesChunk)
171 lapply(indices_workers, computeSynchronesChunk)
174 parallel::stopCluster(cl)
181 #' Compute the WER distances between the synchrones curves (in columns), which are
182 #' returned (e.g.) by \code{computeSynchrones()}
184 #' @param synchrones A big.matrix of synchrones, in columns. The series have same
185 #' length as the series in the initial dataset
186 #' @inheritParams claws
188 #' @return A distances matrix of size K1 x K1
191 computeWerDists = function(synchrones, nvoice, nbytes,endian,ncores_clust=1,
192 verbose=FALSE,parll=TRUE)
194 n <- ncol(synchrones)
195 L <- nrow(synchrones)
196 noctave = ceiling(log2(L)) #min power of 2 to cover serie range
198 # Initialize result as a square big.matrix of size 'number of synchrones'
199 Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
201 # Generate n(n-1)/2 pairs for WER distances computations
207 pairs = c(pairs, lapply(V, function(v) c(i,v)))
210 cwt_file = ".cwt.bin"
211 # Compute the synchrones[,index] CWT, and store it in the binary file above
212 computeSaveCWT = function(index)
214 if (parll && !exists(synchrones)) #avoid going here after first call on a worker
216 require("bigmemory", quietly=TRUE)
217 require("Rwave", quietly=TRUE)
218 require("epclust", quietly=TRUE)
219 synchrones <- bigmemory::attach.big.matrix(synchrones_desc)
221 ts <- scale(ts(synchrones[,index]), center=TRUE, scale=FALSE)
222 ts_cwt = Rwave::cwt(ts, noctave, nvoice, w0=2*pi, twoD=TRUE, plot=FALSE)
225 binarize(as.matrix(c(as.double(Re(ts_cwt)),as.double(Im(ts_cwt)))), cwt_file, 1,
231 cl = parallel::makeCluster(ncores_clust)
232 synchrones_desc <- bigmemory::describe(synchrones)
233 Xwer_dist_desc <- bigmemory::describe(Xwer_dist)
234 parallel::clusterExport(cl, varlist=c("parll","synchrones_desc","Xwer_dist_desc",
235 "noctave","nvoice","verbose","getCWT"), envir=environment())
239 cat(paste("--- Precompute and serialize synchrones CWT\n", sep=""))
243 parallel::parLapply(cl, 1:n, computeSaveCWT)
245 lapply(1:n, computeSaveCWT)
247 # Function to retrieve a synchrone CWT from (binary) file
248 getSynchroneCWT = function(index, L)
250 flat_cwt <- getDataInFile(index, cwt_file, nbytes, endian)
251 cwt_length = length(flat_cwt) / 2
252 re_part = as.matrix(flat_cwt[1:cwt_length], nrow=L)
253 im_part = as.matrix(flat_cwt[(cwt_length+1):(2*cwt_length)], nrow=L)
254 re_part + 1i * im_part
257 # Compute distance between columns i and j in synchrones
258 computeDistanceIJ = function(pair)
262 # parallel workers start with an empty environment
263 require("bigmemory", quietly=TRUE)
264 require("epclust", quietly=TRUE)
265 synchrones <- bigmemory::attach.big.matrix(synchrones_desc)
266 Xwer_dist <- bigmemory::attach.big.matrix(Xwer_dist_desc)
269 i = pair[1] ; j = pair[2]
270 if (verbose && j==i+1 && !parll)
271 cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep=""))
273 # Compute CWT of columns i and j in synchrones
275 cwt_i <- getSynchroneCWT(i, L)
276 cwt_j <- getSynchroneCWT(j, L)
278 # Compute the ratio of integrals formula 5.6 for WER^2
279 # in https://arxiv.org/abs/1101.4744v2 ยง5.3
280 num <- filterMA(Mod(cwt_i * Conj(cwt_j)))
281 WX <- filterMA(Mod(cwt_i * Conj(cwt_i)))
282 WY <- filterMA(Mod(cwt_j * Conj(cwt_j)))
283 wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY))
285 Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2))
286 Xwer_dist[j,i] <- Xwer_dist[i,j]
291 cat(paste("--- Compute WER distances\n", sep=""))
295 parallel::parLapply(cl, pairs, computeDistanceIJ)
297 lapply(pairs, computeDistanceIJ)
300 parallel::stopCluster(cl)
305 Xwer_dist[,] #~small matrix K1 x K1
308 # Helper function to divide indices into balanced sets
309 # If max == TRUE, sets sizes cannot exceed nb_per_set
310 .spreadIndices = function(indices, nb_per_set, max=FALSE)
313 nb_workers = floor( L / nb_per_set )
314 rem = L %% nb_per_set
315 if (nb_workers == 0 || (nb_workers==1 && rem==0))
317 # L <= nb_per_set, simple case
318 indices_workers = list(indices)
322 indices_workers = lapply( seq_len(nb_workers), function(i)
323 indices[(nb_per_set*(i-1)+1):(nb_per_set*i)] )
327 # Sets are not so well balanced, but size is supposed to be critical
328 return ( c( indices_workers, if (rem>0) list((L-rem+1):L) else NULL ) )
331 # Spread the remaining load among the workers
332 rem = L %% nb_per_set
335 index = rem%%nb_workers + 1
336 indices_workers[[index]] = c(indices_workers[[index]], indices[L-rem+1])