'update'
[talweg.git] / pkg / R / computeForecast.R
1 #' Compute forecast
2 #'
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.
6 #'
7 #' @param data Object of class Data, output of \code{getData()}.
8 #' @param indices Indices where to forecast (the day after); integers relative to the
9 #' beginning of data, or (convertible to) Date objects.
10 #' @param forecaster Name of the main forecaster; more details: ?F_<forecastername>
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 #' }
17 #' @param pjump Function to predict the jump at the interface between two days;
18 #' more details: ?J_<functionname>
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.
27 #' @param memory Data depth (in days) to be used for prediction.
28 #' @param horizon Last time step to predict.
29 #' @param ncores Number of cores for parallel execution (1 to disable).
30 #' @param ... Additional parameters for the forecasting models.
31 #'
32 #' @return An object of class Forecast
33 #'
34 #' @examples
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")
37 #' data <- getData(ts_data, exo_data, limit=200)
38 #' pred <- computeForecast(data, 100:130, "Persistence", "Zero",
39 #' predict_from=8, memory=50, horizon=12, ncores=1)
40 #' \dontrun{
41 #' #Sketch for real-time mode:
42 #' data <- Data$new()
43 #' forecaster <- MyForecaster$new(myJumpPredictFunc)
44 #' repeat {
45 #' # As soon as daily predictions are available:
46 #' data$append(
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:
58 #' pred <- forecaster$predictSerie(data, data$getSize(), ...)
59 #' #do_something_with_pred
60 #' } }
61 #' @export
62 computeForecast = function(data, indices, forecaster, pjump, predict_from,
63 memory=Inf, horizon=length(data$getSerie(1)), ncores=3, ...)
64 {
65 # (basic) Arguments sanity checks
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)")
69 if (hasArg("opera") && !list(...)$opera && memory < Inf)
70 memory <- Inf #finite memory in training mode makes no sense
71 horizon = as.integer(horizon)[1]
72 if (horizon<=predict_from || horizon>length(data$getSerie(1)))
73 stop("Horizon in [predict_from+1,24] (hours)")
74 integer_indices = sapply(indices, function(i) dateIndexToInteger(i,data))
75 if (any(integer_indices<=0 | integer_indices>data$getSize()))
76 stop("Indices out of range")
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")
81
82 pred = Forecast$new( sapply(indices, function(i) integerIndexToDate(i,data)) )
83 forecaster_class_name = getFromNamespace(
84 paste(forecaster,"Forecaster",sep=""), "talweg")
85
86 if (!is.null(pjump))
87 pjump <- getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg")
88 forecaster = forecaster_class_name$new(pjump)
89
90 computeOneForecast <- function(i)
91 {
92 list(
93 "forecast" = forecaster$predictSerie(data,i,memory,predict_from,horizon,...),
94 "params" = forecaster$getParameters(),
95 "index" = i )
96 }
97
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
104 # TODO: find a way to fill pred in //...
105 for (i in seq_along(integer_indices))
106 {
107 pred$append(
108 forecast = p[[i]]$forecast,
109 params = p[[i]]$params,
110 index_in_data = p[[i]]$index
111 )
112 }
113 pred
114 }