X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2Fplot.R;fp=pkg%2FR%2Fplot.R;h=0000000000000000000000000000000000000000;hb=e64fcae579a969d17da0c05f02a8c67a626b9102;hp=eb0c81afe7c49a2772c35404ef7108e87dadde91;hpb=1c74a8eca74347b98d6cf4d0e53c02687dda7422;p=talweg.git diff --git a/pkg/R/plot.R b/pkg/R/plot.R deleted file mode 100644 index eb0c81a..0000000 --- a/pkg/R/plot.R +++ /dev/null @@ -1,250 +0,0 @@ -#' Plot curves -#' -#' Plot a range of curves in data -#' -#' @param data Object of class Data -#' @param indices Range of indices (integers or dates) -#' -#' @export -plotCurves <- function(data, indices=seq_len(data$getSize())) -{ - series = data$getSeries(indices) - yrange = quantile(series, probs=c(0.025,0.975), na.rm=TRUE) - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) - for (i in seq_along(indices)) - { - plot(series[,i], type="l", ylim=yrange, - xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10","")) - if (i < length(indices)) - par(new=TRUE) - } -} - -#' Plot error -#' -#' Draw error graphs, potentially from several runs of \code{computeForecast} -#' -#' @param err Error as returned by \code{computeError} -#' @param cols Colors for each error (default: 1,2,3,...) -#' -#' @seealso \code{\link{plotCurves}}, \code{\link{plotPredReal}}, -#' \code{\link{plotSimils}}, \code{\link{plotFbox}}, -#' \code{\link{computeFilaments}, }\code{\link{plotFilamentsBox}}, \code{\link{plotRelVar}} -#' -#' @export -plotError <- function(err, cols=seq_along(err)) -{ - if (!is.null(err$abs)) - err = list(err) - par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) - L = length(err) - yrange = range( sapply(1:L, function(i) ( err[[i]]$abs$day ) ), na.rm=TRUE ) - for (i in seq_len(L)) - { - plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), - ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) - if (i < L) - par(new=TRUE) - } - yrange = range( sapply(1:L, function(i) ( err[[i]]$abs$indices ) ), na.rm=TRUE ) - for (i in seq_len(L)) - { - plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), - ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) - if (i < L) - par(new=TRUE) - } - yrange = range( sapply(1:L, function(i) ( err[[i]]$MAPE$day ) ), na.rm=TRUE ) - for (i in seq_len(L)) - { - plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), - ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) - if (i < L) - par(new=TRUE) - } - yrange = range( sapply(1:L, function(i) ( err[[i]]$MAPE$indices ) ), na.rm=TRUE ) - for (i in seq_len(L)) - { - plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), - ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) - if (i < L) - par(new=TRUE) - } -} - -#' Plot measured / predicted -#' -#' Plot measured curve (in black) and predicted curve (in blue) -#' -#' @param data Object return by \code{getData} -#' @param pred Object as returned by \code{computeForecast} -#' @param index Index in forecasts (integer or date) -#' -#' @export -plotPredReal <- function(data, pred, index) -{ - horizon = length(pred$getSerie(1)) - measure = data$getSerie( pred$getIndexInData(index)+1 )[1:horizon] - prediction = pred$getSerie(index) - yrange = range(measure, prediction) - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=3) - plot(measure, type="l", ylim=yrange, xlab="Temps (en heures)", ylab="PM10") - par(new=TRUE) - plot(prediction, type="l", col="#0000FF", ylim=yrange, xlab="", ylab="") -} - -#' Plot similarities -#' -#' Plot histogram of similarities (weights) -#' -#' @param pred Object as returned by \code{computeForecast} -#' @param index Index in forecasts (integer or date) -#' -#' @export -plotSimils <- function(pred, index) -{ - weights = pred$getParams(index)$weights - if (is.null(weights)) - stop("plotSimils only works on 'Neighbors' forecasts") - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) - hist(pred$getParams(index)$weights, nclass=20, xlab="Poids", ylab="Effectif") -} - -#' Functional boxplot -#' -#' Draw the functional boxplot on the left, and bivariate plot on the right -#' -#' @param data Object return by \code{getData} -#' @param indices integer or date indices to process -#' @param plot_bivariate Should the bivariate plot appear? -#' -#' @export -plotFbox <- function(data, indices=seq_len(data$getSize())) -{ - if (!requireNamespace("rainbow", quietly=TRUE)) - stop("Functional boxplot requires the rainbow package") - - series_matrix = data$getSeries(indices) - # Remove series with NAs - no_NAs_indices = sapply( 1:ncol(series_matrix), - function(i) all(!is.na(series_matrix[,i])) ) - series_matrix = series_matrix[,no_NAs_indices] - - series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) - rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Temps (heures)", ylab="PM10", - plotlegend=FALSE, lwd=2) - rainbow::fboxplot(series_fds, "bivariate", "hdr", plotlegend=FALSE) -} - -#' Compute filaments -#' -#' Get similar days in the past, as black as distances are small -#' -#' @param data Object as returned by \code{getData} -#' @param pred Object of class Forecast -#' @param index Index in forecast (integer or date) -#' @param limit Number of neighbors to consider -#' @param plot Should the result be plotted? -#' -#' @return A list with -#' \itemize{ -#' \item index : index of the current serie ('today') -#' \item neighb_indices : indices of its neighbors -#' \item colors : colors of neighbors curves (shades of gray) -#' } -#' -#' @export -computeFilaments <- function(data, pred, index, limit=60, plot=TRUE) -{ - ref_serie = data$getCenteredSerie( pred$getIndexInData(index) ) - if (any(is.na(ref_serie))) - stop("computeFilaments requires a serie without NAs") - - # Compute colors for each neighbor (from darkest to lightest) - sorted_dists = sort(-log(pred$getParams(index)$weights), index.return=TRUE) - nn = min(limit, length(sorted_dists$x)) - min_dist = min(sorted_dists$x[1:nn]) - max_dist = max(sorted_dists$x[1:nn]) - color_values = floor(19.5*(sorted_dists$x[1:nn]-min_dist)/(max_dist-min_dist)) + 1 - colors = gray.colors(20,0.1,0.9)[color_values] #TODO: 20 == magic number - - if (plot) - { - # Complete series with (past and present) tomorrows - ref_serie = c(ref_serie, data$getCenteredSerie( pred$getIndexInData(index)+1 )) - centered_series = rbind( - data$getCenteredSeries( pred$getParams(index)$indices ), - data$getCenteredSeries( pred$getParams(index)$indices+1 ) ) - yrange = range( ref_serie, quantile(centered_series, probs=c(0.025,0.975), na.rm=TRUE) ) - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2) - for (i in nn:1) - { - plot(centered_series[,sorted_dists$ix[i]], ylim=yrange, type="l", col=colors[i], - xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10 centré","")) - par(new=TRUE) - } - # Also plot ref curve, in red - plot(ref_serie, ylim=yrange, type="l", col="#FF0000", xlab="", ylab="") - abline(v=24, lty=2, col=colors()[56], lwd=1) - } - - list( - "index"=pred$getIndexInData(index), - "neighb_indices"=pred$getParams(index)$indices[sorted_dists$ix[1:nn]], - "colors"=colors) -} - -#' Functional boxplot on filaments -#' -#' Draw the functional boxplot on filaments obtained by \code{computeFilaments} -#' -#' @param data Object return by \code{getData} -#' @param fil Output of \code{computeFilaments} -#' -#' @export -plotFilamentsBox = function(data, fil, ...) -{ - if (!requireNamespace("rainbow", quietly=TRUE)) - stop("Functional boxplot requires the rainbow package") - - series_matrix = rbind( - data$getSeries(fil$neighb_indices), data$getSeries(fil$neighb_indices+1) ) - series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix) - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) - rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Temps (heures)", ylab="PM10", - plotlegend=FALSE, lwd=2) - - # "Magic" found at http://stackoverflow.com/questions/13842560/get-xlim-from-a-plot-in-r - usr <- par("usr") - yr <- (usr[4] - usr[3]) / 27 - par(new=TRUE) - plot(c(data$getSerie(fil$index),data$getSerie(fil$index+1)), type="l", lwd=2, lty=2, - ylim=c(usr[3] + yr, usr[4] - yr), xlab="", ylab="") - abline(v=24, lty=2, col=colors()[56]) -} - -#' Plot relative conditional variability / absolute variability -#' -#' Draw the relative conditional variability / absolute variability based on filaments -#' obtained by \code{computeFilaments} -#' -#' @param data Object return by \code{getData} -#' @param fil Output of \code{computeFilaments} -#' -#' @export -plotRelVar = function(data, fil, ...) -{ - ref_var = c( apply(data$getSeries(fil$neighb_indices),1,sd), - apply(data$getSeries(fil$neighb_indices+1),1,sd) ) - fdays = getNoNA2(data, 1, fil$index-1) - global_var = c( apply(data$getSeries(fdays),1,sd), apply(data$getSeries(fdays+1),1,sd) ) - - yrange = range(ref_var, global_var) - par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5) - plot(ref_var, type="l", col=1, lwd=3, ylim=yrange, - xlab="Temps (heures)", ylab="Écart-type") - par(new=TRUE) - plot(global_var, type="l", col=2, lwd=3, ylim=yrange, xlab="", ylab="") - abline(v=24, lty=2, col=colors()[56]) -}