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