###############################################################################
## 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
###############################################################################

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

## Create directory for storing intermediate results
if(!file.exists("results")) {
    dir.create("results")
}

######################################################################

## Build databases needed for analysis
## For copollutant models (from `seasonalCode' package)
source(file.path("data", "multipollutant.R"))
buildDB(multipollutant, "multipoll")  

######################################################################
## Initial data analysis with smooth periodic functions

## Compute smoothly varying seasonal pollution effects
rm(list = ls(all = TRUE))

registerDB(NULL)

cityname <- dget(file.path("data", "cityList.R"))

extractors <- list(coefficients = function(x) summary(x)$coefficients,
                   cov = vcov)

method <- "periodic"
df.Season <- 1
pollutants <- c("pm10tmean", "l1pm10tmean", "l2pm10tmean")
lags <- as.character(0:2)
results <- vector("list", length = length(lags))
names(results) <- paste("lag", 0:2, sep = "")

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

for(k in seq(along = pollutants)) {
    poll <- pollutants[k]
    lag <- lags[k]        
    lag.results <- vector("list", length = length(cityname))
    names(lag.results) <- cityname

    for(l in seq(along = cityname)) {
        cat(k, method, cityname[l], "\n", sep = " ")
        citydata <- readCity(cityname[l])
        lag.results[[l]] <- try({
          fitCitySeason(data = citydata, pollutant = poll, cause = "death", 
                        season = method, df.Season = df.Season, 
                        extractors = extractors)
        })
    }
    results[[k]] <- lag.results
}
outfile <- paste("seasonal", method, "lag", "012", "dfSeas",
                 df.Season, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)


######################################################################
## Fit multi-pollutant models

rm(list = ls(all = TRUE))

registerDB("multipoll")

cityList <- listDBCities()
extractors <- list(
                   coefficients = function(x) { summary(x)$coefficients },
                   cov = vcov
                   )
cause <- "death"
method <- "stepfun"
poll <- "pm10"
copoll <- "no2"
pollutant <- c(paste("l1", poll, "tmean", sep = ""), paste(copoll, "tmean", sep=""))

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

results <- cityApply(fitCitySeason1, cityList, verbose = TRUE,
                     pollutant = pollutant, cause = cause, season = method,
                     extractors = extractors)
outfile <- paste("seasonal", copoll, method, poll, "lag", 1, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)

copoll <- "o3"
pollutant <- c(paste("l1", poll, "tmean", sep = ""), paste(copoll, "tmean", sep=""))

results <- cityApply(fitCitySeason1, cityList, verbose = TRUE,
                     pollutant = pollutant, cause = cause, season = method,
                     extractors = extractors)
outfile <- paste("seasonal", copoll, method, poll, "lag", 1, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)

copoll <- "so2"
pollutant <- c(paste("l1", poll, "tmean", sep = ""), paste(copoll, "tmean", sep=""))

results <- cityApply(fitCitySeason1, cityList, verbose = TRUE,
                     pollutant = pollutant, cause = cause, season = method,
                     extractors = extractors)
outfile <- paste("seasonal", copoll, method, poll, "lag", 1, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)

## run pm10 only on restricted city list
pollutant <- "l1pm10tmean"

results <- cityApply(fitCitySeason1, cityList, verbose = TRUE,
                     pollutant = pollutant, cause = cause, season = method,
                     extractors = extractors)
outfile <- paste("seasonal", "none", method, poll, "lag", 1, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)

######################################################################
## Compute separate seasonal effects
rm(list = ls(all = TRUE))

registerDB(NULL)

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

method <- "factor2"
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(i, method, cityname, "\n", sep = " ")
        citydata <- readCity(cityname)
        lag.results[[l]] <- try({
            fitCitySeason(data = citydata, pollutant = pollutant, cause = "death",
                          season = method, extractors = extractors)
        })
    }
    results[[i]] <- lag.results
}
outfile <- paste("seasonal", method, "lag", "012", poll, "rda", sep = ".")
save(results, file = file.path("results", outfile), compress = TRUE)