From 638f27f4296727aff62b56643beb9f42aa5b57ef Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Wed, 26 Apr 2017 16:55:41 +0200 Subject: [PATCH] adapt Bruno method into package, add 'operational' mode --- pkg/R/F_Average.R | 18 +++ pkg/R/F_Neighbors.R | 131 +++++++++++++-------- pkg/R/computeForecast.R | 2 + pkg/R/utils.R | 30 +++-- pkg/tests/testthat/test-Forecaster.R | 4 +- pkg/tests/testthat/test-computeFilaments.R | 42 +++---- 6 files changed, 143 insertions(+), 84 deletions(-) diff --git a/pkg/R/F_Average.R b/pkg/R/F_Average.R index 6cd2d6e..bee1974 100644 --- a/pkg/R/F_Average.R +++ b/pkg/R/F_Average.R @@ -23,6 +23,7 @@ AverageForecaster = R6::R6Class("AverageForecaster", first_day = max(1, today-memory) index <- today nb_no_na_series = 0 + opera = ifelse(hasArg("opera"), list(...)$opera, FALSE) repeat { index = index - 7 @@ -35,6 +36,23 @@ AverageForecaster = R6::R6Class("AverageForecaster", nb_no_na_series = nb_no_na_series + 1 } } + if (!opera) + { + # The same, in the future + index <- today + repeat + { + index = index + 7 + if (index > data$getSize()) + break + serie_on_horizon = data$getCenteredSerie(index)[predict_from:horizon] + if (!any(is.na(serie_on_horizon))) + { + avg = avg + serie_on_horizon + nb_no_na_series = nb_no_na_series + 1 + } + } + } avg / nb_no_na_series } ) diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index f140b0b..02536eb 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -51,21 +51,26 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", return (NA) } - # Determine indices of no-NAs days preceded by no-NAs yerstedays - tdays = .getNoNA2(data, max(today-memory,2), today-1) - # Get optional args local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season? simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo" + opera = ifelse(hasArg("opera"), list(...)$opera, FALSE) #operational mode? + + # Determine indices of no-NAs days preceded by no-NAs yerstedays + tdays = .getNoNA2(data, max(today-memory,2), ifelse(opera,today-1,data$getSize())) + if (!opera) + tdays = setdiff(tdays, today) #always exclude current day + + # Shortcut if window is known if (hasArg("window")) { - return ( private$.predictShapeAux(data, - tdays, today, predict_from, horizon, local, list(...)$window, simtype, TRUE) ) + return ( private$.predictShapeAux(data, tdays, today, predict_from, horizon, + local, list(...)$window, simtype, opera, TRUE) ) } # Indices of similar days for cross-validation; TODO: 20 = magic number cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, - days_in=tdays) + days_in=tdays, operational=opera) # Optimize h : h |--> sum of prediction errors on last N "similar" days errorOnLastNdays = function(window, simtype) @@ -76,7 +81,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", { # mix_strategy is never used here (simtype != "mix"), therefore left blank prediction = private$.predictShapeAux(data, tdays, cv_days[i], predict_from, - horizon, local, window, simtype, FALSE) + horizon, local, window, simtype, opera, FALSE) if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 @@ -110,27 +115,27 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", 1 return( private$.predictShapeAux(data, tdays, today, predict_from, horizon, local, - best_window, simtype, TRUE) ) + best_window, simtype, opera, TRUE) ) } ), private = list( - # Precondition: "today" is full (no NAs) + # Precondition: "yersteday until predict_from-1" is full (no NAs) .predictShapeAux = function(data, tdays, today, predict_from, horizon, local, window, - simtype, final_call) + simtype, opera, final_call) { - tdays_cut = tdays[ tdays <= today-1 ] - if (length(tdays_cut) <= 1) + tdays_cut = tdays[ tdays != today ] + if (length(tdays_cut) == 0) return (NA) if (local) { # TODO: 60 == magic number tdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE, - days_in=tdays_cut) - if (length(tdays) <= 1) - return (NA) - # TODO: 10, 12 == magic numbers - tdays = .getConstrainedNeighbs(today,data,tdays,min_neighbs=10,max_neighbs=12) + days_in=tdays_cut, operational=opera) +# if (length(tdays) <= 1) +# return (NA) + # TODO: 10 == magic number + tdays = .getConstrainedNeighbs(today, data, tdays, min_neighbs=10) if (length(tdays) == 1) { if (final_call) @@ -150,14 +155,19 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # Compute endogen similarities using given window window_endo = ifelse(simtype=="mix", window[1], window) - # Distances from last observed day to days in the past - lastSerie = c( data$getSerie(today-1), - data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()] ) - distances2 = sapply(tdays, function(i) { - delta = lastSerie - c(data$getSerie(i-1), - data$getSerie(i)[if (predict_from>=2) 1:(predict_from-1) else c()]) - sqrt(mean(delta^2)) - }) + # Distances from last observed day to selected days in the past + distances2 <- .computeDistsEndo(data, today, tdays, predict_from) + + if (local) + { + max_neighbs = 12 #TODO: 12 = arbitrary number + if (length(distances2) > max_neighbs) + { + ordering <- order(distances2) + tdays <- tdays[ ordering[1:max_neighbs] ] + distances2 <- distances2[ ordering[1:max_neighbs] ] + } + } simils_endo <- .computeSimils(distances2, window_endo) } @@ -167,24 +177,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # Compute exogen similarities using given window window_exo = ifelse(simtype=="mix", window[2], window) - M = matrix( ncol=1+length(tdays), nrow=1+length(data$getExo(1)) ) - M[,1] = c( data$getLevelHat(today), as.double(data$getExoHat(today)) ) - for (i in seq_along(tdays)) - M[,i+1] = c( data$getLevel(tdays[i]), as.double(data$getExo(tdays[i])) ) - - sigma = cov(t(M)) #NOTE: robust covariance is way too slow - # TODO: 10 == magic number; more robust way == det, or always ginv() - sigma_inv = - if (length(tdays) > 10) - solve(sigma) - else - MASS::ginv(sigma) - - # Distances from last observed day to days in the past - distances2 = sapply(seq_along(tdays), function(i) { - delta = M[,1] - M[,i+1] - delta %*% sigma_inv %*% delta - }) + distances2 <- .computeDistsExo(data, today, tdays) simils_exo <- .computeSimils(distances2, window_exo) } @@ -237,12 +230,13 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # @param min_neighbs Minimum number of points in a neighborhood # @param max_neighbs Maximum number of points in a neighborhood # -.getConstrainedNeighbs = function(today, data, tdays, min_neighbs=10, max_neighbs=12) +.getConstrainedNeighbs = function(today, data, tdays, min_neighbs=10) { levelToday = data$getLevelHat(today) - levelYersteday = data$getLevel(today-1) +# levelYersteday = data$getLevel(today-1) distances = sapply(tdays, function(i) { - sqrt((data$getLevel(i-1)-levelYersteday)^2 + (data$getLevel(i)-levelToday)^2) +# sqrt((data$getLevel(i-1)-levelYersteday)^2 + (data$getLevel(i)-levelToday)^2) + abs(data$getLevel(i)-levelToday) }) #TODO: 1, +1, +3 : magic numbers dist_thresh = 1 @@ -256,12 +250,12 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", dist_thresh = dist_thresh + ifelse(dist_thresh>1,3,1) } tdays = tdays[same_pollution] - max_neighbs = 12 - if (nb_neighbs > max_neighbs) - { - # Keep only max_neighbs closest neighbors - tdays = tdays[ order(distances[same_pollution])[1:max_neighbs] ] - } +# max_neighbs = 12 +# if (nb_neighbs > max_neighbs) +# { +# # Keep only max_neighbs closest neighbors +# tdays = tdays[ order(distances[same_pollution])[1:max_neighbs] ] +# } tdays } @@ -282,3 +276,36 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", } exp(-distances2/(sd_dist*window^2)) } + +.computeDistsEndo <- function(data, today, tdays, predict_from) +{ + lastSerie = c( data$getSerie(today-1), + data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()] ) + sapply(tdays, function(i) { + delta = lastSerie - c(data$getSerie(i-1), + data$getSerie(i)[if (predict_from>=2) 1:(predict_from-1) else c()]) + sqrt(mean(delta^2)) + }) +} + +.computeDistsExo <- function(data, today, tdays) +{ + M = matrix( ncol=1+length(tdays), nrow=1+length(data$getExo(1)) ) + M[,1] = c( data$getLevelHat(today), as.double(data$getExoHat(today)) ) + for (i in seq_along(tdays)) + M[,i+1] = c( data$getLevel(tdays[i]), as.double(data$getExo(tdays[i])) ) + + sigma = cov(t(M)) #NOTE: robust covariance is way too slow + # TODO: 10 == magic number; more robust way == det, or always ginv() + sigma_inv = + if (length(tdays) > 10) + solve(sigma) + else + MASS::ginv(sigma) + + # Distances from last observed day to days in the past + sapply(seq_along(tdays), function(i) { + delta = M[,1] - M[,i+1] + delta %*% sigma_inv %*% delta + }) +} diff --git a/pkg/R/computeForecast.R b/pkg/R/computeForecast.R index ca8badd..ef46dd3 100644 --- a/pkg/R/computeForecast.R +++ b/pkg/R/computeForecast.R @@ -56,6 +56,8 @@ computeForecast = function(data, indices, forecaster, pjump, predict_from, predict_from = as.integer(predict_from)[1] if (! predict_from %in% 1:length(data$getSerie(1))) stop("predict_from in [1,24] (hours)") + if (hasArg("opera") && !list(...)$opera && memory < Inf) + memory <- Inf #finite memory in training mode makes no sense horizon = as.integer(horizon)[1] if (horizon<=predict_from || horizon>length(data$getSerie(1))) stop("Horizon in [predict_from+1,24] (hours)") diff --git a/pkg/R/utils.R b/pkg/R/utils.R index a4efd61..bb76996 100644 --- a/pkg/R/utils.R +++ b/pkg/R/utils.R @@ -63,7 +63,8 @@ integerIndexToDate = function(index, data) #' @param days_in Optional set to intersect with results (NULL to discard) #' #' @export -getSimilarDaysIndices = function(index, data, limit, same_season, days_in=NULL) +getSimilarDaysIndices = function(index, data, limit, same_season, + days_in=NULL, operational=TRUE) { index = dateIndexToInteger(index, data) @@ -73,15 +74,30 @@ getSimilarDaysIndices = function(index, data, limit, same_season, days_in=NULL) day_ref = dt_ref$wday #1=monday, ..., 6=saturday, 0=sunday month_ref = as.POSIXlt(data$getTime(index)[1])$mon+1 #month in 1...12 i = index - 1 - while (i >= 1 && length(days) < limit) + if (!operational) + j = index + 1 + while (length(days) < min( limit, ifelse(is.null(days_in),Inf,length(days_in)) )) { - dt = as.POSIXlt(data$getTime(i)[1]) - if ((is.null(days_in) || i %in% days_in) && .isSameDay(dt$wday, day_ref)) + if (i >= 1) { - if (!same_season || .isSameSeason(dt$mon+1, month_ref)) - days = c(days, i) + dt = as.POSIXlt(data$getTime(i)[1]) + if ((is.null(days_in) || i %in% days_in) && .isSameDay(dt$wday, day_ref)) + { + if (!same_season || .isSameSeason(dt$mon+1, month_ref)) + days = c(days, i) + } + i = i - 1 + } + if (!operational && j <= data$getSize()) + { + dt = as.POSIXlt(data$getTime(j)[1]) + if ((is.null(days_in) || j %in% days_in) && .isSameDay(dt$wday, day_ref)) + { + if (!same_season || .isSameSeason(dt$mon+1, month_ref)) + days = c(days, j) + } + j = j + 1 } - i = i - 1 } return ( days ) } diff --git a/pkg/tests/testthat/test-Forecaster.R b/pkg/tests/testthat/test-Forecaster.R index 78e387a..3f5cf9c 100644 --- a/pkg/tests/testthat/test-Forecaster.R +++ b/pkg/tests/testthat/test-Forecaster.R @@ -97,12 +97,12 @@ test_that("Neighbors method behave as expected", { #Situation A: +Zero; correct if jump, wrong otherwise pred00 = computeForecast(data_p, indices, "Neighbors", "Zero", 1, Inf, 24, - simtype="mix", local=FALSE) + simtype="mix", local=FALSE, window=c(1,1)) for (i in 1:7) expect_equal(pred00$getForecast(i), rep(pred_order[i],24)) pred13 = computeForecast(data_p, indices, "Persistence", "Zero", 14, Inf, 24, - simtype="mix", local=FALSE) + simtype="mix", local=FALSE, window=c(1,1)) for (i in 1:7) expect_equal(pred13$getForecast(i), rep(i,24) ) diff --git a/pkg/tests/testthat/test-computeFilaments.R b/pkg/tests/testthat/test-computeFilaments.R index 6169a77..0c58c69 100644 --- a/pkg/tests/testthat/test-computeFilaments.R +++ b/pkg/tests/testthat/test-computeFilaments.R @@ -4,55 +4,51 @@ test_that("output is as expected on simulated series", { data = getDataTest(150) - - - -#TODO: debug - - - - # index 144-1 == 143 : serie type 2 - pred = computeForecast(data, 143, "Neighbors", "Zero", predict_from=8, - horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, h_window=1) + # index 144 : serie type 3, yersteday type 2 + pred = computeForecast(data, 144, "Neighbors", "Zero", predict_from=1, + horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, window=1, opera=TRUE) f = computeFilaments(data, pred, 1, 8, limit=60, plot=FALSE) - # Expected output: 50-3-10 series of type 2, then 23 series of type 3 (closest next) + # Expected output: 50-3-10 series of type 2+1 = 3, + # then 23 series of type 3+1 %% 3 = 1 (3 = closest next) expect_identical(length(f$neighb_indices), as.integer(60)) expect_identical(length(f$colors), as.integer(60)) - expect_equal(f$index, 143) - expect_true(all(I(f$neighb_indices) >= 2)) + expect_equal(f$index, 144) + expect_true(all(I(f$neighb_indices) != 2)) for (i in 1:37) { - expect_equal(I(f$neighb_indices[i]), 2) + expect_equal(I(f$neighb_indices[i]), 3) expect_match(f$colors[i], f$colors[1]) } for (i in 38:60) { - expect_equal(I(f$neighb_indices[i]), 3) + expect_equal(I(f$neighb_indices[i]), 1) expect_match(f$colors[i], f$colors[38]) } expect_match(f$colors[1], "#1*") expect_match(f$colors[38], "#E*") - # index 143-1 == 142 : serie type 1 - pred = computeForecast(data, 143, "Neighbors", "Zero", predict_from=8, - horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, h_window=1) + # index 143 : serie type 2 + pred = computeForecast(data, 143, "Neighbors", "Zero", predict_from=1, + horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, window=1, opera=TRUE) f = computeFilaments(data, pred, 1, 8, limit=50, plot=FALSE) - # Expected output: 50-10-3 series of type 1, then 13 series of type 3 (closest next) - # NOTE: -10 because only past days with no-NAs tomorrow => exclude type 1 in [60,90[ + # Expected output: 50-10-3 series of type 1+1=2, + # then 13 series of type 3+1 %% 3 = 1 (closest next) + # NOTE: -10 because only past tomorrows with no-NAs yerstedays + # => exclude type 2 in [60,90[ expect_identical(length(f$neighb_indices), as.integer(50)) expect_identical(length(f$colors), as.integer(50)) expect_equal(f$index, 143) - expect_true(all(I(f$neighb_indices) != 2)) + expect_true(all(I(f$neighb_indices) <= 2)) for (i in 1:37) { - expect_equal(I(f$neighb_indices[i]), 1) + expect_equal(I(f$neighb_indices[i]), 2) expect_match(f$colors[i], f$colors[1]) } for (i in 38:50) { - expect_equal(I(f$neighb_indices[i]), 3) + expect_equal(I(f$neighb_indices[i]), 1) expect_match(f$colors[i], f$colors[38]) } expect_match(f$colors[1], "#1*") -- 2.44.0