draft final form of current package
[epclust.git] / epclust / R / clustering.R
CommitLineData
4bcfdbee
BA
1#' @name clustering
2#' @rdname clustering
492cd9e7 3#' @aliases clusteringTask1 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
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)
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)
29NULL
30
31#' @rdname clustering
32#' @export
492cd9e7
BA
33clusteringTask1 = function(
34 indices, getContribs, K1, nb_series_per_chunk, ncores_clust=1, verbose=FALSE, parll=TRUE)
5c652979 35{
492cd9e7 36 if (verbose)
e161499b 37 cat(paste("*** Clustering task 1 on ",length(indices)," lines\n", sep=""))
4bcfdbee 38
492cd9e7 39 if (parll)
7b13d0c2 40 {
492cd9e7
BA
41 cl = parallel::makeCluster(ncores_clust)
42 parallel::clusterExport(cl, varlist=c("getContribs","K1","verbose"), envir=environment())
7b13d0c2 43 }
492cd9e7
BA
44 while (length(indices) > K1)
45 {
46 indices_workers = .spreadIndices(indices, nb_series_per_chunk)
e161499b
BA
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 }
492cd9e7
BA
61 }
62 if (parll)
63 parallel::stopCluster(cl)
64
56857861 65 indices #medoids
5c652979
BA
66}
67
4bcfdbee
BA
68#' @rdname clustering
69#' @export
a174b8ea
BA
70clusteringTask2 = function(medoids, K2, getRefSeries, nb_ref_curves,
71 nb_series_per_chunk, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE)
5c652979 72{
e161499b
BA
73 if (verbose)
74 cat(paste("*** Clustering task 2 on ",nrow(medoids)," lines\n", sep=""))
75
bf5c0844
BA
76 if (nrow(medoids) <= K2)
77 return (medoids)
492cd9e7
BA
78 synchrones = computeSynchrones(medoids,
79 getRefSeries, nb_ref_curves, nb_series_per_chunk, ncores_clust, verbose, parll)
a174b8ea 80 distances = computeWerDists(synchrones, nbytes, endian, ncores_clust, verbose, parll)
777c4b02 81 medoids[ computeClusters2(distances,K2,verbose), ]
5c652979
BA
82}
83
bf5c0844
BA
84#' @rdname clustering
85#' @export
e161499b
BA
86computeClusters1 = function(contribs, K1, verbose=FALSE)
87{
88 if (verbose)
89 cat(paste(" computeClusters1() on ",nrow(contribs)," lines\n", sep=""))
bf5c0844 90 cluster::pam(contribs, K1, diss=FALSE)$id.med
e161499b 91}
bf5c0844
BA
92
93#' @rdname clustering
94#' @export
e161499b
BA
95computeClusters2 = function(distances, K2, verbose=FALSE)
96{
97 if (verbose)
98 cat(paste(" computeClusters2() on ",nrow(distances)," lines\n", sep=""))
bf5c0844 99 cluster::pam(distances, K2, diss=TRUE)$id.med
e161499b 100}
bf5c0844 101
4bcfdbee
BA
102#' computeSynchrones
103#'
104#' Compute the synchrones curves (sum of clusters elements) from a matrix of medoids,
105#' using L2 distances.
106#'
24ed5d83 107#' @param medoids big.matrix of medoids (curves of same length as initial series)
4bcfdbee
BA
108#' @param getRefSeries Function to retrieve initial series (e.g. in stage 2 after series
109#' have been replaced by stage-1 medoids)
492cd9e7 110#' @param nb_ref_curves How many reference series? (This number is known at this stage)
4bcfdbee
BA
111#' @inheritParams claws
112#'
24ed5d83
BA
113#' @return A big.matrix of size K1 x L where L = data_length
114#'
4bcfdbee 115#' @export
492cd9e7
BA
116computeSynchrones = function(medoids, getRefSeries,
117 nb_ref_curves, nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE)
e205f218 118{
e161499b
BA
119 if (verbose)
120 cat(paste("--- Compute synchrones\n", sep=""))
121
492cd9e7 122 computeSynchronesChunk = function(indices)
3eef8d3d 123 {
363ae134
BA
124 if (parll)
125 {
126 require("bigmemory", quietly=TRUE)
6ad3f3fd 127 requireNamespace("synchronicity", quietly=TRUE)
363ae134
BA
128 require("epclust", quietly=TRUE)
129 synchrones <- bigmemory::attach.big.matrix(synchrones_desc)
6ad3f3fd 130 counts <- bigmemory::attach.big.matrix(counts_desc)
363ae134
BA
131 medoids <- bigmemory::attach.big.matrix(medoids_desc)
132 m <- synchronicity::attach.mutex(m_desc)
133 }
134
6ad3f3fd
BA
135 ref_series = getRefSeries(indices)
136 nb_series = nrow(ref_series)
137
363ae134 138 #get medoids indices for this chunk of series
2c14dbea 139 mi = computeMedoidsIndices(medoids@address, ref_series)
e161499b
BA
140
141 for (i in seq_len(nb_series))
56857861 142 {
492cd9e7
BA
143 if (parll)
144 synchronicity::lock(m)
6ad3f3fd
BA
145 synchrones[ mi[i], ] = synchrones[ mi[i], ] + ref_series[i,]
146 counts[ mi[i] ] = counts[ mi[i] ] + 1 #TODO: remove counts?
492cd9e7
BA
147 if (parll)
148 synchronicity::unlock(m)
149 }
150 }
151
e161499b 152 K = nrow(medoids) ; L = ncol(medoids)
492cd9e7 153 # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in //
24ed5d83 154 # TODO: if size > RAM (not our case), use file-backed big.matrix
e161499b
BA
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)
24ed5d83 157 # synchronicity is only for Linux & MacOS; on Windows: run sequentially
492cd9e7
BA
158 parll = (requireNamespace("synchronicity",quietly=TRUE)
159 && parll && Sys.info()['sysname'] != "Windows")
160 if (parll)
363ae134 161 {
492cd9e7 162 m <- synchronicity::boost.mutex()
363ae134
BA
163 m_desc <- synchronicity::describe(m)
164 synchrones_desc = bigmemory::describe(synchrones)
6ad3f3fd 165 counts_desc = bigmemory::describe(counts)
363ae134 166 medoids_desc = bigmemory::describe(medoids)
24ed5d83 167 cl = parallel::makeCluster(ncores_clust)
6ad3f3fd
BA
168 parallel::clusterExport(cl, varlist=c("synchrones_desc","counts_desc","counts",
169 "verbose","m_desc","medoids_desc","getRefSeries"), envir=environment())
24ed5d83
BA
170 }
171
492cd9e7 172 indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk)
c45fd663 173 ignored <-
492cd9e7 174 if (parll)
e161499b 175 parallel::parLapply(cl, indices_workers, computeSynchronesChunk)
492cd9e7 176 else
c45fd663 177 lapply(indices_workers, computeSynchronesChunk)
492cd9e7 178
24ed5d83
BA
179 if (parll)
180 parallel::stopCluster(cl)
181
182 #TODO: can we avoid this loop? ( synchrones = sweep(synchrones, 1, counts, '/') )
492cd9e7 183 for (i in seq_len(K))
24ed5d83 184 synchrones[i,] = synchrones[i,] / counts[i,1]
3eef8d3d 185 #NOTE: odds for some clusters to be empty? (when series already come from stage 2)
8702eb86 186 # ...maybe; but let's hope resulting K1' be still quite bigger than K2
24ed5d83
BA
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,]
e205f218 192}
1c6f223e 193
4bcfdbee
BA
194#' computeWerDists
195#'
196#' Compute the WER distances between the synchrones curves (in rows), which are
197#' returned (e.g.) by \code{computeSynchrones()}
198#'
24ed5d83
BA
199#' @param synchrones A big.matrix of synchrones, in rows. The series have same length
200#' as the series in the initial dataset
492cd9e7 201#' @inheritParams claws
4bcfdbee 202#'
777c4b02 203#' @return A matrix of size K1 x K1
24ed5d83 204#'
4bcfdbee 205#' @export
a174b8ea 206computeWerDists = function(synchrones, nbytes,endian,ncores_clust=1,verbose=FALSE,parll=TRUE)
d03c0621 207{
e161499b
BA
208 if (verbose)
209 cat(paste("--- Compute WER dists\n", sep=""))
24ed5d83 210
4bcfdbee
BA
211 n <- nrow(synchrones)
212 delta <- ncol(synchrones)
db6fc17d 213 #TODO: automatic tune of all these parameters ? (for other users)
d03c0621 214 nvoice <- 4
4bcfdbee 215 # noctave = 2^13 = 8192 half hours ~ 180 days ; ~log2(ncol(synchrones))
d7d55bc1
BA
216 noctave = 13
217 # 4 here represent 2^5 = 32 half-hours ~ 1 day
db6fc17d 218 #NOTE: default scalevector == 2^(0:(noctave * nvoice) / nvoice) * s0 (?)
24ed5d83 219 scalevector <- 2^(4:(noctave * nvoice) / nvoice + 1)
db6fc17d 220 #condition: ( log2(s0*w0/(2*pi)) - 1 ) * nvoice + 1.5 >= 1
a174b8ea
BA
221 s0 = 2
222 w0 = 2*pi
db6fc17d
BA
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
e161499b 227 Xwer_dist <- bigmemory::big.matrix(nrow=n, ncol=n, type="double")
e161499b 228
4204e877
BA
229 cwt_file = ".epclust_bin/cwt"
230 #TODO: args, nb_per_chunk, nbytes, endian
231
e161499b
BA
232 # Generate n(n-1)/2 pairs for WER distances computations
233 pairs = list()
4204e877
BA
234 V = seq_len(n)
235 for (i in 1:n)
e161499b
BA
236 {
237 V = V[-1]
4204e877
BA
238 pairs = c(pairs, lapply(V, function(v) c(i,v)))
239 }
a174b8ea 240
4204e877
BA
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)
a174b8ea 252 binarize(c(as.double(Re(res)),as.double(Im(res))), cwt_file, ncol(res), ",", nbytes, endian)
4204e877
BA
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 ...
a174b8ea
BA
274 res <- getDataInFile(c(2*index-1,2*index), cwt_file, nbytes, endian)
275 ###############TODO:
4204e877 276 }
e161499b 277
777c4b02 278 # Distance between rows i and j
e161499b
BA
279 computeDistancesIJ = function(pair)
280 {
2c14dbea 281 if (parll)
363ae134 282 {
2c14dbea
BA
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
e161499b
BA
289 i = pair[1] ; j = pair[2]
290 if (verbose && j==i+1)
291 cat(paste(" Distances (",i,",",j,"), (",i,",",j+1,") ...\n", sep=""))
4204e877
BA
292 cwt_i <- getCWT(i)
293 cwt_j <- getCWT(j)
2c14dbea 294
363ae134
BA
295 num <- epclustFilter(Mod(cwt_i * Conj(cwt_j)))
296 WX <- epclustFilter(Mod(cwt_i * Conj(cwt_i)))
4204e877 297 WY <- epclustFilter(Mod(cwt_j * Conj(cwt_j)))
e161499b 298 wer2 <- sum(colSums(num)^2) / sum(colSums(WX) * colSums(WY))
2c14dbea 299 Xwer_dist[i,j] <- sqrt(delta * ncol(cwt_i) * max(1 - wer2, 0.)) #FIXME: wer2 should be < 1
e161499b
BA
300 Xwer_dist[j,i] <- Xwer_dist[i,j]
301 Xwer_dist[i,i] = 0.
302 }
303
e161499b 304 ignored <-
492cd9e7 305 if (parll)
e161499b 306 parallel::parLapply(cl, pairs, computeDistancesIJ)
492cd9e7 307 else
e161499b 308 lapply(pairs, computeDistancesIJ)
492cd9e7
BA
309
310 if (parll)
311 parallel::stopCluster(cl)
6ad3f3fd 312
492cd9e7 313 Xwer_dist[n,n] = 0.
777c4b02
BA
314 distances <- Xwer_dist[,]
315 rm(Xwer_dist) ; gc()
316 distances #~small matrix K1 x K1
492cd9e7
BA
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
d03c0621 340 }
1c6f223e 341 }
492cd9e7 342 indices_workers
1c6f223e 343}