Commit | Line | Data |
---|---|---|
e030a6e3 BA |
1 | #' @title plot curves |
2 | #' | |
3 | #' @description Plot a range of curves in data | |
4 | #' | |
5 | #' @param data Object of class Data | |
6 | #' @param indices Range of indices (integers or dates) | |
7 | #' | |
8 | #' @export | |
9 | plotCurves <- function(data, indices) | |
10 | { | |
fa8078f9 | 11 | yrange = quantile( range( sapply( indices, function(i) { |
e030a6e3 BA |
12 | serie = c(data$getCenteredSerie(i)) |
13 | if (!all(is.na(serie))) | |
14 | range(serie, na.rm=TRUE) | |
15 | c() | |
fa8078f9 | 16 | }) ), probs=c(0.05,0.95) ) |
e030a6e3 BA |
17 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
18 | for (i in seq_along(indices)) | |
19 | { | |
20 | plot(data$getSerie(indices[i]), type="l", ylim=yrange, | |
21 | xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10","")) | |
22 | if (ii < length(indices)) | |
23 | par(new=TRUE) | |
24 | } | |
25 | } | |
26 | ||
3d69ff21 BA |
27 | #' @title plot measured / predicted |
28 | #' | |
29 | #' @description Plot measured curve (in black) and predicted curve (in red) | |
30 | #' | |
31 | #' @param data Object return by \code{getData} | |
32 | #' @param pred Object as returned by \code{getForecast} | |
33 | #' @param index Index in forecasts | |
34 | #' | |
35 | #' @export | |
36 | plotPredReal <- function(data, pred, index) | |
37 | { | |
38 | horizon = length(pred$getSerie(1)) | |
4e95ec8f | 39 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=3) |
3d69ff21 BA |
40 | measure = data$getSerie(pred$getIndexInData(index)+1)[1:horizon] |
41 | yrange = range( pred$getSerie(index), measure ) | |
4e95ec8f | 42 | plot(measure, type="l", ylim=yrange, xlab="Temps (en heures)", ylab="PM10") |
3d69ff21 | 43 | par(new=TRUE) |
4e95ec8f | 44 | plot(pred$getSerie(index), type="l", col="#0000FF", ylim=yrange, xlab="", ylab="") |
3d69ff21 BA |
45 | } |
46 | ||
fa8078f9 | 47 | #' @title Compute filaments |
3d69ff21 | 48 | #' |
fa8078f9 | 49 | #' @description Get similar days in the past + "past tomorrow", as black as distances are small |
3d69ff21 BA |
50 | #' |
51 | #' @param data Object as returned by \code{getData} | |
52 | #' @param index Index in data | |
53 | #' @param limit Number of neighbors to consider | |
fa8078f9 | 54 | #' @param plot Should the result be plotted? |
3d69ff21 BA |
55 | #' |
56 | #' @export | |
fa8078f9 | 57 | computeFilaments <- function(data, index, limit=60, plot=TRUE) |
3d69ff21 BA |
58 | { |
59 | index = dateIndexToInteger(index, data) | |
60 | ref_serie = data$getCenteredSerie(index) | |
61 | if (any(is.na(ref_serie))) | |
fa8078f9 | 62 | stop("computeFilaments requires a serie without NAs") |
3d69ff21 BA |
63 | L = length(ref_serie) |
64 | first_day = ifelse(length(data$getCenteredSerie(1)<L), 2, 1) | |
65 | distances = sapply(first_day:(index-1), function(i) { | |
66 | sqrt( sum( (ref_serie - data$getCenteredSerie(i))^2 ) / L ) | |
67 | }) | |
68 | # HACK to suppress NA effect while keeping indexation | |
69 | distances[is.na(distances)] = max(distances,na.rm=TRUE) + 1 | |
70 | indices = sort(distances, index.return=TRUE)$ix[1:min(limit,index-first_day)] | |
fa8078f9 | 71 | yrange = quantile( range( ref_serie, sapply( indices, function(i) { |
3d69ff21 BA |
72 | index = i - first_day + 1 |
73 | serie = c(data$getCenteredSerie(index), data$getCenteredSerie(index+1)) | |
74 | if (!all(is.na(serie))) | |
e5aa669a | 75 | return (range(serie, na.rm=TRUE)) |
e030a6e3 | 76 | c() |
fa8078f9 | 77 | }) ), probs=c(0.1,0.9) ) |
3d69ff21 | 78 | grays = gray.colors(20, 0.1, 0.9) #TODO: 20 == magic number |
fa8078f9 BA |
79 | color_values = floor( 20.5 * distances[indices] / (1+max(distances[indices])) ) |
80 | plot_order = sort(color_values, index.return=TRUE)$ix | |
81 | colors = c(grays[ color_values[plot_order] ], "#FF0000") | |
82 | if (plot) | |
3d69ff21 | 83 | { |
fa8078f9 BA |
84 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) |
85 | for ( i in c(plot_order,length(indices)+1) ) | |
86 | { | |
87 | ind = ifelse(i<=length(indices), indices[i] - first_day + 1, index) | |
88 | plot(c(data$getCenteredSerie(ind),data$getCenteredSerie(ind+1)), | |
89 | ylim=yrange, type="l", col=colors[i], | |
90 | xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10 centré","")) | |
91 | if (i <= length(indices)) | |
92 | par(new=TRUE) | |
93 | } | |
3d69ff21 | 94 | } |
fa8078f9 | 95 | list("indices"=c(indices[plot_order]-first_day+1,index), "colors"=colors) |
3d69ff21 BA |
96 | } |
97 | ||
98 | #' @title Plot similarities | |
99 | #' | |
100 | #' @description Plot histogram of similarities (weights) | |
101 | #' | |
102 | #' @param pred Object as returned by \code{getForecast} | |
103 | #' @param index Index in forecasts (not in data) | |
104 | #' | |
105 | #' @export | |
106 | plotSimils <- function(pred, index) | |
107 | { | |
108 | weights = pred$getParams(index)$weights | |
109 | if (is.null(weights)) | |
110 | stop("plotSimils only works on 'Neighbors' forecasts") | |
4e95ec8f BA |
111 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) |
112 | hist(pred$getParams(index)$weights, nclass=20, xlab="Poids", ylab="Effectif") | |
3d69ff21 BA |
113 | } |
114 | ||
115 | #' @title Plot error | |
116 | #' | |
117 | #' @description Draw error graphs, potentially from several runs of \code{getForecast} | |
118 | #' | |
119 | #' @param err Error as returned by \code{getError} | |
09cf9c19 | 120 | #' @param cols Colors for each error (default: 1,2,3,...) |
3d69ff21 BA |
121 | #' |
122 | #' @seealso \code{\link{plotPredReal}}, \code{\link{plotFilaments}}, \code{\link{plotSimils}} | |
123 | #' \code{\link{plotFbox}} | |
124 | #' | |
125 | #' @export | |
09cf9c19 | 126 | plotError <- function(err, cols=seq_along(err)) |
3d69ff21 | 127 | { |
09cf9c19 BA |
128 | if (!is.null(err$abs)) |
129 | err = list(err) | |
4e95ec8f | 130 | par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) |
3d69ff21 BA |
131 | L = length(err) |
132 | yrange = range( sapply(1:L, function(index) ( err[[index]]$abs$day ) ), na.rm=TRUE ) | |
133 | for (i in seq_len(L)) | |
134 | { | |
135 | plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), | |
09cf9c19 | 136 | ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
137 | if (i < L) |
138 | par(new=TRUE) | |
139 | } | |
140 | yrange = range( sapply(1:L, function(index) ( err[[index]]$abs$indices ) ), na.rm=TRUE ) | |
141 | for (i in seq_len(L)) | |
142 | { | |
143 | plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), | |
09cf9c19 | 144 | ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
145 | if (i < L) |
146 | par(new=TRUE) | |
147 | } | |
148 | yrange = range( sapply(1:L, function(index) ( err[[index]]$MAPE$day ) ), na.rm=TRUE ) | |
149 | for (i in seq_len(L)) | |
150 | { | |
151 | plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), | |
09cf9c19 | 152 | ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
153 | if (i < L) |
154 | par(new=TRUE) | |
155 | } | |
156 | yrange = range( sapply(1:L, function(index) ( err[[index]]$MAPE$indices ) ), na.rm=TRUE ) | |
157 | for (i in seq_len(L)) | |
158 | { | |
159 | plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), | |
09cf9c19 | 160 | ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) |
3d69ff21 BA |
161 | if (i < L) |
162 | par(new=TRUE) | |
163 | } | |
164 | } | |
165 | ||
166 | #' @title Functional boxplot | |
167 | #' | |
168 | #' @description Draw the functional boxplot on the left, and bivariate plot on the right | |
169 | #' | |
170 | #' @param data Object return by \code{getData} | |
171 | #' @param fiter Optional filter: return TRUE on indices to process | |
fa8078f9 | 172 | #' @param plot_bivariate Should the bivariate plot appear? |
3d69ff21 BA |
173 | #' |
174 | #' @export | |
fa8078f9 | 175 | plotFbox <- function(data, filter=function(index) TRUE, plot_bivariate=TRUE) |
3d69ff21 BA |
176 | { |
177 | if (!requireNamespace("rainbow", quietly=TRUE)) | |
178 | stop("Functional boxplot requires the rainbow package") | |
179 | ||
180 | start_index = 1 | |
181 | end_index = data$getSize() | |
182 | if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2))) | |
183 | { | |
184 | # Shifted start (7am, or 1pm, or...) | |
185 | start_index = 2 | |
186 | end_index = data$getSize() - 1 | |
187 | } | |
188 | ||
189 | series_matrix = sapply(start_index:end_index, function(index) { | |
190 | as.matrix(data$getSerie(index)) | |
191 | }) | |
192 | # Remove NAs. + filter TODO: merge with previous step: only one pass required... | |
193 | nas_indices = seq_len(ncol(series_matrix))[ sapply( 1:ncol(series_matrix), | |
194 | function(index) ( !filter(index) || any(is.na(series_matrix[,index])) ) ) ] | |
195 | series_matrix = series_matrix[,-nas_indices] | |
196 | ||
197 | series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) | |
fa8078f9 BA |
198 | if (plot_bivariate) |
199 | par(mfrow=c(1,2)) | |
200 | par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) | |
3d69ff21 BA |
201 | rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Temps (heures)", ylab="PM10", |
202 | plotlegend=FALSE, lwd=2) | |
fa8078f9 BA |
203 | if (plot_bivariate) |
204 | rainbow::fboxplot(series_fds, "bivariate", "hdr", plotlegend=FALSE) | |
205 | } | |
206 | ||
207 | #' @title Functional boxplot on filaments | |
208 | #' | |
209 | #' @description Draw the functional boxplot on filaments obtained by \code{computeFilaments} | |
210 | #' | |
211 | #' @param data Object return by \code{getData} | |
212 | #' @param indices Indices as output by \code{computeFilaments} | |
213 | #' | |
214 | #' @export | |
215 | plotFilamentsBox = function(data, indices, ...) | |
216 | { | |
217 | past_neighbs_indices = head(indices,-1) | |
218 | plotFbox(data, function(i) i %in% past_neighbs_indices, plot_bivariate=FALSE) | |
219 | par(new=TRUE) | |
220 | # "Magic" found at http://stackoverflow.com/questions/13842560/get-xlim-from-a-plot-in-r | |
221 | usr <- par("usr") | |
222 | yr <- (usr[4] - usr[3]) / 27 | |
223 | plot(data$getSerie(tail(indices,1)), type="l", lwd=2, lty=2, | |
224 | ylim=c(usr[3] + yr, usr[4] - yr), xlab="", ylab="") | |
3d69ff21 | 225 | } |