ac4ea8ddc40567b72d84c240743fbc38d4e57971
[epclust.git] / epclust / R / main.R
1 #' @title Cluster power curves with PAM in parallel
2 #'
3 #' @description Groups electricity power curves (or any series of similar nature) by applying PAM
4 #' algorithm in parallel to chunks of size \code{nb_series_per_chunk}
5 #'
6 #' @param data Access to the data, which can be of one of the three following types:
7 #' \itemize{
8 #' \item data.frame: each line contains its ID in the first cell, and all values after
9 #' \item connection: any R connection object (e.g. a file) providing lines as described above
10 #' \item function: a custom way to retrieve the curves; it has two arguments: the ranks to be
11 #' retrieved, and the IDs - at least one of them must be present (priority: ranks).
12 #' }
13 #' @param K1 Number of super-consumers to be found after stage 1 (K1 << N)
14 #' @param K2 Number of clusters to be found after stage 2 (K2 << K1)
15 #' @param ntasks Number of tasks (parallel iterations to obtain K1 medoids); default: 1.
16 #' Note: ntasks << N, so that N is "roughly divisible" by N (number of series)
17 #' @param nb_series_per_chunk (Maximum) number of series in each group, inside a task
18 #' @param min_series_per_chunk Minimum number of series in each group
19 #' @param wf Wavelet transform filter; see ?wt.filter. Default: haar
20 #' @param WER "end" to apply stage 2 after stage 1 has iterated and finished, or "mix"
21 #' to apply it after every stage 1
22 #' @param ncores_tasks "MPI" number of parallel tasks (1 to disable: sequential tasks)
23 #' @param ncores_clust "OpenMP" number of parallel clusterings in one task
24 #' @param random Randomize chunks repartition
25 #' @param ... Other arguments to be passed to \code{data} function
26 #'
27 #' @return A data.frame of the final medoids curves (identifiers + values)
28 #'
29 #' @examples
30 #' getData = function(start, n) {
31 #' con = dbConnect(drv = RSQLite::SQLite(), dbname = "mydata.sqlite")
32 #' df = dbGetQuery(con, paste(
33 #' "SELECT * FROM times_values GROUP BY id OFFSET ",start,
34 #' "LIMIT ", n, " ORDER BY date", sep=""))
35 #' return (df)
36 #' }
37 #' #TODO: 3 examples, data.frame / binary file / DB sqLite
38 #' + sampleCurves : wavBootstrap de package wmtsa
39 #' cl = epclust(getData, K1=200, K2=15, ntasks=1000, nb_series_per_chunk=5000, WER="mix")
40 #' @export
41 epclust = function(series,K1,K2,ntasks=1,nb_series_per_chunk=50*K1,min_series_per_chunk=5*K1,
42 wf="haar",WER="end",ncores_tasks=1,ncores_clust=4,random=TRUE,...)
43 {
44 # Check/transform arguments
45 bin_dir = "epclust.bin/"
46 dir.create(bin_dir, showWarnings=FALSE, mode="0755")
47 if (!is.function(series))
48 {
49 series_file = paste(bin_dir,"data",sep="")
50 unlink(series_file)
51 }
52 if (is.matrix(series))
53 serialize(series, series_file)
54 else if (!is.function(series))
55 {
56 tryCatch(
57 {
58 if (is.character(series))
59 series_con = file(series, open="r")
60 else if (!isOpen(series))
61 {
62 open(series)
63 series_con = series
64 }
65 serialize(series_con, series_file)
66 close(series_con)
67 },
68 error=function(e) "series should be a data.frame, a function or a valid connection"
69 )
70 }
71 if (!is.function(series))
72 series = function(indices) getDataInFile(indices, series_file)
73 getSeries = series
74
75 K1 = toInteger(K1, function(x) x>=2)
76 K2 = toInteger(K2, function(x) x>=2)
77 ntasks = toInteger(ntasks, function(x) x>=1)
78 nb_series_per_chunk = toInteger(nb_series_per_chunk, function(x) x>=K1)
79 min_series_per_chunk = toInteger(K1, function(x) x>=K1 && x<=nb_series_per_chunk)
80 ncores_tasks = toInteger(ncores_tasks, function(x) x>=1)
81 ncores_clust = toInteger(ncores_clust, function(x) x>=1)
82 if (WER!="end" && WER!="mix")
83 stop("WER takes values in {'end','mix'}")
84
85 # Serialize all wavelets coefficients (+ IDs) onto a file
86 coefs_file = paste(bin_dir,"coefs",sep="")
87 unlink(coefs_file)
88 index = 1
89 nb_curves = 0
90 repeat
91 {
92 series = getSeries((index-1)+seq_len(nb_series_per_chunk))
93 if (is.null(series))
94 break
95 coeffs_chunk = curvesToCoeffs(series, wf)
96 serialize(coeffs_chunk, coefs_file)
97 index = index + nb_series_per_chunk
98 nb_curves = nb_curves + nrow(coeffs_chunk)
99 }
100 getCoefs = function(indices) getDataInFile(indices, coefs_file)
101 ######TODO: if DB, array rank --> ID at first retrieval, when computing coeffs; so:: NO use of IDs !
102
103 if (nb_curves < min_series_per_chunk)
104 stop("Not enough data: less rows than min_series_per_chunk!")
105 nb_series_per_task = round(nb_curves / ntasks)
106 if (nb_series_per_task < min_series_per_chunk)
107 stop("Too many tasks: less series in one task than min_series_per_chunk!")
108
109 # Cluster coefficients in parallel (by nb_series_per_chunk)
110 indices = if (random) sample(nb_curves) else seq_len(nb_curves)
111 indices_tasks = lapply(seq_len(ntasks), function(i) {
112 upper_bound = ifelse( i<ntasks, min(nb_series_per_task*i,nb_curves), nb_curves )
113 indices[((i-1)*nb_series_per_task+1):upper_bound]
114 })
115 cl_tasks = parallel::makeCluster(ncores_tasks)
116 parallel::clusterExport(cl_tasks,
117 varlist=c("getSeries","getCoefs","K1","K2","WER","nb_series_per_chunk","ncores_clust"),
118 envir=environment())
119 #1000*K1 (or K2) indices (or NOTHING--> series on file)
120 indices = parallel::parLapply(cl_tasks, indices_tasks, clusteringTask)
121 parallel::stopCluster(cl_tasks)
122
123 #Now series must be retrieved from synchrones_file, and have no ID
124 getSeries = function(indices, ids) getDataInFile(indices, synchrones_file)
125
126 # Run step2 on resulting indices or series (from file)
127 computeClusters2(indices=if (WER=="end") indices else NULL, K2, to_file=FALSE)
128 }