| 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) |
| 29 | NULL |
| 30 | |
| 31 | #' @rdname clustering |
| 32 | #' @export |
| 33 | clusteringTask1 = 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 |
| 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(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 K1 x L where L = data_length |
| 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? |
| 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 |
| 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 | } |