# this function is a copy of the "ssm" function from the package "sspir"
# with the exception that lines 
ssm2 <- function (formula, family = gaussian, data = list(), subset = NULL, 
    fit = TRUE, phi = NULL, m0 = NULL, C0 = NULL, Fmat = NULL, 
    Gmat = NULL, Vmat = NULL, Wmat = NULL) 
{
    noTvar <- function(e) {
        if (is.call(e)) 
            if (e[[1]] == as.name("tvar")) 
                return(e[[2]])
            else {
                for (i in 2:length(e)) e[[i]] <- noTvar(e[[i]])
                return(e)
            }
        else e
    }
    removeSpace <- function(string) gsub("[[:space:]]", "", string)
    makeSS <- function(x, y, tvar, tvarNames, ...) {
        dt = rep(1, nrow(x))
        p <- ncol(x)
        ntvar <- length(unique(tvar[tvar != 0]))
        phi <- c(1, rep(1, ntvar))
        names(phi) <- c("epsilon", tvarNames)
        res <- SS(y = ts(matrix(y, ncol = 1), start = start(y), 
            end = end(y), frequency = frequency(y)), x = list(x = x, 
            dt = dt, tvar = tvar), phi = phi, Fmat = function(tt, 
            x, phi) {
            x$x[tt, ]
        }, Gmat = function(tt, x, phi) {
            diag(p)
        }, Vmat = function(tt, x, phi) {
            matrix(phi["epsilon"], 1, 1)
        }, Wmat = function(tt, x, phi) {
            W <- matrix(0, p, p)
            diag(W)[tvar != 0] <- phi[-1][tvar]
            return(W)
        }, C0 = diag(p) * 1e+06, m0 = matrix(0, 1, p))
        return(res)
    }
    cl <- match.call()
    if (is.character(family)) 
        family <- get(family, mode = "function", envir = parent.frame())
    if (is.function(family)) 
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("`family' not recognized")
    }
    if (missing(data)) 
        data <- environment(formula)
    mf <- match.call(expand.dots = FALSE)
    formula.notvar <- noTvar(formula)
    mf$formula <- formula.notvar
    mf$family <- mf$start <- NULL
    mf$drop.unused.levels <- TRUE
    mf$na.action <- "na.pass"
    mf$fit <- NULL
    mf$phi <- NULL
    mf$m0 <- NULL
    mf$C0 <- NULL
    mf$Fmat <- NULL
    mf$Gmat <- NULL
    mf$Wmat <- NULL
    mf$Vmat <- NULL
    mf[[1]] <- as.name("model.frame")
    y <- eval(mf[[2]][[2]], data)
    time <- time(ts(start = start(y), end = end(y), frequency = frequency(y)))
    assign("time", time, env = .GlobalEnv)
    mf <- eval(mf, .GlobalEnv)
    mt <- attr(mf, "terms")
    xvars <- as.character(attr(mt, "variables"))[-1]
    if ((yvar <- attr(mt, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(mf[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    y <- model.response(mf, "numeric")
    x <- if (!is.empty.model(mt)) 
        model.matrix(mt, mf)
    else matrix(, NROW(y), 0)
    term.labels <- unname(sapply(attr(terms(formula.notvar), 
        "term.labels"), removeSpace))
    pAssign <- attr(x, "assign")
    if (attr(terms(formula), "intercept")) {
        pAssign <- pAssign + 1
        pAssign[1] <- 1
        term.labels <- c("(Intercept)", term.labels)
    }
    tvar.terms <- terms(formula, specials = "tvar", keep.order = TRUE)
    idx <- attr(tvar.terms, "specials")$tvar
    if (attr(tvar.terms, "response")) 
        idx <- idx - 1
    if (attr(tvar.terms, "intercept")) 
        if (!any(attr(tvar.terms, "term.labels") == "tvar(1)")) 
            idx <- idx + 1
    tvar <- attr(terms(formula.notvar, keep.order = TRUE), "term.labels")
    if (attr(tvar.terms, "intercept")) 
        tvar <- c("(Intercept)", tvar)
    tvar <- tvar[idx]
    tvar <- unname(sapply(tvar, removeSpace))
    tvarAssign <- match(pAssign, sort(match(tvar, term.labels)))
    tvarAssign[is.na(tvarAssign)] <- 0
    if (any(unlist(lapply(term.labels, substr, 1, 8)) == "polytime")) {
        d <- sum(unlist(lapply(colnames(x), substr, 1, 8)) == 
            "polytime")
        nm <- unlist(lapply(colnames(x), substr, 1, 8))
        polytime.idx <- (1:ncol(x))[nm == "polytime"]
        Pmat <- function(deg) {
            res <- matrix(0, deg + 1, deg + 1)
            diag(res) <- 1
            for (i in 1:deg) {
                res[1:i, i + 1] <- 1/factorial(i:1)
            }
            res
        }
        if (any(unlist(lapply(tvar, substr, 1, 8)) == "polytime")) {
            polnum <- pmatch("polytime", tvar)
            nytvarAssign <- tvarAssign
            nytvarAssign[tvarAssign > polnum] <- tvarAssign[tvarAssign > 
                polnum] + d - 1
            nytvarAssign[tvarAssign == polnum] <- polnum:(polnum + 
                d - 1)
            tvarAssign <- nytvarAssign
            if (polnum > 1) 
                nytvar <- tvar[1:(polnum - 1)]
            else nytvar <- c()
            if (polnum < length(tvar)) 
                nytvar <- c(nytvar, colnames(x)[polytime.idx], 
                  tvar[(polnum + 1):length(tvar)])
            tvar <- nytvar
        }
    }
    if (any(unlist(lapply(term.labels, substr, 1, 9)) == "sumseason")) {
        period <- 1 + sum(unlist(lapply(colnames(x), substr, 
            1, 9)) == "sumseason")
        Cmat <- function(ncol) rbind(-1, cbind(diag(ncol - 1), 
            0))
        nm <- unlist(lapply(colnames(x), substr, 1, 9))
        sumseason.idx <- (1:ncol(x))[nm == "sumseason"]
    }
    if (is.empty.model(mt)) {
        x <- NULL
        z <- list(coefficients = numeric(0), residuals = y, fitted.values = 0 * 
            y, rank = 0, df.residual = length(y))
    }
    else {
#-------------------------
# begin code modification
#-------------------------
# original code in ssm:
#        if (family$family == "binomial") {
#            ntotal <- y[1, ] + y[2, ]
#            y <- y[1, ]
#        }
#        else {
#            ntotal <- NA
#
# revised code to work on cols of a matrix rather than rows:
        if (family$family == "binomial") {
            ntotal <- y[,1] + y[,2]
            y <- y[,1]
        }
        else {
            ntotal <- NA
        }
#-----------------------
# end code modification
#-----------------------
        z <- list()
        z$ss <- makeSS(x, y, w, tvar = tvarAssign, tvarNames = tvar)
        if (any(unlist(lapply(term.labels, substr, 1, 9)) == 
            "sumseason")) {
            oldG <- z$ss$Gmat
            newG <- function(tt, x, phi) {
                Gmat <- oldG(tt, x, phi)
                Gmat[sumseason.idx, sumseason.idx] <- Cmat(period - 
                  1)
                return(Gmat)
            }
            z$ss$Gmat <- newG
            if (any(unlist(lapply(tvar, substr, 1, 9)) == "sumseason")) {
                oldW <- z$ss$Wmat
                newW <- function(tt, x, phi) {
                  Wmat <- oldW(tt, x, phi)
                  diag(Wmat)[sumseason.idx[-1]] <- 0
                  return(Wmat)
                }
                z$ss$Wmat <- newW
            }
        }
        if (any(unlist(lapply(term.labels, substr, 1, 8)) == 
            "polytime")) {
            oldG2 <- z$ss$Gmat
            newG2 <- function(tt, x, phi) {
                Gmat <- oldG2(tt, x, phi)
                Gmat[polytime.idx, polytime.idx] <- Pmat(d - 
                  1)
                return(Gmat)
            }
            z$ss$Gmat <- newG2
        }
    }
    class(z) <- c("ssm", "lm")
    z$contrasts <- attr(x, "contrasts")
    z$xlevels <- xlev
    z$call <- cl
    z$terms <- mt
    z$data <- data
    z$ss$family <- family
    z$ss$ntotal <- ntotal
    if (!is.null(phi)) 
        z$ss$phi[1:length(phi)] <- phi
    if (!is.null(m0)) 
        z$ss$m0 <- m0
    if (!is.null(C0)) 
        z$ss$C0 <- C0
    if (!is.null(Fmat)) 
        z$ss$Fmat <- Fmat
    if (!is.null(Gmat)) 
        z$ss$Gmat <- Gmat
    if (!is.null(Vmat)) 
        z$ss$Vmat <- Vmat
    if (!is.null(Wmat)) 
        z$ss$Wmat <- Wmat
    if (fit) 
        z$ss <- kfs(z)
    z
}

