+
+# 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 tdays Current set of "second 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, tdays, min_neighbs=10, max_neighbs=12)
+{
+ levelToday = data$getLevelHat(today)
+ levelYersteday = data$getLevel(today-1)
+ distances = sapply(tdays, function(i) {
+ sqrt((data$getLevel(i-1)-levelYersteday)^2 + (data$getLevel(i)-levelToday)^2)
+ })
+ #TODO: 1, +1, +3 : magic numbers
+ dist_thresh = 1
+ min_neighbs = min(min_neighbs,length(tdays))
+ repeat
+ {
+ same_pollution = (distances <= dist_thresh)
+ nb_neighbs = sum(same_pollution)
+ if (nb_neighbs >= min_neighbs) #will eventually happen
+ break
+ 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] ]
+ }
+ tdays
+}
+
+# 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))
+}