From f17665c7d3da672163779da686d9f4d1ebad31f9 Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Fri, 24 Feb 2017 21:23:27 +0100 Subject: [PATCH] on the way to R6 class + remove truncated days (simplifications) --- data/scripts/augment_meteo.R | 5 +- pkg/R/Data.R | 32 ++--- pkg/R/F_Neighbors.R | 147 ++++++--------------- pkg/R/getData.R | 23 ++-- reports/report_2017-03-01.7h_average.ipynb | 2 +- 5 files changed, 72 insertions(+), 137 deletions(-) diff --git a/data/scripts/augment_meteo.R b/data/scripts/augment_meteo.R index 11649f3..c762fe8 100644 --- a/data/scripts/augment_meteo.R +++ b/data/scripts/augment_meteo.R @@ -8,7 +8,10 @@ meteo_df$Week = 0 meteo_df$Pollution = -1 #Need to load and aggregate PM10 by days: use getData() from package -data = getData(..., predict_at=0) #TODO: +ts_data = system.file("extdata","pm10_mesures_H_loc.csv",package="talweg") +exo_data = system.file("extdata","meteo_extra_noNAs.csv",package="talweg") +data = getData(ts_data, exo_data, input_tz = "Europe/Paris", + working_tz="Europe/Paris", predict_at=0) for (i in 1:nrow(meteo_df)) { diff --git a/pkg/R/Data.R b/pkg/R/Data.R index d4609f2..4e16805 100644 --- a/pkg/R/Data.R +++ b/pkg/R/Data.R @@ -1,30 +1,27 @@ -#' @title Data +#' Data #' -#' @description Data encapsulation +#' Data encapsulation #' #' @field data List of #' \itemize{ #' \item time: vector of times #' \item serie: centered series #' \item level: corresponding levels +#' \item exo: exogenous variables #' \item exo_hat: predicted exogenous variables -#' \item exo_Dm1: List of measured exogenous variables at day minus 1 #' } #' -#' @exportClass Data -#' @export Data -Data = setRefClass( - Class = "Data", - - fields = list( - data = "list" - ), - - methods = list( +#' @docType class +#' @importFrom R6 R6Class +#' +#' @export +Data = R6Class("Data", + public = list( + data = "list", initialize = function(...) { "Initialize empty Data object" - +#TODO: continue from here callSuper(...) }, getSize = function() @@ -37,12 +34,7 @@ Data = setRefClass( { "'Standard' horizon, from t+1 to midnight" - L1 = length(data[[1]]$serie) - L2 = length(data[[2]]$serie) - if (L1 < L2) - L2 - L1 - else - L1 + 24 - as.POSIXlt( data[[1]]$time[1] )$hour + 1 }, append = function(new_time, new_serie, new_level, new_exo_hat, new_exo) { diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index ffb6d37..43a6a13 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -18,16 +18,16 @@ NeighborsForecaster = setRefClass( # (re)initialize computed parameters params <<- list("weights"=NA, "indices"=NA, "window"=NA) - first_day = max(today - memory, 1) - # The first day is generally not complete: - if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2))) - first_day = 2 + # Get optional args + simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix") #or "endo", or "exo" + kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan" + if (hasArg(h_window)) + return (.predictShapeAux(fdays,today,horizon,list(...)$h_window,kernel,simtype,TRUE)) - # Predict only on (almost) non-NAs days + # HACK for test reports: complete some days with a few NAs, for nicer graphics nas_in_serie = is.na(data$getSerie(today)) if (any(nas_in_serie)) { - #TODO: better define "repairing" conditions (and method) if (sum(nas_in_serie) >= length(nas_in_serie) / 2) return (NA) for (i in seq_along(nas_in_serie)) @@ -52,88 +52,55 @@ NeighborsForecaster = setRefClass( } # Determine indices of no-NAs days followed by no-NAs tomorrows - fdays_indices = c() - for (i in first_day:(today-1)) - { - if ( !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) ) - fdays_indices = c(fdays_indices, i) - } - - #GET OPTIONAL PARAMS - # Similarity computed with exogenous variables ? endogenous ? both ? ("exo","endo","mix") - simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix") - simthresh = ifelse(hasArg("simthresh"), list(...)$simthresh, 0.) - kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan" - mix_strategy = ifelse(hasArg("mix_strategy"), list(...)$mix_strategy, "mult") #or "neighb" - same_season = ifelse(hasArg("same_season"), list(...)$same_season, FALSE) - if (hasArg(h_window)) - return (.predictShapeAux(fdays_indices, today, horizon, list(...)$h_window, kernel, - simtype, simthresh, mix_strategy, TRUE)) - #END GET + first_day = max(today - memory, 1) + fdays = (first_day:(today-1))[ sapply(first_day:(today-1), function(i) { + !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) + }) ] - # Indices for cross-validation; TODO: 45 = magic number - indices = getSimilarDaysIndices(today, limit=45, same_season=same_season) - if (tail(indices,1) == 1) - indices = head(indices,-1) + # Indices of similar days for cross-validation; TODO: 45 = magic number + sdays = getSimilarDaysIndices(today, limit=45, same_season=FALSE) # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days errorOnLastNdays = function(h, kernel, simtype) { error = 0 nb_jours = 0 - for (i in indices) + for (i in intersect(fdays,sdays)) { - # NOTE: predict only on non-NAs days followed by non-NAs (TODO:) - if (!any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1)))) + # mix_strategy is never used here (simtype != "mix"), therefore left blank + prediction = .predictShapeAux(fdays, i, horizon, h, kernel, simtype, FALSE) + if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 - # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = .predictShapeAux(fdays_indices, i, horizon, h, kernel, simtype, - simthresh, "", FALSE) - if (!is.na(prediction[1])) - error = error + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2) + error = error + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2) } } return (error / nb_jours) } - h_best_exo = 1. - if (simtype != "endo" && !(simtype=="mix" && mix_strategy=="neighb")) - { - h_best_exo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel, - simtype="exo")$minimum - } + if (simtype != "endo") + h_best_exo = optimize(errorOnLastNdays, c(0,10), kernel=kernel, simtype="exo")$minimum if (simtype != "exo") - { - h_best_endo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel, - simtype="endo")$minimum - } + h_best_endo = optimize(errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum if (simtype == "endo") - { - return (.predictShapeAux(fdays_indices, today, horizon, h_best_endo, kernel, "endo", - simthresh, "", TRUE)) - } + return (.predictShapeAux(fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) if (simtype == "exo") - { - return (.predictShapeAux(fdays_indices, today, horizon, h_best_exo, kernel, "exo", - simthresh, "", TRUE)) - } + return (.predictShapeAux(fdays, today, horizon, h_best_exo, kernel, "exo", TRUE)) if (simtype == "mix") { - return (.predictShapeAux(fdays_indices, today, horizon, c(h_best_endo,h_best_exo), - kernel, "mix", simthresh, mix_strategy, TRUE)) + h_best_mix = c(h_best_endo,h_best_exo) + return (.predictShapeAux(fdays, today, horizon, h_best_mix, kernel, "mix", TRUE)) } }, # Precondition: "today" is full (no NAs) - .predictShapeAux = function(fdays_indices, today, horizon, h, kernel, simtype, simthresh, - mix_strategy, final_call) + .predictShapeAux = function(fdays, today, horizon, h, kernel, simtype, final_call) { dat = data$data #HACK: faster this way... - fdays_indices = fdays_indices[fdays_indices < today] + fdays = fdays[ fdays < today ] # TODO: 3 = magic number - if (length(fdays_indices) < 3) + if (length(fdays) < 3) return (NA) if (simtype != "exo") @@ -141,10 +108,10 @@ NeighborsForecaster = setRefClass( h_endo = ifelse(simtype=="mix", h[1], h) # Distances from last observed day to days in the past - distances2 = rep(NA, length(fdays_indices)) - for (i in seq_along(fdays_indices)) + distances2 = rep(NA, length(fdays)) + for (i in seq_along(fdays)) { - delta = dat[[today]]$serie - dat[[ fdays_indices[i] ]]$serie + delta = dat[[today]]$serie - dat[[ fdays[i] ]]$serie # Require at least half of non-NA common values to compute the distance if (sum(is.na(delta)) <= 0) #length(delta)/2) distances2[i] = mean(delta^2) #, na.rm=TRUE) @@ -167,13 +134,10 @@ NeighborsForecaster = setRefClass( { h_exo = ifelse(simtype=="mix", h[2], h) - M = matrix( nrow=1+length(fdays_indices), ncol=1+length(dat[[today]]$exo) ) + M = matrix( nrow=1+length(fdays), ncol=1+length(dat[[today]]$exo) ) M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo) ) - for (i in seq_along(fdays_indices)) - { - M[i+1,] = c( dat[[ fdays_indices[i] ]]$level, - as.double(dat[[ fdays_indices[i] ]]$exo) ) - } + for (i in seq_along(fdays)) + M[i+1,] = c( dat[[ fdays[i] ]]$level, as.double(dat[[ fdays[i] ]]$exo) ) sigma = cov(M) #NOTE: robust covariance is way too slow sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed? @@ -188,51 +152,22 @@ NeighborsForecaster = setRefClass( sd_dist = sd(distances2) simils_exo = - if (kernel=="Gauss") { + if (kernel=="Gauss") exp(-distances2/(sd_dist*h_exo^2)) - } else { #Epanechnikov + else { #Epanechnikov u = 1 - distances2/(sd_dist*h_exo^2) u[abs(u)>1] = 0. u } } - if (simtype=="mix") - { - if (mix_strategy == "neighb") - { - #Only (60) most similar days according to exogen variables are kept into consideration - #TODO: 60 = magic number - keep_indices = sort(simils_exo, index.return=TRUE)$ix[1:(min(60,length(simils_exo)))] - simils_endo[-keep_indices] = 0. - } - else #mix_strategy == "mult" - simils_endo = simils_endo * simils_exo - } - similarities = - if (simtype != "exo") { - simils_endo - } else { + if (simtype == "exo") simils_exo - } - - if (simthresh > 0.) - { - max_sim = max(similarities) - # Set to 0 all similarities s where s / max_sim < simthresh, but keep at least 60 - ordering = sort(similarities / max_sim, index.return=TRUE) - if (ordering[60] < simthresh) - { - similarities[ ordering$ix[ - (1:60) ] ] = 0. - } else - { - limit = 61 - while (limit < length(similarities) && ordering[limit] >= simthresh) - limit = limit + 1 - similarities[ ordering$ix[ - 1:limit] ] = 0. - } - } + else if (simtype == "endo") + simils_endo + else #mix + simils_endo * simils_exo prediction = rep(0, horizon) for (i in seq_along(fdays_indices)) @@ -248,7 +183,7 @@ NeighborsForecaster = setRefClass( h_endo } else if (simtype=="exo") { h_exo - } else { + } else { #mix c(h_endo,h_exo) } } diff --git a/pkg/R/getData.R b/pkg/R/getData.R index 8d1a6fa..da4b459 100644 --- a/pkg/R/getData.R +++ b/pkg/R/getData.R @@ -14,6 +14,7 @@ #' see \code{strptime}) #' @param working_tz Timezone to work with ("GMT" or e.g. "Europe/Paris") #' @param predict_at When does the prediction take place ? Integer, in hours. Default: 0 +#' @param limit Number of days to extract (default: Inf, for "all") #' #' @return An object of class Data #' @@ -61,7 +62,7 @@ getData = function(ts_data, exo_data, input_tz="GMT", date_format="%d/%m/%Y %H:% line = 1 #index in PM10 file (24 lines for 1 cell) nb_lines = nrow(ts_df) nb_exos = ( ncol(exo_df) - 1 ) / 2 - data = list() #new("Data") + data = Data$new() i = 1 #index of a cell in data while (line <= nb_lines) { @@ -78,18 +79,22 @@ getData = function(ts_data, exo_data, input_tz="GMT", date_format="%d/%m/%Y %H:% break } - # NOTE: if predict_at does not cut days at midnight, exogenous vars need to be shifted - exo_hat = as.data.frame( exo_df[ - ifelse(predict_at>0,max(1,i-1),i) , (1+nb_exos+1):(1+2*nb_exos) ] ) - exo = as.data.frame( exo_df[ ifelse(predict_at>0,max(1,i-1),i) , 2:(1+nb_exos) ] ) + exo = as.data.frame( exo_df[i,2:(1+nb_exos)] ) + exo_hat = as.data.frame( exo_df[i,(1+nb_exos+1):(1+2*nb_exos)] ) level = mean(serie, na.rm=TRUE) centered_serie = serie - level - #data$append(time, centered_serie, level, exo_hat, exo_Jm1) #too slow; TODO: use R6 class - data[[length(data)+1]] = list("time"=time, "serie"=centered_serie, "level"=level, - "exo_hat"=exo_hat, "exo"=exo) + data$append(time, centered_serie, level, exo, exo_hat) if (i >= limit) break i = i + 1 } - new("Data",data=data) + if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(1))) + data$removeFirst() + if (length(data$getCenteredSerie( data$getSize() )) < + length(data$getCenteredSerie( data$getSize()-1 ))) + { + data$removeLast() + } + + data } diff --git a/reports/report_2017-03-01.7h_average.ipynb b/reports/report_2017-03-01.7h_average.ipynb index 1776673..795307c 100644 --- a/reports/report_2017-03-01.7h_average.ipynb +++ b/reports/report_2017-03-01.7h_average.ipynb @@ -591,7 +591,7 @@ "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", - "version": "3.3.2" + "version": "3.2.3" } }, "nbformat": 4, -- 2.44.0