'update'
[talweg.git] / pkg / R / computeForecast.R
CommitLineData
af3b84f4 1#' Compute forecast
3d69ff21 2#'
4f3fdbb8
BA
3#' Predict time-series curves ("today" from predict_from to horizon) at the selected days
4#' indices ("today" from 1am to predict_from-1). This function just runs a loop over all
5#' requested indices, and stores the individual forecasts into a Forecast object.
c4c329f6 6#'
102bcfda 7#' @param data Object of class Data, output of \code{getData()}.
2057c793 8#' @param indices Indices where to forecast (the day after); integers relative to the
e169b5d5
BA
9#' beginning of data, or (convertible to) Date objects.
10#' @param forecaster Name of the main forecaster; more details: ?F_<forecastername>
4f3fdbb8
BA
11#' \itemize{
12#' \item Persistence : use last (similar) day
13#' \item Neighbors : weighted similar days
14#' \item Average : average curve of all same day-in-week
15#' \item Zero : just output 0 (benchmarking purpose)
16#' }
e169b5d5
BA
17#' @param pjump Function to predict the jump at the interface between two days;
18#' more details: ?J_<functionname>
4f3fdbb8
BA
19#' \itemize{
20#' \item Persistence : use last (similar) day
21#' \item Neighbors: re-use the weights from F_Neighbors
22#' \item Zero: just output 0 (no adjustment)
23#' }
24#' If pjump=NULL, then no adjustment is performed (output of \code{predictShape()} is
25#' used directly).
26#' @param predict_from First time step to predict.
e169b5d5 27#' @param memory Data depth (in days) to be used for prediction.
4f3fdbb8 28#' @param horizon Last time step to predict.
e169b5d5
BA
29#' @param ncores Number of cores for parallel execution (1 to disable).
30#' @param ... Additional parameters for the forecasting models.
3d69ff21 31#'
a66a84b5 32#' @return An object of class Forecast
3d69ff21
BA
33#'
34#' @examples
e169b5d5
BA
35#' ts_data <- system.file("extdata","pm10_mesures_H_loc.csv",package="talweg")
36#' exo_data <- system.file("extdata","meteo_extra_noNAs.csv",package="talweg")
4f3fdbb8 37#' data <- getData(ts_data, exo_data, limit=200)
3ddf1c12 38#' pred <- computeForecast(data, 100:130, "Persistence", "Zero",
4f3fdbb8
BA
39#' predict_from=8, memory=50, horizon=12, ncores=1)
40#' \dontrun{
41#' #Sketch for real-time mode:
e169b5d5 42#' data <- Data$new()
e169b5d5 43#' forecaster <- MyForecaster$new(myJumpPredictFunc)
3d69ff21 44#' repeat {
4f3fdbb8 45#' # As soon as daily predictions are available:
c1be9898 46#' data$append(
4f3fdbb8
BA
47#' level_hat=predicted_level,
48#' exo_hat=predicted_exogenous)
49#' # When a day ends:
50#' data$append(
51#' level=observed_level,
52#' exo=observed_exogenous)
53#' # And, at every hour:
54#' data$append(
55#' time=current_hour,
56#' value=current_PM10)
57#' # Finally, a bit before predict_from hour:
102bcfda 58#' pred <- forecaster$predictSerie(data, data$getSize(), ...)
3d69ff21 59#' #do_something_with_pred
4f3fdbb8 60#' } }
3d69ff21 61#' @export
d2ab47a7
BA
62computeForecast = function(data, indices, forecaster, pjump, predict_from,
63 memory=Inf, horizon=length(data$getSerie(1)), ncores=3, ...)
3d69ff21 64{
e030a6e3 65 # (basic) Arguments sanity checks
d2ab47a7
BA
66 predict_from = as.integer(predict_from)[1]
67 if (! predict_from %in% 1:length(data$getSerie(1)))
68 stop("predict_from in [1,24] (hours)")
638f27f4
BA
69 if (hasArg("opera") && !list(...)$opera && memory < Inf)
70 memory <- Inf #finite memory in training mode makes no sense
3d69ff21 71 horizon = as.integer(horizon)[1]
d2ab47a7
BA
72 if (horizon<=predict_from || horizon>length(data$getSerie(1)))
73 stop("Horizon in [predict_from+1,24] (hours)")
98e958ca 74 integer_indices = sapply(indices, function(i) dateIndexToInteger(i,data))
a66a84b5 75 if (any(integer_indices<=0 | integer_indices>data$getSize()))
3d69ff21 76 stop("Indices out of range")
4f3fdbb8
BA
77 if (!is.character(forecaster))
78 stop("forecaster (name): character")
79 if (!is.null(pjump) && !is.character(pjump))
80 stop("pjump (function): character or NULL")
3d69ff21 81
98e958ca 82 pred = Forecast$new( sapply(indices, function(i) integerIndexToDate(i,data)) )
72b9c501
BA
83 forecaster_class_name = getFromNamespace(
84 paste(forecaster,"Forecaster",sep=""), "talweg")
4f3fdbb8
BA
85
86 if (!is.null(pjump))
87 pjump <- getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg")
88 forecaster = forecaster_class_name$new(pjump)
5e838b3e 89
1e8327df 90 computeOneForecast <- function(i)
a866acb3 91 {
1e8327df
BA
92 list(
93 "forecast" = forecaster$predictSerie(data,i,memory,predict_from,horizon,...),
94 "params" = forecaster$getParameters(),
95 "index" = i )
a866acb3 96 }
5e838b3e 97
1e8327df
BA
98 p <-
99 if (ncores > 1 && requireNamespace("parallel",quietly=TRUE))
100 parallel::mclapply(integer_indices, computeOneForecast, mc.cores=ncores)
101 else
102 lapply(integer_indices, computeOneForecast)
103
ee8b1b4e
BA
104 # TODO: find a way to fill pred in //...
105 for (i in seq_along(integer_indices))
106 {
107 pred$append(
72b9c501
BA
108 forecast = p[[i]]$forecast,
109 params = p[[i]]$params,
110 index_in_data = p[[i]]$index
ee8b1b4e
BA
111 )
112 }
25b75559 113 pred
3d69ff21 114}