various fixes (especially in plotFbox and Neighbors shape predict
[talweg.git] / pkg / R / F_Neighbors.R
CommitLineData
e030a6e3 1#' @include Forecaster.R
3d69ff21 2#'
e030a6e3 3#' @title Neighbors Forecaster
3d69ff21
BA
4#'
5#' @description Predict tomorrow as a weighted combination of "futures of the past" days.
e030a6e3
BA
6#' Inherits \code{\link{Forecaster}}
7NeighborsForecaster = setRefClass(
8 Class = "NeighborsForecaster",
9 contains = "Forecaster",
3d69ff21
BA
10
11 methods = list(
12 initialize = function(...)
13 {
14 callSuper(...)
15 },
e030a6e3 16 predictShape = function(today, memory, horizon, ...)
3d69ff21
BA
17 {
18 # (re)initialize computed parameters
19 params <<- list("weights"=NA, "indices"=NA, "window"=NA)
20
21 first_day = max(today - memory, 1)
22 # The first day is generally not complete:
23 if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2)))
24 first_day = 2
25
09cf9c19
BA
26 # Predict only on (almost) non-NAs days
27 nas_in_serie = is.na(data$getSerie(today))
28 if (any(nas_in_serie))
29 {
30 #TODO: better define "repairing" conditions (and method)
31 if (sum(nas_in_serie) >= length(nas_in_serie) / 2)
32 return (NA)
33 for (i in seq_along(nas_in_serie))
34 {
35 if (nas_in_serie[i])
36 {
37 #look left
38 left = i-1
39 while (left>=1 && nas_in_serie[left])
40 left = left-1
41 #look right
42 right = i+1
43 while (right<=length(nas_in_serie) && nas_in_serie[right])
44 right = right+1
45 #HACK: modify by-reference Data object...
46 data$data[[today]]$serie[i] <<-
47 if (left==0) data$data[[today]]$serie[right]
48 else if (right==0) data$data[[today]]$serie[left]
49 else (data$data[[today]]$serie[left] + data$data[[today]]$serie[right]) / 2.
50 }
51 }
52 }
3d69ff21
BA
53
54 # Determine indices of no-NAs days followed by no-NAs tomorrows
55 fdays_indices = c()
56 for (i in first_day:(today-1))
57 {
58 if ( !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) )
59 fdays_indices = c(fdays_indices, i)
60 }
61
62 #GET OPTIONAL PARAMS
63 # Similarity computed with exogenous variables ? endogenous ? both ? ("exo","endo","mix")
99f83c9a 64 simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix")
3d69ff21 65 simthresh = ifelse(hasArg("simthresh"), list(...)$simthresh, 0.)
99f83c9a
BA
66 kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan"
67 mix_strategy = ifelse(hasArg("mix_strategy"), list(...)$mix_strategy, "mult") #or "neighb"
68 same_season = ifelse(hasArg("same_season"), list(...)$same_season, FALSE)
3d69ff21 69 if (hasArg(h_window))
e030a6e3 70 return (.predictShapeAux(fdays_indices, today, horizon, list(...)$h_window, kernel,
99f83c9a 71 simtype, simthresh, mix_strategy, TRUE))
3d69ff21
BA
72 #END GET
73
74 # Indices for cross-validation; TODO: 45 = magic number
75 indices = getSimilarDaysIndices(today, limit=45, same_season=same_season)
99f83c9a
BA
76 if (tail(indices,1) == 1)
77 indices = head(indices,-1)
3d69ff21
BA
78
79 # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days
80 errorOnLastNdays = function(h, kernel, simtype)
81 {
82 error = 0
83 nb_jours = 0
84 for (i in indices)
85 {
86 # NOTE: predict only on non-NAs days followed by non-NAs (TODO:)
87 if (!any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))))
88 {
89 nb_jours = nb_jours + 1
90 # mix_strategy is never used here (simtype != "mix"), therefore left blank
e030a6e3
BA
91 prediction = .predictShapeAux(fdays_indices, i, horizon, h, kernel, simtype,
92 simthresh, "", FALSE)
3d69ff21
BA
93 if (!is.na(prediction[1]))
94 error = error + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2)
95 }
96 }
97 return (error / nb_jours)
98 }
99
100 h_best_exo = 1.
101 if (simtype != "endo" && !(simtype=="mix" && mix_strategy=="neighb"))
102 {
103 h_best_exo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel,
104 simtype="exo")$minimum
105 }
106 if (simtype != "exo")
107 {
108 h_best_endo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel,
109 simtype="endo")$minimum
110 }
111
112 if (simtype == "endo")
113 {
e030a6e3 114 return (.predictShapeAux(fdays_indices, today, horizon, h_best_endo, kernel, "endo",
3d69ff21
BA
115 simthresh, "", TRUE))
116 }
117 if (simtype == "exo")
118 {
e030a6e3 119 return (.predictShapeAux(fdays_indices, today, horizon, h_best_exo, kernel, "exo",
3d69ff21
BA
120 simthresh, "", TRUE))
121 }
122 if (simtype == "mix")
123 {
e030a6e3
BA
124 return (.predictShapeAux(fdays_indices, today, horizon, c(h_best_endo,h_best_exo),
125 kernel, "mix", simthresh, mix_strategy, TRUE))
3d69ff21
BA
126 }
127 },
128 # Precondition: "today" is full (no NAs)
e030a6e3 129 .predictShapeAux = function(fdays_indices, today, horizon, h, kernel, simtype, simthresh,
3d69ff21
BA
130 mix_strategy, final_call)
131 {
132 dat = data$data #HACK: faster this way...
133
134 fdays_indices = fdays_indices[fdays_indices < today]
135 # TODO: 3 = magic number
136 if (length(fdays_indices) < 3)
137 return (NA)
138
139 if (simtype != "exo")
140 {
141 h_endo = ifelse(simtype=="mix", h[1], h)
142
143 # Distances from last observed day to days in the past
144 distances2 = rep(NA, length(fdays_indices))
145 for (i in seq_along(fdays_indices))
146 {
147 delta = dat[[today]]$serie - dat[[ fdays_indices[i] ]]$serie
148 # Require at least half of non-NA common values to compute the distance
149 if (sum(is.na(delta)) <= 0) #length(delta)/2)
150 distances2[i] = mean(delta^2) #, na.rm=TRUE)
151 }
152
153 sd_dist = sd(distances2)
99f83c9a
BA
154 if (sd_dist < .Machine$double.eps)
155 sd_dist = 1 #mostly for tests... FIXME:
3d69ff21 156 simils_endo =
99f83c9a 157 if (kernel=="Gauss")
3d69ff21 158 exp(-distances2/(sd_dist*h_endo^2))
99f83c9a 159 else { #Epanechnikov
3d69ff21
BA
160 u = 1 - distances2/(sd_dist*h_endo^2)
161 u[abs(u)>1] = 0.
162 u
163 }
164 }
165
166 if (simtype != "endo")
167 {
168 h_exo = ifelse(simtype=="mix", h[2], h)
169
dea7ff86
BA
170 M = matrix( nrow=1+length(fdays_indices), ncol=1+length(dat[[today]]$exo) )
171 M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo) )
3d69ff21
BA
172 for (i in seq_along(fdays_indices))
173 {
174 M[i+1,] = c( dat[[ fdays_indices[i] ]]$level,
dea7ff86 175 as.double(dat[[ fdays_indices[i] ]]$exo) )
3d69ff21
BA
176 }
177
178 sigma = cov(M) #NOTE: robust covariance is way too slow
613a986f 179 sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed?
3d69ff21
BA
180
181 # Distances from last observed day to days in the past
182 distances2 = rep(NA, nrow(M)-1)
183 for (i in 2:nrow(M))
184 {
185 delta = M[1,] - M[i,]
186 distances2[i-1] = delta %*% sigma_inv %*% delta
187 }
188
189 sd_dist = sd(distances2)
190 simils_exo =
191 if (kernel=="Gauss") {
192 exp(-distances2/(sd_dist*h_exo^2))
193 } else { #Epanechnikov
194 u = 1 - distances2/(sd_dist*h_exo^2)
195 u[abs(u)>1] = 0.
196 u
197 }
198 }
199
200 if (simtype=="mix")
201 {
202 if (mix_strategy == "neighb")
203 {
204 #Only (60) most similar days according to exogen variables are kept into consideration
205 #TODO: 60 = magic number
206 keep_indices = sort(simils_exo, index.return=TRUE)$ix[1:(min(60,length(simils_exo)))]
207 simils_endo[-keep_indices] = 0.
3d69ff21 208 }
99f83c9a
BA
209 else #mix_strategy == "mult"
210 simils_endo = simils_endo * simils_exo
3d69ff21
BA
211 }
212
213 similarities =
214 if (simtype != "exo") {
215 simils_endo
216 } else {
217 simils_exo
218 }
219
220 if (simthresh > 0.)
221 {
222 max_sim = max(similarities)
223 # Set to 0 all similarities s where s / max_sim < simthresh, but keep at least 60
224 ordering = sort(similarities / max_sim, index.return=TRUE)
225 if (ordering[60] < simthresh)
226 {
227 similarities[ ordering$ix[ - (1:60) ] ] = 0.
228 } else
229 {
230 limit = 61
231 while (limit < length(similarities) && ordering[limit] >= simthresh)
232 limit = limit + 1
233 similarities[ ordering$ix[ - 1:limit] ] = 0.
234 }
235 }
236
237 prediction = rep(0, horizon)
238 for (i in seq_along(fdays_indices))
239 prediction = prediction + similarities[i] * dat[[ fdays_indices[i]+1 ]]$serie[1:horizon]
3d69ff21 240 prediction = prediction / sum(similarities, na.rm=TRUE)
99f83c9a 241
3d69ff21
BA
242 if (final_call)
243 {
244 params$weights <<- similarities
245 params$indices <<- fdays_indices
246 params$window <<-
247 if (simtype=="endo") {
248 h_endo
249 } else if (simtype=="exo") {
250 h_exo
251 } else {
252 c(h_endo,h_exo)
253 }
254 }
99f83c9a 255
3d69ff21
BA
256 return (prediction)
257 }
258 )
259)