###################################################
### chunk number 8: periodicLags012
###################################################
## Periodic Seasonal by Region (all lags at once)
rm(list = ls())

source(file.path("data", "utils.R"))

exclude <- c("anch", "hono")

## Load the results
load(file.path("results", "seasonal.periodic.lag.012.dfSeas.1.rda"))
results <- lapply(results, function(r) r[setdiff(names(r), exclude)])

## Get region indicators for each city
data(cities)
cityRegion <- with(cities, region[match(names(results[[1]]), city)])
regionNames <- sort(unique(cityRegion))
regionInd <- 1 * outer(cityRegion, regionNames, "==")
colnames(regionInd) <- regionNames

ndays <- 365
x <- 1:ndays
B <- harmonic(x, 1, ndays, intercept = TRUE)

bc <- lapply(results, extractBetaCov, pollutant = "pm10")
set.seed(100302)
gg <- lapply(bc, function(x) {
    initTLNise()
    tlnise(x$beta, x$cov, regionInd, prnt = FALSE, intercept = FALSE,
           seed = sample(100000, 1))
})
pooledRegion <- lapply(gg, function(g) {
    matrix(g$gamma[,1], byrow = TRUE, ncol = ncol(regionInd),
           nrow = nrow(g$theta))
})
pooledRegionSD <- lapply(gg, function(g) {
    matrix(g$gamma[,2], byrow = TRUE, ncol = ncol(regionInd),
           nrow = nrow(g$theta))
})
b <- lapply(pooledRegion, function(p) B %*% p)
zero <- matrix(0, nrow = ndays, ncol = 7 * 3)
curveRegion <- lapply(seq(along = pooledRegion), function(j) {
    pr <- pooledRegion[[j]]
    lapply(1:ncol(pr), function(i) {
        B <- zero
        ## This is weird
        B[, seq(i, 7 * 3, 7)] <- harmonic(x, 1, ndays, intercept = TRUE)
        s <- sqrt(diag(B %*% gg[[j]]$Dgamma %*% t(B)))
        cbind(b[[j]][,i], b[[j]][,i] - 2*s, b[[j]][,i] + 2*s)
    })
})

## Overall estimate
set.seed(10321342)
gg <- lapply(bc, function(b) {
    initTLNise()
    tlnise(b$beta, b$cov, rep(1, NROW(b$beta)), seed = sample(100000, 1),
           prnt = FALSE)
})
V <- lapply(gg, function(g) diag(B %*% g$Dgamma %*% t(B)))
b <- lapply(gg, function(g) B %*% g$gamma[,1])
curve <- lapply(seq(along = b), function(i) {
    cbind(b[[i]], b[[i]] - 2*sqrt(V[[i]]), b[[i]] + 2*sqrt(V[[i]]))
})

curveTotal <- curveRegion
for(i in seq(along = curveTotal)) {
    curveTotal[[i]] <- append(curveTotal[[i]], list(curve[[i]]))
}

regionFullNames <- c("Industrial Midwest", "Northeast", "Northwest",
                     "Southern California", "Southeast", "Southwest",
                     "Upper Midwest", "All Regions")

y0 <- do.call("rbind", curveTotal[[1]])[,1] * 1000
y1 <- do.call("rbind", curveTotal[[2]])[,1] * 1000
conf1 <- do.call("rbind", curveTotal[[2]])[, 2:3] * 1000
conf1 <- rbind(conf1, conf1, conf1)
y2 <- do.call("rbind", curveTotal[[3]])[,1] * 1000
rng <- range(unlist(curveTotal)) * 1000 + c(-1,1) * .07
xx <- rep(x, length(regionFullNames))
f <- ordered(rep(regionFullNames, each = 365), levels = regionFullNames)

## trellis.device(dev = x11, height = 5, width = 8, color = FALSE)

l <- trellis.par.get("superpose.line")
l$lty <- c(2, 1, 4)
l$lwd <- c(1, 1.5, 1)
trellis.par.set("superpose.line", l)

p <- xyplot(y0 + y1 + y2 ~ xx | f, allow.multiple = TRUE, type = "l",
            ylim = rng, layout = c(4, 2), as.table = TRUE,
            panel = function(x, y, subscripts, groups, ...) {
                s <- subscripts[1:365]
                xx <- x[1:365]
                llines(xx, conf1[s, 1], lty = 1, lwd = 1.5)
                llines(xx, conf1[s, 2], lty = 1, lwd = 1.5)
                panel.superpose(x, y, subscripts, groups, ...)
                panel.abline(h = 0, lty = 5)
            }, 
            xlab = "Day in year",
            strip = strip.custom(bg = 0),
            scales = list(alternating = 1, tck = c(0.8, 0), x = list(cex=.8)),
            ylab = expression(paste("% increase in mortality with ",
                "10-", mu, "g/", m^3, " increase in ", PM[10]))
            )
print(p)