| 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 | #' first clustering algorithm on a contributions matrix, while the latter clusters |
| 15 | #' a set of series inside one task (~nb_items_clust1) |
| 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 | #' @inheritParams computeSynchrones |
| 21 | #' @inheritParams claws |
| 22 | #' |
| 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) |
| 26 | NULL |
| 27 | |
| 28 | #' @rdname clustering |
| 29 | #' @export |
| 30 | clusteringTask1 = function(indices, getContribs, K1, algoClust1, nb_items_clust1, |
| 31 | ncores_clust=1, verbose=FALSE, parll=TRUE) |
| 32 | { |
| 33 | if (parll) |
| 34 | { |
| 35 | cl = parallel::makeCluster(ncores_clust, outfile = "") |
| 36 | parallel::clusterExport(cl, c("getContribs","K1","verbose"), envir=environment()) |
| 37 | } |
| 38 | # Iterate clustering algorithm 1 until K1 medoids are found |
| 39 | while (length(indices) > K1) |
| 40 | { |
| 41 | # Balance tasks by splitting the indices set - as evenly as possible |
| 42 | indices_workers = .spreadIndices(indices, nb_items_clust1) |
| 43 | if (verbose) |
| 44 | cat(paste("*** [iterated] Clustering task 1 on ",length(indices)," series\n", sep="")) |
| 45 | indices <- |
| 46 | if (parll) |
| 47 | { |
| 48 | unlist( parallel::parLapply(cl, indices_workers, function(inds) { |
| 49 | require("epclust", quietly=TRUE) |
| 50 | inds[ algoClust1(getContribs(inds), K1) ] |
| 51 | }) ) |
| 52 | } |
| 53 | else |
| 54 | { |
| 55 | unlist( lapply(indices_workers, function(inds) |
| 56 | inds[ algoClust1(getContribs(inds), K1) ] |
| 57 | ) ) |
| 58 | } |
| 59 | } |
| 60 | if (parll) |
| 61 | parallel::stopCluster(cl) |
| 62 | |
| 63 | indices #medoids |
| 64 | } |
| 65 | |
| 66 | #' @rdname clustering |
| 67 | #' @export |
| 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) |
| 70 | { |
| 71 | if (verbose) |
| 72 | cat(paste("*** Clustering task 2 on ",ncol(medoids)," synchrones\n", sep="")) |
| 73 | |
| 74 | if (ncol(medoids) <= K2) |
| 75 | return (medoids) |
| 76 | |
| 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) |
| 81 | |
| 82 | # B) Compute the WER distances (Wavelets Extended coefficient of deteRmination) |
| 83 | distances = computeWerDists( |
| 84 | synchrones, nvoice, nbytes, endian, ncores_clust, verbose, parll) |
| 85 | |
| 86 | # C) Apply clustering algorithm 2 on the WER distances matrix |
| 87 | if (verbose) |
| 88 | cat(paste("*** algoClust2() on ",nrow(distances)," items\n", sep="")) |
| 89 | medoids[ ,algoClust2(distances,K2) ] |
| 90 | } |
| 91 | |
| 92 | #' computeSynchrones |
| 93 | #' |
| 94 | #' Compute the synchrones curves (sum of clusters elements) from a matrix of medoids, |
| 95 | #' using euclidian distance. |
| 96 | #' |
| 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 |
| 102 | #' |
| 103 | #' @return A big.matrix of size L x K1 where L = length of a serie |
| 104 | #' |
| 105 | #' @export |
| 106 | computeSynchrones = function(medoids, getRefSeries, nb_ref_curves, |
| 107 | nb_series_per_chunk, ncores_clust=1,verbose=FALSE,parll=TRUE) |
| 108 | { |
| 109 | # Synchrones computation is embarassingly parallel: compute it by chunks of series |
| 110 | computeSynchronesChunk = function(indices) |
| 111 | { |
| 112 | if (parll) |
| 113 | { |
| 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) |
| 121 | } |
| 122 | |
| 123 | # Obtain a chunk of reference series |
| 124 | ref_series = getRefSeries(indices) |
| 125 | nb_series = ncol(ref_series) |
| 126 | |
| 127 | # Get medoids indices for this chunk of series |
| 128 | mi = computeMedoidsIndices(medoids@address, ref_series) |
| 129 | |
| 130 | # Update synchrones using mi above |
| 131 | for (i in seq_len(nb_series)) |
| 132 | { |
| 133 | if (parll) |
| 134 | synchronicity::lock(m) #locking required because several writes at the same time |
| 135 | synchrones[, mi[i] ] = synchrones[, mi[i] ] + ref_series[,i] |
| 136 | if (parll) |
| 137 | synchronicity::unlock(m) |
| 138 | } |
| 139 | NULL |
| 140 | } |
| 141 | |
| 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") |
| 148 | if (parll) |
| 149 | { |
| 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")) |
| 159 | } |
| 160 | |
| 161 | if (verbose) |
| 162 | cat(paste("--- Compute ",K," synchrones with ",nb_ref_curves," series\n", sep="")) |
| 163 | |
| 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) |
| 167 | ignored <- |
| 168 | if (parll) |
| 169 | parallel::parLapply(cl, indices_workers, computeSynchronesChunk) |
| 170 | else |
| 171 | lapply(indices_workers, computeSynchronesChunk) |
| 172 | |
| 173 | if (parll) |
| 174 | parallel::stopCluster(cl) |
| 175 | |
| 176 | return (synchrones) |
| 177 | } |
| 178 | |
| 179 | #' computeWerDists |
| 180 | #' |
| 181 | #' Compute the WER distances between the synchrones curves (in columns), which are |
| 182 | #' returned (e.g.) by \code{computeSynchrones()} |
| 183 | #' |
| 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 |
| 187 | #' |
| 188 | #' @return A distances matrix of size K1 x K1 |
| 189 | #' |
| 190 | #' @export |
| 191 | computeWerDists = function(synchrones, nvoice, nbytes,endian,ncores_clust=1, |
| 192 | verbose=FALSE,parll=TRUE) |
| 193 | { |
| 194 | n <- ncol(synchrones) |
| 195 | L <- nrow(synchrones) |
| 196 | noctave = ceiling(log2(L)) #min power of 2 to cover serie range |
| 197 | |
| 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") |
| 200 | |
| 201 | # Generate n(n-1)/2 pairs for WER distances computations |
| 202 | pairs = list() |
| 203 | V = seq_len(n) |
| 204 | for (i in 1:n) |
| 205 | { |
| 206 | V = V[-1] |
| 207 | pairs = c(pairs, lapply(V, function(v) c(i,v))) |
| 208 | } |
| 209 | |
| 210 | cwt_file = ".cwt.bin" |
| 211 | # Compute the synchrones[,index] CWT, and store it in the binary file above |
| 212 | computeSaveCWT = function(index) |
| 213 | { |
| 214 | if (parll && !exists(synchrones)) #avoid going here after first call on a worker |
| 215 | { |
| 216 | require("bigmemory", quietly=TRUE) |
| 217 | require("Rwave", quietly=TRUE) |
| 218 | require("epclust", quietly=TRUE) |
| 219 | synchrones <- bigmemory::attach.big.matrix(synchrones_desc) |
| 220 | } |
| 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) |
| 223 | |
| 224 | # Serialization |
| 225 | binarize(as.matrix(c(as.double(Re(ts_cwt)),as.double(Im(ts_cwt)))), cwt_file, 1, |
| 226 | ",", nbytes, endian) |
| 227 | } |
| 228 | |
| 229 | if (parll) |
| 230 | { |
| 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()) |
| 236 | } |
| 237 | |
| 238 | if (verbose) |
| 239 | cat(paste("--- Precompute and serialize synchrones CWT\n", sep="")) |
| 240 | |
| 241 | ignored <- |
| 242 | if (parll) |
| 243 | parallel::parLapply(cl, 1:n, computeSaveCWT) |
| 244 | else |
| 245 | lapply(1:n, computeSaveCWT) |
| 246 | |
| 247 | # Function to retrieve a synchrone CWT from (binary) file |
| 248 | getSynchroneCWT = function(index, L) |
| 249 | { |
| 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 |
| 255 | } |
| 256 | |
| 257 | # Compute distance between columns i and j in synchrones |
| 258 | computeDistanceIJ = function(pair) |
| 259 | { |
| 260 | if (parll) |
| 261 | { |
| 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) |
| 267 | } |
| 268 | |
| 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="")) |
| 272 | |
| 273 | # Compute CWT of columns i and j in synchrones |
| 274 | L = nrow(synchrones) |
| 275 | cwt_i <- getSynchroneCWT(i, L) |
| 276 | cwt_j <- getSynchroneCWT(j, L) |
| 277 | |
| 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)) |
| 284 | |
| 285 | Xwer_dist[i,j] <- sqrt(L * ncol(cwt_i) * (1 - wer2)) |
| 286 | Xwer_dist[j,i] <- Xwer_dist[i,j] |
| 287 | Xwer_dist[i,i] <- 0. |
| 288 | } |
| 289 | |
| 290 | if (verbose) |
| 291 | cat(paste("--- Compute WER distances\n", sep="")) |
| 292 | |
| 293 | ignored <- |
| 294 | if (parll) |
| 295 | parallel::parLapply(cl, pairs, computeDistanceIJ) |
| 296 | else |
| 297 | lapply(pairs, computeDistanceIJ) |
| 298 | |
| 299 | if (parll) |
| 300 | parallel::stopCluster(cl) |
| 301 | |
| 302 | unlink(cwt_file) |
| 303 | |
| 304 | Xwer_dist[n,n] = 0. |
| 305 | Xwer_dist[,] #~small matrix K1 x K1 |
| 306 | } |
| 307 | |
| 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) |
| 311 | { |
| 312 | L = length(indices) |
| 313 | nb_workers = floor( L / nb_per_set ) |
| 314 | rem = L %% nb_per_set |
| 315 | if (nb_workers == 0 || (nb_workers==1 && rem==0)) |
| 316 | { |
| 317 | # L <= nb_per_set, simple case |
| 318 | indices_workers = list(indices) |
| 319 | } |
| 320 | else |
| 321 | { |
| 322 | indices_workers = lapply( seq_len(nb_workers), function(i) |
| 323 | indices[(nb_per_set*(i-1)+1):(nb_per_set*i)] ) |
| 324 | |
| 325 | if (max) |
| 326 | { |
| 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 ) ) |
| 329 | } |
| 330 | |
| 331 | # Spread the remaining load among the workers |
| 332 | rem = L %% nb_per_set |
| 333 | while (rem > 0) |
| 334 | { |
| 335 | index = rem%%nb_workers + 1 |
| 336 | indices_workers[[index]] = c(indices_workers[[index]], indices[L-rem+1]) |
| 337 | rem = rem - 1 |
| 338 | } |
| 339 | } |
| 340 | indices_workers |
| 341 | } |