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

## Setup database for multipollutant analysis

multipollutant <- function(dataframe) {
    cityName <- as.character(dataframe$city[1])
    pollNames <- c("l1pm10tmean", "so2tmean", "o3tmean", "no2tmean")

    use <- complete.cases(dataframe[, pollNames])
    n <- sum(use)
    
    if(n < 90) {
        cat("excluded by # observations\t")
        return(NULL)
    }
    if(cityName %in% c("gdrp", "cinc", "clev", "evan", "indi", "milw",
                       "prov", "stlo")) {
        cat("excluded by hitlist\t")
        return(NULL)
    }
    
    Age2Ind <- as.numeric(dataframe[, "agecat"] == 2)
    Age3Ind <- as.numeric(dataframe[, "agecat"] == 3)
    n <- as.numeric(table(dataframe[, "agecat"])[1])
    yearday <- strptime(as.character(dataframe[1:n, "date"]), "%Y%m%d")$yday
    SeasonTime <- yearday + 1
    indlist <- NMMAPSdata:::genSeasonInd(dataframe)
    seas.ind <- with(indlist,
                     factor(winter * 1 + spring * 2 + summer * 3 + fall * 4,
                            labels = c("Winter", "Spring", "Summer", "Fall")))
    Season <- rep(seas.ind, 3)
    tmpd <- dataframe[, "tmpd"]
    rm7tmpd <- filter(tmpd, filter = c(0, rep(1/7, 7)), sides = 1, 
                      circular = FALSE)
    is.na(dataframe[, "death"]) <- as.logical(dataframe[, "markdeath"])
    is.na(dataframe[, "cvd"]) <- as.logical(dataframe[, "markcvd"])
    is.na(dataframe[, "resp"]) <- as.logical(dataframe[, "markresp"])
    cvdresp <- dataframe[, "cvd"] + dataframe[, "resp"]
    dataframe[, "dow"] <- as.factor(dataframe[, "dow"])
    dataframe[, "agecat"] <- as.factor(dataframe[, "agecat"])
    varList <- c("cvd", "death", "resp", "tmpd", "rmtmpd", "dptp", 
                 "rmdptp", "time", "agecat", "dow", "pm10tmean",
                 "so2tmean", "o3tmean", "no2tmean", 
                 paste("l", 1:7, "pm10tmean", sep = ""))
    data.frame(dataframe[, varList],
               cvdresp = cvdresp,
               SeasonTime = SeasonTime, 
               Season = Season,
               rm7tmpd = as.numeric(rm7tmpd),
               Age2Ind = Age2Ind,
               Age3Ind = Age3Ind)
}