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 }