reaster <- function(fixed, random, pred, fam, varvar, idvar, root,
    famlist = fam.default(), origin, data, effects, sigma, response)
    UseMethod("reaster")

reaster.default <- function(fixed, random, pred, fam, varvar, idvar, root,
    famlist = fam.default(), origin, data, effects, sigma, response)
{
    stopifnot(is.matrix(fixed))
    stopifnot(is.numeric(fixed))
    stopifnot(is.finite(fixed))

    if (! is.list(random))
        random <- list(random)
    for (irandom in seq(along = random)) {
        r <- random[[irandom]]
        if (! is.matrix(r))
            stop("random must be matrix or list of matrices")
        if (! is.numeric(r))
            stop("random must be numeric matrix or list of such")
        if (! all(is.finite(r)))
            stop("random must be all finite numeric matrix or list of such")
        if (nrow(r) != nrow(fixed))
            stop("fixed and random effect model matrices with different row dimensions")
    }

    stopifnot(is.vector(response))
    stopifnot(is.numeric(response))
    stopifnot(is.finite(response))
    if(length(response) != nrow(fixed))
        stop("length of response vector != row dimension of model matrices")

    stopifnot(is.vector(idvar))
    stopifnot(is.numeric(idvar))
    nind <- length(unique(idvar))

    stopifnot(is.factor(varvar))
    stopifnot(is.character(levels(varvar)))
    nnode <- nlevels(varvar)

    if (nind * nnode != length(response))
        stop("idvar or varvar wrong number of levels")

    varvarmat <- matrix(as.character(varvar), nind, nnode)
    idvarmat <- matrix(as.vector(idvar), nind, nnode)
    foo <- apply(varvarmat, 2, function(x) length(unique(x)))
    bar <- apply(idvarmat, 1, function(x) length(unique(x)))
    if (! (all(foo == 1) & all(bar == 1)))
        stop("idvar or varvar wrong structure")

    stopifnot(is.vector(pred))
    stopifnot(is.numeric(pred))
    stopifnot(pred %in% seq(0, nlevels(varvar)))

    stopifnot(is.vector(fam))
    stopifnot(is.numeric(fam))
    stopifnot(fam %in% seq(1, length(famlist)))

    stopifnot(is.vector(root))
    stopifnot(length(root) == length(response))
    stopifnot(is.numeric(root))
    stopifnot(is.finite(root))
    stopifnot(root > 0)

    if (! missing(origin)) {
        stopifnot(is.vector(origin))
        stopifnot(length(origin) == length(response))
        stopifnot(is.numeric(origin))
        stopifnot(is.finite(origin))
    }

    if (! missing(sigma)) {
        stopifnot(is.vector(sigma))
        stopifnot(is.numeric(sigma))
        stopifnot(is.finite(sigma))
        stopifnot(length(sigma) == length(as.list(random)))
    }

    ##### fit fixed effect aster model (even if effects supplied)
    ncoef <- ncol(fixed)
    modmat <- array(fixed, dim = c(nind, nnode, ncoef))
    dimnames(modmat) <- list(NULL, NULL, colnames(fixed))
    y <- matrix(response, nind, nnode)
    asterargs <- list(x = y, root = matrix(root, nind, nnode),
        pred = pred, fam = fam, modmat = modmat, famlist = famlist)
    if (! missing(origin))
        asterargs$origin <- matrix(origin, nind, nnode)
    aout <- try(do.call(aster, asterargs))
    if (inherits(aout, "try-error") || (! aout$converged)) {
        asterargs$method <- "nlm"
        aout <- try(do.call(aster, asterargs))
        if (inherits(aout, "try-error") || (! aout$converged))
            stop("cannot get convergence in fixed effects model")
    }

    ##### figure out what columns to drop from fixed effects matrix
    ikeep <- match(names(aout$coefficients), colnames(fixed))
    dropped <- colnames(fixed)[- ikeep]
    nfix.orig <- ncol(fixed)
    fixed <- fixed[ , ikeep]
    nfix <- ncol(fixed)
    idx.fixed <- seq(1, nfix)

    nrand <- sapply(random, ncol)
    blurfle <- names(random)
    if (is.null(blurfle)) {
        blurfle <- paste("R", seq(along = random), sep = "")
    } else {
        if (length(unique(blurfle)) != length(blurfle))
            stop("duplicate names for argument random")
    }
    splitter <- factor(rep(blurfle, times = nrand), levels = blurfle)
    idx.random <- nfix + seq(1, sum(nrand))

    if (! missing(effects)) {
        stopifnot(is.vector(effects))
        stopifnot(is.numeric(effects))
        stopifnot(is.finite(effects))
        if (length(dropped) != 0)
            if (length(effects) == nfix.orig + sum(nrand)) {
                warning("shortening effects to agree with columns dropped from fixed")
                effects.fixed <- effects[1:nfix.orig]
                effects.random <- effects[- (1:nfix.orig)]
                effects.fixed <- effects.fixed[ikeep]
                effects <- c(effects.fixed, effects.random)
            }
        if (length(effects) != ncol(fixed) + sum(sapply(random, ncol)))
            stop("length(effects) not sum of column dimensions of fixed and random effects model matrices")
    }

    if (missing(sigma) | missing(effects)) {
        sigma.start <- rep(1, length(nrand))
        eff.start <- c(aout$coefficients, rep(0, sum(nrand)))
        tout <- try(trust(objfun = penmlogl, eff.start,
            rinit = 1, rmax = 10,
            sigma = sigma.start, fixed = fixed, random = random,
            obj = aout, iterlim = 1000), silent = TRUE)
        if (inherits(tout, "try-error") || (! tout$converged))
            stop("cannot get convergence in penalized likelihood step")
        ### FIX ME ###
        ### need to try another optimization method when trust fails

        eff <- tout$argument * tout$scale
        eff.split <- split(eff[idx.random], splitter)
        sigma.start <- sapply(eff.split, function(x) sqrt(mean(x^2)))
    } else {
        sigma.start <- sigma
        eff.start <- effects
    }

    cache <- new.env(parent = emptyenv())
    oout <- suppressWarnings(try(optim(sigma.start, pickle, parm = eff.start,
        fixed = fixed, random = random, obj = aout, y = y, cache = cache),
        silent = TRUE))
    if (inherits(oout, "try-error") || oout$convergence != 0)
        stop("step 1 part 1 (optim Nelder-Mead with pickle) failed")

    sigma.mle <- oout$par
    trout <- trust(objfun = penmlogl, cache$parm, rinit = 1, rmax = 10,
        sigma = sigma.mle, fixed = fixed, random = random, obj = aout, y = y)
    if (! trout$converged)
        stop("step 1 part 2 (trust with penmlogl) failed")
    parm.mle <- trout$argument
    zwz.mle <- makezwz(sigma.mle, parm = parm.mle,
        fixed = fixed, random = random, obj = aout, y = y)
    save.iter <- NULL

    # now iterate pickle3 and trust
    foo <- function(alphaceesigma)
        pickle3(alphaceesigma, fixed = fixed, random = random,
            obj = aout, y = y, zwz = zwz.mle, deriv = 2)
    repeat {
        tout <- trust(foo, c(parm.mle, sigma.mle), rinit = 1, rmax = 10,
            iterlim = 1000)
        if (! tout$converged)
            stop("step 2 (trust with pickle3) failed")
        save.iter <- c(save.iter, tout$iterations)
        sigma.old <- as.vector(sigma.mle)
        sigma.mle <- tout$argument[nfix + sum(nrand) + seq(along = nrand)]
        parm.mle <- tout$argument[seq(1, nfix + sum(nrand))]
        zwz.mle <- makezwz(sigma.mle, parm.mle,
            fixed = fixed, random = random, obj = aout, y = y)
        if (isTRUE(all.equal(sigma.old, as.vector(sigma.mle)))) break
    }

    alpha.mle <- parm.mle[1:nfix]
    c.mle <- parm.mle[- (1:nfix)]

    # fix up negative sigma

    idx <- rep(seq(along = sigma.mle), times = nrand)
    for (k in seq(along = sigma.mle))
        if (sigma.mle[k] < 0) {
            sigma.mle[k] <- (- sigma.mle[k])
            eek <- (idx == k)
            c.mle[eek] <- (- c.mle[eek])
    }

    a.mle <- rep(sigma.mle, times = nrand)
    b.mle <- a.mle * c.mle

    if (! is.null(colnames(fixed)))
        names(alpha.mle) <- colnames(fixed)
    if (! is.null(colnames(random)))
        names(sigma.mle) <- names(random)
    if (all(! is.null(sapply(random, colnames)))) {
        names(b.mle) <- Reduce(c, Map(colnames, random))
        names(c.mle) <- names(b.mle)
    }

    if (missing(origin)) origin <- NULL
    result <- list(obj = aout, fixed = fixed, random = random,
        dropped = dropped, sigma = sigma.mle, nu = sigma.mle^2,
        c = c.mle, b = b.mle, alpha = alpha.mle, zwz = zwz.mle,
        response = response, origin = origin,
        iterations = save.iter, counts = oout$counts)
    class(result) <- "reaster"
    return(result)
}

