1 ## File : 06_predictions.r
2 ## Description : Calculer les synchrones pour chaque groupe
3 ## obtenu par le clustering.
7 if(Sys.info()[4] == "mojarrita"){
8 setwd("~/Documents/projects/2014_EDF-Orsay-Lyon2/codes/")
10 setwd("~/2014_EDF-Orsay-Lyon2/codes/")
13 ## 1. Read auxiliar data files ####
15 identifiants <- read.table("identifs.txt")[ ,1]
16 dates0 <- read.table("datesall.txt")[, 1]
17 dates09 <- as.character(dates0[grep("2009", dates0)])
18 dates10 <- as.character(dates0[grep("2010", dates0)])
19 dates <- c(dates09, dates10)
20 rm(dates0, dates09, dates10)
22 n <- length(identifiants)
25 # groups for the predictor
26 gr_calen <- read.table("calendar_transition_groups-1996-2011.txt")
27 i_start_calen <- which(rownames(gr_calen) == "2009-01-01")
28 i_end_calen <- which(rownames(gr_calen) == "2010-12-31")
29 gr <- gr_calen$gr[i_start_calen:i_end_calen]
30 gr <- gsub("Thu", "TWT", gr)
31 gr <- gsub("Wed", "TWT", gr)
32 gr <- gsub("Tue", "TWT", gr)
34 days <- seq.Date(from = as.Date("2009/01/01"),
35 to = as.Date("2010/12/31"),
38 #synchros09 <- as.matrix(read.table("~/tmp/2009_synchros200.txt"))
39 #synchros10 <- as.matrix(read.table("~/tmp/2010_synchros200.txt"))
40 synchros09 <- as.matrix(read.table("~/tmp/2009_synchros200RC.txt"))
41 synchros10 <- as.matrix(read.table("~/tmp/2010_synchros200RC.txt"))
43 synchros09 <- t(synchros09) # series must be on rows
44 synchros10 <- t(synchros10) # series must be on rows
46 ## 2. Imputation of missing values ####
47 # (done by linear interpolation)
48 nas <- which(is.na(synchros09)[, 1]) # some 1/1/2009 are missing
49 synchros09[nas, 1] <- rowMeans(synchros09[nas, 2:4])
51 imput09 <- synchros09[, 4180:4181] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2)
52 imput10 <- synchros10[, 4132:4133] %*% matrix(c(2/3, 1/3, 1/3, 2/3), 2)
54 synchros09 <- cbind(synchros09[, 1:4180], imput09, synchros09[, 4181:17518])
55 synchros10 <- cbind(synchros10[, 1:4132], imput10, synchros10[, 4133:17518])
58 rm(imput09, imput10, nas)
59 rm(gr_calen, i_start_calen, i_end_calen)
61 # Set start and end timepoints for the prediction
62 i_00 <- which(days == "2009-01-01")
63 n_00 <- which(days == "2009-12-31")
64 n_01 <- which(days == "2010-12-31")
67 ## 3. Construct predictors ####
69 synchros <- (cbind( synchros09, synchros10 ))[-201, ]
71 #clustfactors <- dir("~/tmp/", pattern = 'clustfactorAC')
72 #clustfactors <- dir("~/tmp/", pattern = 'clustfactorRC')
73 clustfactors <- dir("~/tmp/", pattern = 'clustfactorWER')
75 maxK <- length(clustfactors) - 1
77 performK <- data.frame(K = rep(NA_real_, maxK),
78 mape = rep(NA_real_, maxK),
79 rmse = rep(NA_real_, maxK))
82 for(cf in seq_along(clustfactors)) {
83 #clustfactor <- read.table(file = "~/tmp/clustfactor3.txt")
84 clustfactor <- read.table(file = clustfactors[cf])
86 K <- nrow(unique(clustfactor))
90 SC <- matrix(0, ncol = ncol(synchros), nrow = K)
92 clustk <- which(clustfactor == k)
93 if(length(clustk) > 0) {
94 if(length(clustk) > 1) {
95 SCk <- colSums(synchros[which(clustfactor == k), ])
97 SCk <- synchros[which(clustfactor == k), ]
99 SC[k, ] <- SC[k, ] + SCk
104 mat_Load <- matrix(colSums(SC), ncol = 48, byrow = TRUE)
105 mat_ks <- lapply(1:K, function(k) matrix(SC[k, ], ncol = 48, byrow = TRUE))
108 ## 2. Transform to wkdata (performs DWT) ##
110 wk_Load <- wkdata(X = c(t(mat_Load)),
111 gr = gr[1:nrow(mat_Load)],
114 wks <- lapply(mat_ks, function(mat_k)
115 wkdata(X = c(t(mat_k)), gr = gr[1:nrow(mat_Load)], p = 48))
117 ## 3. Begin of big loop over ##
118 grid <- n_00:(n_01 - 1) # these are the reference days
119 #h_Load <- h_k1 <- h_k2 <- h_k3 <- double(length(grid))
121 # output matrix of predicted days
122 # names of the predicted days (reference days + 1)
123 pred_Load <- matrix(NA_real_, ncol = 48, nrow = length(grid))
124 preds <- lapply(1:K, function(k) pred_Load)
127 for(q in seq_along(grid)[-c(121, 122, 227, 228)]) { #seq_along(grid)[-c(46:47)]
128 if((q %/% 10) == (q / 10)) cat(paste(" iteration: ", q, "\n"))
131 select_Load <- select.wkdata(wk_Load, i_00:Q)
132 aux_Load <- wavkerfun(obj = select_Load, kerneltype = "gaussian",
133 EPS = .Machine$double.eps)
134 pred_Load[q, ] <- predict.wavkerfun(obj = aux_Load, cent = "DIFF")$X
136 selects <- lapply(wks, function(k) select.wkdata(k, i_00:Q))
137 auxs <- lapply(selects,
138 function(k) wavkerfun(obj = k,
139 kerneltype = "gaussian",
140 EPS = .Machine$double.eps))
142 preds[[k]][q, ] <- predict.wavkerfun(obj = auxs[[k]], cent = "DIFF")$X
147 #rm(wks, selects, auxs, aux_Load, mat_ks, q, Q, wk_Load)
149 # Compute the prediction of synchrone by the sum of individuals
150 pred_Load_sum <- Reduce('+', preds)
152 ## 4. Analisis of the predictions ####
153 resids_Load <- pred_Load - mat_Load[grid + 1, ] # i.e. reference days + 1
154 resids_Load_sum <- pred_Load_sum - mat_Load[grid + 1, ] # i.e. reference days + 1
156 mape_Load <- 100 * rowMeans(abs(resids_Load / mat_Load[grid + 1, ]))
157 rmse_Load <- sqrt(rowMeans(resids_Load^2))
159 mape_Load_sum <- 100 * rowMeans(abs(resids_Load_sum / mat_Load[grid + 1, ]))
160 rmse_Load_sum <- sqrt(rowMeans(resids_Load_sum^2))
162 performK[cf, 2] <- mean(mape_Load, na.rm = TRUE)
163 performK[cf, 3] <- mean(mape_Load_sum, na.rm = TRUE)
164 #print(mean(mape_Load, na.rm = TRUE))
165 #print(mean(mape_Load_sum, na.rm = TRUE))
170 performKwer <- performK
172 #performKrc <- read.table("~/Documents/projects/2014_EDF-Orsay-Lyon2/res/perform200RC.txt")
173 #write.table(file = '~/Documents/projects/2014_EDF-Orsay-Lyon2/res/perform200RC.txt',
175 write.table(file = '~/Documents/projects/2014_EDF-Orsay-Lyon2/res/perform200WER_3.txt',
178 #performK <- data.frame(K = performKrc$K,
179 # rc = performKrc$rmse,
180 # wer = performKwer$rmse)
183 performK <- performK[order(performK$K), ]
186 #plot(performK$K, performK$rc, type = 'l', ylim = c(1.7,2.2),
188 #lines(performK$K, performK$wer, col = 3, lwd = 2)
189 #abline(h = mean(performKrc$mape), lwd = 2)
190 #legend("topright", c('Sync 32K', 'Cl. RC', "Cl. WER"), lwd = 2, col = 1:3)