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 verbose TRUE to print basic traces (runs beginnings)
  31 #' @param ... Additional parameters for the forecasting models.
  32 #'
  33 #' @return An object of class Forecast
  34 #'
  35 #' @examples
  36 #' ts_data <- system.file("extdata","pm10_mesures_H_loc.csv",package="talweg")
  37 #' exo_data <- system.file("extdata","meteo_extra_noNAs.csv",package="talweg")
  38 #' data <- getData(ts_data, exo_data, limit=200)
  39 #' pred <- computeForecast(data, 100:130, "Persistence", "Zero",
  40 #'   predict_from=8, memory=50, horizon=12, ncores=1)
  41 #' \dontrun{
  42 #' #Sketch for real-time mode:
  43 #' data <- Data$new()
  44 #' forecaster <- MyForecaster$new(myJumpPredictFunc)
  45 #' repeat {
  46 #'   # As soon as daily predictions are available:
  47 #'   data$append(
  48 #'     level_hat=predicted_level,
  49 #'     exo_hat=predicted_exogenous)
  50 #'   # When a day ends:
  51 #'   data$append(
  52 #'     level=observed_level,
  53 #'     exo=observed_exogenous)
  54 #'   # And, at every hour:
  55 #'   data$append(
  56 #'     time=current_hour,
  57 #'     value=current_PM10)
  58 #'   # Finally, a bit before predict_from hour:
  59 #'   pred <- forecaster$predictSerie(data, data$getSize(), ...)
  60 #'   #do_something_with_pred
  61 #' } }
  62 #' @export
  63 computeForecast = function(data, indices, forecaster, pjump, predict_from,
  64     memory=Inf, horizon=length(data$getSerie(1)), ncores=3, verbose=FALSE, ...)
  65 {
  66     # (basic) Arguments sanity checks
  67     predict_from = as.integer(predict_from)[1]
  68     if (! predict_from %in% 1:length(data$getSerie(1)))
  69         stop("predict_from in [1,24] (hours)")
  70     if (hasArg("opera") && !list(...)$opera && memory < Inf)
  71         memory <- Inf #finite memory in training mode makes no sense
  72     horizon = as.integer(horizon)[1]
  73     if (horizon<=predict_from || horizon>length(data$getSerie(1)))
  74         stop("Horizon in [predict_from+1,24] (hours)")
  75     integer_indices = sapply(indices, function(i) dateIndexToInteger(i,data))
  76     if (any(integer_indices<=0 | integer_indices>data$getSize()))
  77         stop("Indices out of range")
  78     if (!is.character(forecaster))
  79         stop("forecaster (name): character")
  80     if (!is.null(pjump) && !is.character(pjump))
  81         stop("pjump (function): character or NULL")
  82
  83     pred = Forecast$new( sapply(indices, function(i) integerIndexToDate(i,data)) )
  84     forecaster_class_name = getFromNamespace(
  85         paste(forecaster,"Forecaster",sep=""), "talweg")
  86
  87     if (!is.null(pjump))
  88         pjump <- getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg")
  89     forecaster = forecaster_class_name$new(pjump)
  90
  91     computeOneForecast <- function(i)
  92     {
  93         if (verbose)
  94             print(paste("Index",i))
  95         list(
  96             "forecast" = forecaster$predictSerie(data,i,memory,predict_from,horizon,...),
  97             "params" = forecaster$getParameters(),
  98             "index" = i )
  99     }
 100
 101     p <-
 102         if (ncores > 1 && requireNamespace("parallel",quietly=TRUE))
 103             parallel::mclapply(integer_indices, computeOneForecast, mc.cores=ncores)
 104         else
 105             lapply(integer_indices, computeOneForecast)
 106
 107     # TODO: find a way to fill pred in //...
 108     for (i in seq_along(integer_indices))
 109     {
 110         pred$append(
 111             forecast = p[[i]]$forecast,
 112             params = p[[i]]$params,
 113             index_in_data = p[[i]]$index
 114         )
 115     }
 116     pred
 117 }