Permalink
Browse files

Rewrite vforecast

  • Loading branch information...
1 parent 719170f commit 5f7a42546caec9681910b2a9b4040809e5e33fdd @eodus committed Apr 27, 2012
Showing with 53 additions and 39 deletions.
  1. +53 −39 R/forecast.R
View
92 R/forecast.R
@@ -132,31 +132,34 @@ apply.lrf <- function(F, lrf, len = 1) {
}
-"vforecast.1d-ssa" <- function(this, groups, len = 1,
- ...) {
- L <- this$window
- K <- this$length - L + 1;
- L.s <- min(L, K);
- N <- K + L - 1 + len + L - 1;
- N.res <- K + L - 1 + len;
-
- dv <- c(1:(L.s-1), rep(L.s, N-2*L.s+2), (L.s-1):1);
+"vforecast.1d-ssa" <- function(this, groups, len = 1,
+ ...,
+ cache = TRUE) {
+ #TODO Implement convolution via FFTW in C-code
convolve.open <- function(F, G) {
NF <- length(F)
NG <- length(G)
-
+
NN <- nextn(NF+NG-1)
ZFZ <- c(rep(0, NG-1), F, rep(0, NN-(NF+NG-1)));
GZ <- c(G, rep(0, NN-NG));
-
+
res <- fft(fft(ZFZ)*Conj(fft(GZ)), inverse = TRUE)/NN;
-
+
Re(res)[1:(NF+NG-1)];
}
+ L <- this$window;
+ K <- this$length - L + 1;
+ L.s <- min(L, K);
+ N <- K + L - 1 + len + L - 1;
+ N.res <- K + L - 1 + len;
+
+ dv <- c(1:(L.s-1), rep(L.s, N-2*L.s+2), (L.s-1):1);
+
if (missing(groups))
- groups <- as.list(1:min(nlambda(this), nu(this)))
+ groups <- as.list(1:min(nlambda(this), nu(this)));
# Determine the upper bound of desired eigentriples
desired <- max(unlist(groups));
@@ -168,38 +171,49 @@ apply.lrf <- function(F, lrf, len = 1) {
lambda <- .get(this, "lambda");
U <- .get(this, "U");
- V <- if (nv(this) >= desired) .get(this, "V") else NULL
-
- out <- list()
+ if (nv(this) >= desired) {
+ V <- .get(this, "V");
+ } else {
+ if (cache){
+ # FIXME We can evaluate less FV here
+ V <- calc.v(this, idx = 1:desired);
+ .set(this, "V", V);
+ } else {
+ V <- NULL;
+ }
+ }
+
+ out <- list();
for (i in seq_along(groups)) {
- ET <- unique(groups[[i]])
-
- Vl <- matrix(NA, N.res, length(ET));
- Uet <- U[ , ET, drop=FALSE];
-
- Vl[1:K, ] <- (if (is.null(V))
- calc.v(this, idx = ET)
- else V[ , ET]) %*% diag(lambda[ET], nrow = length(ET))
-
- U.head <- Uet[-L, , drop=FALSE];
- U.tail <- Uet[-1, , drop=FALSE];
+ group <- unique(groups[[i]]);
+
+ Uet <- U[, group, drop = FALSE];
+ Vet <- if (is.null(V)) calc.v(this, idx = group) else V[, group, drop = FALSE];
+
+ Z <- rbind(t(lambda[group] * t(Vet)), matrix(NA, len + L - 1, length(group)));
+
+ U.head <- Uet[-L, , drop = FALSE];
+ U.tail <- Uet[-1, , drop = FALSE];
+
+ #TODO Get rid of solve
P <- solve(t(U.head) %*% U.head, t(U.head) %*% U.tail);
-
- for (j in (K+1):(K+len+L-1)) {
- Vl[j, ] <- P %*% Vl[j-1, ];
+
+ for(j in (K + 1):(K + len + L - 1)) {
+ Z[j, ] <- P %*% Z[j - 1, ];
}
-
- res <- rep(0, N);
- for (j in 1:length(ET)) {
- res <- res + convolve.open(Vl[ , j], rev(Uet[ , j]));
+
+ res <- double(N);
+ for(j in seq_along(group)) {
+ res <- res + convolve.open(Z[ , j], rev(Uet[ , j]));
}
-
+
out[[i]] <- (res/dv)[1:N.res];
- };
-
+ }
+
names(out) <- paste(sep = "", "F", 1:length(groups));
-
- out
+
+ # Forecasted series can be pretty huge...
+ invisible(out);
}
rforecast.ssa <- function(x, groups, len = 1,

0 comments on commit 5f7a425

Please sign in to comment.