reaster.formula <- function(fixed, random, pred, fam, varvar, idvar, root,
    famlist = fam.default(), origin, data, effects, sigma, response)
{
    stopifnot(inherits(fixed, "formula"))
    if (! is.list(random))
        random <- list(random)
    if (!  all(sapply(random, function(x) inherits(x, "formula"))))
        stop("random must be formula or list of formulas")
    stopifnot(missing(data) || is.data.frame(data))

    save.fixed <- fixed
    save.random <- random

    oldopt <- options(na.action = na.fail)
    on.exit(options(oldopt))

    ##### stuff copied from glm.R and not understood #####
    ##### see also http://developer.r-project.org/model-fitting-functions.txt

    call <- match.call()
    if(missing(data))
        data <- environment(fixed)
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("fixed", "data", "varvar", "idvar", "root"), names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    fred <- as.list(mf)
    names(fred)[2] <- "formula"
    mf <- as.call(fred)
    mf <- eval.parent(mf)

    mt <- attr(mf, "terms")
    if (missing(response))
        response <- model.response(mf, "numeric")
    if (is.empty.model(mt)) {
        stop("empty model")
    } else {
        fixed <- model.matrix(mt, mf)
    }
    rownames(fixed) <- NULL

    varvar <- mf[["(varvar)"]]
    idvar <- mf[["(idvar)"]]
    root <- mf[["(root)"]]

    ##### end of stuff copied from glm.R and not understood #####

    for (irandom in seq(along = random)) {
        f <- random[[irandom]]
        foo <- lm(f, data = data, method = "model.frame")
        bar <- model.matrix(f, data = foo)
        rownames(bar) <- NULL
        random[[irandom]] <- bar
    }

    result <- NextMethod("reaster", response = response)

    result$formula <- list(fixed = save.fixed, random = save.random)
    result$call <- call
    class(result) <- c("reaster.formula", "reaster")
    return(result)
}

