###############################################################################
## Fit seasonal models
## Copyright (C) 2004, Roger D. Peng <rpeng@jhsph.edu> and Aidan McDermott
##     
## 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
###############################################################################

## NOTE: Before the function `fitCitySeason' can be used, the data
## must be preprocessed using the `seasonal' function in the
## `NMMAPSdata' package.  Type `?seasonal' after loading the
## `NMMAPSdata' package into R for more information.  `fitCitySeason'
## will not work correctly with the full/raw database.


## Fit a single city seasonal model using glm() and ns()
library(splines)

## Note: Just running `fitCitySeason(city)', where `city' is the data
## frame for a particular city, should fit the usual NMMAPS model,
## except with seasonal main effects.

fitCitySeason <- function(data, pollutant = "l1pm10tmean", cause = "death",
                          season = c("none", "smooth", "stepfun"),
                          tempModel = c("default", "rm7", "tempInt"),
                          dfyr.Time = 7, pdfyr.time = 0.15, df.Temp = 6,
                          df.Dew = 3, df.Season = 1, extractors = NULL) {
    season <- match.arg(season)
    tempModel <- match.arg(tempModel)

    ## Specify temperature model
    temp.info <- setupTemp(data, df.Temp, tempModel)
    data <- temp.info[["adj.data"]]
    temp.f <- temp.info[["temp.f"]]
        
    ## Need to modify degrees of freedom based on missingness of data
    sub <- data[, c("time", "agecat", "tmpd", "rmtmpd", "dptp", "rmdptp",
                    temp.info[["addedVars"]], cause, pollutant)]
    subset <- complete.cases(sub)
    df.Time <- round( numdf(subset, dfyr.Time) )
    df.time <- round( df.Time * pdfyr.time )

    ## Don't setup smooth function of time where there are incomplete cases
    is.na(data[, "time"]) <- !subset

    ## Setup other formulas
    covariates.f <- paste(cause, "~ dow + agecat + Season")
    weather.f <- paste(c(paste("ns(dptp,", df.Dew, ")"),
                         paste("ns(rmdptp,", df.Dew, ")")), 
                       collapse = " + ")

    ## Setup smooth functions of time
    time.f <- paste(c(paste("ns(time,", df.Time, ")"),
                      paste("I(ns(time,", df.time, ")*Age2Ind)"),
                      paste("I(ns(time,", df.time, ")*Age3Ind)")),
                    collapse = "+")
    
    poll.f <- switch(season,
                     none = paste(pollutant, collapse = " + "),
                     smooth = {
                         nfreq <- df.Season
                         paste("periodicBasis(SeasonTime,", nfreq,
                               ",365, intercept = TRUE):", pollutant,
                               sep = "", collapse = "+")
                     },
                     stepfun = {
                         paste(paste("Season", pollutant, sep = ":"),
                               collapse = "+")
                     }
                     )
    form.str <- paste(covariates.f, time.f, temp.f, weather.f,
                      poll.f, sep = " + ")
    modelFormula <- as.formula(form.str)

    ## Fit the model!
    fit <- glm(modelFormula, family = quasipoisson, data = data,
               control = glm.control(epsilon = 1e-10, maxit = 1000),
               na.action = na.omit)

    rval <- if(is.null(extractors))
        fit
    else 
        lapply(extractors, function(f) f(fit))
    invisible(rval)
}

periodicBasis <- function(x, nfreq, period, intercept = FALSE) {
    pi <- base::pi
    stopifnot(nfreq > 0)
    x <- as.numeric(x)
    N <- seq(0, nfreq - 1)
    k <- 2^N * 2 * pi / period
    M <- outer(x, k)
    sinM <- apply(M, 2, sin)
    cosM <- apply(M, 2, cos)

    if(!intercept) 
        cbind(sinM, cosM)
    else 
        cbind(1, sinM, cosM)
}

setupTemp <- function(dataframe, df.Temp, tempModel) {
    default.temp.f <- paste(c(paste("ns(tmpd,", df.Temp, ")"),
                              paste("ns(rmtmpd,", df.Temp, ")")),
                            collapse = " + ")
    orig.namelist <- names(dataframe)

    temp.f <- switch(tempModel,
                     default = default.temp.f,
                     rm7 = paste(default.temp.f, "ns(rm7tmpd, 3)", sep = "+"),
                     tempInt = {
                         paste("ns(tmpd,", df.Temp, "):", 
                               "ns(rmtmpd,", df.Temp, ")", sep = "")
                     })
    list(adj.data = dataframe, temp.f = temp.f,
         addedVars = setdiff(names(dataframe), orig.namelist))
}


## -----------------------------------------------------------------------
## Some support functions used in model fitting
## -----------------------------------------------------------------------

##--------------------------------------------
numdf <- function(usedata,num=5){
    n <- length(usedata)
    use <- usedata[1:(n/3)]        
    ll  <- round(length(use)/12)

    ## this is to eliminate the warning message the length of use is
    ## not a multiple of 12
    usenew <- use[1:(ll*12)]
    
    mat <- matrix(usenew, ncol=12,byrow=T)
    m   <- sum(ceiling(apply(mat,1,sum,na.rm=T)/12)) ##-365.25/12   
    df  <- round(12*m/365.25*num)   
    max(1,df)
}
##--------------------------------------------