Commit | Line | Data |
---|---|---|
98e958ca | 1 | #' Plot curves |
e030a6e3 | 2 | #' |
af3b84f4 | 3 | #' Plot a range of curves in data |
e030a6e3 BA |
4 | #' |
5 | #' @param data Object of class Data | |
6 | #' @param indices Range of indices (integers or dates) | |
7 | #' | |
8 | #' @export | |
1e20780e | 9 | plotCurves <- function(data, indices=seq_len(data$getSize())) |
e030a6e3 | 10 | { |
98e958ca | 11 | series = data$getSeries(indices) |
a5a3a294 | 12 | yrange = quantile(series, probs=c(0.025,0.975), na.rm=TRUE) |
e030a6e3 BA |
13 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
14 | for (i in seq_along(indices)) | |
15 | { | |
98e958ca | 16 | plot(series[,i], type="l", ylim=yrange, |
e030a6e3 | 17 | xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10","")) |
1e20780e | 18 | if (i < length(indices)) |
e030a6e3 BA |
19 | par(new=TRUE) |
20 | } | |
21 | } | |
22 | ||
af3b84f4 | 23 | #' Plot error |
3d69ff21 | 24 | #' |
af3b84f4 | 25 | #' Draw error graphs, potentially from several runs of \code{computeForecast} |
3d69ff21 | 26 | #' |
99f83c9a | 27 | #' @param err Error as returned by \code{computeError} |
09cf9c19 | 28 | #' @param cols Colors for each error (default: 1,2,3,...) |
3d69ff21 | 29 | #' |
98e958ca BA |
30 | #' @seealso \code{\link{plotCurves}}, \code{\link{plotPredReal}}, |
31 | #' \code{\link{plotSimils}}, \code{\link{plotFbox}}, | |
32 | #' \code{\link{computeFilaments}, }\code{\link{plotFilamentsBox}}, \code{\link{plotRelVar}} | |
3d69ff21 BA |
33 | #' |
34 | #' @export | |
09cf9c19 | 35 | plotError <- function(err, cols=seq_along(err)) |
3d69ff21 | 36 | { |
09cf9c19 BA |
37 | if (!is.null(err$abs)) |
38 | err = list(err) | |
4e95ec8f | 39 | par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) |
3d69ff21 | 40 | L = length(err) |
98e958ca | 41 | yrange = range( sapply(1:L, function(i) ( err[[i]]$abs$day ) ), na.rm=TRUE ) |
3d69ff21 BA |
42 | for (i in seq_len(L)) |
43 | { | |
44 | plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), | |
09cf9c19 | 45 | ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
46 | if (i < L) |
47 | par(new=TRUE) | |
48 | } | |
98e958ca | 49 | yrange = range( sapply(1:L, function(i) ( err[[i]]$abs$indices ) ), na.rm=TRUE ) |
3d69ff21 BA |
50 | for (i in seq_len(L)) |
51 | { | |
52 | plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), | |
09cf9c19 | 53 | ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
54 | if (i < L) |
55 | par(new=TRUE) | |
56 | } | |
98e958ca | 57 | yrange = range( sapply(1:L, function(i) ( err[[i]]$MAPE$day ) ), na.rm=TRUE ) |
3d69ff21 BA |
58 | for (i in seq_len(L)) |
59 | { | |
60 | plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), | |
09cf9c19 | 61 | ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
62 | if (i < L) |
63 | par(new=TRUE) | |
64 | } | |
98e958ca | 65 | yrange = range( sapply(1:L, function(i) ( err[[i]]$MAPE$indices ) ), na.rm=TRUE ) |
3d69ff21 BA |
66 | for (i in seq_len(L)) |
67 | { | |
68 | plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), | |
09cf9c19 | 69 | ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
70 | if (i < L) |
71 | par(new=TRUE) | |
72 | } | |
73 | } | |
74 | ||
98e958ca BA |
75 | #' Plot measured / predicted |
76 | #' | |
77 | #' Plot measured curve (in black) and predicted curve (in blue) | |
78 | #' | |
79 | #' @param data Object return by \code{getData} | |
80 | #' @param pred Object as returned by \code{computeForecast} | |
81 | #' @param index Index in forecasts (integer or date) | |
82 | #' | |
83 | #' @export | |
84 | plotPredReal <- function(data, pred, index) | |
85 | { | |
86 | horizon = length(pred$getSerie(1)) | |
87 | measure = data$getSerie( pred$getIndexInData(index)+1 )[1:horizon] | |
88 | prediction = pred$getSerie(index) | |
89 | yrange = range(measure, prediction) | |
90 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=3) | |
91 | plot(measure, type="l", ylim=yrange, xlab="Temps (en heures)", ylab="PM10") | |
92 | par(new=TRUE) | |
93 | plot(prediction, type="l", col="#0000FF", ylim=yrange, xlab="", ylab="") | |
94 | } | |
95 | ||
96 | #' Plot similarities | |
97 | #' | |
98 | #' Plot histogram of similarities (weights) | |
99 | #' | |
100 | #' @param pred Object as returned by \code{computeForecast} | |
101 | #' @param index Index in forecasts (integer or date) | |
102 | #' | |
103 | #' @export | |
104 | plotSimils <- function(pred, index) | |
105 | { | |
106 | weights = pred$getParams(index)$weights | |
107 | if (is.null(weights)) | |
108 | stop("plotSimils only works on 'Neighbors' forecasts") | |
109 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
110 | hist(pred$getParams(index)$weights, nclass=20, xlab="Poids", ylab="Effectif") | |
111 | } | |
112 | ||
af3b84f4 | 113 | #' Functional boxplot |
3d69ff21 | 114 | #' |
af3b84f4 | 115 | #' Draw the functional boxplot on the left, and bivariate plot on the right |
3d69ff21 BA |
116 | #' |
117 | #' @param data Object return by \code{getData} | |
98e958ca | 118 | #' @param indices integer or date indices to process |
fa8078f9 | 119 | #' @param plot_bivariate Should the bivariate plot appear? |
3d69ff21 BA |
120 | #' |
121 | #' @export | |
98e958ca | 122 | plotFbox <- function(data, indices=seq_len(data$getSize())) |
3d69ff21 BA |
123 | { |
124 | if (!requireNamespace("rainbow", quietly=TRUE)) | |
125 | stop("Functional boxplot requires the rainbow package") | |
126 | ||
98e958ca BA |
127 | series_matrix = data$getSeries(indices) |
128 | # Remove series with NAs | |
af3b84f4 BA |
129 | no_NAs_indices = sapply( 1:ncol(series_matrix), |
130 | function(i) all(!is.na(series_matrix[,i])) ) | |
99f83c9a | 131 | series_matrix = series_matrix[,no_NAs_indices] |
3d69ff21 BA |
132 | |
133 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) | |
fa8078f9 | 134 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
3d69ff21 BA |
135 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Temps (heures)", ylab="PM10", |
136 | plotlegend=FALSE, lwd=2) | |
98e958ca BA |
137 | rainbow::fboxplot(series_fds, "bivariate", "hdr", plotlegend=FALSE) |
138 | } | |
139 | ||
140 | #' Compute filaments | |
141 | #' | |
142 | #' Get similar days in the past, as black as distances are small | |
143 | #' | |
144 | #' @param data Object as returned by \code{getData} | |
145 | #' @param index Index in data (integer or date) | |
146 | #' @param limit Number of neighbors to consider | |
147 | #' @param plot Should the result be plotted? | |
148 | #' | |
149 | #' @return A list with | |
150 | #' \itemize{ | |
151 | #' \item index : index of the current serie ('today') | |
152 | #' \item neighb_indices : indices of its neighbors | |
153 | #' \item colors : colors of neighbors curves (shades of gray) | |
154 | #' } | |
155 | #' | |
156 | #' @export | |
157 | computeFilaments <- function(data, index, limit=60, plot=TRUE) | |
158 | { | |
159 | ref_serie = data$getCenteredSerie(index) | |
160 | if (any(is.na(ref_serie))) | |
161 | stop("computeFilaments requires a serie without NAs") | |
162 | ||
163 | # Determine indices of no-NAs days followed by no-NAs tomorrows | |
164 | fdays = getNoNA2(data, 1, dateIndexToInteger(index,data)-1) | |
165 | # Series + tomorrows in columns, ref_serie first | |
166 | centered_series = data$getCenteredSeries(fdays) | |
167 | ||
168 | # Obtain neighbors (closest for euclidian norm) | |
169 | L = length(ref_serie) | |
170 | distances = sqrt( colSums( (centered_series - ref_serie)^2 / L ) ) | |
171 | sorted_distances = sort(distances, index.return=TRUE) | |
172 | ||
173 | # Compute colors for each neighbor (from darkest to lightest) | |
174 | nn = min(limit, length(distances)) | |
175 | min_dist = min(sorted_distances$x[1:nn]) | |
176 | max_dist = max(sorted_distances$x[1:nn]) | |
177 | color_values = floor( 19.5 * (sorted_distances$x[1:nn]-min_dist) / (max_dist-min_dist) ) + 1 | |
178 | colors = gray.colors(20,0.1,0.9)[color_values] #TODO: 20 == magic number | |
179 | ||
180 | if (plot) | |
181 | { | |
182 | # Complete series with (past and present) tomorrows | |
183 | ref_serie = c(ref_serie,data$getCenteredSerie(index+1)) | |
184 | centered_series = rbind( centered_series, data$getCenteredSeries(fdays+1) ) | |
a5a3a294 | 185 | yrange = quantile(cbind(ref_serie,centered_series), probs=c(0.025,0.975), na.rm=TRUE) |
98e958ca BA |
186 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) |
187 | for (i in nn:1) | |
188 | { | |
189 | plot(centered_series[,sorted_distances$ix[i]], ylim=yrange, type="l", col=colors[i], | |
190 | xlab=ifelse(i==nn,"Temps (en heures)",""), ylab=ifelse(i==nn,"PM10 centré","")) | |
191 | par(new=TRUE) | |
192 | } | |
193 | # Also plot ref curve, in red | |
194 | plot(ref_serie, ylim=yrange, type="l", col="#FF0000", xlab="", ylab="") | |
6d50a76f | 195 | abline(v=24, lty=2, col=colors()[56], lwd=1) |
98e958ca BA |
196 | } |
197 | ||
198 | list("index"=index,"neighb_indices"=fdays[sorted_distances$ix[1:nn]],"colors"=colors) | |
fa8078f9 BA |
199 | } |
200 | ||
af3b84f4 | 201 | #' Functional boxplot on filaments |
fa8078f9 | 202 | #' |
af3b84f4 | 203 | #' Draw the functional boxplot on filaments obtained by \code{computeFilaments} |
fa8078f9 BA |
204 | #' |
205 | #' @param data Object return by \code{getData} | |
98e958ca | 206 | #' @param fil Output of \code{computeFilaments} |
fa8078f9 BA |
207 | #' |
208 | #' @export | |
98e958ca | 209 | plotFilamentsBox = function(data, fil, ...) |
fa8078f9 | 210 | { |
98e958ca BA |
211 | if (!requireNamespace("rainbow", quietly=TRUE)) |
212 | stop("Functional boxplot requires the rainbow package") | |
213 | ||
214 | series_matrix = rbind( | |
215 | data$getSeries(fil$neighb_indices), data$getSeries(fil$neighb_indices+1) ) | |
216 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) | |
217 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
218 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Temps (heures)", ylab="PM10", | |
219 | plotlegend=FALSE, lwd=2) | |
220 | ||
fa8078f9 BA |
221 | # "Magic" found at http://stackoverflow.com/questions/13842560/get-xlim-from-a-plot-in-r |
222 | usr <- par("usr") | |
223 | yr <- (usr[4] - usr[3]) / 27 | |
98e958ca | 224 | par(new=TRUE) |
a5a3a294 | 225 | plot(c(data$getSerie(fil$index),data$getSerie(fil$index+1)), type="l", lwd=2, lty=2, |
fa8078f9 | 226 | ylim=c(usr[3] + yr, usr[4] - yr), xlab="", ylab="") |
d5fe1f30 | 227 | abline(v=24, lty=2, col=colors()[56]) |
3d69ff21 | 228 | } |
16b1c049 | 229 | |
af3b84f4 | 230 | #' Plot relative conditional variability / absolute variability |
16b1c049 | 231 | #' |
af3b84f4 BA |
232 | #' Draw the relative conditional variability / absolute variability based on filaments |
233 | #' obtained by \code{computeFilaments} | |
16b1c049 BA |
234 | #' |
235 | #' @param data Object return by \code{getData} | |
98e958ca | 236 | #' @param fil Output of \code{computeFilaments} |
16b1c049 BA |
237 | #' |
238 | #' @export | |
98e958ca | 239 | plotRelVar = function(data, fil, ...) |
16b1c049 | 240 | { |
98e958ca BA |
241 | ref_var = c( apply(data$getSeries(fil$neighb_indices),1,sd), |
242 | apply(data$getSeries(fil$neighb_indices+1),1,sd) ) | |
243 | fdays = getNoNA2(data, 1, fil$index-1) | |
244 | global_var = c( apply(data$getSeries(fdays),1,sd), apply(data$getSeries(fdays+1),1,sd) ) | |
16b1c049 | 245 | |
af3b84f4 | 246 | yrange = range(ref_var, global_var) |
16b1c049 | 247 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
af3b84f4 BA |
248 | plot(ref_var, type="l", col=1, lwd=3, ylim=yrange, |
249 | xlab="Temps (heures)", ylab="Écart-type") | |
16b1c049 | 250 | par(new=TRUE) |
98e958ca | 251 | plot(global_var, type="l", col=2, lwd=3, ylim=yrange, xlab="", ylab="") |
16b1c049 BA |
252 | abline(v=24, lty=2, col=colors()[56]) |
253 | } |