pkg/R/plot.R

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