Commit | Line | Data |
---|---|---|
ad642dc6 BA |
1 | ## File : 00_plots-energycon.r |
2 | ## Description : Using the full data matrix, we extract handy features to | |
3 | ## cluster. | |
4 | ||
5 | rm(list = ls()) | |
6 | ||
7 | library(Rwave) # CWT | |
8 | #library(cluster) # pam | |
9 | #library(flexclust) # kcca | |
10 | source("aux.r") # auxiliary clustering functions | |
11 | source("sowas-superseded.r") # auxiliary CWT functions | |
12 | source('http://eric.univ-lyon2.fr/~jcugliari/codes/functional-clustering.r') | |
13 | setwd("~/recherche/03_projects/2014_EDF-Orsay-Lyon2/codes/") | |
14 | ||
15 | ||
16 | ## 1. Read auxiliar data files #### | |
17 | ||
18 | #identifiants <- read.table("identifs.txt")[ ,1] | |
19 | dates0 <- read.table("datesall.txt")[, 1] | |
20 | dates <- dates0[grep("2009", dates0)] | |
21 | #rm(dates0) | |
22 | ||
23 | #n <- length(identifiants) | |
24 | p <- length(dates) | |
25 | ||
26 | #blocks <- c(rep(6500, 3), 5511) | |
27 | ||
28 | # table( substr(dates, 11, 15) ) # Sunlight time saving produces an | |
29 | # unbalanced number of time points | |
30 | # per time stepa across the year | |
31 | ||
32 | ||
33 | ## 2. Process the large file #### | |
34 | ||
35 | # if(exists("con")) close(con) | |
36 | # con <- file("~/tmp/2009_full.txt") # Establish a connection to the file | |
37 | # open(con, "r") # Open the connection | |
38 | # | |
39 | # nb <- 4 | |
40 | # actual <- readLines(con = con, n = nb )[-3] | |
41 | # auxmat <- matrix(unlist(strsplit(actual, " ")), ncol = p + 1, byrow = TRUE) | |
42 | # rm(actual) | |
43 | # | |
44 | # datamat <- t(apply(auxmat[, -1], 1, as.numeric)) | |
45 | # rownames(datamat) <- substr(auxmat[, 1], 2, 7) | |
46 | # rm(auxmat) | |
47 | # | |
48 | # datamat <- datamat[, 1:(48 * 14)] | |
49 | # p <- ncol(datamat) | |
50 | # | |
51 | # auxDWT <- t(apply(datamat, 1, toDWT)) | |
52 | # matcontrib <- t(apply(auxDWT, 1, contrib)) | |
53 | # rm(auxDWT) | |
54 | # | |
55 | # close(con); rm(con) # close connection to the file | |
56 | ||
57 | ||
58 | load("~/data/Irlande/Data_CER_clean/SME.RData") | |
59 | ||
60 | smp <- c(4, 79, 126) | |
61 | ||
62 | datamat <- t(SME[1:(48 * 7), smp]) | |
63 | p <- ncol(datamat) | |
64 | ||
65 | auxDWT <- t(apply(datamat, 1, toDWT)) | |
66 | matcontrib <- t(apply(auxDWT, 1, contrib)) | |
67 | rm(auxDWT) | |
68 | ||
69 | matplot(t(datamat), type = "l", lty = 1, | |
70 | col = 1:8, lwd = 2) | |
71 | ||
72 | pdf("~/courbes.pdf") | |
73 | #par(mai = c(1, 1, 0.8, 0.6), mfcol = c(4, 2), cex = 2) | |
74 | for(courbe in 1:ncol(SME)) { | |
75 | plot(SME[1:(48 * 7), courbe], main = paste(courbe), | |
76 | xlab = "", ylab = "", type = "l") | |
77 | } | |
78 | dev.off() | |
79 | ||
80 | ||
81 | ||
82 | ||
83 | ## Plots for ENERGYCON full article #### | |
84 | op <- par() | |
85 | ||
86 | col <- c("grey", "black", "black") | |
87 | ## Curves | |
88 | pdf('~/curves.pdf', width = 12) | |
89 | par(mai = c(1.3, 1.3, 0.1, 0.1), cex = 1.2) | |
90 | matplot(t(datamat), type = "l", lty = c(1, 1, 2), | |
91 | lwd = 2, col = col, # c(1, 2, 4), | |
92 | ylab = "Load", xlab = "Time (1/2 hours)") | |
93 | legend("top", c("Cust. A", "Cust. B", "Cust. C"), | |
94 | col = col, lty = c(1, 1, 2), ncol = 3, lwd = 2) | |
95 | # matplot(scale(t(datamat), scale = F), | |
96 | # type = "l", lty = 1, col = c(1, 2, 4), lwd = 2, | |
97 | # ylab = "Load", xlab = "Time (1/2 hours)") | |
98 | dev.off() | |
99 | ||
100 | normi <- function(x) x / max(x) | |
101 | normi2 <- function(x) (x - min(x)) / (max(x) - min(x)) | |
102 | ||
103 | matplot(apply(datamat, 1, normi2), | |
104 | type = "l", lty = 1, col = c(1, 2, 4), lwd = 2, | |
105 | ylab = "Load", xlab = "Time (1/2 hours)") | |
106 | ||
107 | ||
108 | cont_dist <- dist(scale(matcontrib)) | |
109 | cont_dist <- cont_dist / max(cont_dist) | |
110 | ||
111 | ||
112 | delta <- p | |
113 | n <- nrow(datamat) | |
114 | ||
115 | ## _.a CWT -- Filtering the lowest freqs (>6m) #### | |
116 | nvoice <- 4 | |
117 | # noctave4 = 2^13 = 8192 half hours ~ 180 days | |
118 | noctave4 <- adjust.noctave(N = delta, dt = 1, s0 = 2, | |
119 | tw = 0, noctave = 10) | |
120 | # 10 here represent 2^10 = 1024 half-hours ~ 2 weeks day | |
121 | scalevector4 <- 2^(4:(noctave4 * nvoice) / nvoice) * 2 | |
122 | lscvect4 <- length(scalevector4) | |
123 | lscvect <- lscvect4 # i should clean my code: werFam demands a lscvect | |
124 | Xcwt4 <- toCWT(datamat, noctave = noctave4, dt = 1, | |
125 | scalevector = scalevector4, | |
126 | lt = delta, smooth = FALSE, | |
127 | nvoice = nvoice) # observations node with CWT | |
128 | ||
129 | Xcwt2 <- matrix(NA_complex_, nrow= n, ncol= 2 + length(c(Xcwt4[, ,1]))) | |
130 | ||
131 | ||
132 | for(i in 1:n) | |
133 | Xcwt2[i,] <- c(delta, lscvect, Xcwt4[,,i] / max(Mod(Xcwt4[,,i])) ) | |
134 | ||
135 | rm(conso, Xcwt4); gc() | |
136 | ||
137 | lscvect <- 41 ## very very nasty: toCWT changes scalevector to 41 (instead of | |
138 | ## the original length--37--) | |
139 | ## _.b WER^2 distances ######## | |
140 | Xwer_dist <- matrix(0.0, n, n) | |
141 | for(i in 1:(n - 1)){ | |
142 | mat1 <- vect2mat(Xcwt2[i,]) | |
143 | for(j in (i + 1):n){ | |
144 | mat2 <- vect2mat(Xcwt2[j,]) | |
145 | num <- Mod(mat1 * Conj(mat2)) | |
146 | WX <- Mod(mat1 * Conj(mat1)) | |
147 | WY <- Mod(mat2 * Conj(mat2)) | |
148 | smsmnum <- smCWT(num, scalevector = scalevector4) | |
149 | smsmWX <- smCWT(WX, scalevector = scalevector4) | |
150 | smsmWY <- smCWT(WY, scalevector = scalevector4) | |
151 | wer2 <- sum(colSums(smsmnum)^2) / | |
152 | sum( sum(colSums(smsmWX) * colSums(smsmWY)) ) | |
153 | Xwer_dist[i, j] <- sqrt(delta * lscvect * (1 - wer2)) | |
154 | Xwer_dist[j, i] <- Xwer_dist[i, j] | |
155 | } | |
156 | } | |
157 | diag(Xwer_dist) <- numeric(n) | |
158 | ||
159 | ||
160 | Xwer_dist <- Xwer_dist / max(Xwer_dist) | |
161 | ||
162 | ||
163 | pdf("~/cmdscale.pdf", width = 14) | |
164 | #layout(matrix(1:2, 1, 2)) | |
165 | par(mai = c(1, 1, 0.8, 0.6), mfcol = c(1, 2), cex = 2) | |
166 | plot(cmdscale(cont_dist), pch = c("A", "B", "C"),#c(15, 17, 19), | |
167 | main = "RC based distance", | |
168 | #col = c(1, 2, 4), | |
169 | xlim = c(-1, 1), ylim = c(-1, 1), | |
170 | xlab = "", ylab = "", asp = 1) | |
171 | abline(h= 0, lwd = 2); abline(v = 0, lwd = 2); grid(lty = 1) | |
172 | ||
173 | plot(cmdscale(Xwer_dist), pch = c("A", "B", "C"),#c(15, 17, 19), | |
174 | main = "WER distance", | |
175 | #col = c(1, 2, 4), | |
176 | xlim = c(-1, 1), ylim = c(-1, 1), | |
177 | xlab = "", ylab = "", asp = 1) | |
178 | abline(h= 0, lwd = 2); abline(v = 0, lwd = 2); grid(lty = 1) | |
179 | dev.off() | |
180 | ||
181 | ||
182 | ||
183 | ## Effectives distribution | |
184 | #load('../res/clfit500.Rdata') | |
185 | pdf("~/distro500.pdf", width = 12) | |
186 | plot(sort(table(clfit$clustering), decreasing = TRUE), | |
187 | type ="s", ylab = "Effectives", xlab = "Class", lwd = 2) | |
188 | abline(v = 200) | |
189 | #dev.off() | |
190 |