Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
221 lines (206 sloc) 7.18 KB
# Combination approach
# Author: Alan Lee; largely improved by Earo Wang
UpdateC <- function(c.list) {
k <- length(c.list)
div <- 1L
nvec <- numeric(k)
comb.vec <- NULL
for (i in 1:k) {
m <- c.list[[i]][[2L]]
cc <- c.list[[i]][[1L]]
nvec[i] <- dim(cc)[1L]
div <- div - sum(m * (cc %*% m))
comb.vec <- c(comb.vec, m)
}
d <- sum(comb.vec)
div <- div + d
sum.nvec <- sum(nvec)
c.star <- matrix(, sum.nvec, sum.nvec)
hi <- cumsum(nvec)
lo <- cumsum(c(1L, nvec[-length(nvec)]))
for (i in 1:k) {
cm1 <- as.vector(c.list[[i]][[1L]] %*% c.list[[i]][[2L]])
row.range <- lo[i]:hi[i]
for (j in 1:k) {
cm2 <- as.vector(c.list[[j]][[1L]] %*% c.list[[j]][[2]])
cinsert <- outer(1 - cm1, 1 - cm2)/div
if (i == j) {
cinsert <- c.list[[i]][[1L]] + cinsert
}
col.range <- lo[j]:hi[j]
c.star[row.range, col.range] <- cinsert
}
}
return(list(C = c.star, nvec = comb.vec))
}
CombineH <- function(fcasts, nodes) {
class(fcasts) <- "matrix" # drop "ts" object to process faster
# Split fcasts to a list
levels <- cumsum(Mnodes(nodes))
l.levels <- length(levels)
flist <- lapply(2L:l.levels, function(x) {
fcasts[, seq(levels[x - 1L] + 1L, levels[x]), drop = FALSE]
})
flist <- c(list(fcasts[, 1L, drop = FALSE]), flist)
rm(fcasts)
# Start with the last level
lenl <- length(levels)
lenn <- length(nodes)
last.nodes <- nodes[[lenn]]
last.len <- length(last.nodes)
cmat <- lapply(1:last.len, function(x)
list(matrix(1/(last.nodes[x] + 1)), last.nodes[x]))
idx <- c(0, cumsum(last.nodes))
smat <- lapply(1L:last.len, function(x)
flist[[lenl]][, (idx[x] + 1L):idx[x + 1L], drop = FALSE]
+ flist[[lenl - 1L]][, x])
if (lenn == 1L) { # A simple hierarchy with only 2 levels
cmat <- UpdateC(cmat[1L])$C
smat <- apply(smat[[1L]], 2, function(x) x + flist[[1L]])
sums <- rowsum(t(smat), rep(1L:last.len, last.nodes))
comb <- smat - rep(cmat %*% sums, last.nodes)
} else { # more than 2 levels
# Recursively update C matrix from L - 1 to 1
for (i in 1L:(lenn - 1L)) {
newn <- nodes[[lenn - i]]
newl <- length(newn)
new.cmat <- vector(length = newl, mode = "list")
new.smat <- vector(length = newl, mode = "list")
idx <- c(0L, cumsum(newn))
for (j in 1L:newl) {
new.cmat[[j]] <- UpdateC(cmat[(idx[j] + 1L):idx[j + 1L]])
sblock <- smat[(idx[j] + 1L):idx[j + 1L]]
sblock <- do.call("cbind", sblock)
new.smat[[j]] <- flist[[lenl - i - 1L]][, j] + sblock
}
cmat <- new.cmat
smat <- new.smat
}
cmat <- cmat[[1L]]$C
comb <- t(apply(smat[[1L]], 1, function(x) {
sums <- rowsum(x, rep(1L:last.len, last.nodes))
return(x - rep(cmat %*% sums, last.nodes))
}))
}
colnames(comb) <- NULL
return(comb)
}
# Combination with weights
UpdateCw <- function(c.list, d1.vec, d0) {
l.c <- length(c.list)
comb.vec <- NULL
nvec <- numeric(l.c)
div <- d0
for (i in 1L:l.c) {
m <- c.list[[i]][[2L]]
cmat <- c.list[[i]][[1L]]
d <- d1.vec[m]
nvec[i] <- length(m)
div <- div + sum(d) - sum(d * (cmat %*% d))
comb.vec <- c(comb.vec, m)
}
len.comb <- length(comb.vec)
c.star <- matrix(, nrow = len.comb, ncol = len.comb)
hi <- cumsum(nvec)
lo <- cumsum(c(1L, nvec[-length(nvec)]))
for (i in 1L:l.c) {
di <- d1.vec[c.list[[i]][[2L]]]
cd1 <- as.vector(c.list[[i]][[1L]] %*% di)
row.range <- lo[i]:hi[i]
for (j in 1L:l.c) {
col.range <- lo[j]:hi[j]
dj <- d1.vec[c.list[[j]][[2L]]]
cd2 <- as.vector(c.list[[j]][[1L]] %*% dj)
cinsert <- outer(1L - cd1, 1L - cd2)/div
if (i == j) {
cinsert <- c.list[[i]][[1L]] + cinsert
}
c.star[row.range, col.range] <- cinsert
}
}
return(list(cmat = c.star, m = comb.vec))
}
CombineHw <- function(fcasts, nodes, weights) {
class(fcasts) <- "matrix" # drop "ts" object to process faster
h <- nrow(fcasts)
# Split fcasts to a list
levels <- cumsum(Mnodes(nodes))
l.levels <- length(levels)
flist <- lapply(2L:l.levels, function(x) {
fcasts[, seq(levels[x - 1L] + 1L, levels[x]), drop = FALSE]
})
flist <- c(list(fcasts[, 1L, drop = FALSE]), flist)
rm(fcasts)
# Split weights to a list
wlist <- lapply(2L:l.levels, function(x) {
weights[seq(levels[x - 1L] + 1L, levels[x])]
})
wlist <- c(list(weights[1L]), wlist)
# Start with the last level
lenl <- length(levels)
lenn <- length(nodes)
last.nodes <- nodes[[lenn]]
last.len <- length(last.nodes)
lastg <- rep(1L:last.len, last.nodes)
dlist <- split(1/wlist[[l.levels]], lastg)
d1vec <- sapply(dlist, sum)
d0 <- 1/wlist[[l.levels - 1L]]
cmat <- lapply(1:last.len, function(x)
list(matrix(1L/(d0[x] + sum(dlist[[x]]))), x))
idx <- c(0, cumsum(last.nodes))
smat <- lapply(1L:last.len, function(x) {
yy <- sweep(flist[[lenl]][, (idx[x] + 1L):idx[x + 1L],
drop = FALSE], 2,
wlist[[lenl]][(idx[x] + 1L):idx[x + 1L]], "*")
tmp.yy <- sweep(flist[[lenl - 1L]][, x, drop = FALSE], 2,
wlist[[lenl - 1L]][x], "*")
return(apply(yy, 2, function(x) x + tmp.yy))
})
if (lenn == 1L) { # A simple hierarchy with only 2 levels
cmat <- UpdateCw(cmat[1L], d1vec, d0)$cmat
dvec <- unlist(dlist)
smat <- apply(smat[[1L]], 2, function(x) x + flist[[1L]] * wlist[[1L]])
sums <- rowsum(t(smat), rep(1L:last.len, last.nodes))
comb <- (smat - rep(cmat %*% sums, last.nodes)) * dvec
} else { # more than 2 levels
# Recursively update C matrix from L - 1 to 1
for (i in 1L:(lenn - 1L)) {
d0 <- 1/wlist[[lenl - i - 1L]]
newn <- nodes[[lenn - i]]
newl <- length(newn)
new.cmat <- vector(length = newl, mode = "list")
new.smat <- vector(length = newl, mode = "list")
idx <- c(0L, cumsum(newn))
for (j in 1L:newl) {
new.cmat[[j]] <- UpdateCw(cmat[(idx[j] + 1L):idx[j + 1L]], d1vec,
d0[j])
tmpw <- wlist[[lenl - i - 1L]][j]
tmps <- sweep(flist[[lenl - i - 1L]][, j, drop = FALSE], 2,
tmpw, "*")
sblock <- smat[(idx[j] + 1L):idx[j + 1L]]
if (h == 1) {
sblock <- unlist(sblock)
new.smat[[j]] <- sblock + tmps
} else {
sblock <- do.call("cbind", sblock)
new.smat[[j]] <- apply(sblock, 2, function(x) x + tmps)
}
}
cmat <- new.cmat
smat <- new.smat
}
cmat <- cmat[[1L]]$cmat
dvec <- unlist(dlist)
if (h == 1) {
sums <- rowsum(smat[[1L]] * dvec, rep(1L:last.len, last.nodes))
comb <- (smat[[1L]] - rep(cmat %*% sums, last.nodes)) * dvec
} else {
comb <- t(apply(smat[[1L]], 1, function(x) {
sums <- rowsum(x * dvec, rep(1L:last.len, last.nodes))
return((x - rep(cmat %*% sums, last.nodes)) * dvec)
}))
}
}
colnames(comb) <- NULL
return(comb)
}