From f17665c7d3da672163779da686d9f4d1ebad31f9 Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
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