on the way to R6 class + remove truncated days (simplifications)

diff --git a/data/scripts/augment_meteo.R b/data/scripts/augment_meteo.R
index 11649f3..c762fe8 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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,