###############################################################################
## Fit seasonal models
## Copyright (C) 2004, Roger D. Peng <rpeng@jhsph.edu>
##     
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
###############################################################################


######################################################################
## Compute city specific estimates without seasonality
library(NMMAPSdata)
library(splines)
library(tsModelSpec)

rm(list = ls(all = TRUE))

registerDB(NULL)

cc <- dget(file.path("data", "cityList.R"))
extractors <- list(coefficients = function(x) summary(x)$coefficients,
                   cov = vcov)
results <- vector("list", length = length(cc))

poll <- "pm10"
pollutants <- paste(c("", "l1", "l2"), paste(poll, "tmean", sep=""), sep="")
results <- vector("list", length = length(pollutants))
names(results) <- paste("lag", 0:2, sep = "")

source(file.path("data", "model-fitting.R"))

for(i in seq(along = pollutants)) {
    pollutant <- pollutants[i]
    lag.results <- vector("list", length = length(cc))
    names(lag.results) <- cc
    
    for(l in seq(along = cc)) {
        cityname <- cc[l]
        cat(cityname, "\n", sep = " ")
        citydata <- readCity(cityname)
        lag.results[[l]] <- try({
            fitCitySeason(data = citydata, pollutant = pollutant, cause = "death",
                          season = "none", extractors = extractors)
        })
    }
    results[[i]] <- lag.results
}
outfile <- paste("city-specific-est", poll, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)


######################################################################
## Compute seasonally varying estimates using orthogonal sine/cosine basis

library(NMMAPSdata)
library(splines)
library(tsModelSpec)

rm(list = ls(all = TRUE))

registerDB(NULL)

cityList <- local({
    d <- dget(file.path("data", "cityList.R"))
    d[-match(c("anch", "hono"), d)]
})
extractors <- list(
                   coefficients = function(x) { summary(x)$coefficients },
                   cov = vcov
                   )
cause <- "death"
method <- "periodic2"
poll <- "pm10"
lags <- c(0, 1, 2)
results <- vector("list", length = length(lags))
names(results) <- paste("Lag", lags, sep = "")

source(file.path("data", "model-fitting.R"))

for(i in seq(along = lags)) {
    lag <- lags[i]
    cat(lag, " ")
    pollutant <- paste("Lag(", poll, "tmean, ", lag, ")", sep = "")

    lagresults <- cityApply(fitCitySeason, cityList = cityList, verbose = TRUE,
                            hookFun = extractors, pollutant = pollutant, 
                            cause = cause, season = method)
    names(lagresults) <- cityList
    results[[i]] <- lagresults
}

outfile <- paste("seasonal", "periodic-ortho", "lag", paste(lags, collapse = ""),
                 poll, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)

####################################################################
## Sensitivity analysis wrt degrees of freedom in smooth function of time

library(NMMAPSdata)
library(splines)
library(tsModelSpec)

rm(list = ls(all = TRUE))

cc <- dget(file.path("data", "cityList.R"))
extractors <- list(
                   coefficients = function(x) { summary(x)$coefficients },
                   cov = vcov
                   )
cause <- "death"
method <- "periodic"
poll <- "pm10"
lags <- c(0, 1, 2)
df.Season <- 1
dfVec <- seq(3, 13, 2)
results <- vector("list", length = length(lags))
names(results) <- paste("Lag", lags, sep = "")

postProcess <- function(x) { 
    covm <- x$cov
    n <- rownames(covm)
    i <- grep("pm10", n, fixed = TRUE)
    x$cov <- covm[i, i]
    cf <- x$coefficients
    n <- rownames(cf)
    i <- grep("pm10", n, fixed = TRUE)
    x$coefficients <- cf[i, ]
    x
}

source(file.path("data", "model-fitting.R"))

for(i in seq(along = lags)) {
    lag <- lags[i]
    cat(lag, " ")
    pollutant <- paste("Lag(", poll, "tmean, ", lag, ")", sep = "")

    lagresults <- lapply(cc, function(city) { 
        cat(city, " ") 
        cityresults <- multiDFFit(dfVec, city, pollutant = pollutant,
                                  cause = cause, season = method,
                                  df.Season = df.Season,
                                  extractors = extractors)
        names(cityresults) <- paste("df", dfVec, sep = "")
        lapply(cityresults, postProcess)
    })
    names(lagresults) <- cc
    results[[i]] <- lagresults
}

outfile <- paste("seasonal", method, "sens", "lag", paste(lags, collapse = ""),
                 "df", df.Season, poll, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)