+ computeSynchronesChunk = function(indices)
+ {
+ if (parll)
+ {
+ require("bigmemory", quietly=TRUE)
+ requireNamespace("synchronicity", quietly=TRUE)
+ require("epclust", quietly=TRUE)
+ synchrones <- bigmemory::attach.big.matrix(synchrones_desc)
+ if (sync_mean)
+ counts <- bigmemory::attach.big.matrix(counts_desc)
+ medoids <- bigmemory::attach.big.matrix(medoids_desc)
+ m <- synchronicity::attach.mutex(m_desc)
+ }
+
+ ref_series = getRefSeries(indices)
+ nb_series = ncol(ref_series)
+
+ # Get medoids indices for this chunk of series
+ mi = computeMedoidsIndices(medoids@address, ref_series)
+
+ for (i in seq_len(nb_series))
+ {
+ if (parll)
+ synchronicity::lock(m)
+ synchrones[, mi[i] ] = synchrones[, mi[i] ] + ref_series[,i]
+ if (sync_mean)
+ counts[ mi[i] ] = counts[ mi[i] ] + 1
+ if (parll)
+ synchronicity::unlock(m)
+ }
+ }
+
+ K = ncol(medoids) ; L = nrow(medoids)
+ # Use bigmemory (shared==TRUE by default) + synchronicity to fill synchrones in //
+ # TODO: if size > RAM (not our case), use file-backed big.matrix
+ synchrones = bigmemory::big.matrix(nrow=L, ncol=K, type="double", init=0.)
+ if (sync_mean)
+ counts = bigmemory::big.matrix(nrow=K, ncol=1, type="double", init=0)
+ # synchronicity is only for Linux & MacOS; on Windows: run sequentially
+ parll = (requireNamespace("synchronicity",quietly=TRUE)
+ && parll && Sys.info()['sysname'] != "Windows")
+ if (parll)
+ {
+ m <- synchronicity::boost.mutex()
+ m_desc <- synchronicity::describe(m)
+ synchrones_desc = bigmemory::describe(synchrones)
+ if (sync_mean)
+ counts_desc = bigmemory::describe(counts)
+ medoids_desc = bigmemory::describe(medoids)
+ cl = parallel::makeCluster(ncores_clust)
+ varlist=c("synchrones_desc","sync_mean","m_desc","medoids_desc","getRefSeries")
+ if (sync_mean)
+ varlist = c(varlist, "counts_desc")
+ parallel::clusterExport(cl, varlist, envir=environment())
+ }
+
+ if (verbose)
+ {
+ if (verbose)
+ cat(paste("--- Compute ",K," synchrones with ",nb_ref_curves," series\n", sep=""))
+ }
+ indices_workers = .spreadIndices(seq_len(nb_ref_curves), nb_series_per_chunk)
+ ignored <-
+ if (parll)
+ parallel::parLapply(cl, indices_workers, computeSynchronesChunk)
+ else
+ lapply(indices_workers, computeSynchronesChunk)
+
+ if (parll)
+ parallel::stopCluster(cl)
+
+ if (!sync_mean)
+ return (synchrones)
+
+ #TODO: can we avoid this loop? ( synchrones = sweep(synchrones, 2, counts, '/') )
+ for (i in seq_len(K))
+ synchrones[,i] = synchrones[,i] / counts[i]
+ #NOTE: odds for some clusters to be empty? (when series already come from stage 2)
+ # ...maybe; but let's hope resulting K1' be still quite bigger than K2
+ noNA_rows = sapply(seq_len(K), function(i) all(!is.nan(synchrones[,i])))
+ if (all(noNA_rows))
+ return (synchrones)
+ # Else: some clusters are empty, need to slice synchrones
+ bigmemory::as.big.matrix(synchrones[,noNA_rows])