0f895bd16483795302623b53ea983f3ab9c2cd92
[talweg.git] / 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 }