'update'
[talweg.git] / pkg / R / plot.R
... / ...
CommitLineData
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
9plotCurves <- 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
30plotError <- 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
76plotPredReal <- 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
102plotSimils <- 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
122plotFbox <- 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
158computeFilaments <- 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
216plotFilamentsBox = 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
247plotRelVar = 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}