summary.reaster <- function(object, standard.deviation = TRUE, ...)
{
    stopifnot(inherits(object, "reaster"))
    stopifnot(inherits(object$obj, "aster"))
    stopifnot(is.logical(standard.deviation))
    stopifnot(length(standard.deviation) == 1)

    alpha <- object$alpha
    sigma <- object$sigma
    nu <- object$nu
    cee <- object$c
    bee <- object$b
    zwz <- object$zwz
    fixed <- object$fixed
    random <- object$random
    if (is.matrix(random))
        random <- list(random)
    nfix <- ncol(fixed)
    nrand <- sapply(random, ncol)

    obj <- object$obj
    y <- matrix(object$response, nrow = nrow(obj$x), ncol = ncol(obj$x))
    origin <- object$origin

    if (is.null(origin)) {
        iz <- is.zero(c(alpha, bee, nu), fixed, random, obj, y, zwz = zwz)
    } else {
        iz <- is.zero(c(alpha, bee, nu), fixed, random, obj, y,
            origin = origin, zwz = zwz)
    }
    nu[iz] <- 0
    sigma[iz] <- 0
    izbee <- rep(iz, times = nrand)
    bee[izbee] <- 0

    has.se <- c(rep(TRUE, length(alpha)), ! iz)
    if (all(iz)) {
        se.alpha <- sqrt(diag(solve(obj$fisher)))
        se.bee <- rep(NA_real_, length(bee))
        se.sigma <- rep(NA_real_, length(sigma))
        se.nu <- rep(NA_real_, length(nu))
    } else {
        subrandom <- random[! iz]
        subnu <- nu[! iz]
        subbee <- bee[! izbee]
        subzwz <- zwz[! izbee, , drop = FALSE]
        subzwz <- subzwz[ , ! izbee, drop = FALSE]
        if (is.null(origin)) {
            qout <- quickle(c(alpha, subnu), subbee, fixed, subrandom,
                obj, y, zwz = subzwz, deriv = 2)
        } else {
            qout <- quickle(c(alpha, subnu), subbee, fixed, subrandom,
                obj, y, origin = origin, zwz = subzwz, deriv = 2)
        }
        subfish <- qout$hessian
        eout <- eigen(subfish, symmetric = TRUE, only.values = TRUE)
        goodfish <- min(eout$values) > 0
        if (goodfish) {
            subfish.inv <- solve(subfish)
            se.subparm <- sqrt(diag(subfish.inv))
        } else {
            warning(paste("estimated Fisher information matrix not positive",
               "definite, making all standard errors infinite"))
            se.subparm <- rep(Inf, nrow(subfish))
        }
        se.parm <- rep(NA_real_, length(has.se))
        se.parm[has.se] <- se.subparm
        se.alpha <- se.parm[seq(along = alpha)]
        se.nu <- se.parm[- seq(along = alpha)]
        se.sigma <- se.nu / (2 * sigma)
    }

    foo <- alpha
    foo <- cbind(foo, se.alpha)
    foo <- cbind(foo, foo[ , 1] / foo[ , 2])
    foo <- cbind(foo, 2 * pnorm(- abs(foo[ , 3])))
    rownames(foo) <- colnames(object$fixed)
    colnames(foo) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")

    bar <- sigma
    bar <- cbind(bar, se.sigma)
    bar <- cbind(bar, bar[ , 1] / bar[ , 2])
    bar <- cbind(bar, pnorm(- abs(bar[ , 3])))
    blurfle <- names(object$random)
    if (is.null(blurfle))
        blurfle <- paste("sigma", seq(along = object$random), sep = "")
    rownames(bar) <- blurfle
    colnames(bar) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)/2")

    baz <- nu
    baz <- cbind(baz, se.nu)
    baz <- cbind(baz, baz[ , 1] / baz[ , 2])
    baz <- cbind(baz, pnorm(- abs(baz[ , 3])))
    blurfle <- names(object$random)
    if (is.null(blurfle))
        blurfle <- paste("nu", seq(along = object$random), sep = "")
    rownames(baz) <- blurfle
    colnames(baz) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)/2")

    return(structure(list(alpha = foo, sigma = bar, nu = baz,
        object = object, standard.deviation = standard.deviation,
        fisher = subfish), class = "summary.reaster"))
}

print.summary.reaster <-
    function (x, digits = max(3, getOption("digits") - 3),
        signif.stars = getOption("show.signif.stars"), ...)
{
    if (! is.null(x$object$call)) {
        cat("\nCall:\n")
        cat(paste(deparse(x$object$call), sep="\n", collapse="\n"),
            "\n\n", sep="")
    }

    cat("\nFixed Effects:\n")
    printCoefmat(x$alpha, digits = digits,
        signif.stars = signif.stars, na.print = "NA", ...)

    if (x$standard.deviation) {
        cat("\nSquare Roots of Variance Components (P-values are one-tailed):\n")
        printCoefmat(x$sigma, digits = digits,
            signif.stars = signif.stars, na.print = "NA", ...)
    } else {
        cat("\nVariance Components (P-values are one-tailed):\n")
        printCoefmat(x$nu, digits = digits,
            signif.stars = signif.stars, na.print = "NA", ...)
    }

    if (! is.null(x$bee)) {
        cat("\nRandom Effects:\n")
        printCoefmat(x$bee, digits = digits,
            signif.stars = signif.stars, na.print = "NA", ...)
    }

    return(invisible(x))
}