Commit | Line | Data |
---|---|---|
98e958ca | 1 | #' Plot curves |
e030a6e3 | 2 | #' |
102bcfda | 3 | #' Plot a range of curves in data. |
e030a6e3 | 4 | #' |
102bcfda | 5 | #' @inheritParams computeError |
e030a6e3 BA |
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 | 13 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
4f3fdbb8 | 14 | matplot(series, type="l", ylim=yrange, xlab="Time (hours)", ylab="PM10") |
e030a6e3 BA |
15 | } |
16 | ||
af3b84f4 | 17 | #' Plot error |
3d69ff21 | 18 | #' |
102bcfda | 19 | #' Draw error graphs, potentially from several runs of \code{computeForecast()}. |
3d69ff21 | 20 | #' |
102bcfda | 21 | #' @param err Error as returned by \code{computeError()} |
09cf9c19 | 22 | #' @param cols Colors for each error (default: 1,2,3,...) |
81cb8b9e BA |
23 | <<<<<<< HEAD |
24 | ======= | |
3a38473a | 25 | #' @param agg Aggregation level ("day", "week" or "month") |
81cb8b9e | 26 | >>>>>>> 3a38473a435221cc15b6215b0d6677cd370dc2d6 |
3d69ff21 | 27 | #' |
98e958ca | 28 | #' @seealso \code{\link{plotCurves}}, \code{\link{plotPredReal}}, |
72b9c501 BA |
29 | #' \code{\link{plotSimils}}, \code{\link{plotFbox}}, \code{\link{computeFilaments}}, |
30 | #' \code{\link{plotFilamentsBox}}, \code{\link{plotRelVar}} | |
3d69ff21 BA |
31 | #' |
32 | #' @export | |
aa5397f1 | 33 | plotError <- function(err, cols=seq_along(err), agg="day") |
3d69ff21 | 34 | { |
09cf9c19 BA |
35 | if (!is.null(err$abs)) |
36 | err = list(err) | |
10886062 | 37 | par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
3d69ff21 | 38 | L = length(err) |
9b9bb2d4 | 39 | |
aa5397f1 BA |
40 | yrange = range( sapply(1:L, function(i) err[[i]]$abs$day), na.rm=TRUE ) |
41 | matplot(sapply( seq_len(L), function(i) err[[i]]$abs$day ), type="l", | |
42 | xlab="Time (hours)", ylab="Mean |y - y_hat|", ylim=yrange, col=cols, lwd=2, lty=1) | |
9b9bb2d4 | 43 | |
aa5397f1 BA |
44 | agg_curves <- sapply( seq_len(L), function(i) { |
45 | curve <- err[[i]]$abs$indices | |
46 | delta <- if (agg=="day") 1 else if (agg=="week") 7 else if (agg=="month") 30 | |
47 | vapply( seq(1,length(curve),delta), function(i) { | |
48 | mean(curve[i:(i+delta-1)], na.rm=TRUE) | |
49 | }, vector("double",1), USE.NAMES=FALSE ) | |
50 | }) | |
51 | yrange = range(agg_curves, na.rm=TRUE) | |
52 | matplot(agg_curves, type="l", xlab=paste("Time (",agg,"s)", sep=""), | |
53 | ylab="Mean |y - y_hat|", ylim=yrange, col=cols, lwd=2, lty=1) | |
9b9bb2d4 | 54 | |
aa5397f1 BA |
55 | yrange = range( sapply(1:L, function(i) err[[i]]$MAPE$day), na.rm=TRUE ) |
56 | matplot(sapply( seq_len(L), function(i) err[[i]]$MAPE$day ), type="l", | |
57 | xlab="Time (hours)", ylab="Mean MAPE", ylim=yrange, col=cols, lwd=2, lty=1) | |
9b9bb2d4 | 58 | |
aa5397f1 BA |
59 | agg_curves <- sapply( seq_len(L), function(i) { |
60 | curve <- err[[i]]$MAPE$indices | |
61 | delta <- if (agg=="day") 1 else if (agg=="week") 7 else if (agg=="month") 30 | |
62 | vapply( seq(1,length(curve),delta), function(i) { | |
63 | mean(curve[i:(i+delta-1)], na.rm=TRUE) | |
64 | }, vector("double",1), USE.NAMES=FALSE ) | |
65 | }) | |
66 | yrange = range(agg_curves, na.rm=TRUE) | |
67 | matplot(agg_curves, type="l", xlab=paste("Time (",agg,"s)", sep=""), | |
68 | ylab="Mean MAPE", ylim=yrange, col=cols, lwd=2, lty=1) | |
3d69ff21 BA |
69 | } |
70 | ||
98e958ca BA |
71 | #' Plot measured / predicted |
72 | #' | |
102bcfda | 73 | #' Plot measured curve (in black) and predicted curve (in blue). |
98e958ca | 74 | #' |
102bcfda | 75 | #' @inheritParams computeError |
98e958ca BA |
76 | #' @param index Index in forecasts (integer or date) |
77 | #' | |
78 | #' @export | |
79 | plotPredReal <- function(data, pred, index) | |
80 | { | |
72b9c501 | 81 | prediction = pred$getForecast(index) |
9b9bb2d4 BA |
82 | measure = data$getSerie( pred$getIndexInData(index) )[1:length(pred$getForecast(1))] |
83 | ||
84 | # Remove the common part, where prediction == measure | |
2e0ef04b BA |
85 | dot_mark <- ifelse(prediction[1]==measure[1], |
86 | which.max(seq_along(prediction)[prediction==measure]), 0) | |
9b9bb2d4 BA |
87 | prediction = prediction[(dot_mark+1):length(prediction)] |
88 | measure = measure[(dot_mark+1):length(measure)] | |
89 | ||
98e958ca BA |
90 | yrange = range(measure, prediction) |
91 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=3) | |
4e25de2c | 92 | plot(measure, type="l", ylim=yrange, xlab="Time (hours)", ylab="PM10") |
98e958ca BA |
93 | par(new=TRUE) |
94 | plot(prediction, type="l", col="#0000FF", ylim=yrange, xlab="", ylab="") | |
95 | } | |
96 | ||
97 | #' Plot similarities | |
98 | #' | |
102bcfda | 99 | #' Plot histogram of similarities (weights), for 'Neighbors' method. |
98e958ca | 100 | #' |
102bcfda | 101 | #' @inheritParams computeError |
98e958ca BA |
102 | #' @param index Index in forecasts (integer or date) |
103 | #' | |
104 | #' @export | |
105 | plotSimils <- function(pred, index) | |
106 | { | |
107 | weights = pred$getParams(index)$weights | |
108 | if (is.null(weights)) | |
109 | stop("plotSimils only works on 'Neighbors' forecasts") | |
9b9bb2d4 BA |
110 | par(mfrow=c(1,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
111 | small_weights = weights[ weights < 1/length(weights) ] | |
112 | large_weights = weights[ weights >= 1/length(weights) ] | |
113 | hist(small_weights, nclass=25, main="", xlab="Weight < 1/N", ylab="Count") | |
114 | hist(large_weights, nclass=25, main="", xlab="Weight >= 1/N", ylab="Count") | |
98e958ca BA |
115 | } |
116 | ||
af3b84f4 | 117 | #' Functional boxplot |
3d69ff21 | 118 | #' |
102bcfda | 119 | #' Draw the functional boxplot on the left, and bivariate plot on the right. |
3d69ff21 | 120 | #' |
102bcfda BA |
121 | #' @inheritParams computeError |
122 | #' @inheritParams plotCurves | |
3d69ff21 BA |
123 | #' |
124 | #' @export | |
98e958ca | 125 | plotFbox <- function(data, indices=seq_len(data$getSize())) |
3d69ff21 BA |
126 | { |
127 | if (!requireNamespace("rainbow", quietly=TRUE)) | |
128 | stop("Functional boxplot requires the rainbow package") | |
129 | ||
98e958ca BA |
130 | series_matrix = data$getSeries(indices) |
131 | # Remove series with NAs | |
af3b84f4 BA |
132 | no_NAs_indices = sapply( 1:ncol(series_matrix), |
133 | function(i) all(!is.na(series_matrix[,i])) ) | |
99f83c9a | 134 | series_matrix = series_matrix[,no_NAs_indices] |
3d69ff21 BA |
135 | |
136 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) | |
fa8078f9 | 137 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
4e25de2c | 138 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Time (hours)", ylab="PM10", |
3d69ff21 | 139 | plotlegend=FALSE, lwd=2) |
98e958ca BA |
140 | rainbow::fboxplot(series_fds, "bivariate", "hdr", plotlegend=FALSE) |
141 | } | |
142 | ||
143 | #' Compute filaments | |
144 | #' | |
102bcfda BA |
145 | #' Obtain similar days in the past, and (optionally) plot them -- as black as distances |
146 | #' are small. | |
98e958ca | 147 | #' |
102bcfda | 148 | #' @inheritParams computeError |
8f84543c | 149 | #' @param index Index in forecast (integer or date) |
98e958ca BA |
150 | #' @param limit Number of neighbors to consider |
151 | #' @param plot Should the result be plotted? | |
152 | #' | |
153 | #' @return A list with | |
154 | #' \itemize{ | |
155 | #' \item index : index of the current serie ('today') | |
156 | #' \item neighb_indices : indices of its neighbors | |
157 | #' \item colors : colors of neighbors curves (shades of gray) | |
158 | #' } | |
159 | #' | |
160 | #' @export | |
9b9bb2d4 | 161 | computeFilaments <- function(data, pred, index, limit=60, plot=TRUE) |
98e958ca | 162 | { |
9b9bb2d4 BA |
163 | weights <- pred$getParams(index)$weights |
164 | if (is.null(weights) || is.na(pred$getParams(index)$weights[1])) | |
98e958ca BA |
165 | stop("computeFilaments requires a serie without NAs") |
166 | ||
9b9bb2d4 BA |
167 | nn <- min(limit, length(weights)) |
168 | sorted_dists = sort(-log(weights), index.return=TRUE) | |
169 | # Compute colors for each neighbor (from darkest to lightest), if weights differ | |
170 | if ( any( weights != weights[1] ) ) | |
171 | { | |
172 | min_dist = min(sorted_dists$x[1:nn]) | |
173 | max_dist = max(sorted_dists$x[1:nn]) | |
174 | color_values = floor(19.5*(sorted_dists$x[1:nn]-min_dist)/(max_dist-min_dist)) + 1 | |
175 | colors = gray.colors(20,0.1,0.9)[color_values] #TODO: 20 == magic number | |
176 | } | |
177 | else | |
178 | colors <- rep(colors()[17], length(weights)) | |
98e958ca BA |
179 | |
180 | if (plot) | |
181 | { | |
182 | # Complete series with (past and present) tomorrows | |
d2ab47a7 BA |
183 | ref_serie = c( data$getCenteredSerie( pred$getIndexInData(index)-1 ), |
184 | data$getCenteredSerie( pred$getIndexInData(index) ) ) | |
8f84543c | 185 | centered_series = rbind( |
9e0f25f6 BA |
186 | data$getCenteredSeries( pred$getParams(index)$indices-1 ), |
187 | data$getCenteredSeries( pred$getParams(index)$indices ) ) | |
72b9c501 BA |
188 | yrange = range( ref_serie, |
189 | quantile(centered_series, probs=c(0.025,0.975), na.rm=TRUE) ) | |
98e958ca BA |
190 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) |
191 | for (i in nn:1) | |
192 | { | |
8f84543c | 193 | plot(centered_series[,sorted_dists$ix[i]], ylim=yrange, type="l", col=colors[i], |
4e25de2c | 194 | xlab=ifelse(i==1,"Time (hours)",""), ylab=ifelse(i==1,"Centered PM10","")) |
98e958ca BA |
195 | par(new=TRUE) |
196 | } | |
197 | # Also plot ref curve, in red | |
198 | plot(ref_serie, ylim=yrange, type="l", col="#FF0000", xlab="", ylab="") | |
9b9bb2d4 BA |
199 | dot_mark <- 0.5 + which.max( pred$getForecast(1) == |
200 | data$getSerie( pred$getIndexInData(1) )[1:length(pred$getForecast(1))] ) | |
201 | abline(v=24+dot_mark, lty=2, col=colors()[56], lwd=1) | |
98e958ca BA |
202 | } |
203 | ||
8f84543c BA |
204 | list( |
205 | "index"=pred$getIndexInData(index), | |
206 | "neighb_indices"=pred$getParams(index)$indices[sorted_dists$ix[1:nn]], | |
207 | "colors"=colors) | |
fa8078f9 BA |
208 | } |
209 | ||
af3b84f4 | 210 | #' Functional boxplot on filaments |
fa8078f9 | 211 | #' |
102bcfda | 212 | #' Draw the functional boxplot on filaments obtained by \code{computeFilaments()}. |
fa8078f9 | 213 | #' |
102bcfda | 214 | #' @inheritParams computeError |
98e958ca | 215 | #' @param fil Output of \code{computeFilaments} |
9b9bb2d4 | 216 | #' @param predict_from First predicted time step |
fa8078f9 BA |
217 | #' |
218 | #' @export | |
d2ab47a7 | 219 | plotFilamentsBox = function(data, fil, predict_from) |
fa8078f9 | 220 | { |
98e958ca BA |
221 | if (!requireNamespace("rainbow", quietly=TRUE)) |
222 | stop("Functional boxplot requires the rainbow package") | |
223 | ||
224 | series_matrix = rbind( | |
3fd7377d | 225 | data$getSeries(fil$neighb_indices-1), data$getSeries(fil$neighb_indices) ) |
98e958ca | 226 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) |
3fd7377d | 227 | |
98e958ca | 228 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
4e25de2c | 229 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Time (hours)", ylab="PM10", |
98e958ca BA |
230 | plotlegend=FALSE, lwd=2) |
231 | ||
72b9c501 | 232 | # "Magic": http://stackoverflow.com/questions/13842560/get-xlim-from-a-plot-in-r |
fa8078f9 BA |
233 | usr <- par("usr") |
234 | yr <- (usr[4] - usr[3]) / 27 | |
98e958ca | 235 | par(new=TRUE) |
9b9bb2d4 | 236 | plot(c(data$getSerie(fil$index-1),data$getSerie(fil$index)), type="l", lwd=2, lty=2, |
fa8078f9 | 237 | ylim=c(usr[3] + yr, usr[4] - yr), xlab="", ylab="") |
d2ab47a7 | 238 | abline(v=24+predict_from-0.5, lty=2, col=colors()[56]) |
3d69ff21 | 239 | } |
16b1c049 | 240 | |
af3b84f4 | 241 | #' Plot relative conditional variability / absolute variability |
16b1c049 | 242 | #' |
af3b84f4 | 243 | #' Draw the relative conditional variability / absolute variability based on filaments |
102bcfda | 244 | #' obtained by \code{computeFilaments()}. |
16b1c049 | 245 | #' |
102bcfda BA |
246 | #' @inheritParams computeError |
247 | #' @inheritParams plotFilamentsBox | |
16b1c049 BA |
248 | #' |
249 | #' @export | |
d2ab47a7 | 250 | plotRelVar = function(data, fil, predict_from) |
16b1c049 | 251 | { |
3fd7377d BA |
252 | ref_var = c( apply(data$getSeries(fil$neighb_indices-1),1,sd), |
253 | apply(data$getSeries(fil$neighb_indices),1,sd) ) | |
cf3bb001 | 254 | tdays = .getNoNA2(data, 2, fil$index) |
72b9c501 | 255 | global_var = c( |
cf3bb001 BA |
256 | apply(data$getSeries(tdays-1),1,sd), |
257 | apply(data$getSeries(tdays),1,sd) ) | |
16b1c049 | 258 | |
af3b84f4 | 259 | yrange = range(ref_var, global_var) |
16b1c049 | 260 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
af3b84f4 | 261 | plot(ref_var, type="l", col=1, lwd=3, ylim=yrange, |
4e25de2c | 262 | xlab="Time (hours)", ylab="Standard deviation") |
16b1c049 | 263 | par(new=TRUE) |
98e958ca | 264 | plot(global_var, type="l", col=2, lwd=3, ylim=yrange, xlab="", ylab="") |
d2ab47a7 | 265 | abline(v=24+predict_from-0.5, lty=2, col=colors()[56]) |
16b1c049 | 266 | } |