X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors.R;h=51b57304347c66faeb90528fa15428a5ec3b726b;hb=689aa1d3c58e8b1e9fdfa9b895ca38c7228f56cc;hp=12595d998f82d4883069f4c02a32d458b21aaa70;hpb=102bcfda4afbb5cfee885cbee0f55545624168fd;p=talweg.git

diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R
index 12595d9..51b5730 100644
--- a/pkg/R/F_Neighbors.R
+++ b/pkg/R/F_Neighbors.R
@@ -28,9 +28,11 @@
 #'     obtain the final prediction.
 #' }
 #'
+#' @usage f <- NeighborsForecaster$new(pjump)
+#'
 #' @docType class
 #' @format R6 class, inherits Forecaster
-#' @alias F_Neighbors
+#' @aliases F_Neighbors
 #'
 NeighborsForecaster = R6::R6Class("NeighborsForecaster",
 	inherit = Forecaster,
@@ -46,7 +48,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster",
 				return (NA)
 
 			# Determine indices of no-NAs days followed by no-NAs tomorrows
-			fdays = getNoNA2(data, max(today-memory,1), today-1)
+			fdays = .getNoNA2(data, max(today-memory,1), today-1)
 
 			# Get optional args
 			local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season?
@@ -118,33 +120,14 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster",
 
 			if (local)
 			{
-				# Neighbors: days in "same season"; TODO: 60 == magic number...
+				# TODO: 60 == magic number
 				fdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE,
 					days_in=fdays_cut)
 				if (length(fdays) <= 1)
 					return (NA)
-				levelToday = data$getLevel(today)
-				distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday))
-				#TODO: 2, 10, 3, 12 magic numbers here...
-				dist_thresh = 2
-				min_neighbs = min(10,length(fdays))
-				repeat
-				{
-					same_pollution = (distances <= dist_thresh)
-					nb_neighbs = sum(same_pollution)
-					if (nb_neighbs >= min_neighbs) #will eventually happen
-						break
-					dist_thresh = dist_thresh + 3
-				}
-				fdays = fdays[same_pollution]
-				max_neighbs = 12
-				if (nb_neighbs > max_neighbs)
-				{
-					# Keep only max_neighbs closest neighbors
-					fdays = fdays[
-						sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ]
-				}
-				if (length(fdays) == 1) #the other extreme...
+				# TODO: 10, 12 == magic numbers
+				fdays = .getConstrainedNeighbs(today,data,fdays,min_neighbs=10,max_neighbs=12)
+				if (length(fdays) == 1)
 				{
 					if (final_call)
 					{
@@ -170,13 +153,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster",
 					mean(delta^2)
 				})
 
-				sd_dist = sd(distances2)
-				if (sd_dist < .25 * sqrt(.Machine$double.eps))
-				{
-#					warning("All computed distances are very close: stdev too small")
-					sd_dist = 1 #mostly for tests... FIXME:
-				}
-				simils_endo = exp(-distances2/(sd_dist*window_endo^2))
+				simils_endo <- .computeSimils(distances2, window_endo)
 			}
 
 			if (simtype == "exo" || simtype == "mix")
@@ -203,13 +180,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster",
 					delta %*% sigma_inv %*% delta
 				})
 
-				sd_dist = sd(distances2)
-				if (sd_dist < .25 * sqrt(.Machine$double.eps))
-				{
-#					warning("All computed distances are very close: stdev too small")
-					sd_dist = 1 #mostly for tests... FIXME:
-				}
-				simils_exo = exp(-distances2/(sd_dist*window_exo^2))
+				simils_exo <- .computeSimils(distances2, window_exo)
 			}
 
 			similarities =
@@ -246,3 +217,57 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster",
 		}
 	)
 )
+
+# getConstrainedNeighbs
+#
+# Get indices of neighbors of similar pollution level (among same season + day type).
+#
+# @param today Index of current day
+# @param data Object of class Data
+# @param fdays Current set of "first days" (no-NA pairs)
+# @param min_neighbs Minimum number of points in a neighborhood
+# @param max_neighbs Maximum number of points in a neighborhood
+#
+.getConstrainedNeighbs = function(today, data, fdays, min_neighbs=10, max_neighbs=12)
+{
+	levelToday = data$getLevel(today)
+	distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday))
+	#TODO: 2, +3 : magic numbers
+	dist_thresh = 2
+	min_neighbs = min(min_neighbs,length(fdays))
+	repeat
+	{
+		same_pollution = (distances <= dist_thresh)
+		nb_neighbs = sum(same_pollution)
+		if (nb_neighbs >= min_neighbs) #will eventually happen
+			break
+		dist_thresh = dist_thresh + 3
+	}
+	fdays = fdays[same_pollution]
+	max_neighbs = 12
+	if (nb_neighbs > max_neighbs)
+	{
+		# Keep only max_neighbs closest neighbors
+		fdays = fdays[
+			sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ]
+	}
+	fdays
+}
+
+# compute similarities
+#
+# Apply the gaussian kernel on computed squared distances.
+#
+# @param distances2 Squared distances
+# @param window Window parameter for the kernel
+#
+.computeSimils <- function(distances2, window)
+{
+	sd_dist = sd(distances2)
+	if (sd_dist < .25 * sqrt(.Machine$double.eps))
+	{
+#		warning("All computed distances are very close: stdev too small")
+		sd_dist = 1 #mostly for tests... FIXME:
+	}
+	exp(-distances2/(sd_dist*window^2))
+}