| 1 | #' CLAWS: CLustering with wAvelets and Wer distanceS |
| 2 | #' |
| 3 | #' Cluster electricity power curves (or any series of similar nature) by applying a |
| 4 | #' two stage procedure in parallel (see details). |
| 5 | #' Input series must be sampled on the same time grid, no missing values. |
| 6 | #' |
| 7 | #' @details Summary of the function execution flow: |
| 8 | #' \enumerate{ |
| 9 | #' \item Compute and serialize all contributions, obtained through discrete wavelet |
| 10 | #' decomposition (see Antoniadis & al. [2013]) |
| 11 | #' \item Divide series into \code{ntasks} groups to process in parallel. In each task: |
| 12 | #' \enumerate{ |
| 13 | #' \item iterate the first clustering algorithm on its aggregated outputs, |
| 14 | #' on inputs of size \code{nb_items_clust} |
| 15 | #' \item optionally, if WER=="mix": |
| 16 | #' a) compute the K1 synchrones curves, |
| 17 | #' b) compute WER distances (K1xK1 matrix) between synchrones and |
| 18 | #' c) apply the second clustering algorithm |
| 19 | #' } |
| 20 | #' \item Launch a final task on the aggregated outputs of all previous tasks: |
| 21 | #' in the case WER=="end" this task takes indices in input, otherwise |
| 22 | #' (medoid) curves |
| 23 | #' } |
| 24 | #' The main argument -- \code{getSeries} -- has a quite misleading name, since it can be |
| 25 | #' either a [big.]matrix, a CSV file, a connection or a user function to retrieve |
| 26 | #' series; the name was chosen because all types of arguments are converted to a function. |
| 27 | #' When \code{getSeries} is given as a function, it must take a single argument, |
| 28 | #' 'indices', integer vector equal to the indices of the curves to retrieve; |
| 29 | #' see SQLite example. The nature and role of other arguments should be clear |
| 30 | #' |
| 31 | #' @param getSeries Access to the (time-)series, which can be of one of the three |
| 32 | #' following types: |
| 33 | #' \itemize{ |
| 34 | #' \item [big.]matrix: each column contains the (time-ordered) values of one time-serie |
| 35 | #' \item connection: any R connection object providing lines as described above |
| 36 | #' \item character: name of a CSV file containing series in rows (no header) |
| 37 | #' \item function: a custom way to retrieve the curves; it has only one argument: |
| 38 | #' the indices of the series to be retrieved. See SQLite example |
| 39 | #' } |
| 40 | #' @param K1 Number of clusters to be found after stage 1 (K1 << N [number of series]) |
| 41 | #' @param K2 Number of clusters to be found after stage 2 (K2 << K1) |
| 42 | #' @param nb_per_chunk (Maximum) number of items to retrieve in one batch, for both types of |
| 43 | #' retrieval: resp. series and contribution; in a vector of size 2 |
| 44 | #' @param nb_items_clust1 (Maximum) number of items in input of the clustering algorithm |
| 45 | #' for stage 1 |
| 46 | #' @param wav_filt Wavelet transform filter; see ?wavelets::wt.filter |
| 47 | #' @param contrib_type Type of contribution: "relative", "logit" or "absolute" (any prefix) |
| 48 | #' @param WER "end" to apply stage 2 after stage 1 has fully iterated, or "mix" to apply |
| 49 | #' stage 2 at the end of each task |
| 50 | #' @param random TRUE (default) for random chunks repartition |
| 51 | #' @param ntasks Number of tasks (parallel iterations to obtain K1 [if WER=="end"] |
| 52 | #' or K2 [if WER=="mix"] medoids); default: 1. |
| 53 | #' Note: ntasks << N (number of series), so that N is "roughly divisible" by ntasks |
| 54 | #' @param ncores_tasks Number of parallel tasks (1 to disable: sequential tasks) |
| 55 | #' @param ncores_clust Number of parallel clusterings in one task (4 should be a minimum) |
| 56 | #' @param sep Separator in CSV input file (if any provided) |
| 57 | #' @param nbytes Number of bytes to serialize a floating-point number; 4 or 8 |
| 58 | #' @param endian Endianness for (de)serialization ("little" or "big") |
| 59 | #' @param verbose Level of verbosity (0/FALSE for nothing or 1/TRUE for all; devel stage) |
| 60 | #' @param parll TRUE to fully parallelize; otherwise run sequentially (debug, comparison) |
| 61 | #' |
| 62 | #' @return A matrix of the final K2 medoids curves, in columns |
| 63 | #' |
| 64 | #' @references Clustering functional data using Wavelets [2013]; |
| 65 | #' A. Antoniadis, X. Brossat, J. Cugliari & J.-M. Poggi. |
| 66 | #' Inter. J. of Wavelets, Multiresolution and Information Procesing, |
| 67 | #' vol. 11, No 1, pp.1-30. doi:10.1142/S0219691313500033 |
| 68 | #' |
| 69 | #' @examples |
| 70 | #' \dontrun{ |
| 71 | #' # WER distances computations are too long for CRAN (for now) |
| 72 | #' |
| 73 | #' # Random series around cos(x,2x,3x)/sin(x,2x,3x) |
| 74 | #' x = seq(0,500,0.05) |
| 75 | #' L = length(x) #10001 |
| 76 | #' ref_series = matrix( c(cos(x),cos(2*x),cos(3*x),sin(x),sin(2*x),sin(3*x)), ncol=6 ) |
| 77 | #' library(wmtsa) |
| 78 | #' series = do.call( cbind, lapply( 1:6, function(i) |
| 79 | #' do.call(cbind, wmtsa::wavBootstrap(ref_series[i,], n.realization=400)) ) ) |
| 80 | #' #dim(series) #c(2400,10001) |
| 81 | #' medoids_ascii = claws(series, K1=60, K2=6, nb_per_chunk=c(200,500), verbose=TRUE) |
| 82 | #' |
| 83 | #' # Same example, from CSV file |
| 84 | #' csv_file = "/tmp/epclust_series.csv" |
| 85 | #' write.table(series, csv_file, sep=",", row.names=FALSE, col.names=FALSE) |
| 86 | #' medoids_csv = claws(csv_file, K1=60, K2=6, nb_per_chunk=c(200,500)) |
| 87 | #' |
| 88 | #' # Same example, from binary file |
| 89 | #' bin_file <- "/tmp/epclust_series.bin" |
| 90 | #' nbytes <- 8 |
| 91 | #' endian <- "little" |
| 92 | #' binarize(csv_file, bin_file, 500, nbytes, endian) |
| 93 | #' getSeries <- function(indices) getDataInFile(indices, bin_file, nbytes, endian) |
| 94 | #' medoids_bin <- claws(getSeries, K1=60, K2=6, nb_per_chunk=c(200,500)) |
| 95 | #' unlink(csv_file) |
| 96 | #' unlink(bin_file) |
| 97 | #' |
| 98 | #' # Same example, from SQLite database |
| 99 | #' library(DBI) |
| 100 | #' series_db <- dbConnect(RSQLite::SQLite(), "file::memory:") |
| 101 | #' # Prepare data.frame in DB-format |
| 102 | #' n <- nrow(series) |
| 103 | #' time_values <- data.frame( |
| 104 | #' id = rep(1:n,each=L), |
| 105 | #' time = rep( as.POSIXct(1800*(0:n),"GMT",origin="2001-01-01"), L ), |
| 106 | #' value = as.double(t(series)) ) |
| 107 | #' dbWriteTable(series_db, "times_values", times_values) |
| 108 | #' # Fill associative array, map index to identifier |
| 109 | #' indexToID_inDB <- as.character( |
| 110 | #' dbGetQuery(series_db, 'SELECT DISTINCT id FROM time_values')[,"id"] ) |
| 111 | #' serie_length <- as.integer( dbGetQuery(series_db, |
| 112 | #' paste("SELECT COUNT * FROM time_values WHERE id == ",indexToID_inDB[1],sep="")) ) |
| 113 | #' getSeries <- function(indices) { |
| 114 | #' request <- "SELECT id,value FROM times_values WHERE id in (" |
| 115 | #' for (i in indices) |
| 116 | #' request <- paste(request, indexToID_inDB[i], ",", sep="") |
| 117 | #' request <- paste(request, ")", sep="") |
| 118 | #' df_series <- dbGetQuery(series_db, request) |
| 119 | #' as.matrix(df_series[,"value"], nrow=serie_length) |
| 120 | #' } |
| 121 | #' medoids_db = claws(getSeries, K1=60, K2=6, nb_per_chunk=c(200,500)) |
| 122 | #' dbDisconnect(series_db) |
| 123 | #' |
| 124 | #' # All computed medoids should be the same: |
| 125 | #' digest::sha1(medoids_ascii) |
| 126 | #' digest::sha1(medoids_csv) |
| 127 | #' digest::sha1(medoids_bin) |
| 128 | #' digest::sha1(medoids_db) |
| 129 | #' } |
| 130 | #' @export |
| 131 | claws <- function(getSeries, K1, K2, |
| 132 | nb_per_chunk,nb_items_clust1=7*K1 #volumes of data |
| 133 | wav_filt="d8",contrib_type="absolute", #stage 1 |
| 134 | WER="end", #stage 2 |
| 135 | random=TRUE, #randomize series order? |
| 136 | ntasks=1, ncores_tasks=1, ncores_clust=4, #parallelism |
| 137 | sep=",", #ASCII input separator |
| 138 | nbytes=4, endian=.Platform$endian, #serialization (write,read) |
| 139 | verbose=FALSE, parll=TRUE) |
| 140 | { |
| 141 | # Check/transform arguments |
| 142 | if (!is.matrix(getSeries) && !bigmemory::is.big.matrix(getSeries) |
| 143 | && !is.function(getSeries) |
| 144 | && !methods::is(getSeries,"connection") && !is.character(getSeries)) |
| 145 | { |
| 146 | stop("'getSeries': [big]matrix, function, file or valid connection (no NA)") |
| 147 | } |
| 148 | K1 <- .toInteger(K1, function(x) x>=2) |
| 149 | K2 <- .toInteger(K2, function(x) x>=2) |
| 150 | if (!is.numeric(nb_per_chunk) || length(nb_per_chunk)!=2) |
| 151 | stop("'nb_per_chunk': numeric, size 2") |
| 152 | nb_per_chunk[1] <- .toInteger(nb_per_chunk[1], function(x) x>=1) |
| 153 | # A batch of contributions should have at least as many elements as a batch of series, |
| 154 | # because it always contains much less values |
| 155 | nb_per_chunk[2] <- max(.toInteger(nb_per_chunk[2],function(x) x>=1), nb_per_chunk[1]) |
| 156 | nb_items_clust1 <- .toInteger(nb_items_clust1, function(x) x>K1) |
| 157 | random <- .toLogical(random) |
| 158 | tryCatch |
| 159 | ( |
| 160 | {ignored <- wavelets::wt.filter(wav_filt)}, |
| 161 | error = function(e) stop("Invalid wavelet filter; see ?wavelets::wt.filter") |
| 162 | ) |
| 163 | ctypes = c("relative","absolute","logit") |
| 164 | contrib_type = ctypes[ pmatch(contrib_type,ctypes) ] |
| 165 | if (is.na(contrib_type)) |
| 166 | stop("'contrib_type' in {'relative','absolute','logit'}") |
| 167 | if (WER!="end" && WER!="mix") |
| 168 | stop("'WER': in {'end','mix'}") |
| 169 | random <- .toLogical(random) |
| 170 | ntasks <- .toInteger(ntasks, function(x) x>=1) |
| 171 | ncores_tasks <- .toInteger(ncores_tasks, function(x) x>=1) |
| 172 | ncores_clust <- .toInteger(ncores_clust, function(x) x>=1) |
| 173 | if (!is.character(sep)) |
| 174 | stop("'sep': character") |
| 175 | nbytes <- .toInteger(nbytes, function(x) x==4 || x==8) |
| 176 | verbose <- .toLogical(verbose) |
| 177 | parll <- .toLogical(parll) |
| 178 | |
| 179 | # Serialize series if required, to always use a function |
| 180 | bin_dir <- ".epclust_bin/" |
| 181 | dir.create(bin_dir, showWarnings=FALSE, mode="0755") |
| 182 | if (!is.function(getSeries)) |
| 183 | { |
| 184 | if (verbose) |
| 185 | cat("...Serialize time-series\n") |
| 186 | series_file = paste(bin_dir,"data",sep="") ; unlink(series_file) |
| 187 | binarize(getSeries, series_file, nb_series_per_chunk, sep, nbytes, endian) |
| 188 | getSeries = function(inds) getDataInFile(inds, series_file, nbytes, endian) |
| 189 | } |
| 190 | |
| 191 | # Serialize all computed wavelets contributions into a file |
| 192 | contribs_file = paste(bin_dir,"contribs",sep="") ; unlink(contribs_file) |
| 193 | index = 1 |
| 194 | nb_curves = 0 |
| 195 | if (verbose) |
| 196 | cat("...Compute contributions and serialize them\n") |
| 197 | nb_curves = binarizeTransform(getSeries, |
| 198 | function(series) curvesToContribs(series, wf, ctype), |
| 199 | contribs_file, nb_series_per_chunk, nbytes, endian) |
| 200 | getContribs = function(indices) getDataInFile(indices, contribs_file, nbytes, endian) |
| 201 | |
| 202 | if (nb_curves < K2) |
| 203 | stop("Not enough data: less series than final number of clusters") |
| 204 | nb_series_per_task = round(nb_curves / ntasks) |
| 205 | if (nb_series_per_task < K2) |
| 206 | stop("Too many tasks: less series in one task than final number of clusters") |
| 207 | |
| 208 | runTwoStepClustering = function(inds) |
| 209 | { |
| 210 | if (parll && ntasks>1) |
| 211 | require("epclust", quietly=TRUE) |
| 212 | indices_medoids = clusteringTask1( |
| 213 | inds, getContribs, K1, nb_series_per_chunk, ncores_clust, verbose, parll) |
| 214 | if (WER=="mix") |
| 215 | { |
| 216 | if (parll && ntasks>1) |
| 217 | require("bigmemory", quietly=TRUE) |
| 218 | medoids1 = bigmemory::as.big.matrix( getSeries(indices_medoids) ) |
| 219 | medoids2 = clusteringTask2(medoids1, K2, getSeries, nb_curves, nb_series_per_chunk, |
| 220 | nbytes, endian, ncores_clust, verbose, parll) |
| 221 | binarize(medoids2, synchrones_file, nb_series_per_chunk, sep, nbytes, endian) |
| 222 | return (vector("integer",0)) |
| 223 | } |
| 224 | indices_medoids |
| 225 | } |
| 226 | |
| 227 | # Cluster contributions in parallel (by nb_series_per_chunk) |
| 228 | indices_all = if (random) sample(nb_curves) else seq_len(nb_curves) |
| 229 | indices_tasks = lapply(seq_len(ntasks), function(i) { |
| 230 | upper_bound = ifelse( i<ntasks, min(nb_series_per_task*i,nb_curves), nb_curves ) |
| 231 | indices_all[((i-1)*nb_series_per_task+1):upper_bound] |
| 232 | }) |
| 233 | if (verbose) |
| 234 | { |
| 235 | message = paste("...Run ",ntasks," x stage 1", sep="") |
| 236 | if (WER=="mix") |
| 237 | message = paste(message," + stage 2", sep="") |
| 238 | cat(paste(message,"\n", sep="")) |
| 239 | } |
| 240 | if (WER=="mix") |
| 241 | {synchrones_file = paste(bin_dir,"synchrones",sep="") ; unlink(synchrones_file)} |
| 242 | if (parll && ntasks>1) |
| 243 | { |
| 244 | cl = parallel::makeCluster(ncores_tasks, outfile="") |
| 245 | varlist = c("getSeries","getContribs","K1","K2","verbose","parll", |
| 246 | "nb_series_per_chunk","ntasks","ncores_clust","sep","nbytes","endian") |
| 247 | if (WER=="mix") |
| 248 | varlist = c(varlist, "synchrones_file") |
| 249 | parallel::clusterExport(cl, varlist=varlist, envir = environment()) |
| 250 | } |
| 251 | |
| 252 | # 1000*K1 indices [if WER=="end"], or empty vector [if WER=="mix"] --> series on file |
| 253 | indices <- |
| 254 | if (parll && ntasks>1) |
| 255 | unlist( parallel::parLapply(cl, indices_tasks, runTwoStepClustering) ) |
| 256 | else |
| 257 | unlist( lapply(indices_tasks, runTwoStepClustering) ) |
| 258 | if (parll && ntasks>1) |
| 259 | parallel::stopCluster(cl) |
| 260 | |
| 261 | getRefSeries = getSeries |
| 262 | if (WER=="mix") |
| 263 | { |
| 264 | indices = seq_len(ntasks*K2) |
| 265 | #Now series must be retrieved from synchrones_file |
| 266 | getSeries = function(inds) getDataInFile(inds, synchrones_file, nbytes, endian) |
| 267 | #Contributions must be re-computed |
| 268 | unlink(contribs_file) |
| 269 | index = 1 |
| 270 | if (verbose) |
| 271 | cat("...Serialize contributions computed on synchrones\n") |
| 272 | ignored = binarizeTransform(getSeries, |
| 273 | function(series) curvesToContribs(series, wf, ctype), |
| 274 | contribs_file, nb_series_per_chunk, nbytes, endian) |
| 275 | } |
| 276 | |
| 277 | # Run step2 on resulting indices or series (from file) |
| 278 | if (verbose) |
| 279 | cat("...Run final // stage 1 + stage 2\n") |
| 280 | indices_medoids = clusteringTask1( |
| 281 | indices, getContribs, K1, nb_series_per_chunk, ncores_tasks*ncores_clust, verbose, parll) |
| 282 | medoids1 = bigmemory::as.big.matrix( getSeries(indices_medoids) ) |
| 283 | medoids2 = clusteringTask2(medoids1, K2, getRefSeries, nb_curves, nb_series_per_chunk, |
| 284 | nbytes, endian, ncores_tasks*ncores_clust, verbose, parll) |
| 285 | |
| 286 | # Cleanup |
| 287 | unlink(bin_dir, recursive=TRUE) |
| 288 | |
| 289 | medoids2[,] |
| 290 | } |
| 291 | |
| 292 | #' curvesToContribs |
| 293 | #' |
| 294 | #' Compute the discrete wavelet coefficients for each series, and aggregate them in |
| 295 | #' energy contribution across scales as described in https://arxiv.org/abs/1101.4744v2 |
| 296 | #' |
| 297 | #' @param series [big.]matrix of series (in columns), of size L x n |
| 298 | #' @inheritParams claws |
| 299 | #' |
| 300 | #' @return A [big.]matrix of size log(L) x n containing contributions in columns |
| 301 | #' |
| 302 | #' @export |
| 303 | curvesToContribs = function(series, wav_filt, contrib_type) |
| 304 | { |
| 305 | L = nrow(series) |
| 306 | D = ceiling( log2(L) ) |
| 307 | nb_sample_points = 2^D |
| 308 | apply(series, 2, function(x) { |
| 309 | interpolated_curve = spline(1:L, x, n=nb_sample_points)$y |
| 310 | W = wavelets::dwt(interpolated_curve, filter=wf, D)@W |
| 311 | nrj = rev( sapply( W, function(v) ( sqrt( sum(v^2) ) ) ) ) |
| 312 | if (contrib_type!="absolute") |
| 313 | nrj = nrj / sum(nrj) |
| 314 | if (contrib_type=="logit") |
| 315 | nrj = - log(1 - nrj) |
| 316 | nrj |
| 317 | }) |
| 318 | } |
| 319 | |
| 320 | # Check integer arguments with functional conditions |
| 321 | .toInteger <- function(x, condition) |
| 322 | { |
| 323 | errWarn <- function(ignored) |
| 324 | paste("Cannot convert argument' ",substitute(x),"' to integer", sep="") |
| 325 | if (!is.integer(x)) |
| 326 | tryCatch({x = as.integer(x)[1]; if (is.na(x)) stop()}, |
| 327 | warning = errWarn, error = errWarn) |
| 328 | if (!condition(x)) |
| 329 | { |
| 330 | stop(paste("Argument '",substitute(x), |
| 331 | "' does not verify condition ",body(condition), sep="")) |
| 332 | } |
| 333 | x |
| 334 | } |
| 335 | |
| 336 | # Check logical arguments |
| 337 | .toLogical <- function(x) |
| 338 | { |
| 339 | errWarn <- function(ignored) |
| 340 | paste("Cannot convert argument' ",substitute(x),"' to logical", sep="") |
| 341 | if (!is.logical(x)) |
| 342 | tryCatch({x = as.logical(x)[1]; if (is.na(x)) stop()}, |
| 343 | warning = errWarn, error = errWarn) |
| 344 | x |
| 345 | } |