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