.packageName <- "gss"
cdssden <- ## Evaluate conditional density
function (object,x,cond,int=NULL) {
    if (class(object)!="ssden")
        stop("gss error in cdssden: not a ssden object")
    if (nrow(cond)!=1)
        stop("gss error in cdssden: condition has to be a single point")
    xnames <- NULL
    for (i in colnames(object$mf))
        if (all(i!=colnames(cond))) xnames <- c(xnames,i)
    if (any(length(xnames)==c(0,ncol(object$mf))))
        stop("gss error in cdssden: not a conditional density")
    if (length(xnames)==1&is.vector(x)) {
        x <- data.frame(x)
        colnames(x) <- xnames
    }
    if (!all(sort(xnames)==sort(colnames(x))))
        stop("gss error in cdssden: mismatched variable names")
    ## Calculate normalizing constant
    if (is.null(int)) {
        if (length(xnames)==1) {
            ## Gauss-Legendre quadrature
            mn <- min(object$domain[,xnames])
            mx <- max(object$domain[,xnames])
            quad <- gauss.quad(200,c(mn,mx))
            xmesh <- data.frame(quad$pt)
            colnames(xmesh) <- xnames
        }
        else {
            ## Smolyak cubature
            domain <- object$domain[,colnames(x)]
            code <- c(15,14,13)
            quad <- smolyak.quad(ncol(x),code[ncol(x)-1])
            for (i in 1:ncol(x)) {
                wk <- x[,i]
                jk <- ssden(~wk,domain=data.frame(wk=domain[,i]),alpha=2,
                            id.basis=object$id.basis)
                quad$pt[,i] <- qssden(jk,quad$pt[,i])
                quad$wt <- quad$wt/dssden(jk,quad$pt[,i])
            }
            jk <- wk <- NULL
            xmesh <- data.frame(quad$pt)
            colnames(xmesh) <- colnames(x)
        }
        xx <- cond[rep(1,nrow(xmesh)),,drop=FALSE]
        int <- sum(dssden(object,cbind(xmesh,xx))*quad$wt)
    }
    ## Return value
    xx <- cond[rep(1,nrow(x)),,drop=FALSE]
    list(pdf=dssden(object,cbind(x,xx))/int,int=int)
}

cpssden <- ## Compute cdf for univariate conditional density
function(object,q,cond,int=NULL) {
    if (class(object)!="ssden")
        stop("gss error in cpssden: not a ssden object")
    xnames <- NULL
    for (i in colnames(object$mf))
        if (all(i!=colnames(cond))) xnames <- c(xnames,i)
    if (length(xnames)!=1)
        stop("gss error in cpssden: not a 1-D conditional density")
    mn <- min(object$domain[,xnames])
    mx <- max(object$domain[,xnames])
    if (is.null(int)) int <- cdssden(object,mn,cond)$int
    order.q <- rank(q)
    p <- q <- sort(q)
    q.dup <- duplicated(q)
    p[q<=mn] <- 0
    p[q>=mx] <- 1
    kk <- (1:length(q))[q>mn&q<mx]
    for (i in kk) {
        if (q.dup[i]) {
            p[i] <- p.dup
            next
        }
        nqd.l <- max(20,ceiling((q[i]-mn)/(mx-mn)*200))
        qd.l <- gauss.quad(nqd.l,c(mn,q[i]))
        p.l <- sum(cdssden(object,qd.l$pt,cond,int)$pdf*qd.l$wt)
        nqd.u <- max(20,ceiling((mx-q[i])/(mx-mn)*200))
        qd.u <- gauss.quad(nqd.u,c(q[i],mx))
        p.u <- sum(cdssden(object,qd.u$pt,cond,int)$pdf*qd.u$wt)
        p[i] <- p.dup <- p.l/(p.l+p.u)
    }
    p[order.q]
}

cqssden <- ## Compute quantiles for univariate conditional density
function(object,p,cond,int=NULL) {
    if (class(object)!="ssden")
        stop("gss error in cqssden: not a ssden object")
    xnames <- NULL
    for (i in colnames(object$mf))
        if (all(i!=colnames(cond))) xnames <- c(xnames,i)
    if (length(xnames)!=1)
        stop("gss error in cqssden: not a 1-D conditional density")
    mn <- min(object$domain[,xnames])
    mx <- max(object$domain[,xnames])
    if (is.null(int)) int <- cdssden(object,mn,cond)$int
    order.p <- rank(p)
    q <- p <- sort(p)
    p.dup <- duplicated(p)
    q[p<=0] <- mn
    q[p>=1] <- mx
    kk <- (1:length(p))[p>0&p<1]
    quad <- gauss.quad(200,c(mn,mx))
    q.wk <- quad$pt
    p.wk <- cumsum(quad$wt*cdssden(object,q.wk,cond,int)$pdf)
    for (i in kk) {
        if (p.dup[i]) {
            q[i] <- q.dup
            next
        }
        j <- which.min(abs(p[i]-p.wk))
        q0 <- q.wk[j]
        p0 <- cpssden(object,q0,cond,int)
        if (p0==p[i]) {
            q[i] <- q0
            next
        }
        if (p0<p[i]) {
            q.l <- q0
            p.l <- p0
            while (p0<p[i]) {
                j <- j+1
                q0 <- ifelse(is.null(q.wk[j]),mx,q.wk[j])
                p0 <- cpssden(object,q0,cond,int)
            }
            q.u <- q0
            p.u <- p0
        }
        else {
            q.u <- q0
            p.u <- p0
            while (p0>p[i]) {
                j <- j-1
                q0 <- ifelse(is.null(q.wk[j]),mn,q.wk[j])
                p0 <- cpssden(object,q0,cond,int)
            }
            q.l <- q0
            p.l <- p0
        }
        while (abs(p0-p[i])>1e-10) {
            q0 <- q.l+(p[i]-p.l)/(p.u-p.l)*(q.u-q.l)
            p0 <- cpssden(object,q0,cond,int)
            if (p0>p[i]) {
                q.u <- q0
                p.u <- p0
            }
            else {
                q.l <- q0
                p.l <- p0
            }
        }
        q[i] <- q.dup <- q0
    }
    q[order.p]
}
## Fit single smoothing parameter density
sspdsty <- function(s,r,q,cnt,qd.s,qd.r,qd.wt,prec,maxiter,alpha)
{
    nxi <- dim(r)[1]
    nobs <- dim(r)[2]
    nqd <- length(qd.wt)
    if (!is.null(s)) nnull <- dim(s)[1]
    else nnull <- 0
    nxis <- nxi+nnull
    if (is.null(cnt)) cnt <- 0
    ## cv function
    cv <- function(lambda) {
        fit <- .Fortran("dnewton",
                        cd=as.double(cd), as.integer(nxis),
                        as.double(10^(lambda+theta)*q), as.integer(nxi),
                        as.double(rbind(10^theta*r,s)), as.integer(nobs),
                        as.integer(sum(cnt)), as.integer(cnt),
                        as.double(t(rbind(10^theta*qd.r,qd.s))), as.integer(nqd),
                        as.double(qd.wt),
                        as.double(prec), as.integer(maxiter),
                        as.double(.Machine$double.eps),
                        wk=double(2*(nqd+nobs)+nxis*(nxis+4)+max(nxis,3)),
                        info=integer(1),PACKAGE="gss")
        if (fit$info==1) stop("gss error in ssden: Newton iteration diverges")
        if (fit$info==2) warning("gss warning in ssden: Newton iteration fails to converge")
        assign("cd",fit$cd,inherit=TRUE)
        assign("int",fit$wk[3],inherit=TRUE)
        cv <- alpha*fit$wk[2]-fit$wk[1]
        alpha.wk <- max(0,log.la0-lambda-5)*(3-alpha) + alpha
        alpha.wk <- min(alpha.wk,3)
        adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*fit$wk[2],0)
        cv+adj
    }
    ## initialization
    mu.r <- apply(qd.wt*t(qd.r),2,sum)/sum(qd.wt)
    v.r <- apply(qd.wt*t(qd.r^2),2,sum)/sum(qd.wt)
    mu.s <- apply(qd.wt*t(qd.s),2,sum)/sum(qd.wt)
    v.s <- apply(qd.wt*t(qd.s^2),2,sum)/sum(qd.wt)
    if (is.null(s)) theta <- 0
    else theta <- log10(sum(v.s-mu.s^2)/nnull/sum(v.r-mu.r^2)*nxi) / 2
    log.la0 <- log10(sum(v.r-mu.r^2)/sum(diag(q))) + theta
    ## lambda search
    cd <- rep(0,nxi+nnull)
    int <- NULL
    la <- log.la0
    repeat {
        mn <- la-1
        mx <- la+1
        zz <- nlm0(cv,c(mn,mx))
        if (min(zz$est-mn,mx-zz$est)>=1e-3) break
        else la <- zz$est
    }
    ## return
    jk1 <- cv(zz$est)
    c <- cd[1:nxi]
    if (nnull) d <- cd[nxi+(1:nnull)]
    else d <- NULL
    list(lambda=zz$est,theta=theta,c=c,d=d,int=int,cv=zz$min)
}

## Fit multiple smoothing parameter density
mspdsty <- function(s,r,q,cnt,qd.s,qd.r,qd.wt,prec,maxiter,alpha)
{
    nxi <- dim(r)[1]
    nobs <- dim(r)[2]
    nqd <- length(qd.wt)
    nq <- dim(q)[3]
    if (!is.null(s)) nnull <- dim(s)[1]
    else nnull <- 0
    nxis <- nxi+nnull
    if (is.null(cnt)) cnt <- 0
    ## cv function
    cv <- function(theta) {
        r.wk <- q.wk <- qd.r.wk <- 0
        for (i in 1:nq) {
            r.wk <- r.wk + 10^theta[i]*r[,,i]
            q.wk <- q.wk + 10^theta[i]*q[,,i]
            qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i]
        }
        fit <- .Fortran("dnewton",
                        cd=as.double(cd), as.integer(nxis),
                        as.double(10^lambda*q.wk), as.integer(nxi),
                        as.double(rbind(r.wk,s)), as.integer(nobs),
                        as.integer(sum(cnt)), as.integer(cnt),
                        as.double(t(rbind(qd.r.wk,qd.s))), as.integer(nqd),
                        as.double(qd.wt),
                        as.double(prec), as.integer(maxiter),
                        as.double(.Machine$double.eps),
                        wk=double(2*(nqd+nobs)+nxis*(nxis+4)+max(nxis,3)),
                        info=integer(1),PACKAGE="gss")
        if (fit$info==1) stop("gss error in ssden: Newton iteration diverges")
        if (fit$info==2) warning("gss warning in ssden: Newton iteration fails to converge")
        assign("cd",fit$cd,inherit=TRUE)
        assign("int",fit$wk[3],inherit=TRUE)
        cv <- alpha*fit$wk[2]-fit$wk[1]
        alpha.wk <- max(0,theta-log.th0-5)*(3-alpha) + alpha
        alpha.wk <- min(alpha.wk,3)
        adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*fit$wk[2],0)
        cv+adj
    }
    cv.wk <- function(theta) cv.scale*cv(theta)+cv.shift
    ## initialization
    theta <- -log10(apply(q,3,function(x)sum(diag(x))))
    r.wk <- q.wk <- qd.r.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
        qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i]
    }
    ## theta adjustment
    z <- sspdsty(s,r.wk,q.wk,cnt,qd.s,qd.r.wk,qd.wt,prec,maxiter,alpha)
    theta <- theta + z$theta
    for (i in 1:nq) {
        theta[i] <- 2*theta[i] + log10(t(z$c)%*%q[,,i]%*%z$c)
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
        qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i]
    }
    mu <- apply(qd.wt*t(qd.r.wk),2,sum)/sum(qd.wt)
    v <- apply(qd.wt*t(qd.r.wk^2),2,sum)/sum(qd.wt)
    log.la0 <- log10(sum(v-mu^2)/sum(diag(q.wk)))
    log.th0 <- theta-log.la0
    ## lambda search
    z <- sspdsty(s,r.wk,q.wk,cnt,qd.s,qd.r.wk,qd.wt,prec,maxiter,alpha)
    lambda <- z$lambda
    log.th0 <- log.th0 + z$lambda
    theta <- theta + z$theta
    cd <- c(z$c,z$d)
    int <- z$int
    ## theta search
    counter <- 0
    ## scale and shift cv
    tmp <- abs(cv(theta))
    cv.scale <- 1
    cv.shift <- 0
    if (tmp<1&tmp>10^(-4)) {
        cv.scale <- 10/tmp
        cv.shift <- 0
    }
    if (tmp<10^(-4)) {
        cv.scale <- 10^2
        cv.shift <- 10
    }
    repeat {
        zz <- nlm(cv.wk,theta,stepmax=1,ndigit=7)
        if (zz$code<=3)  break
        theta <- zz$est        
        counter <- counter + 1
        if (counter>=5) {
            warning("gss warning in ssden: CV iteration fails to converge")
            break
        }
    }
    ## return
    jk1 <- cv(zz$est)
    c <- cd[1:nxi]
    if (nnull) d <- cd[nxi+(1:nnull)]
    else d <- NULL
    list(lambda=lambda,theta=zz$est,c=c,d=d,int=int,cv=zz$min)
}
dssden <- ## Evaluate density estimate
function (object,x) {
    if (class(object)!="ssden") stop("gss error in dssden: not a ssden object")
    if (dim(object$mf)[2]==1&is.vector(x)) {
        x <- data.frame(x)
        colnames(x) <- colnames(object$mf)
    }
    s <- NULL
    r <- matrix(0,dim(x)[1],length(object$id.basis))
    nq <- 0
    for (label in object$terms$labels) {
        xx <- object$mf[object$id.basis,object$terms[[label]]$vlist]
        x.new <- x[,object$terms[[label]]$vlist]
        nphi <- object$terms[[label]]$nphi
        nrk <-  object$terms[[label]]$nrk
        if (nphi) {
            phi <-  object$terms[[label]]$phi
            for (i in 1:nphi) {
                s <- cbind(s,phi$fun(x.new,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            rk <- object$terms[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq + 1
                r <- r + 10^object$theta[nq]*rk$fun(x.new,xx,nu=i,env=rk$env,out=TRUE)
            }
        }
    }
    as.vector(exp(cbind(s,r)%*%c(object$d,object$c))/object$int)
}

pssden <- ## Compute cdf for univariate density estimate
function(object,q) {
    if (class(object)!="ssden") stop("gss error in pssden: not a ssden object")
    if (dim(object$mf)[2]!=1) stop("gss error in pssden: not a 1-D density")
    mn <- min(object$domain)
    mx <- max(object$domain)
    order.q <- rank(q)
    p <- q <- sort(q)
    q.dup <- duplicated(q)
    p[q<=mn] <- 0
    p[q>=mx] <- 1
    kk <- (1:length(q))[q>mn&q<mx]
    for (i in kk) {
        if (q.dup[i]) {
            p[i] <- p.dup
            next
        }
        nqd.l <- max(20,ceiling((q[i]-mn)/(mx-mn)*200))
        qd.l <- gauss.quad(nqd.l,c(mn,q[i]))
        p.l <- sum(dssden(object,qd.l$pt)*qd.l$wt)
        nqd.u <- max(20,ceiling((mx-q[i])/(mx-mn)*200))
        qd.u <- gauss.quad(nqd.u,c(q[i],mx))
        p.u <- sum(dssden(object,qd.u$pt)*qd.u$wt)
        p[i] <- p.dup <- p.l/(p.l+p.u)
    }
    p[order.q]
}

qssden <- ## Compute quantiles for univariate density estimate
function(object,p) {
    if (class(object)!="ssden") stop("gss error in qssden: not a ssden object")
    if (dim(object$mf)[2]!=1) stop("gss error in qssden: not a 1-D density")
    mn <- min(object$domain)
    mx <- max(object$domain)
    order.p <- rank(p)
    q <- p <- sort(p)
    p.dup <- duplicated(p)
    q[p<=0] <- mn
    q[p>=1] <- mx
    kk <- (1:length(p))[p>0&p<1]
    q.wk <- object$quad$pt[,1]
    p.wk <- cumsum(object$quad$wt*dssden(object,q.wk))
    for (i in kk) {
        if (p.dup[i]) {
            q[i] <- q.dup
            next
        }
        j <- which.min(abs(p[i]-p.wk))
        q0 <- q.wk[j]
        p0 <- pssden(object,q0)
        if (p0==p[i]) {
            q[i] <- q0
            next
        }
        if (p0<p[i]) {
            q.l <- q0
            p.l <- p0
            while (p0<p[i]) {
                j <- j+1
                q0 <- ifelse(is.null(q.wk[j]),mx,q.wk[j])
                p0 <- pssden(object,q0)
            }
            q.u <- q0
            p.u <- p0
        }
        else {
            q.u <- q0
            p.u <- p0
            while (p0>p[i]) {
                j <- j-1
                q0 <- ifelse(is.null(q.wk[j]),mn,q.wk[j])
                p0 <- pssden(object,q0)
            }
            q.l <- q0
            p.l <- p0
        }
        while (abs(p0-p[i])>1e-10) {
            q0 <- q.l+(p[i]-p.l)/(p.u-p.l)*(q.u-q.l)
            p0 <- pssden(object,q0)
            if (p0>p[i]) {
                q.u <- q0
                p.u <- p0
            }
            else {
                q.l <- q0
                p.l <- p0
            }
        }
        q[i] <- q.dup <- q0
    }
    q[order.p]
}
##%%%%%%%%%%  Binomial Family %%%%%%%%%%

## Make pseudo data for logistic regression
mkdata.binomial <- function(y,eta,wt,offset)
{
    if (is.vector(y)) y <- as.matrix(y)
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    if (dim(y)[2]==1) {
        if ((max(y)>1)|(min(y)<0))
            stop("gss error: binomial responses should be between 0 and 1")
    }
    else {
        if (min(y)<0)
            stop("gss error: paired binomial response should be nonnegative")
        wt <- wt * (y[,1]+y[,2])
        y <- y[,1]/(y[,1]+y[,2])
    }
    p <- 1-1/(1+exp(eta))
    u <- p - y
    w <- p*(1-p)
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt)
}

## Calculate deviance residuals for logistic regression
dev.resid.binomial <- function(y,eta,wt)
{
    if (is.vector(y)) y <- as.matrix(y)
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (dim(y)[2]>1) {
        wt <- wt * (y[,1]+y[,2])
        y <- y[,1]/(y[,1]+y[,2])
    }
    p <- 1-1/(1+exp(eta))
    as.vector(2*wt*(y*log(ifelse(y==0,1,y/p))
                    +(1-y)*log(ifelse(y==1,1,(1-y)/(1-p)))))
}

## Calculate null deviance for logistic regression
dev.null.binomial <- function(y,wt,offset)
{
    if (is.vector(y)) y <- as.matrix(y)
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (dim(y)[2]>1) {
        wt <- wt * (y[,1]+y[,2])
        y <- y[,1]/(y[,1]+y[,2])
    }
    p <- sum(wt*y)/sum(wt)
    if (!is.null(offset)) {
        eta <- log(p/(1-p)) - mean(offset)
        repeat {
            p <- 1-1/(1+exp(eta+offset))
            u <- p - y
            w <- p*(1-p)
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
    }
    sum(2*wt*(y*log(ifelse(y==0,1,y/p))
              +(1-y)*log(ifelse(y==1,1,(1-y)/(1-p)))))
}


##%%%%%%%%%%  Poisson Family %%%%%%%%%%

## Make pseudo data for Poisson regression
mkdata.poisson <- function(y,eta,wt,offset)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    if (is.null(offset)) offset <- rep(0,length(y))
    if (min(y)<0)
        stop("gss error: Poisson response should be nonnegative")
    lambda <- exp(eta)
    u <- lambda - y
    w <- lambda
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt)
}

## Calculate deviance residuals for Poisson regression
dev.resid.poisson <- function(y,eta,wt)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    lambda <- exp(eta)
    as.vector(2*wt*(y*log(ifelse(y==0,1,y/lambda))-(y-lambda)))
}

## Calculate null deviance for Poisson regression
dev.null.poisson <- function(y,wt,offset)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    lambda <- sum(wt*y)/sum(wt)
    if (!is.null(offset)) {
        eta <- log(lambda) - mean(offset)
        repeat {
            lambda <- exp(eta+offset)
            u <- lambda - y
            w <- lambda
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
    }
    sum(2*wt*(y*log(ifelse(y==0,1,y/lambda))-(y-lambda)))
}


##%%%%%%%%%%  Gamma Family %%%%%%%%%%

## Make pseudo data for Gamma regression
mkdata.Gamma <- function(y,eta,wt,offset)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    if (is.null(offset)) offset <- rep(0,length(y))
    if (min(y)<=0)
        stop("gss error: gamma responses should be positive")
    mu <- exp(eta)
    u <- 1-y/mu
    w <- y/mu
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt)
}

## Calculate deviance residuals for Gamma regression
dev.resid.Gamma <- function(y,eta,wt)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    mu <- exp(eta)
    as.vector(2*wt*(-log(y/mu)+(y-mu)/mu))
}

## Calculate null deviance for Gamma regression
dev.null.Gamma <-
function(y,wt,offset) {
  if (is.null(wt)) wt <- rep(1,length(y))
  mu <- sum(wt*y)/sum(wt)
  if (!is.null(offset)) {
    eta <- log(mu)-mean(offset)
    repeat {
      mu <- exp(eta+offset)
      u <- 1-y/mu
      w <- y/mu
      eta.new <- eta-sum(wt*u)/sum(wt*w)
      if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
      eta <- eta.new    
    }
  }
  sum(2*wt*(-log(y/mu)+(y-mu)/mu))
}


##%%%%%%%%%%  Inverse Gaussian Family %%%%%%%%%%

## Make pseudo data for IG regression
mkdata.inverse.gaussian <- function(y,eta,wt,offset)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    if (is.null(offset)) offset <- rep(0,length(y))
    if (min(y)<=0)
        stop("gss error: inverse gaussian responses should be positive")
    mu <- exp(eta)
    u <- (1-y/mu)/mu
    w <- 1/mu
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt)
}

## Calculate deviance residuals for IG regression
dev.resid.inverse.gaussian <- function(y,eta,wt)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    mu <- exp(eta)
    as.vector(wt*((y-mu)^2/(y*mu^2)))
}

## Calculate null deviance for IG regression
dev.null.inverse.gaussian <- function(y,wt,offset)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    mu <- sum(wt*y)/sum(wt)
    if (!is.null(offset)) {
        eta <- log(mu)-mean(offset)
        repeat {
            mu <- exp(eta+offset)
            u <- (1-y/mu)/mu
            w <- 1/mu
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
    }
    sum(wt*((y-mu)^2/(y*mu^2)))
}


##%%%%%%%%%%  Negative Binomial Family %%%%%%%%%%

## Make pseudo data for NB regression
mkdata.nbinomial <- function(y,eta,wt,offset,nu)
{
    if (is.vector(y)) y <- as.matrix(y)
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    if (dim(y)[2]==2) {
        if (min(y[,1])<0)
            stop("gss error: negative binomial response should be nonnegative")
        if (min(y[,2])<=0)
            stop("gss error: negative binomial size should be positive")
        p <- 1-1/(1+exp(eta))
        u <- (y[,1]+y[,2])*p-y[,2]
        w <- (y[,1]+y[,2])*p*(1-p)
        ywk <- eta-u/w-offset
        wt <- w*wt
        list(ywk=ywk,wt=wt)
    }
    else {
        if (min(y)<0)
            stop("gss error: negative binomial response should be nonnegative")
        p <- 1-1/(1+exp(eta))
        if (is.null(nu)) log.nu <- log(mean(y*exp(eta)))
        else log.nu <- log(nu)
        repeat {
            nu <- exp(log.nu)
            ua <- sum(digamma(y+nu)-digamma(nu)+log(p))*nu
            wa <- sum(trigamma(y+nu)-trigamma(nu))*nu*nu+ua
            log.nu.new <- log.nu - ua/wa
            if (abs(log.nu-log.nu.new)/(1+abs(log.nu))<1e-7) break
            log.nu <- log.nu.new
        }
        u <- (y+nu)*p-nu
        w <- (y+nu)*p*(1-p)
        ywk <- eta-u/w-offset
        wt <- w*wt
        list(ywk=ywk,wt=wt,nu=nu)
    }
}

## Calculate deviance residuals for NB regression
dev.resid.nbinomial <- function(y,eta,wt)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    p <- 1-1/(1+exp(eta))
    as.vector(2*wt*(y[,1]*log(ifelse(y[,1]==0,1,y[,1]/(y[,1]+y[,2])/(1-p)))
                    +y[,2]*log(y[,2]/(y[,1]+y[,2])/p)))
}

## Calculate null deviance for NB regression
dev.null.nbinomial <- function(y,wt,offset)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    p <- sum(wt*y[,2])/sum(wt*y)
    if (!is.null(offset)) {
        eta <- log(p/(1-p)) - mean(offset)
        repeat {
            p <- 1-1/(1+exp(eta+offset))
            u <- (y[,1]+y[,2])*p-y[,2]
            w <- (y[,1]+y[,2])*p*(1-p)
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
    }
    sum(2*wt*(y[,1]*log(ifelse(y[,1]==0,1,y[,1]/(y[,1]+y[,2])/(1-p)))
              +y[,2]*log(y[,2]/(y[,1]+y[,2])/p)))
}
##%%%%%%%%%%  Binomial Family %%%%%%%%%%

## Calculate CV score for binomial regression
cv.binomial <- function(y,eta,wt,hat,alpha)
{
    if (is.vector(y)) y <- as.matrix(y)
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (dim(y)[2]==1) {
        if ((max(y)>1)|(min(y)<0))
            stop("gss error: binomial responses should be between 0 and 1")
        m <- rep(1,dim(y)[1])
    }
    else {
        if (min(y)<0)
            stop("gss error: paired binomial response should be nonnegative")
        m <- y[,1]+y[,2]
        y <- y[,1]/m
    }
    wtt <- wt * m
    p <- 1-1/(1+exp(eta))
    w <- p*(1-p)
    lkhd <- -sum(wtt*(y*eta+log(1-p)))/sum(wtt)
    aux1 <- sum(hat/w)/(sum(wtt)-sum(hat))
    aux2 <- sum(wtt*y*(1-p))/sum(wtt)
    list(score=lkhd+abs(alpha)*aux1*aux2,varht=1,w=as.vector(wtt*w))
}


##%%%%%%%%%%  Poisson Family %%%%%%%%%%

## Calculate CV score for Poisson regression
cv.poisson <- function(y,eta,wt,hat,alpha,sr,q)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    if (min(y)<0)
        stop("gss error: Poisson response should be nonnegative")
    nxi <- ncol(q)
    nn <- ncol(sr)
    nnull <- nn-nxi
    lambda <- exp(eta)
    w <- as.vector(lambda)
    lkhd <- -sum(wt*(y*eta-lambda))/sum(wt*y)
    ## matrix H
    mu <- apply(wt*w*sr,2,sum)/sum(wt*w)
    v <- t(sr)%*%(wt*w*sr)/sum(wt*w)-outer(mu,mu)
    v[(nnull+1):nn,(nnull+1):nn] <- v[(nnull+1):nn,(nnull+1):nn]+q/sum(wt*y)
    ## Cholesky decomposition of H
    z <- .Fortran("dchdc",
                  v=as.double(v), as.integer(nn), as.integer(nn),
                  double(nn), jpvt=as.integer(rep(0,nn)),
                  as.integer(1), rkv=integer(1),
                  PACKAGE="base")[c("v","jpvt","rkv")]
    v <- matrix(z$v,nn,nn)
    rkv <- z$rkv
    while (v[rkv,rkv]<v[1,1]*sqrt(.Machine$double.eps)) rkv <- rkv-1
    if (rkv<nn) v[(rkv+1):nn,(rkv+1):nn] <- diag(v[1,1],nn-rkv)
    ## trace
    mu <- apply(wt*y*sr,2,sum)/sum(wt*y)
    sr <- sqrt(wt*y)*t(t(sr)-mu)
    sr <- backsolve(v,t(sr[,z$jpvt]),tran=TRUE)
    aux1 <- sum(sr^2)
    aux2 <- 1/sum(wt*y)/(sum(wt*y)-1)
    list(score=lkhd+abs(alpha)*aux1*aux2,varht=1,w=as.vector(wt*w))
}


##%%%%%%%%%%  Gamma Family %%%%%%%%%%

## Calculate CV score for Gamma regression
cv.Gamma <- function(y,eta,wt,hat,rss,alpha)
{
    if (is.null(wt)) wt <- rep(1,length(y))
    if (min(y)<=0)
        stop("gss error: gamma responses should be positive")
    mu <- exp(eta)
    u <- 1-y/mu
    w <- y/mu
    lkhd <- sum(wt*(y/mu+eta))/sum(wt)
    aux1 <- sum(hat/w)/(sum(wt)-sum(hat))
    aux2 <- -sum(wt*y*u/mu)/sum(wt)
    list(score=lkhd+alpha*aux1*aux2,varht=rss/(1-mean(hat)),w=as.vector(wt*w))
}


##%%%%%%%%%%  Inverse Gaussian Family %%%%%%%%%%
############  THIS DOES NOT WORK  ##############
## Calculate CV score for inverse gaussian regression
#cv.inverse.gaussian <- function(y,eta,wt,hat,rss,alpha)
#{
#    if (is.null(wt)) wt <- rep(1,length(y))
#    if (min(y)<=0)
#        stop("gss error: inverse gaussian responses should be positive")
#    mu <- exp(eta)
#    u <- (1-y/mu)/mu
#    w <- 1/mu
#    eta1 <- eta+hat/(1-hat)*u/w
#    mu1 <- exp(eta1)
#    lkhd <- sum(wt*((y/mu/2-1)/mu))/sum(wt)
#    aux1 <- sum(hat/w)/(sum(wt)-sum(hat))
#    aux2 <- -sum(wt*y*u/mu/mu)/sum(wt)
#    aux3 <- -sum(wt*y*(1/2/mu/mu-1/2/mu1/mu1))/sum(wt)
#    list(score=lkhd+alpha*aux3,varht=rss/(1-mean(hat)),w=as.vector(wt*w))
#}
############  THIS DOES NOT WORK  ##############


##%%%%%%%%%%  Negative Binomial Family %%%%%%%%%%

## Calculate CV score for NB regression
cv.nbinomial <- function(y,eta,wt,hat,alpha)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    if (min(y[,1])<0)
        stop("gss error: negative binomial response should be nonnegative")
    if (min(y[,2])<=0)
        stop("gss error: negative binomial size should be positive")
    p <- 1-1/(1+exp(eta))
    u <- (y[,1]+y[,2])*p-y[,2]
    w <- (y[,1]+y[,2])*p*(1-p)
    lkhd <- sum(wt*(-(y[,1]+y[,2])*log(1-p)-y[,2]*eta))/sum(wt)
    lkhd <- lkhd+sum(wt*(lgamma(y[,2])-lgamma(y[,1]+y[,2])))/sum(wt)
    aux1 <- sum(hat/w)/(sum(wt)-sum(hat))
    aux2 <- sum(wt*y[,1]*u*p)/sum(wt)
    list(score=lkhd+alpha*aux1*aux2,varht=1,w=as.vector(wt*w))
}


##%%%%%%%%%%  Weibull Family %%%%%%%%%%

## Calculate CV score for Weibull regression
cv.weibull <- function(y,eta,wt,hat,nu,alpha)
{
    if (is.vector(y)) stop("gss error: missing censoring indicator")
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (any(zz<0)|any(zz>=xx))
        stop("gss error: inconsistent life time data")
    u <- nu*(delta-(xx^nu-zz^nu)*exp(-nu*eta))
    w <- nu^2*(xx^nu-zz^nu)*exp(-nu*eta)
    lkhd <- sum(wt*((xx^nu-zz^nu)*exp(-nu*eta)-delta*(nu*(log(xx)-eta)+log(nu))))/sum(wt)
    aux1 <- sum(hat/w)/(sum(wt)-sum(hat))
    aux2 <- sum(wt*nu*delta*u)/sum(wt)
    list(score=lkhd+alpha*aux1*aux2,varht=1,w=as.vector(wt*w))
}


##%%%%%%%%%%  Log Normal Family %%%%%%%%%%

## Calculate CV score for log normal regression
cv.lognorm <- function(y,eta,wt,hat,nu,alpha)
{
    if (is.vector(y)) stop("gss error: missing censoring indicator")
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (any(zz<0)|any(zz>=xx))
        stop("gss error: inconsistent life time data")
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    s.xx <- ifelse(xx<7,dnorm(xx)/(1-pnorm(xx)),xx+.15)
    s.zz <- ifelse(zz<7,dnorm(zz)/(1-pnorm(zz)),zz+.15)
    s.xx <- pmax(s.xx,s.zz)
    u <- nu*(delta*(s.xx-xx)-(s.xx-s.zz))
    w <- (s.xx^2/2-xx*s.xx+xx^2/2+log(s.xx)+log(2*pi)/2)
    w <- nu^2*(w-ifelse(s.zz==0,0,(s.zz^2/2-zz*s.zz+zz^2/2+log(s.zz)+log(2*pi)/2)))
    w <- ifelse(w<1e-6,1e-6,w)
    aux1 <- sum(hat/w)/(sum(wt)-sum(hat))
    aux2 <- sum(wt*nu*delta*(s.xx-xx)*u)/sum(wt)
    s.xx <- ifelse(xx<7,log(1-pnorm(xx)),-xx^2/2-log(xx+.15)-log(2*pi)/2)
    s.zz <- ifelse(zz<7,log(1-pnorm(zz)),-zz^2/2-log(zz+.15)-log(2*pi)/2)
    s.xx <- pmin(s.xx,s.zz)
    lkhd <- sum(wt*(delta*(xx^2/2+s.xx-log(nu))+s.zz-s.xx))/sum(wt)
    list(score=lkhd+alpha*aux1*aux2,varht=1,w=as.vector(wt*w))
}


##%%%%%%%%%%  Log Logistic Family %%%%%%%%%%

## Calculate CV score for log logistic regression
cv.loglogis <- function(y,eta,wt,hat,nu,alpha)
{
    if (is.vector(y)) stop("gss error: missing censoring indicator")
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (any(zz<0)|any(zz>=xx))
        stop("gss error: inconsistent life time data")
    xx <- 1/(1+exp(nu*(log(xx)-eta)))
    zz <- 1/(1+exp(nu*(log(zz)-eta)))
    u <- nu*(delta*xx-(zz-xx))
    w <- nu^2/2*(zz^2-xx^2)
    lkhd <- sum(wt*(delta*(-log(1-xx)-log(nu))+log(zz)-log(xx)))/sum(wt)
    aux1 <- sum(hat/w)/(sum(wt)-sum(hat))
    aux2 <- sum(wt*nu*delta*xx*u)/sum(wt)
    list(score=lkhd+alpha*aux1*aux2,varht=1,w=as.vector(wt*w))
}
##%%%%%%%%%%  Binomial Family %%%%%%%%%%
y0.binomial <- function(y,eta0,wt)
{
    if (is.matrix(y)) wt <- wt * (y[,1]+y[,2])
    p <- plogis(eta0)
    list(p=p,eta=eta0,wt=wt)
}
proj0.binomial <- function(y0,eta,offset)
{
    if (is.null(offset)) offset <- rep(0,length(eta))
    p <- plogis(eta)
    u <- p - y0$p
    w <- p*(1-p)
    ywk <- eta-u/w-offset
    wt <- w*y0$wt
    kl <- mean(y0$wt*(y0$p*(y0$eta-eta)+log((1-y0$p)/(1-p))))
    list(ywk=ywk,wt=wt,kl=kl,u=wt*u)
}
kl.binomial <- function(eta0,eta1,wt)
{
    p0 <- plogis(eta0)
    p1 <- plogis(eta1)
    mean(wt*(p0*(eta0-eta1)+log((1-p0)/(1-p1))))
}
cfit.binomial <- function(y,wt,offset)
{
    if (is.vector(y)) y <- as.matrix(y)
    if (dim(y)[2]>1) {
        wt <- wt * (y[,1]+y[,2])
        y <- y[,1]/(y[,1]+y[,2])
    }
    p <- sum(wt*y)/sum(wt)
    if (is.null(offset)) eta <- rep(qlogis(p),length(y))
    else {
        eta <- qlogis(p)-mean(offset)
        repeat {
            p <- plogis(eta+offset)
            u <- p - y
            w <- p*(1-p)
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
        eta <- eta + offset
    }
    eta
}


##%%%%%%%%%%  Poisson Family %%%%%%%%%%
y0.poisson <- function(eta0)
{
    lambda <- exp(eta0)
    list(lambda=lambda,eta=eta0)
}
proj0.poisson <- function(y0,eta,wt,offset)
{
    if (is.null(offset)) offset <- rep(0,length(eta))
    lambda <- exp(eta)
    u <- lambda - y0$lambda
    w <- lambda
    ywk <- eta-u/w-offset
    kl <- mean(wt*(y0$lambda*(y0$eta-eta)-y0$lambda+lambda))
    wt <- w*wt
    list(ywk=ywk,wt=wt,kl=kl,u=wt*u)
}
kl.poisson <- function(eta0,eta1,wt)
{
    lambda0 <- exp(eta0)
    lambda1 <- exp(eta1)
    mean(wt*(lambda0*(eta0-eta1)-lambda0+lambda1))
}
cfit.poisson <- function(y,wt,offset)
{
    lambda <- sum(wt*y)/sum(wt)
    if (is.null(offset)) eta <- rep(log(lambda),length(y))
    else {
        eta <- log(lambda) - mean(offset)
        repeat {
            lambda <- exp(eta+offset)
            u <- lambda - y
            w <- lambda
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
        eta <- eta + offset
    }
    eta
}


##%%%%%%%%%%  Gamma Family %%%%%%%%%%
y0.Gamma <- function(eta0)
{
    mu <- exp(eta0)
    list(mu=mu)
}
proj0.Gamma <- function(y0,eta,wt,offset)
{
    if (is.null(offset)) offset <- rep(0,length(y))
    mu <- exp(eta)
    u <- 1-y0$mu/mu
    w <- y0$mu/mu
    ywk <- eta-u/w-offset
    kl <- mean(wt*(y0$mu*(-1/y0$mu+1/mu)+log(mu/y0$mu)))
    wt <- w*wt
    list(ywk=ywk,wt=wt,kl=kl,u=wt*u)
}
kl.Gamma <- function(eta0,eta1,wt)
{
    mu0 <- exp(eta0)
    mu1 <- exp(eta1)
    mean(wt*(mu0*(-1/mu0+1/mu1)+log(mu1/mu0)))
}
cfit.Gamma <- function(y,wt,offset)
{
    mu <- sum(wt*y)/sum(wt)
    if (is.null(offset)) eta <- rep(log(mu),length(y))
    else {
        eta <- log(mu)-mean(offset)
        repeat {
            mu <- exp(eta+offset)
            u <- 1-y/mu
            w <- y/mu
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
        eta <- eta + offset
    }
    eta
}


##%%%%%%%%%%  Negative Binomial Family %%%%%%%%%%
y0.nbinomial <- function(y,eta0,nu)
{
    if (!is.vector(y)) {
        nu <- y[,2]
        y <- y[,1]
    }
    mu <- nu*exp(-eta0)
    list(y=y,nu=nu,mu=mu,eta=eta0)
}
proj0.nbinomial <- function(y0,eta,wt,offset)
{
    if (is.null(offset)) offset <- rep(0,length(eta))
    p <- plogis(eta)
    u <- (y0$mu+y0$nu)*p-y0$nu
    w <- (y0$mu+y0$nu)*p*(1-p)
    ywk <- eta-u/w-offset
    kl <- mean(wt*((y0$nu+y0$mu)*log((1+exp(eta))/(1+exp(y0$eta)))
                   +y0$nu*(y0$eta-eta)))
    wt <- w*wt
    list(ywk=ywk,wt=wt,kl=kl,u=wt*u)
}
kl.nbinomial <- function(eta0,eta1,wt,nu)
{
    mu0 <- nu*exp(-eta0)
    mean(wt*((nu+mu0)*log((1+exp(eta1))/(1+exp(eta0)))+nu*(eta0-eta1)))
}
cfit.nbinomial <- function(y,wt,offset,nu)
{
    if (!is.vector(y)) {
        nu <- y[,2]
        y <- y[,1]
    }
    p <- sum(wt*nu)/sum(wt*(y+nu))
    if (is.null(offset)) eta <- rep(qlogis(p),length(y))
    else {
        eta <- qlogis(p)-mean(offset)
        repeat {
            p <- 1-1/(1+exp(eta+offset))
            u <- (y+nu)*p-nu
            w <- (y+nu)*p*(1-p)
            eta.new <- eta-sum(wt*u)/sum(wt*w)
            if (abs(eta-eta.new)/(1+abs(eta))<1e-7) break
            eta <- eta.new    
        }
        eta <- eta + offset
    }
    eta
}


##%%%%%%%%%%  Weibull Family %%%%%%%%%%
y0.weibull <- function(y,eta0,nu)
{
    xx <- y[,1]
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    lam <- exp(-nu*eta0)
    list(lam=lam,eta=eta0,int=(xx^nu-zz^nu))
}
proj0.weibull <- function(y0,eta,wt,offset,nu)
{
    if (is.null(offset)) offset <- rep(0,length(eta))
    u <- nu*(y0$lam-exp(-nu*eta))
    w <- nu*nu*exp(-nu*eta)
    ywk <- eta-u/w-offset
    kl <- mean(wt*y0$int*(y0$lam*nu*(eta-y0$eta)+exp(-nu*eta)-y0$lam))
    u <- y0$int*u
    w <- y0$int*w
    wt <- w*wt
    list(ywk=ywk,wt=wt,kl=kl,u=wt*u)
}
kl.weibull <- function(eta0,eta1,wt,nu,int)
{
    lam0 <- exp(-nu*eta0)
    lam1 <- exp(-nu*eta1)
    mean(wt*int*(lam0*nu*(eta1-eta0)+lam1-lam0))
}
cfit.weibull <- function(y,wt,offset,nu)
{
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (is.null(offset)) offset <- rep(0,length(xx))
    eta <- log(sum(wt*(xx^nu-zz^nu)*exp(-nu*offset))/sum(wt*delta))/nu
    eta + offset    
}


##%%%%%%%%%%  Lognorm Family %%%%%%%%%%
y0.lognorm <- function(y,eta0,nu)
{
    xx <- y[,1]
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    quad <- gauss.quad(50,c(0,1))
    list(eta=eta0,xx=xx,zz=zz,q.pt=quad$pt,q.wt=quad$wt)
}
proj0.lognorm <- function(y0,eta,wt,offset,nu)
{
    if (is.null(offset)) offset <- rep(0,length(eta))
    u <- NULL
    kl <- 0
    for (i in 1:length(eta)) {
        q.pt <- y0$q.pt*(y0$xx[i]-y0$zz[i])+y0$zz[i]
        q.wt <- y0$q.wt*(y0$xx[i]-y0$zz[i])
        z0 <- nu*(log(q.pt)-y0$eta[i])
        z1 <- nu*(log(q.pt)-eta[i])
        lam0 <- ifelse(z0<7,dnorm(z0)/(1-pnorm(z0)),z0+.15)
        lam1 <- ifelse(z1<7,dnorm(z1)/(1-pnorm(z1)),z1+.15)
        u <- c(u,nu*nu*sum(q.wt*(lam0-lam1)*(lam1-z1)/q.pt))
        kl <- kl + nu*sum(q.wt*(lam0*log(lam0/lam1)+lam1-lam0)/q.pt)
    }
    xx <- nu*(log(y0$xx)-eta)
    zz <- nu*(log(y0$zz)-eta)
    s.xx <- ifelse(xx<7,dnorm(xx)/(1-pnorm(xx)),xx+.15)
    s.zz <- ifelse(zz<7,dnorm(zz)/(1-pnorm(zz)),zz+.15)
    s.xx <- pmax(s.xx,s.zz)
    w <- (s.xx^2/2-xx*s.xx+xx^2/2+log(s.xx)+log(2*pi)/2)
    w <- nu^2*(w-ifelse(s.zz==0,0,(s.zz^2/2-zz*s.zz+zz^2/2+log(s.zz)+log(2*pi)/2)))
    w <- pmax(w,1e-6)
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt,kl=kl/length(eta),u=wt*u)
}
kl.lognorm <- function(eta0,eta1,wt,nu,y0)
{
    kl <- 0
    for (i in 1:length(eta0)) {
        q.pt <- y0$q.pt*(y0$xx[i]-y0$zz[i])+y0$zz[i]
        q.wt <- y0$q.wt*(y0$xx[i]-y0$zz[i])
        z0 <- nu*(log(q.pt)-eta0[i])
        z1 <- nu*(log(q.pt)-eta1[i])
        lam0 <- ifelse(z0<7,dnorm(z0)/(1-pnorm(z0)),z0+.15)
        lam1 <- ifelse(z1<7,dnorm(z1)/(1-pnorm(z1)),z1+.15)
        kl <- kl + nu*sum(q.wt*(lam0*log(lam0/lam1)+lam1-lam0)/q.pt)
    }
    kl/length(eta0)
}
cfit.lognorm <- function(y,wt,offset,nu)
{
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (is.null(offset)) offset <- rep(0,length(xx))
    lkhd <- function(eta) {
        eta <- eta + offset
        xx.wk <- nu*(log(xx)-eta)
        zz.wk <- nu*(log(zz)-eta)
        -sum(wt*(delta*(-xx.wk^2/2-log(1-pnorm(xx.wk)))
                 +log((1-pnorm(xx.wk))/(1-pnorm(zz.wk)))))
    }
    nlm(lkhd,mean(log(xx)-offset),stepmax=1)$est + offset
}


##%%%%%%%%%%  Loglogis Family %%%%%%%%%%
y0.loglogis <- function(y,eta0,nu)
{
    xx <- y[,1]
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    quad <- gauss.quad(50,c(0,1))
    list(eta=eta0,xx=xx,zz=zz,q.pt=quad$pt,q.wt=quad$wt)
}
proj0.loglogis <- function(y0,eta,wt,offset,nu)
{
    if (is.null(offset)) offset <- rep(0,length(eta))
    e0 <- exp(-nu*y0$eta)
    e1 <- exp(-nu*eta)
    kl <- sum(log((1+y0$xx^nu*e1)*(1+y0$zz^nu*e0)/(1+y0$zz^nu*e1)/(1+y0$xx^nu*e0))
              +nu*(eta-y0$eta)*log((1+y0$xx^nu*e0)/(1+y0$zz^nu*e0)))
    xx <- 1/(1+y0$xx^nu*e1)
    zz <- 1/(1+y0$zz^nu*e1)
    u <- -nu*(zz-xx)
    w <- nu^2/2*(zz^2-xx^2)
    for (i in 1:length(eta)) {
        q.pt <- y0$q.pt*(y0$xx[i]-y0$zz[i])+y0$zz[i]
        q.wt <- y0$q.wt*(y0$xx[i]-y0$zz[i])
        u[i] <- u[i]+nu^2*sum(q.wt*q.pt^(nu-1)*e0[i]
                              /(1+q.pt^nu*e0[i])/(1+q.pt^nu*e1[i]))
        kl <- kl + nu*sum(q.wt*q.pt^(nu-1)*e0[i]/(1+q.pt^nu*e0[i])
                          *log((1+q.pt^nu*e1[i])/(1+q.pt^nu*e0[i])))
    }
    w <- pmax(w,1e-6)
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt,kl=kl/length(eta),u=wt*u)
}
kl.loglogis <- function(eta0,eta1,wt,nu,y0)
{
    e0 <- exp(-nu*eta0)
    e1 <- exp(-nu*eta1)
    kl <- sum(log((1+y0$xx^nu*e1)*(1+y0$zz^nu*e0)/(1+y0$zz^nu*e1)/(1+y0$xx^nu*e0))
              +nu*(eta1-eta0)*log((1+y0$xx^nu*e0)/(1+y0$zz^nu*e0)))
    for (i in 1:length(eta0)) {
        q.pt <- y0$q.pt*(y0$xx[i]-y0$zz[i])+y0$zz[i]
        q.wt <- y0$q.wt*(y0$xx[i]-y0$zz[i])
        kl <- kl + nu*sum(q.wt*q.pt^(nu-1)*e0[i]/(1+q.pt^nu*e0[i])
                          *log((1+q.pt^nu*e1[i])/(1+q.pt^nu*e0[i])))
    }
    kl/length(eta0)
}
cfit.loglogis <- function(y,wt,offset,nu)
{
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (is.null(offset)) offset <- rep(0,length(xx))
    lkhd <- function(eta) {
        eta <- eta + offset
        xx.wk <- nu*(log(xx)-eta)
        zz.wk <- nu*(log(zz)-eta)
        -sum(wt*(delta*(xx.wk-log(1+exp(xx.wk)))
                 -log((1+exp(xx.wk))/(1+exp(zz.wk)))))
    }
    nlm(lkhd,mean(log(xx)-offset),stepmax=1)$est + offset
}
##%%%%%%%%%%  Weibull Family %%%%%%%%%%

## Make pseudo data for Weibull regression
mkdata.weibull <- function(y,eta,wt,offset,nu)
{
    if (is.vector(y)) stop("gss error: missing censoring indicator")
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (any(zz<0)|any(zz>=xx))
        stop("gss error: inconsistent life time data")
    if (nu[[2]]) {
        lkhd <- function(log.nu) {
            nu <- exp(log.nu)
            -sum(wt*(delta*(nu*(log(xx)-eta)+log.nu)
                 -(xx^nu-zz^nu)*exp(-nu*eta)))
        }
        if (is.null(nu[[1]])) nu[[1]] <- 1
        nu[[1]] <- exp(nlm(lkhd,log(nu[[1]]),stepmax=.5)$est)
    }
    u <- nu[[1]]*(delta-(xx^nu[[1]]-zz^nu[[1]])*exp(-nu[[1]]*eta))
    w <- nu[[1]]^2*(xx^nu[[1]]-zz^nu[[1]])*exp(-nu[[1]]*eta)
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt,nu=nu)
}

## Calculate deviance residuals for Weibull regression
dev.resid.weibull <- function(y,eta,wt,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    z <- -delta*nu*(log(xx)-eta)+(xx^nu-zz^nu)*exp(-nu*eta)
    as.numeric(2*wt*(z+delta*(log(xx^nu)-log(xx^nu-zz^nu)-1)))
}

## Calculate null deviance for Weibull regression
dev.null.weibull <- function(y,wt,offset,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (is.null(offset)) offset <- rep(0,length(xx))
    eta <- log(sum(wt*(xx^nu-zz^nu)*exp(-nu*offset))/sum(wt*delta))/nu
    eta <- eta + offset    
    z <- -delta*nu*(log(xx)-eta)+(xx^nu-zz^nu)*exp(-nu*eta)
    sum(2*wt*(z+delta*(log(xx^nu)-log(xx^nu-zz^nu)-1)))
}


##%%%%%%%%%%  Log Normal Family %%%%%%%%%%

## Make pseudo data for log normal regression
mkdata.lognorm <- function(y,eta,wt,offset,nu)
{
    if (is.vector(y)) stop("gss error: missing censoring indicator")
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (any(zz<0)|any(zz>=xx))
        stop("gss error: inconsistent life time data")
    if (nu[[2]]) {
        lkhd <- function(log.nu) {
            nu <- exp(log.nu)
            xx.wk <- nu*(log(xx)-eta)
            zz.wk <- nu*(log(zz)-eta)
            -sum(wt*(delta*(-xx.wk^2/2-log(1-pnorm(xx.wk))+log.nu)
                 +log((1-pnorm(xx.wk))/(1-pnorm(zz.wk)))))
        }
        if (is.null(nu[[1]])) nu[[1]] <- 1
        nu[[1]] <- exp(nlm(lkhd,log(nu[[1]]),stepmax=.5)$est)
    }
    xx <- nu[[1]]*(log(xx)-eta)
    zz <- nu[[1]]*(log(zz)-eta)
    s.xx <- ifelse(xx<7,dnorm(xx)/(1-pnorm(xx)),xx+.15)
    s.zz <- ifelse(zz<7,dnorm(zz)/(1-pnorm(zz)),zz+.15)
    s.xx <- pmax(s.xx,s.zz)
    u <- nu[[1]]*(delta*(s.xx-xx)-(s.xx-s.zz))
    w <- (s.xx^2/2-xx*s.xx+xx^2/2+log(s.xx)+log(2*pi)/2)
    w <- nu[[1]]^2*(w-ifelse(s.zz==0,0,(s.zz^2/2-zz*s.zz+zz^2/2+log(s.zz)+log(2*pi)/2)))
    w <- ifelse(w<1e-6,1e-6,w)
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt,nu=nu)
}

## Calculate deviance residuals for log normal regression
dev.resid.lognorm <- function(y,eta,wt,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    dev0 <- NULL
    for (i in 1:length(xx)) {
        if (!delta[i]|!zz[i]) dev0 <- c(dev0,0)
        else {
            fun.wk <- function(eta) {
                (nu*(log(xx[i])-eta))^2/2+log(1-pnorm(nu*(log(zz[i])-eta)))
            }
            dev0 <- c(dev0,nlm(fun.wk,log(xx[i]),stepmax=1)$min)
        }
    }
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    s.xx <- log(1-pnorm(xx))
    s.zz <- log(1-pnorm(zz))
    z <- -delta*(-xx^2/2-s.xx)-s.xx+s.zz
    as.numeric(2*wt*(z-dev0))
}

dev0.resid.lognorm <- function(y,eta,wt,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    s.xx <- ifelse(xx<7,log(1-pnorm(xx)),-xx^2/2-log(xx+.15)-log(2*pi)/2)
    s.zz <- ifelse(zz<7,log(1-pnorm(zz)),-zz^2/2-log(zz+.15)-log(2*pi)/2)
    s.xx <- pmin(s.xx,s.zz)
    z <- -delta*(-xx^2/2-s.xx)-s.xx+s.zz
    as.numeric(2*wt*z)
}

## Calculate null deviance for log normal regression
dev.null.lognorm <- function(y,wt,offset,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    dev0 <- NULL
    for (i in 1:length(xx)) {
        if (!delta[i]|!zz[i]) dev0 <- c(dev0,0)
        else {
            fun.wk <- function(eta) {
                (nu*(log(xx[i])-eta))^2/2+log(1-pnorm(nu*(log(zz[i])-eta)))
            }
            dev0 <- c(dev0,nlm(fun.wk,log(xx[i]),stepmax=1)$min)
        }
    }
    if (is.null(offset)) offset <- rep(0,length(xx))
    lkhd <- function(eta) {
        eta <- eta + offset
        xx.wk <- nu*(log(xx)-eta)
        zz.wk <- nu*(log(zz)-eta)
        -sum(wt*(delta*(-xx.wk^2/2-log(1-pnorm(xx.wk)))
                 +log((1-pnorm(xx.wk))/(1-pnorm(zz.wk)))))
    }
    eta <- nlm(lkhd,mean(log(xx)-offset),stepmax=1)$est + offset
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    z <- -delta*(-xx^2/2-log(1-pnorm(xx)))-log((1-pnorm(xx))/(1-pnorm(zz)))
    sum(2*wt*(z-dev0))
}


##%%%%%%%%%%  Log Logistic Family %%%%%%%%%%

## Make pseudo data for log logistic regression
mkdata.loglogis <- function(y,eta,wt,offset,nu)
{
    if (is.vector(y)) stop("gss error: missing censoring indicator")
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    if (is.null(offset)) offset <- rep(0,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    if (any(zz<0)|any(zz>=xx))
        stop("gss error: inconsistent life time data")
    if (nu[[2]]) {
        lkhd <- function(log.nu) {
            nu <- exp(log.nu)
            xx.wk <- nu*(log(xx)-eta)
            zz.wk <- nu*(log(zz)-eta)
            -sum(wt*(delta*(xx.wk-log(1+exp(xx.wk))+log.nu)
                 -log((1+exp(xx.wk))/(1+exp(zz.wk)))))
        }
        if (is.null(nu[[1]])) nu[[1]] <- 1
        nu[[1]] <- exp(nlm(lkhd,log(nu[[1]]),stepmax=.5)$est)
    }
    xx <- 1/(1+exp(nu[[1]]*(log(xx)-eta)))
    zz <- 1/(1+exp(nu[[1]]*(log(zz)-eta)))
    u <- nu[[1]]*(delta*xx-(zz-xx))
    w <- nu[[1]]^2/2*(zz^2-xx^2)
    w <- pmax(w,1e-6)
    ywk <- eta-u/w-offset
    wt <- w*wt
    list(ywk=ywk,wt=wt,nu=nu)
}

## Calculate deviance residuals for log logistic regression
dev.resid.loglogis <- function(y,eta,wt,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    dev0 <- NULL
    for (i in 1:length(xx)) {
        if (!delta[i]) dev0 <- c(dev0,0)
        else {
            if (!zz[i]) dev0 <- c(dev0,2*log(2))
            else {
                if ((xx[i]/zz[i])^nu<=2) dev0 <- c(dev0,nu*log(xx[i]/zz[i]))
                else dev0 <- c(dev0,2*log(2)-log(xx[i]^nu/(xx[i]^nu-zz[i]^nu)))
            }
        }
    }
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    z <- -delta*(xx-log(1+exp(xx)))+log((1+exp(xx))/(1+exp(zz)))
    as.numeric(2*wt*(z-dev0))
}

dev0.resid.loglogis <- function(y,eta,wt,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    z <- -delta*(xx-log(1+exp(xx)))+log((1+exp(xx))/(1+exp(zz)))
    as.numeric(2*wt*z)
}

## Calculate null deviance for log logistic regression
dev.null.loglogis <- function(y,wt,offset,nu)
{
    if (is.null(wt)) wt <- rep(1,dim(y)[1])
    xx <- y[,1]
    delta <- as.logical(y[,2])
    if (dim(y)[2]>=3) zz <- y[,3]
    else zz <- rep(0,length(xx))
    dev0 <- NULL
    for (i in 1:length(xx)) {
        if (!delta[i]) dev0 <- c(dev0,0)
        else {
            if (!zz[i]) dev0 <- c(dev0,2*log(2))
            else {
                if ((xx[i]/zz[i])^nu<=2) dev0 <- c(dev0,nu*log(xx[i]/zz[i]))
                else dev0 <- c(dev0,2*log(2)-log(xx[i]^nu/(xx[i]^nu-zz[i]^nu)))
            }
        }
    }
    if (is.null(offset)) offset <- rep(0,length(xx))
    lkhd <- function(eta) {
        eta <- eta + offset
        xx.wk <- nu*(log(xx)-eta)
        zz.wk <- nu*(log(zz)-eta)
        -sum(wt*(delta*(xx.wk-log(1+exp(xx.wk)))
                 -log((1+exp(xx.wk))/(1+exp(zz.wk)))))
    }
    eta <- nlm(lkhd,mean(log(xx)-offset))$est + offset
    xx <- nu*(log(xx)-eta)
    zz <- nu*(log(zz)-eta)
    z <- -delta*(xx-log(1+exp(xx)))+log((1+exp(xx))/(1+exp(zz)))
    sum(2*wt*(z-dev0))
}
## Obtain fitted values from ssanova objects
fitted.ssanova <- function(object,...)
{
    mf <- object$mf
    if (!is.null(object$random)) mf$random <- I(object$random$z)
    predict(object,mf)
}

## Obtain residuals from ssanova objects
residuals.ssanova <- function(object,...)
{
    y <- model.response(object$mf,"numeric")
    as.numeric(y-fitted.ssanova(object))
}

## Obtain fitted values in working scale from gssanova objects
fitted.gssanova <- function(object,...)
{
    as.numeric(object$eta)
}

## Obtain residuals from gssanova objects
residuals.gssanova <- function(object,type="working",...)
{
    y <- model.response(object$mf,"numeric")
    wt <- model.weights(object$mf)
    offset <- NULL
    if ((object$family=="nbinomial")&(!is.null(object$nu))) y <- cbind(y,object$nu)
    dat <- switch(object$family,
                  binomial=mkdata.binomial(y,object$eta,wt,offset),
                  nbinomial=mkdata.nbinomial(y,object$eta,wt,offset,nu),
                  poisson=mkdata.poisson(y,object$eta,wt,offset),
                  inverse.gaussian=mkdata.inverse.gaussian(y,object$eta,wt,offset),
                  Gamma=mkdata.Gamma(y,object$eta,wt,offset),
                  weibull=mkdata.weibull(y,object$eta,wt,offset,list(object$nu,FALSE)),
                  lognorm=mkdata.lognorm(y,object$eta,wt,offset,list(object$nu,FALSE)),
                  loglogis=mkdata.loglogis(y,object$eta,wt,offset,list(object$nu,FALSE)))
    res <- as.numeric(dat$ywk - object$eta)
    if (!is.na(charmatch(type,"deviance"))) {
        dev.resid <- switch(object$family,
                            binomial=dev.resid.binomial(y,object$eta,wt),
                            nbinomial=dev.resid.nbinomial(y,object$eta,wt),
                            poisson=dev.resid.poisson(y,object$eta,wt),
                            inverse.gaussian=dev.resid.inverse.gaussian(y,object$eta,wt),
                            Gamma=dev.resid.Gamma(y,object$eta,wt),
                            weibull=dev.resid.weibull(y,object$eta,wt,object$nu),
                            lognorm=dev.resid.lognorm(y,object$eta,wt,object$nu),
                            loglogis=dev.resid.loglogis(y,object$eta,wt,object$nu))
        res <- sqrt(dev.resid)*sign(res)
    }
    res
}
gauss.quad <- ## Generate Gauss-Legendre quadrature
function(size,interval) {
    if (interval[1]>=interval[2])
        warning("gss warning in gauss.quad: interval limits swapped")
    z <- .Fortran("gaussq",
                  as.integer(1),
                  as.integer(size),
                  as.double(0), as.double(0),
                  as.integer(0),
                  as.double(c(-1,1)), double(size),
                  t=double(size), w=double(size),
                  PACKAGE="gss")
    mn <- min(interval[1:2])
    range <- abs(interval[1]-interval[2])
    pt <- mn+range*(z$t+1)/2
    wt <- range*z$w/2
    list(pt=pt,wt=wt)
}
## Fit Single Smoothing Parameter REGression by Performance-Oriented Iteration
sspregpoi <- function(family,s,q,y,wt,offset,method="u",
                      varht=1,nu,prec=1e-7,maxiter=30)
{
    ## Check inputs
    if (is.vector(s)) s <- as.matrix(s)
    if (!(is.matrix(s)&is.matrix(q)&is.character(method))) {
        stop("gss error in sspregpoi: inputs are of wrong types")
    }
    nobs <- dim(s)[1]
    nnull <- dim(s)[2]
    if (!((dim(s)[1]==nobs)&(dim(q)[1]==nobs)&(dim(q)[2]==nobs)
          &(nobs>=nnull)&(nnull>0))) {
        stop("gss error in sspregpoi: inputs have wrong dimensions")
    }
    ## Set method for smoothing parameter selection
    code <- (1:3)[c("v","m","u")==method]
    if (!length(code)) {
        stop("gss error: unsupported method for smoothing parameter selection")
    }
    eta <- rep(0,nobs)
    nla0 <- log10(mean(abs(diag(q))))
    limnla <- nla0+c(-.5,.5)
    iter <- 0
    if (family=="nbinomial") nu <- NULL
    else nu <- list(nu,is.null(nu))
    repeat {
        iter <- iter+1
        dat <- switch(family,
                      binomial=mkdata.binomial(y,eta,wt,offset),
                      nbinomial=mkdata.nbinomial(y,eta,wt,offset,nu),
                      poisson=mkdata.poisson(y,eta,wt,offset),
                      inverse.gaussian=mkdata.inverse.gaussian(y,eta,wt,offset),
                      Gamma=mkdata.Gamma(y,eta,wt,offset),
                      weibull=mkdata.weibull(y,eta,wt,offset,nu),
                      lognorm=mkdata.lognorm(y,eta,wt,offset,nu),
                      loglogis=mkdata.loglogis(y,eta,wt,offset,nu))
        nu <- dat$nu
        w <- as.vector(sqrt(dat$wt))
        ywk <- w*dat$ywk
        swk <- w*s
        qwk <- w*t(w*q)
        ## Call RKPACK driver DSIDR
        z <- .Fortran("dsidr0",
                      as.integer(code),
                      swk=as.double(swk), as.integer(nobs),
                      as.integer(nobs), as.integer(nnull),
                      as.double(ywk),
                      qwk=as.double(qwk), as.integer(nobs),
                      as.double(0), as.integer(-1), as.double(limnla),
                      nlambda=double(1), score=double(1), varht=as.double(varht),
                      c=double(nobs), d=double(nnull),
                      qraux=double(nnull), jpvt=integer(nnull),
                      double(3*nobs),
                      info=integer(1),PACKAGE="gss")
        ## Check info for error
        if (info<-z$info) {               
            if (info>0)
                stop("gss error in sspregpoi: matrix s is rank deficient")
            if (info==-2)
                stop("gss error in sspregpoi: matrix q is indefinite")
            if (info==-1)
                stop("gss error in sspregpoi: input data have wrong dimensions")
            if (info==-3)
                stop("gss error in sspregpoi: unknown method for smoothing parameter selection.")
        }
        eta.new <- (ywk-10^z$nlambda*z$c)/w
        if (!is.null(offset)) eta.new <- eta.new + offset
        disc <- sum(dat$wt*((eta-eta.new)/(1+abs(eta)))^2)/sum(dat$wt)
        limnla <- pmax(z$nlambda+c(-.5,.5),nla0-5)
        if (disc<prec) break
        if (iter>=maxiter) {
            warning("gss warning in gssanova: performance-oriented iteration fails to converge")
            break
        }
        eta <- eta.new
    }
    ## Return the fit
    if (is.list(nu)) nu <- nu[[1]]
    c(list(method=method,theta=0,w=as.vector(dat$wt),
           eta=as.vector(eta),iter=iter,nu=nu),
      z[c("c","d","nlambda","score","varht","swk","qraux","jpvt","qwk")])
}

## Fit Multiple Smoothing Parameter REGression by Performance-Oriented Iteration
mspregpoi <- function(family,s,q,y,wt,offset,method="u",
                      varht=1,nu,prec=1e-7,maxiter=30)
{
    ## Check inputs
    if (is.vector(s)) s <- as.matrix(s)
    if (!(is.matrix(s)&is.array(q)&(length(dim(q))==3)
          &is.character(method))) {
        stop("gss error in mspregpoi: inputs are of wrong types")
    }
    nobs <- dim(s)[1]
    nnull <- dim(s)[2]
    nq <- dim(q)[3]
    if (!((dim(s)[1]==nobs)&(dim(q)[1]==nobs)&(dim(q)[2]==nobs)
          &(nobs>=nnull)&(nnull>0))) {
        stop("gss error in sspregpoi: inputs have wrong dimensions")
    }
    ## Set method for smoothing parameter selection
    code <- (1:3)[c("v","m","u")==method]
    if (!length(code)) {
        stop("gss error: unsupported method for smoothing parameter selection")
    }
    eta <- rep(0,nobs)
    init <- 0
    theta <- rep(0,nq)
    iter <- 0
    if (family=="nbinomial") nu <- NULL
    else nu <- list(nu,is.null(nu))
    qwk <- array(0,c(nobs,nobs,nq))
    repeat {
        iter <- iter+1
        dat <- switch(family,
                      binomial=mkdata.binomial(y,eta,wt,offset),
                      nbinomial=mkdata.nbinomial(y,eta,wt,offset,nu),
                      poisson=mkdata.poisson(y,eta,wt,offset),
                      inverse.gaussian=mkdata.inverse.gaussian(y,eta,wt,offset),
                      Gamma=mkdata.Gamma(y,eta,wt,offset),
                      weibull=mkdata.weibull(y,eta,wt,offset,nu),
                      lognorm=mkdata.lognorm(y,eta,wt,offset,nu),
                      loglogis=mkdata.loglogis(y,eta,wt,offset,nu))
        nu <- dat$nu
        w <- as.vector(sqrt(dat$wt))
        ywk <- w*dat$ywk
        swk <- w*s
        for (i in 1:nq) qwk[,,i] <- w*t(w*q[,,i])
        ## Call RKPACK driver DMUDR
        z <- .Fortran("dmudr0",
                      as.integer(code),
                      as.double(swk),   # s
                      as.integer(nobs), as.integer(nobs), as.integer(nnull),
                      as.double(qwk),   # q
                      as.integer(nobs), as.integer(nobs), as.integer(nq),
                      as.double(ywk),   # y
                      as.double(0), as.integer(init),
                      as.double(prec), as.integer(maxiter),
                      theta=as.double(theta), nlambda=double(1),
                      score=double(1), varht=as.double(varht),
                      c=double(nobs), d=double(nnull),
                      double(nobs*nobs*(nq+2)),
                      info=integer(1),PACKAGE="gss")[c("theta","nlambda","c","info")]
        ## Check info for error
        if (info<-z$info) {               
            if (info>0)
                stop("gss error in mspreg: matrix s is rank deficient")
            if (info==-2)
                stop("gss error in mspreg: matrix q is indefinite")
            if (info==-1)
                stop("gss error in mspreg: input data have wrong dimensions")
            if (info==-3)
                stop("gss error in mspreg: unknown method for smoothing parameter selection.")
            if (info==-4)
                stop("gss error in mspreg: iteration fails to converge, try to increase maxiter")
            if (info==-5)
                stop("gss error in mspreg: iteration fails to find a reasonable descent direction")
        }
        eta.new <- (ywk-10^z$nlambda*z$c)/w
        if (!is.null(offset)) eta.new <- eta.new + offset
        disc <- sum(dat$wt*((eta-eta.new)/(1+abs(eta)))^2)/sum(dat$wt)
        if (disc<prec) break
        if (iter>=maxiter) {
            warning("gss warning in gssanova: performance-oriented iteration fails to converge")
            break
        }
        init <- 1
        theta <- z$theta
        eta <- eta.new
    }
    qqwk <- 10^z$theta[1]*qwk[,,1]
    for (i in 2:nq) qqwk <- qqwk + 10^z$theta[i]*qwk[,,i]
    ## Call RKPACK driver DSIDR
    z <- .Fortran("dsidr0",
                  as.integer(code),
                  swk=as.double(swk), as.integer(nobs),
                  as.integer(nobs), as.integer(nnull),
                  as.double(ywk),
                  qwk=as.double(qqwk), as.integer(nobs),
                  as.double(0), as.integer(0), double(2),
                  nlambda=double(1), score=double(1), varht=as.double(varht),
                  c=double(nobs), d=double(nnull),
                  qraux=double(nnull), jpvt=integer(nnull),
                  double(3*nobs),
                  info=integer(1),PACKAGE="gss")
    ## Check info for error
    if (info<-z$info) {               
        if (info>0)
            stop("gss error in sspregpoi: matrix s is rank deficient")
        if (info==-2)
            stop("gss error in sspregpoi: matrix q is indefinite")
        if (info==-1)
            stop("gss error in sspregpoi: input data have wrong dimensions")
        if (info==-3)
            stop("gss error in sspregpoi: unknown method for smoothing parameter selection.")
    }
    ## Return the fit
    if (is.list(nu)) nu <- nu[[1]]
    c(list(method=method,theta=theta,w=as.vector(dat$wt),
           eta=as.vector(eta),iter=iter,nu=nu),
      z[c("c","d","nlambda","score","varht","swk","qraux","jpvt","qwk")])
}
## Fit Single Smoothing Parameter Non-Gaussian REGression
sspngreg1 <- function(family,s,r,q,y,wt,offset,alpha,nu,random)
{
    nobs <- nrow(r)
    nxi <- ncol(r)
    if (!is.null(s)) {
        if (is.vector(s)) nnull <- 1
        else nnull <- ncol(s)
    }
    else nnull <- 0
    if (!is.null(random)) nz <- ncol(as.matrix(random$z))
    else nz <- 0
    nxiz <- nxi + nz
    nn <- nxiz + nnull
    ## cv function
    cv <- function(lambda) {
        if (nu[[2]]) {
            la.wk <- lambda[-2]
            nu.wk <- list(exp(lambda[2]),FALSE)
        }
        else {
            la.wk <- lambda
            nu.wk <- nu
        }
        if (is.null(random)) q.wk <- 10^(la.wk+theta)*q
        else {
            q.wk <- matrix(0,nxiz,nxiz)
            q.wk[1:nxi,1:nxi] <- 10^(la.wk[1]+theta)*q
            q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
                10^(2*ran.scal)*random$sigma$fun(la.wk[-1],random$sigma$env)
        }
        alpha.wk <- max(0,log.la0-la.wk[1]-5)*(3-alpha) + alpha
        alpha.wk <- min(alpha.wk,3)
        z <- ngreg1(dc,family,cbind(s,10^theta*r),q.wk,y,wt,offset,nu.wk,alpha.wk)
        assign("dc",z$dc,inherit=TRUE)
        assign("fit",z[c(1:3,5:10)],inherit=TRUE)
        z$score
    }
    cv.wk <- function(lambda) cv.scale*cv(lambda)+cv.shift
    ## initialization
    tmp <- sum(r^2)
    if (is.null(s)) theta <- 0
    else theta <- log10(sum(s^2)/nnull/tmp*nxi) / 2
    log.la0 <- log10(tmp/sum(diag(q))) + theta
    if (!is.null(random)) {
        ran.scal <- theta - log10(sum(random$z^2)/nz/tmp*nxi) / 2
        r <- cbind(r,10^(ran.scal-theta)*random$z)
    }
    else ran.scal <- NULL
    if (nu[[2]]&is.null(nu[[1]])) {
        eta <- rep(0,nobs)
        wk <- switch(family,
                      nbinomial=mkdata.nbinomial(y,eta,wt,offset,NULL),
                      weibull=mkdata.weibull(y,eta,wt,offset,nu),
                      lognorm=mkdata.lognorm(y,eta,wt,offset,nu),
                      loglogis=mkdata.loglogis(y,eta,wt,offset,nu))
        nu[[1]] <- wk$nu[[1]]
    }
    ## lambda search
    dc <- rep(0,nn)
    fit <- NULL
    la <- log.la0
    if (nu[[2]]) la <- c(la, log(nu[[1]]))
    if (!is.null(random)) la <- c(la,random$init)
    if (length(la)-1) {
        counter <- 0
        ## scale and shift cv
        tmp <- abs(cv(la))
        cv.scale <- 1
        cv.shift <- 0
        if (tmp<1&tmp>10^(-4)) {
            cv.scale <- 10/tmp
            cv.shift <- 0
        }
        if (tmp<10^(-4)) {
            cv.scale <- 10^2
            cv.shift <- 10
        }
        repeat {
            zz <- nlm(cv.wk,la,stepmax=1,ndigit=7)
            if (zz$code<=3) break
            la <- zz$est
            counter <- counter + 1
            if (counter>=5) {
                warning("gss warning in ssanova1: iteration for model selection fails to converge")
                break
            }
        }
    }
    else {
        repeat {
            mn <- la-1
            mx <- la+1
            zz <- nlm0(cv,c(mn,mx))
            if (min(zz$est-mn,mx-zz$est)>=1e-3) break
            else la <- zz$est
        }
    }
    ## return
    jk <- cv(zz$est)
    if (nu[[2]]) {
        nu.wk <- exp(zz$est[2])
        zz$est <- zz$est[-2]
    }
    else nu.wk <- NULL
    if (is.null(random)) q.wk <- 10^theta*q
    else {
        q.wk <- matrix(0,nxiz,nxiz)
        q.wk[1:nxi,1:nxi] <- 10^theta*q
        q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
            10^(2*ran.scal-zz$est[1])*random$sigma$fun(zz$est[-1],random$sigma$env)
    }
    zzz <- eigen(q.wk,TRUE)
    rkq <- min(fit$rkv-nnull,sum(zzz$val/zzz$val[1]>sqrt(.Machine$double.eps)))
    val <- zzz$val[1:rkq]
    vec <- zzz$vec[,1:rkq,drop=FALSE]
    qinv <- vec%*%diag(1/val,rkq)%*%t(vec)
    se.aux <- t(fit$w*cbind(s,10^theta*r))%*%(10^theta*r)%*%qinv
    c <- fit$dc[nnull+(1:nxi)]
    if (nnull) d <- fit$dc[1:nnull]
    else d <- NULL
    if (nz) b <- 10^(ran.scal)*fit$dc[nnull+nxi+(1:nz)]
    else b <- NULL
    c(list(theta=theta,ran.scal=ran.scal,c=c,d=d,b=b,nlambda=zz$est[1],
           zeta=zz$est[-1],nu=nu.wk),fit[-1],list(qinv=qinv,se.aux=se.aux))
}

## Fit Multiple Smoothing Parameter Non-Gaussian REGression
mspngreg1 <- function(family,s,r,q,y,wt,offset,alpha,nu,random)
{
    nobs <- nrow(r)
    nxi <- ncol(r)
    if (!is.null(s)) {
        if (is.vector(s)) nnull <- 1
        else nnull <- ncol(s)
    }
    else nnull <- 0
    if (!is.null(random)) nz <-ncol(as.matrix(random$z))
    else nz <- 0
    nxiz <- nxi + nz
    nn <- nxiz + nnull
    nq <- dim(q)[3]
    ## cv function
    cv <- function(theta) {
        if (nu[[2]]) {
            the.wk <- theta[-(nq+1)]
            nu.wk <- list(exp(theta[nq+1]),FALSE)
        }
        else {
            the.wk <- theta
            nu.wk <- nu
        }
        r.wk <- qq.wk <- 0
        for (i in 1:nq) {
            r.wk <- r.wk + 10^the.wk[i]*r[,,i]
            qq.wk <- qq.wk + 10^the.wk[i]*q[,,i]
        }
        if (is.null(random)) q.wk <- 10^nlambda*qq.wk
        else {
            r.wk <- cbind(r.wk,10^(ran.scal)*random$z)
            q.wk <- matrix(0,nxiz,nxiz)
            q.wk[1:nxi,1:nxi] <- 10^nlambda*qq.wk
            q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
                10^(2*ran.scal)*random$sigma$fun(the.wk[-(1:nq)],random$sigma$env)
        }
        alpha.wk <- max(0,the.wk[1:nq]-log.th0-5)*(3-alpha) + alpha
        alpha.wk <- min(alpha.wk,3)
        z <- ngreg1(dc,family,cbind(s,r.wk),q.wk,y,wt,offset,nu.wk,alpha.wk)
        assign("dc",z$dc,inherit=TRUE)
        assign("fit",z[c(1:3,5:10)],inherit=TRUE)
        z$score
    }
    cv.wk <- function(theta) cv.scale*cv(theta)+cv.shift
    ## initialization
    theta <- -log10(apply(q,3,function(x)sum(diag(x))))
    r.wk <- q.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
    }
    ## theta adjustment
    z <- sspngreg1(family,s,r.wk,q.wk,y,wt,offset,alpha,nu,random)
    if (nu[[2]]) nu[[1]] <- z$nu
    theta <- theta + z$theta
    r.wk <- q.wk <- 0
    for (i in 1:nq) {
        theta[i] <- 2*theta[i] + log10(t(z$c)%*%q[,,i]%*%z$c)
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
    }
    log.la0 <- log10(sum(r.wk^2)/sum(diag(q.wk)))
    log.th0 <- theta-log.la0
    ## lambda search
    z <- sspngreg1(family,s,r.wk,q.wk,y,wt,offset,alpha,nu,random)
    if (nu[[2]]) nu[[1]] <- z$nu
    nlambda <- z$nlambda
    log.th0 <- log.th0 + z$lambda
    theta <- theta + z$theta
    if (!is.null(random)) ran.scal <- z$ran.scal
    ## theta search
    dc <- rep(0,nn)
    fit <- NULL
    if (nu[[2]]) theta <- c(theta, log(nu[[1]]))
    if (!is.null(random)) theta <- c(theta,z$zeta)
    counter <- 0
    tmp <- abs(cv(theta))
    cv.scale <- 1
    cv.shift <- 0
    if (tmp<1&tmp>10^(-4)) {
        cv.scale <- 10/tmp
        cv.shift <- 0
    }
    if (tmp<10^(-4)) {
        cv.scale <- 10^2
        cv.shift <- 10
    }
    repeat {
        zz <- nlm(cv.wk,theta,stepmax=1,ndigit=7)
        if (zz$code<=3)  break
        theta <- zz$est        
        counter <- counter + 1
        if (counter>=5) {
            warning("gss warning in gssanova1: iteration for model selection fails to converge")
            break
        }
    }
    ## return
    jk <- cv(zz$est)
    if (nu[[2]]) {
        nu.wk <- exp(zz$est[nq+1])
        zz$est <- zz$est[-(nq+1)]
    }
    else nu.wk <- NULL
    r.wk <- qq.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^zz$est[i]*r[,,i]
        qq.wk <- qq.wk + 10^zz$est[i]*q[,,i]
    }
    if (is.null(random)) q.wk <- qq.wk
    else {
        r.wk <- cbind(r.wk,10^(ran.scal)*random$z)
        q.wk <- matrix(0,nxiz,nxiz)
        q.wk[1:nxi,1:nxi] <- qq.wk
        q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
            10^(2*ran.scal-nlambda)*random$sigma$fun(zz$est[-(1:nq)],random$sigma$env)
    }
    zzz <- eigen(q.wk,TRUE)
    rkq <- min(fit$rkv-nnull,sum(zzz$val/zzz$val[1]>sqrt(.Machine$double.eps)))
    val <- zzz$val[1:rkq]
    vec <- zzz$vec[,1:rkq,drop=FALSE]
    qinv <- vec%*%diag(1/val,rkq)%*%t(vec)
    se.aux <- t(fit$w*cbind(s,r.wk))%*%r.wk%*%qinv
    c <- fit$dc[nnull+(1:nxi)]
    if (nnull) d <- fit$dc[1:nnull]
    else d <- NULL
    if (nz) b <- 10^(ran.scal)*fit$dc[nnull+nxi+(1:nz)]
    else b <- NULL
    c(list(theta=zz$est[1:nq],c=c,d=d,b=b,nlambda=nlambda,zeta=zz$est[-(1:nq)],nu=nu.wk),
      fit[-1],list(qinv=qinv,se.aux=se.aux))
}

## Non-Gaussian regression with fixed smoothing parameters
ngreg1 <- function(dc,family,sr,q,y,wt,offset,nu,alpha)
{
    nobs <- nrow(sr)
    nn <- ncol(sr)
    nxi <- nrow(q)
    nnull <- nn - nxi
    ## initialization
    cc <- dc[nnull+(1:nxi)]
    eta <- sr%*%dc
    if (!is.null(offset)) eta <- eta + offset
    if ((family=="nbinomial")&is.vector(y)) y <- cbind(y,nu[[1]])
    iter <- 0
    flag <- 0
    dev <- switch(family,
                  binomial=dev.resid.binomial(y,eta,wt),
                  nbinomial=dev.resid.nbinomial(y,eta,wt),
                  poisson=dev.resid.poisson(y,eta,wt),
                  Gamma=dev.resid.Gamma(y,eta,wt),
                  weibull=dev.resid.weibull(y,eta,wt,nu[[1]]),
                  lognorm=dev0.resid.lognorm(y,eta,wt,nu[[1]]),
                  loglogis=dev0.resid.loglogis(y,eta,wt,nu[[1]]))
    dev <- sum(dev) + t(cc)%*%q%*%cc
    ## Newton iteration
    repeat {
        iter <- iter+1
        dat <- switch(family,
                      binomial=mkdata.binomial(y,eta,wt,offset),
                      nbinomial=mkdata.nbinomial(y,eta,wt,offset,nu),
                      poisson=mkdata.poisson(y,eta,wt,offset),
                      Gamma=mkdata.Gamma(y,eta,wt,offset),
                      weibull=mkdata.weibull(y,eta,wt,offset,nu),
                      lognorm=mkdata.lognorm(y,eta,wt,offset,nu),
                      loglogis=mkdata.loglogis(y,eta,wt,offset,nu))
        ## weighted least squares fit
        w <- as.vector(sqrt(dat$wt))
        ywk <- w*dat$ywk
        srwk <- w*sr
        if (!is.finite(sum(w,ywk,srwk))) {
            if (flag) stop("gss error in gssanova1: Newton iteration diverges")
            eta <- rep(0,nobs)
            iter <- 0
            flag <- 1
            next
        }
        z <- .Fortran("reg",
                      as.double(srwk), as.integer(nobs), as.integer(nnull),
                      as.double(q), as.integer(nxi), as.double(ywk),
                      as.integer(4),
                      double(1), double(1), double(1), dc=double(nn),
                      as.double(.Machine$double.eps),
                      double(nn*nn), double(nn),
                      as.integer(c(rep(1,nnull),rep(0,nxi))),
                      double(max(nobs,nn)), integer(1), integer(1),
                      PACKAGE="gss")["dc"]
        dc.diff <- z$dc-dc
        adj <- 0
        repeat {
            dc.new <- dc + dc.diff
            cc <- dc.new[nnull+(1:nxi)]
            eta.new <- sr%*%dc.new
            if (!is.null(offset)) eta.new <- eta.new + offset
            dev.new <- switch(family,
                              binomial=dev.resid.binomial(y,eta.new,wt),
                              nbinomial=dev.resid.nbinomial(y,eta.new,wt),
                              poisson=dev.resid.poisson(y,eta.new,wt),
                              Gamma=dev.resid.Gamma(y,eta.new,wt),
                              weibull=dev.resid.weibull(y,eta.new,wt,nu[[1]]),
                              lognorm=dev0.resid.lognorm(y,eta.new,wt,nu[[1]]),
                              loglogis=dev0.resid.loglogis(y,eta.new,wt,nu[[1]]))
            dev.new <- sum(dev.new) + t(cc)%*%q%*%cc
            if (!is.finite(dev.new)) dev.new <- Inf
            if (dev.new-dev<(1+abs(dev))*1e-1) break
            adj <- 1
            dc.diff <- dc.diff/2
        }
        disc <- sum(dat$wt*((eta-eta.new)/(1+abs(eta)))^2)/sum(dat$wt)
        if (!is.finite(disc)) {
            if (flag) stop("gss error in gssanova1: Newton iteration diverges")
            eta <- rep(0,nobs)
            iter <- 0
            flag <- 1
            next
        }
        dc <- dc.new
        eta <- eta.new
        dev <- dev.new
        if (adj) next
        if (disc<1e-7) break
        if (iter<=30) next
        if (!flag) {
            eta <- rep(0,nobs)
            iter <- 0
            flag <- 1
        }
        else {
            warning("gss warning in gssanova1: Newton iteration fails to converge")
            break
        }
    }
    ## calculate cv
    dat <- switch(family,
                  binomial=mkdata.binomial(y,eta,wt,offset),
                  nbinomial=mkdata.nbinomial(y,eta,wt,offset,nu),
                  poisson=mkdata.poisson(y,eta,wt,offset),
                  Gamma=mkdata.Gamma(y,eta,wt,offset),
                  weibull=mkdata.weibull(y,eta,wt,offset,nu),
                  lognorm=mkdata.lognorm(y,eta,wt,offset,nu),
                  loglogis=mkdata.loglogis(y,eta,wt,offset,nu))
    ## weighted least squares fit
    w <- as.vector(sqrt(dat$wt))
    ywk <- w*dat$ywk
    srwk <- w*sr
    z <- .Fortran("reg",
                  as.double(srwk), as.integer(nobs), as.integer(nnull),
                  as.double(q), as.integer(nxi), as.double(ywk),
                  as.integer(5),
                  double(1), double(1), double(1), dc=double(nn),
                  as.double(.Machine$double.eps),
                  chol=double(nn*nn), double(nn),
                  jpvt=as.integer(c(rep(1,nnull),rep(0,nxi))),
                  hat=double(max(nobs+1,nn)), rkv=integer(1), integer(1),
                  PACKAGE="gss")[c("dc","chol","jpvt","hat","rkv")]
    cv <- switch(family,
                 binomial=cv.binomial(y,eta,wt,z$hat[1:nobs],alpha),
                 poisson=cv.poisson(y,eta,wt,z$hat[1:nobs],alpha,sr,q),
                 Gamma=cv.Gamma(y,eta,wt,z$hat[1:nobs],z$hat[nobs+1],alpha),
                 nbinomial=cv.nbinomial(y,eta,wt,z$hat[1:nobs],alpha),
                 weibull=cv.weibull(y,eta,wt,z$hat[1:nobs],nu[[1]],alpha),
                 lognorm=cv.lognorm(y,eta,wt,z$hat[1:nobs],nu[[1]],alpha),
                 loglogis=cv.loglogis(y,eta,wt,z$hat[1:nobs],nu[[1]],alpha))
    c(z,cv,list(eta=eta))
}
## Fit gssanova model
gssanova <- function(formula,family,type="cubic",data=list(),
                     weights,subset,offset,na.action=na.omit,
                     partial=NULL,method=NULL,varht=1,nu=NULL,
                     prec=1e-7,maxiter=30,ext=.05,order=2)
{
    ## Obtain model frame and model terms
    mf <- match.call()
    mf$family <- mf$type <- mf$partial <- NULL
    mf$method <- mf$varht <- mf$nu <- NULL
    mf$prec <- mf$maxiter <- mf$ext <- mf$order <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf,sys.frame(sys.parent()))
    if (type=="cubic") term <- mkterm.cubic(mf,ext)
    if (type=="linear") term <- mkterm.linear(mf,ext)
    if (type=="tp") term <- mkterm.tp(mf,order,mf,1)
    if (is.null(term)) stop("gss error in gssanova: unknown type")
    ## Specify default method
    if (is.null(method)) {
        method <- switch(family,
                         binomial="u",
                         nbinomial="u",
                         poisson="u",
                         inverse.gaussian="v",
                         Gamma="v",
                         weibull="u",
                         lognorm="u",
                         loglogis="u")
    }
    ## Generate s, q, and y
    nobs <- dim(mf)[1]
    s <- q <- NULL
    nq <- 0
    for (label in term$labels) {
        if (label=="1") {
            s <- cbind(s,rep(1,len=nobs))
            next
        }
        x <- mf[,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            phi <- term[[label]]$phi
            for (i in 1:nphi)
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
        }
        if (nrk) {
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq+1
                q <- array(c(q,rk$fun(x,x,nu=i,env=rk$env,out=TRUE)),c(nobs,nobs,nq))
            }
        }
    }
    ## Add the partial term
    if (!is.null(partial)) {
        if (is.vector(partial)) partial <- as.matrix(partial)
        if (dim(partial)[1]!=dim(mf)[1])
            stop("gss error in gssanova: partial data are of wrong size")
        term$labels <- c(term$labels,"partial")
        term$partial <- list(nphi=dim(partial)[2],nrk=0,
                             iphi=ifelse(is.null(s),0,dim(s)[2])+1)
        s <- cbind(s,partial)
        mf$partial <- partial
    }
    if (qr(s)$rank<dim(s)[2])
        stop("gss error in gssanova: fixed effects are linearly dependent")
    y <- model.response(mf,"numeric")
    wt <- model.weights(mf)
    offset <- model.offset(mf)
    if (!is.null(offset)) {
        term$labels <- c(term$labels,"offset")
        term$offset <- list(nphi=0,nrk=0)
    }
    if (!nq) stop("gss error in gssanova: use glm for models with only fixed effects")
    ## Fit the model
    if (nq==1) {
        q <- q[,,1]
        z <- sspregpoi(family,s,q,y,wt,offset,method,varht,nu,prec,maxiter)
    }
    else z <- mspregpoi(family,s,q,y,wt,offset,method,varht,nu,prec,maxiter)
    ## Brief description of model terms
    desc <- NULL
    for (label in term$labels)
        desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
    desc <- rbind(desc,apply(desc,2,sum))
    rownames(desc) <- c(term$labels,"total")
    colnames(desc) <- c("Unpenalized","Penalized")
    ## Return the results
    obj <- c(list(call=match.call(),family=family,mf=mf,terms=term,desc=desc),z)
    class(obj) <- c("gssanova","ssanova")
    obj
}
## Fit gssanova model
gssanova1 <- function(formula,family,type="cubic",data=list(),
                      weights,subset,offset,na.action=na.omit,
                      partial=NULL,alpha=NULL,nu=NULL,
                      id.basis=NULL,nbasis=NULL,seed=NULL,random=NULL,
                      ext=.05,order=2)
{
    if (!(family%in%c("binomial","poisson","Gamma","nbinomial","weibull","lognorm","loglogis")))
        stop("gss error in gssanova1: family not implemented")
    if (is.null(alpha)) {
        alpha <- 1.4
        if (family=="binomial") alpha <- 1
    }
    ## Obtain model frame and model terms
    mf <- match.call()
    mf$family <- mf$type <- mf$partial <- NULL
    mf$method <- mf$varht <- mf$nu <- NULL
    mf$alpha <- mf$id.basis <- mf$nbasis <- mf$seed <- NULL
    mf$random <- mf$ext <- mf$order <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf,sys.frame(sys.parent()))
    wt <- model.weights(mf)
    ## Generate sub-basis
    nobs <- dim(mf)[1]
    if (is.null(id.basis)) {
        if (is.null(nbasis))  nbasis <- max(30,ceiling(10*nobs^(2/9)))
        if (nbasis>=nobs)  nbasis <- nobs
        if (!is.null(seed))  set.seed(seed)
        id.basis <- sample(nobs,nbasis,prob=wt)
    }
    else {
        if (max(id.basis)>nobs|min(id.basis)<1)
            stop("gss error in gssanova1: id.basis out of range")
        nbasis <- length(id.basis)
    }
    ## Generate terms
    if (type=="cubic") term <- mkterm.cubic(mf,ext)
    if (type=="linear") term <- mkterm.linear(mf,ext)
    if (type=="tp") term <- mkterm.tp(mf,order,mf[id.basis,],1)
    if (is.null(term)) stop("gss error in gssanova1: unknown type")
    ## Generate random
    if (!is.null(random)) {
        if (class(random)=="formula") random <- mkran(random,data)
    }
    ## Generate s, r, q, and y
    s <- r <- NULL
    nq <- 0
    for (label in term$labels) {
        if (label=="1") {
            s <- cbind(s,rep(1,len=nobs))
            next
        }
        x <- mf[,term[[label]]$vlist]
        x.basis <- mf[id.basis,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            phi <- term[[label]]$phi
            for (i in 1:nphi)
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
        }
        if (nrk) {
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq+1
                r <- array(c(r,rk$fun(x,x.basis,nu=i,env=rk$env,out=TRUE)),c(nobs,nbasis,nq))
            }
        }
    }
    if (is.null(r))
        stop("gss error in gssanova1: use glm for models with only fixed effects")
    else q <- r[id.basis,,,drop=FALSE]
    ## Add the partial term
    if (!is.null(partial)) {
        if (is.vector(partial)) partial <- as.matrix(partial)
        if (dim(partial)[1]!=dim(mf)[1])
            stop("gss error in gssanova1: partial data are of wrong size")
        term$labels <- c(term$labels,"partial")
        term$partial <- list(nphi=dim(partial)[2],nrk=0,
                             iphi=ifelse(is.null(s),0,dim(s)[2])+1)
        s <- cbind(s,partial)
        mf$partial <- partial
    }
    if (qr(s)$rank<dim(s)[2])
        stop("gss error in gssanova1: fixed effects are linearly dependent")
    ## Prepare the data
    y <- model.response(mf,"numeric")
    offset <- model.offset(mf)
    if (!is.null(offset)) {
        term$labels <- c(term$labels,"offset")
        term$offset <- list(nphi=0,nrk=0)
    }
    nu.wk <- list(NULL,FALSE)
    if ((family=="nbinomial")&is.vector(y)) nu.wk <- list(NULL,TRUE)
    if (family%in%c("weibull","lognorm","loglogis")) {
        if (is.null(nu)) nu.wk <- list(nu,TRUE)
        else nu.wk <- list(nu,FALSE)
    }
    ## Fit the model
    if (nq==1) {
        r <- r[,,1]
        q <- q[,,1]
        z <- sspngreg1(family,s,r,q,y,wt,offset,alpha,nu.wk,random)
    }
    else z <- mspngreg1(family,s,r,q,y,wt,offset,alpha,nu.wk,random)
    ## Brief description of model terms
    desc <- NULL
    for (label in term$labels)
        desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
    desc <- rbind(desc,apply(desc,2,sum))
    rownames(desc) <- c(term$labels,"total")
    colnames(desc) <- c("Unpenalized","Penalized")
    ## Return the results
    obj <- c(list(call=match.call(),family=family,mf=mf,terms=term,desc=desc,
                  alpha=alpha,id.basis=id.basis,random=random),z)
    class(obj) <- c("gssanova1","ssanova1","gssanova","ssanova")
    obj
}
hzdrate.sshzd <- ## Evaluate hazard estimate
function (object,x,se=FALSE) {
    if (class(object)!="sshzd")
        stop("gss error in hzdrate.sshzd: not a sshzd object")
    if (dim(object$mf)[2]==1&is.vector(x)) {
        x <- data.frame(x)
        colnames(x) <- colnames(object$mf)
    }
    s <- NULL
    r <- matrix(0,dim(x)[1],length(object$id.basis))
    nq <- 0
    for (label in object$terms$labels) {
        if (label=="1") {
            s <- cbind(s,rep(1,dim(x)[1]))
            next
        }
        xx <- object$mf[object$id.basis,object$terms[[label]]$vlist]
        x.new <- x[,object$terms[[label]]$vlist]
        nphi <- object$terms[[label]]$nphi
        nrk <-  object$terms[[label]]$nrk
        if (nphi) {
            phi <-  object$terms[[label]]$phi
            for (i in 1:nphi) {
                s <- cbind(s,phi$fun(x.new,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            rk <- object$terms[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq + 1
                r <- r + 10^object$theta[nq]*rk$fun(x.new,xx,nu=i,env=rk$env,out=TRUE)
            }
        }
    }
    rs <- cbind(r,s)
    if (!se) as.vector(exp(rs%*%c(object$c,object$d)))
    else {
        fit <- as.vector(exp(rs%*%c(object$c,object$d)))
        se.fit <- .Fortran("hzdaux2",
                           as.double(object$se.aux$v), as.integer(dim(rs)[2]),
                           as.integer(object$se.aux$jpvt),
                           as.double(t(rs)), as.integer(dim(rs)[1]),
                           se=double(dim(rs)[1]), PACKAGE="gss")[["se"]]
        list(fit=fit,se.fit=se.fit)
    }
}

hzdcurve.sshzd <- ## Evaluate hazard curve for plotting
function (object,time,covariates=NULL,se=FALSE) {
    tname <- object$tname
    xnames <- object$xnames
    if (class(object)!="sshzd")
        stop("gss error in hzdcurve.sshzd: not a sshzd object")
    if (length(xnames)&&(!all(xnames%in%names(covariates))))
        stop("gss error in survexp.sshzd: missing covariates")
    mn <- min(object$domain[[tname]])
    mx <- max(object$domain[[tname]])
    if ((min(time)<mn)|(max(time)>mx))
        stop("gss error in hzdcurve.sshzd: time range over the domain")
    if (length(xnames)) {
        xx <- covariates[,xnames,drop=FALSE]
        xy <- data.frame(matrix(0,length(time),length(xnames)+1))
        names(xy) <- c(tname,xnames)
        xy[,tname] <- time
    }
    else xx <- NULL
    if (!se) {
        if (is.null(xx))
            zz <- hzdrate.sshzd(object,time)
        else {
            zz <- NULL
            for (i in 1:dim(xx)[1]) {
                xy[,xnames] <- xx[rep(i,length(time)),]
                zz <- cbind(zz,hzdrate.sshzd(object,xy))
            }
            zz <- zz[,,drop=TRUE]
        }
        zz
    }
    else {
        if (is.null(xx))
            zz <- hzdrate.sshzd(object,time,TRUE)
        else {
            fit <- se.fit <- NULL
            for (i in 1:dim(xx)[1]) {
                xy[,xnames] <- xx[rep(i,length(time)),]
                wk <- hzdrate.sshzd(object,xy,TRUE)
                fit <- cbind(fit,wk$fit)
                se.fit <- cbind(se.fit,wk$se.fit)
            }
            zz <- list(fit=fit[,,drop=TRUE],se.fit=se.fit[,,drop=TRUE])
        }
        zz
    }
}

survexp.sshzd <- ## Compute expected survival
function(object,time,covariates=NULL,start=0) {
    tname <- object$tname
    xnames <- object$xnames
    ## Check inputs
    if (class(object)!="sshzd")
        stop("gss error in survexp.sshzd: not a sshzd object")
    if (length(xnames)&&(!all(xnames%in%names(covariates))))
        stop("gss error in survexp.sshzd: missing covariates")
    lmt <- cbind(start,time)
    if (any(lmt[,1]>lmt[,2]))
        stop("gss error in survexp.sshzd: start after follow-up time")
    nt <- dim(lmt)[1]
    if (is.null(covariates)) ncov <- 1
    else ncov <- dim(covariates)[1]
    if (length(xnames)&&(nt-1)&&(ncov-1)&&(nt-ncov))
        stop("gss error in survexp.sshzd: size mismatch")
    mn <- min(object$domain[[tname]])
    mx <- max(object$domain[[tname]])
    if ((min(start)<mn)|(max(time)>mx))
        stop("gss error in survexp.sshzd: time range over the domain")
    ## Calculate
    if (is.null(covariates)) {
        zz <- NULL
        for (i in 1:nt) {
            nqd <- max(20,ceiling((lmt[i,2]-lmt[i,1])/(mx-mn)*200))
            quad <- gauss.quad(nqd,lmt[i,])
            zz <- c(zz,sum(quad$wt*hzdrate.sshzd(object,quad$pt)))
        }
    }
    else {
        if (ncov>nt)
            lmt <- matrix(lmt,ncov,2,byrow=TRUE)
        if (ncov<nt)
            covariates <- covariates[rep(1,nt),,drop=FALSE]
        zz <- NULL
        for (i in 1:max(ncov,nt)) {
            nqd <- max(20,ceiling((lmt[i,2]-lmt[i,1])/(mx-mn)*200))
            quad <- gauss.quad(nqd,lmt[i,])
            wk <- covariates[rep(i,nqd),,drop=FALSE]
            wk[[tname]] <- quad$pt
            zz <- c(zz,sum(quad$wt*hzdrate.sshzd(object,wk)))
        }
    }
    exp(-zz)
}
## Make RK for linear splines
mkrk.linear <- function(range)
{
    ## Create the environment
    env <- list(min=min(range), max=max(range))
    ## Create the rk function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ##% Check the inputs
        if (!(is.vector(x)&is.vector(y))) {
            stop("gss error in rk: inputs are of wrong types")
        }
        if ((min(x,y)<env$min)|(max(x,y)>env$max)) {
            stop("gss error in rk: inputs are out of range")
        }
        ##% Scale the inputs
        x <- (x-env$min)/(env$max-env$min)
        y <- (y-env$min)/(env$max-env$min)
        ##% Return the result
        rk <- function(x,y) {
            k1 <- function(x) (x-.5)
            k2 <- function(x) ((x-.5)^2-1/12)/2
            k1(x)*k1(y)+k2(abs(x-y))
        }
        if (outer.prod) outer(x,y,rk)
        else rk(x,y)
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make RK for cubic splines
mkrk.cubic <- function(range)
{
    ## Create the environment
    env <- list(min=min(range), max=max(range))
    ## Create the rk function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ##% Check the inputs
        if (!(is.vector(x)&is.vector(y))) {
            stop("gss error in rk: inputs are of wrong types")
        }
        if ((min(x,y)<env$min)|(max(x,y)>env$max)) {
            stop("gss error in rk: inputs are out of range")
        }
        ##% Scale the inputs
        x <- (x-env$min)/(env$max-env$min)
        y <- (y-env$min)/(env$max-env$min)
        ##% Return the result
        rk <- function(x,y) {
            k2 <- function(x) ((x-.5)^2-1/12)/2
            k4 <- function(x) ((x-.5)^4-(x-.5)^2/2+7/240)/24
            k2(x)*k2(y)-k4(abs(x-y))
        }
        if (outer.prod) outer(x,y,rk)
        else rk(x,y)
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make phi function for cubic splines
mkphi.cubic <- function(range)
{
    ## Create the environment
    env <- list(min=min(range), max=max(range))
    ## Create the phi function
    fun <- function(x,env) {
        ##% Check the input
        if (!is.vector(x)) {
            stop("gss error in phi: inputs are of wrong types")
        }
        if ((min(x)<env$min)|(max(x)>env$max)) {
            stop("gss error in phi: inputs are out of range")
        }
        ##% Return the result
        (x-env$min)/(env$max-env$min)-.5
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}
## Make RK for thin-plate splines
mkrk.tp <- function(dm,order,mesh,weight=1)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    if (xor(is.vector(mesh),dm==1)
        |xor(is.matrix(mesh),dm>=2)) {
        stop("gss error in mkrk.tp: mismatched inputs")
    }
    if ((min(weight)<0)|(max(weight)<=0)) {
        stop("gss error in mkrk.tp: negative weights")
    }
    ## Set weights
    if (is.vector(mesh)) N <- length(mesh)
    else N <- dim(mesh)[1]
    weight <- rep(weight,len=N)
    weight <- sqrt(weight/sum(weight))
    ## Obtain orthonormal basis
    phi.p <- mkphi.tp.p(dm,order)
    nnull <- choose(dm+order-1,dm)
    s <- NULL
    for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env))
    s <- qr(weight*s)
    if (s$rank<nnull) {
        stop("gss error in mkrk.tp: insufficient normalizing mesh for thin-plate spline")
    }
    q <- qr.Q(s)
    r <- qr.R(s)
    ## Set Q^{T}E(|u_{i}-u_{j}|)Q
    rk.p <- mkrk.tp.p(dm,order)
    pep <- weight*t(weight*rk.p$fun(mesh,mesh,rk.p$env,out=TRUE))
    pep <- t(q)%*%pep%*%q
    ## Create the environment
    env <- list(dim=dm,order=order,weight=weight,
                phi.p=phi.p,rk.p=rk.p,q=q,r=r,mesh=mesh,pep=pep)
    ## Create the rk function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ## Check inputs
        if (env$dim==1) {
            if (!(is.vector(x)&is.vector(y))) {
                stop("gss error in rk: inputs are of wrong types")
            }
            nx <- length(x)
            ny <- length(y)
        }
        else {
            if (is.vector(x)) x <- t(as.matrix(x))
            if (env$dim!=dim(x)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
            nx <- dim(x)[1]
            if (is.vector(y)) y <- t(as.matrix(y))
            if (env$dim!=dim(y)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
            ny <- dim(y)[1]
        }
        ## Return the results
        nnull <- choose(env$dim+env$order-1,env$dim)
        if (outer.prod) {
            phix <- phiy <- NULL
            for (nu in 1:nnull) {
                phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env))
                phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env))
            }
            phix <- backsolve(env$r,phix,tr=TRUE)
            phiy <- backsolve(env$r,phiy,tr=TRUE)
            ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE)
            ex <- env$weight*ex
            ex <- t(env$q)%*%ex
            ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE)
            ey <- env$weight*ey
            ey <- t(env$q)%*%ey
            env$rk.p$fun(x,y,env$rk.p$env,out=TRUE)-t(phix)%*%ey-
                t(ex)%*%phiy+t(phix)%*%env$pep%*%phiy
        }
        else {
            N <- max(nx,ny)
            phix <- phiy <- NULL
            for (nu in 1:nnull) {
                phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env))
                phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env))
            }
            phix <- backsolve(env$r,phix,tr=TRUE)
            phix <- matrix(phix,nnull,N)
            phiy <- backsolve(env$r,phiy,tr=TRUE)
            phiy <- matrix(phiy,nnull,N)
            ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE)
            ex <- env$weight*ex
            ex <- t(env$q)%*%ex
            ex <- matrix(ex,nnull,N)
            ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE)
            ey <- env$weight*ey
            ey <- t(env$q)%*%ey
            ey <- matrix(ey,nnull,N)
            fn1 <- function(x,n) x[1:n]%*%x[n+(1:n)]
            fn2 <- function(x,pep,n) t(x[1:n])%*%pep%*%x[n+(1:n)]
            env$rk.p$fun(x,y,env$rk.p$env)-apply(rbind(phix,ey),2,fn1,nnull)-
                apply(rbind(phiy,ex),2,fn1,nnull)+
                    apply(rbind(phix,phiy),2,fn2,env$pep,nnull)
        }
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make phi function for thin-plate splines
mkphi.tp <-  function(dm,order,mesh,weight)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    if (xor(is.vector(mesh),dm==1)
        |xor(is.matrix(mesh),dm>=2)) {
        stop("gss error in mkphi.tp: mismatched inputs")
    }
    if ((min(weight)<0)|(max(weight)<=0)) {
        stop("gss error in mkphi.tp: negative weights")
    }
    ## Set weights
    if (is.vector(mesh)) N <- length(mesh)
    else N <- dim(mesh)[1]
    weight <- rep(weight,len=N)
    weight <- sqrt(weight/sum(weight))
    ## Create the environment
    phi.p <- mkphi.tp.p(dm,order)
    nnull <- choose(dm+order-1,dm)
    s <- NULL
    for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env))
    s <- qr(weight*s)
    if (s$rank<nnull) {
        stop("gss error in mkphi: insufficient normalizing mesh for thin-plate spline")
    }
    r <- qr.R(s)
    env <- list(dim=dm,order=order,phi.p=phi.p,r=r)
    ## Create the phi function
    fun <- function(x,nu,env) {
        nnull <- choose(env$dim+env$order-1,env$dim)
        phix <- NULL
        for(i in 1:nnull)
            phix <- rbind(phix,env$phi.p$fun(x,i,env$phi.p$env))
        t(backsolve(env$r,phix,tr=TRUE))[,nu]
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make pseudo RK for thin-plate splines
mkrk.tp.p <- function(dm,order)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    ## Create the environment
    if (dm%%2) {                    
        theta <- gamma(dm/2-order)/2^(2*order)/pi^(dm/2)/gamma(order)
    }
    else {
        theta <- (-1)^(dm/2+order+1)/2^(2*order-1)/pi^(dm/2)/
            gamma(order)/gamma(order-dm/2+1)
    }
    env <- list(dim=dm,order=order,theta=theta)
    ## Create the rk.p function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ## Check inputs
        if (env$dim==1) {
            if (!(is.vector(x)&is.vector(y))) {
                stop("gss error in rk: inputs are of wrong types")
            }
        }
        else {
            if (is.vector(x)) x <- t(as.matrix(x))
            if (env$dim!=dim(x)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
            if (is.vector(y)) y <- t(as.matrix(y))
            if (env$dim!=dim(y)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
        }
        ## Return the results
        if (outer.prod) {               
            if (env$dim==1) {
                fn1 <- function(x,y) abs(x-y)
                d <- outer(x,y,fn1)
            }
            else {
                fn2 <- function(x,y) sqrt(sum((x-y)^2))
                d <- NULL
                for (i in 1:dim(y)[1]) d <- cbind(d,apply(x,1,fn2,y[i,]))
            }
        }
        else {
            if (env$dim==1) d <- abs(x-y)
            else {
                N <- max(dim(x)[1],dim(y)[1])
                x <- t(matrix(t(x),env$dim,N))
                y <- t(matrix(t(y),env$dim,N))
                fn <- function(x) sqrt(sum(x^2))
                d <- apply(x-y,1,fn)
            }
        }
        power <- 2*env$order-env$dim
        switch(1+env$dim%%2,
               env$theta*d^power*log(ifelse(d>0,d,1)),
               env$theta*d^power)
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make pseudo phi function for thin-plate splines
mkphi.tp.p <- function(dm,order)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    ## Create the environment
    pol.code <- NULL
    for (i in 0:(order^dm-1)) {
        ind <- i; code <- NULL
        for (j in 1:dm) {
            code <- c(code,ind%%order)
            ind <- ind%/%order
        }
        if (sum(code)<order) pol.code <- cbind(pol.code,code)
    }
    env <- list(dim=dm,pol.code=pol.code)
    ## Create the phi function  
    fun <- function(x,nu,env) {
        if (env$dim==1) x <- as.matrix(x)
        if (env$dim!=dim(x)[2]) {
            stop("gss error in phi: inputs are of wrong dimensions")
        }
        apply(t(x)^env$pol.code[,nu],2,prod)
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}
## Make random effects for mixed-effect models
mkran <- function(formula,data)
{
    attach(data)
    ## decipher formula
    form.wk <- terms.formula(formula)[[2]]
    if (!("|"%in%strsplit(deparse(form.wk),'')[[1]]))
        stop("gss error in mkran: missing | in grouping formula")
    term.wk <- strsplit(deparse(form.wk),' \\| ')[[1]]
    ## make matrix Z
    z2.wk <- eval(parse(text=term.wk[2]))
    if (!is.factor(z2.wk))
        stop(paste("gss error in mkran: ", term.wk[2], " should be a factor"))
    z <- NULL
    lvl.z2 <- levels(z2.wk)
    for (i in lvl.z2) z <- cbind(z,as.numeric(z2.wk==i))
    ## make sigma function
    if (term.wk[1]=="1") {
        init <- 0
        env <- length(levels(z2.wk))
        fun <- function(zeta,env) diag(10^(-zeta),env)
        sigma <- list(fun=fun,env=env)
    }
    else {
        z1.wk <- eval(parse(text=term.wk[1]))
        if (!is.factor(z1.wk))
            stop(paste("gss error in mkran: ", term.wk[1], " should be a factor"))
        ind <- lvl.wk <- NULL
        nz <- length(lvl.z2)
        nsig <- length(levels(z1.wk))
        for (i in levels(z1.wk)) {
            zz.wk <- z2.wk[z1.wk==i,drop=TRUE]
            ind <- c(ind,list((1:nz)[lvl.z2%in%levels(zz.wk)]))
            lvl.wk <- c(lvl.wk,levels(zz.wk))
        }
        if (max(table(lvl.wk)>1))
            stop("gss error in mkran: ", term.wk[2], " should be nested under ", term.wk[1])
        init <- rep(0, length(levels(z1.wk)))
        env <- list(size=nz,nsig=nsig,ind=ind)
        fun <- function(zeta,env) {
            wk <- rep(0,env$size)
            for (i in 1:env$nsig) wk[env$ind[[i]]] <- 10^(-zeta[i])
            diag(wk)
        }
        sigma <- list(fun=fun,env=env)
    }
    list(z=z,sigma=sigma,init=init)
}
## Make RK for nominal shrinkage
mkrk.nominal <- function(levels)
{
    k <- length(levels)
    if (k<2) stop("gss error: factor should have at least two levels")
    code <- 1:k
    names(code) <- as.character(levels)
    ## Create the environment
    env <- list(code=code,table=diag(k)-1/k)
    ## Create the rk function
    fun <- function(x, y, env, outer.prod = FALSE) {
        if (!(is.factor(x)&is.factor(y))) {
            stop("gss error in rk: inputs are of wrong types")
        }
        x <- as.numeric(env$code[as.character(x)])
        y <- as.numeric(env$code[as.character(y)])
        if (any(is.na(c(x,y)))) {
            stop("gss error in rk: unknown factor levels")
        }
        if (outer.prod) env$table[x, y]
        else env$table[cbind(x,y)]
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make RK for ordinal shrinkage
mkrk.ordinal <- function(levels)
{
    k <- length(levels)
    if (k<2) stop("gss error: factor should have at least two levels")
    code <- 1:k
    names(code) <- as.character(levels)
    ## penalty matrix
    if (k==2) {
        B <- diag(.25,2)
        B[1,2] <- B[2,1] <- -.25
    }
    else {
        B <- diag(2,k)
        B[1,1] <- B[k,k] <- 1
        diag(B[-1,-k]) <- diag(B[-k,-1]) <- -1
        ## Moore-Penrose inverse
        B <- eigen(B)
        B <- B$vec[,-k] %*% diag(1/B$val[-k]) %*% t(B$vec[,-k])
        tol <- sqrt(.Machine$double.eps)
        B <- ifelse(abs(B)<tol,0,B)
    }
    ## Create the environment
    env <- list(code=code,table=B)
    ## Create the rk function
    fun <- function(x, y, env, outer.prod = FALSE) {
        if (!(is.factor(x)&is.factor(y))) {
            stop("gss error in rk: inputs are of wrong types")
        }
        x <- as.numeric(env$code[as.character(x)])
        y <- as.numeric(env$code[as.character(y)])
        if (any(is.na(c(x,y)))) {
            stop("gss error in rk: unknown factor levels")
        }
        if (outer.prod) env$table[x, y]
        else env$table[cbind(x,y)]
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}
## Make phi and rk for cubic spline model terms
mkterm.cubic <- function(mf,ext)
{
    ## Obtain model terms
    mt <- attr(mf,"terms")
    xvars <- as.character(attr(mt,"variables"))[-1]
    xfacs <- attr(mt,"factors")
    term.labels <- labels(mt)
    if (attr(attr(mf,"terms"),"intercept"))
        term.labels <- c("1",term.labels)
    ## Create the phi and rk functions
    term <- list(labels=term.labels)
    iphi.wk <- 1
    irk.wk <- 1
    for (label in term.labels) {
        iphi <- irk <- phi <- rk <- NULL
        if (label=="1") {
            ## the constant term
            iphi <- iphi.wk
            iphi.wk <- iphi.wk + 1
            term[[label]] <- list(iphi=iphi,nphi=1,nrk=0)
            next
        }
        vlist <- xvars[as.logical(xfacs[,label])]
        x <- mf[,vlist]
        dm <- length(vlist)
        if (dm==1) {
            if (!is.factor(x)) {
                ## numeric variable
                mx <- max(x)
                mn <- min(x)
                range <- mx - mn
                ## phi
                phi.env <- mkphi.cubic(c(mn,mx)+c(-1,1)*ext*range)
                phi.fun <- function(x,nu=1,env) env$fun(x,env$env)
                nphi <- 1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=phi.env)
                ## rk
                rk.env <- mkrk.cubic(c(mn,mx)+c(-1,1)*ext*range)
                rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                    env$fun(x,y,env$env,outer.prod)
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
            else {
                ## factor variable
                if (!is.ordered(x)) fun.env <- mkrk.nominal(levels(x))
                else fun.env <- mkrk.ordinal(levels(x))
                if (nlevels(x)>2) {
                    ## phi
                    nphi <- 0
                    ## rk
                    rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                        env$fun(x,y,env$env,outer.prod)
                    }
                    nrk <- 1
                    irk <- irk.wk
                    irk.wk <- irk.wk + nrk
                    rk <- list(fun=rk.fun,env=fun.env)
                }
                else {
                    ## phi
                    phi.fun <- function(x,nu=1,env) {
                        wk <- as.factor(names(env$env$code)[1])
                        env$fun(x,wk,env$env)
                    }
                    nphi <- 1
                    iphi <- iphi.wk
                    iphi.wk <- iphi.wk + nphi
                    phi <- list(fun=phi.fun,env=fun.env)
                    ## rk
                    nrk <- 0
                }
            }
        }    
        else {
            bin.fac <- n.phi <- phi.list <- rk.list <- NULL
            for (i in 1:dm) {
                if (!is.factor(x[[i]])) {
                    ## numeric variable
                    mx <- max(x[[i]])
                    mn <- min(x[[i]])
                    range <- mx - mn
                    phi.wk <- mkphi.cubic(c(mn,mx)+c(-1,1)*ext*range)
                    rk.wk <- mkrk.cubic(c(mn,mx)+c(-1,1)*ext*range)
                    n.phi <- c(n.phi,1)
                    bin.fac <- c(bin.fac,0)
                }
                else {
                    ## factor variable
                    if (!is.ordered(x[[i]]))
                        rk.wk <- mkrk.nominal(levels(x[[i]]))
                    else rk.wk <- mkrk.ordinal(levels(x[[i]]))
                    phi.wk <- rk.wk
                    n.phi <- c(n.phi,0)
                    bin.fac <- c(bin.fac,!(nlevels(x[[i]])>2))
                }
                phi.list <- c(phi.list,list(phi.wk))
                rk.list <- c(rk.list,list(rk.wk))
            }
            ## phi
            if (sum(n.phi+bin.fac)<dm) nphi <- 0
            else {
                phi.env <- list(dim=dm,n.phi=n.phi,phi=phi.list)
                phi.fun <- function(x,nu=1,env) {
                    z <- 1
                    for (i in 1:env$dim) {
                        if (env$n.phi[i])
                            z <- z * env$phi[[i]]$fun(x[[i]],env$phi[[i]]$env)
                        else {
                            wk <- as.factor(names(env$phi[[i]]$env$code)[1])
                            z <- z * env$phi[[i]]$fun(x[[i]],wk,env$phi[[i]]$env)
                        }
                    }
                    z
                }
                nphi <- 1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=phi.env)
            }
            ## rk
            rk.env <- list(dim=dm,n.phi=n.phi,nphi=nphi,phi=phi.list,rk=rk.list)
            rk.fun <- function(x,y,nu,env,outer.prod=FALSE) {
                div <- env$n.phi + 1
                ind <- nu - 1 + env$nphi
                z <- 1
                for (i in 1:env$dim) {
                    code <- ind%%div[i] + 1
                    ind <- ind%/%div[i]
                    if (code==div[i])
                        z <- z * env$rk[[i]]$fun(x[[i]],y[[i]],
                                                 env$rk[[i]]$env,outer.prod)
                    else {
                        phix <- env$phi[[i]]$fun(x[[i]],env$phi[[i]]$env)
                        phiy <- env$phi[[i]]$fun(y[[i]],env$phi[[i]]$env)
                        if (outer.prod) z <- z * outer(phix,phiy)
                        else z <- z * phix * phiy
                    }
                }
                z
            }
            nrk <- prod(n.phi+1) - nphi
            irk <- irk.wk
            irk.wk <- irk.wk + nrk
            rk <- list(fun=rk.fun,env=rk.env)
        }
        term[[label]] <- list(vlist=vlist,
                              iphi=iphi,nphi=nphi,phi=phi,
                              irk=irk,nrk=nrk,rk=rk)
    }
    term
}
## Make phi and rk for cubic spline model terms
mkterm.cubic1 <- function(mf,range)
{
    ## Obtain model terms
    mt <- attr(mf,"terms")
    xvars <- as.character(attr(mt,"variables"))[-1]
    xfacs <- attr(mt,"factors")
    term.labels <- labels(mt)
    if (attr(attr(mf,"terms"),"intercept"))
        term.labels <- c("1",term.labels)
    ## Create the phi and rk functions
    term <- list(labels=term.labels)
    iphi.wk <- 1
    irk.wk <- 1
    for (label in term.labels) {
        iphi <- irk <- phi <- rk <- NULL
        if (label=="1") {
            ## the constant term
            iphi <- iphi.wk
            iphi.wk <- iphi.wk + 1
            term[[label]] <- list(iphi=iphi,nphi=1,nrk=0)
            next
        }
        vlist <- xvars[as.logical(xfacs[,label])]
        x <- mf[,vlist]
        lmt <- range[,vlist]
        dm <- length(vlist)
        if (dm==1) {
            if (!is.factor(x)) {
                ## numeric variable
                mn <- min(lmt)
                mx <- max(lmt)
                ## phi
                phi.env <- mkphi.cubic(c(mn,mx))
                phi.fun <- function(x,nu=1,env) env$fun(x,env$env)
                nphi <- 1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=phi.env)
                ## rk
                rk.env <- mkrk.cubic(c(mn,mx))
                rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                    env$fun(x,y,env$env,outer.prod)
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
            else {
                ## factor variable
                if (!is.ordered(x)) fun.env <- mkrk.nominal(levels(x))
                else fun.env <- mkrk.ordinal(levels(x))
                if (nlevels(x)>2) {
                    ## phi
                    nphi <- 0
                    ## rk
                    rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                        env$fun(x,y,env$env,outer.prod)
                    }
                    nrk <- 1
                    irk <- irk.wk
                    irk.wk <- irk.wk + nrk
                    rk <- list(fun=rk.fun,env=fun.env)
                }
                else {
                    ## phi
                    phi.fun <- function(x,nu=1,env) {
                        wk <- as.factor(names(env$env$code)[1])
                        env$fun(x,wk,env$env)
                    }
                    nphi <- 1
                    iphi <- iphi.wk
                    iphi.wk <- iphi.wk + nphi
                    phi <- list(fun=phi.fun,env=fun.env)
                    ## rk
                    nrk <- 0
                }
            }
        }    
        else {
            bin.fac <- n.phi <- phi.list <- rk.list <- NULL
            for (i in 1:dm) {
                if (!is.factor(x[[i]])) {
                    ## numeric variable
                    mn <- min(lmt[[i]])
                    mx <- max(lmt[[i]])
                    phi.wk <- mkphi.cubic(c(mn,mx))
                    rk.wk <- mkrk.cubic(c(mn,mx))
                    n.phi <- c(n.phi,1)
                    bin.fac <- c(bin.fac,0)
                }
                else {
                    ## factor variable
                    if (!is.ordered(x[[i]]))
                        rk.wk <- mkrk.nominal(levels(x[[i]]))
                    else rk.wk <- mkrk.ordinal(levels(x[[i]]))
                    phi.wk <- rk.wk
                    n.phi <- c(n.phi,0)
                    bin.fac <- c(bin.fac,!(nlevels(x[[i]])>2))
                }
                phi.list <- c(phi.list,list(phi.wk))
                rk.list <- c(rk.list,list(rk.wk))
            }
            ## phi
            if (sum(n.phi+bin.fac)<dm) nphi <- 0
            else {
                phi.env <- list(dim=dm,n.phi=n.phi,phi=phi.list)
                phi.fun <- function(x,nu=1,env) {
                    z <- 1
                    for (i in 1:env$dim) {
                        if (env$n.phi[i])
                            z <- z * env$phi[[i]]$fun(x[[i]],env$phi[[i]]$env)
                        else {
                            wk <- as.factor(names(env$phi[[i]]$env$code)[1])
                            z <- z * env$phi[[i]]$fun(x[[i]],wk,env$phi[[i]]$env)
                        }
                    }
                    z
                }
                nphi <- 1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=phi.env)
            }
            ## rk
            rk.env <- list(dim=dm,n.phi=n.phi,nphi=nphi,phi=phi.list,rk=rk.list)
            rk.fun <- function(x,y,nu,env,outer.prod=FALSE) {
                div <- env$n.phi + 1
                ind <- nu - 1 + env$nphi
                z <- 1
                for (i in 1:env$dim) {
                    code <- ind%%div[i] + 1
                    ind <- ind%/%div[i]
                    if (code==div[i])
                        z <- z * env$rk[[i]]$fun(x[[i]],y[[i]],
                                                 env$rk[[i]]$env,outer.prod)
                    else {
                        phix <- env$phi[[i]]$fun(x[[i]],env$phi[[i]]$env)
                        phiy <- env$phi[[i]]$fun(y[[i]],env$phi[[i]]$env)
                        if (outer.prod) z <- z * outer(phix,phiy)
                        else z <- z * phix * phiy
                    }
                }
                z
            }
            nrk <- prod(n.phi+1) - nphi
            irk <- irk.wk
            irk.wk <- irk.wk + nrk
            rk <- list(fun=rk.fun,env=rk.env)
        }
        term[[label]] <- list(vlist=vlist,
                              iphi=iphi,nphi=nphi,phi=phi,
                              irk=irk,nrk=nrk,rk=rk)
    }
    term
}
## Make phi and rk for linear spline model terms
mkterm.linear <- function(mf,ext)
{
    ## Obtain model terms
    mt <- attr(mf,"terms")
    xvars <- as.character(attr(mt,"variables"))[-1]
    xfacs <- attr(mt,"factors")
    term.labels <- labels(mt)
    if (attr(attr(mf,"terms"),"intercept"))
        term.labels <- c("1",term.labels)
    ## Create the phi and rk functions
    term <- list(labels=term.labels)
    iphi.wk <- irk.wk <- 1
    for (label in term.labels) {
        iphi <- irk <- phi <- rk <- NULL
        if (label=="1") {
            ## the constant term
            iphi <- iphi.wk
            iphi.wk <- iphi.wk + 1
            term[[label]] <- list(iphi=iphi,nphi=1,nrk=0)
            next
        }
        vlist <- xvars[as.logical(xfacs[,label])]
        x <- mf[,vlist]
        dm <- length(vlist)
        if (dm==1) {
            if (!is.factor(x)) {
                ## numeric variable
                mx <- max(x)
                mn <- min(x)
                range <- mx - mn
                ## phi
                nphi <- 0
                ## rk
                rk.env <- mkrk.linear(c(mn,mx)+c(-1,1)*ext*range)
                rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                    env$fun(x,y,env$env,outer.prod)
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
            else {
                ## factor variable
                if (!is.ordered(x)) fun.env <- mkrk.nominal(levels(x))
                else fun.env <- mkrk.ordinal(levels(x))
                if (nlevels(x)>2) {
                    ## phi
                    nphi <- 0
                    ## rk
                    rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                        env$fun(x,y,env$env,outer.prod)
                    }
                    nrk <- 1
                    irk <- irk.wk
                    irk.wk <- irk.wk + nrk
                    rk <- list(fun=rk.fun,env=fun.env)
                }
                else {
                    ## phi
                    phi.fun <- function(x,nu=1,env) {
                        wk <- as.factor(names(env$env$code)[1])
                        env$fun(x,wk,env$env)
                    }
                    nphi <- 1
                    iphi <- iphi.wk
                    iphi.wk <- iphi.wk + nphi
                    phi <- list(fun=phi.fun,env=fun.env)
                    ## rk
                    nrk <- 0
                }
            }
        }
        else {
            bin.fac <- rk.list <- NULL
            for (i in 1:dm) {
                if (!is.factor(x[[i]])) {
                    ## numeric variable
                    mx <- max(x[[i]])
                    mn <- min(x[[i]])
                    range <- mx - mn
                    rk.wk <- mkrk.linear(c(mn,mx)+c(-1,1)*ext*range)
                    bin.fac <- c(bin.fac,0)
                }
                else {
                    ## factor variable
                    if (!is.ordered(x[[i]]))
                        rk.wk <- mkrk.nominal(levels(x[[i]]))
                    else rk.wk <- mkrk.ordinal(levels(x[[i]]))
                    bin.fac <- c(bin.fac,!(nlevels(x[[i]])>2))
                }
                rk.list <- c(rk.list,list(rk.wk))
            }
            rk.env <- list(dim=dm,rk=rk.list)
            if (sum(bin.fac)==dm) {
                ## phi
                phi.fun <- function(x,nu=1,env) {
                    z <- 1
                    for (i in 1:env$dim) {
                        wk <- as.factor(names(env$rk[[i]]$env$code)[1])
                        z <- z * env$rk[[i]]$fun(x[[i]],wk,env$rk[[i]]$env)
                    }
                    z
                }
                nphi <- 1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=rk.env)
                ## rk
                nrk <- 0
            }
            else {             
                ## phi
                nphi <- 0
                ## rk
                rk.fun <- function(x,y,nu,env,outer.prod=FALSE) {
                    z <- 1
                    for (i in 1:env$dim)
                        z <- z * env$rk[[i]]$fun(x[[i]],y[[i]],
                                                 env$rk[[i]]$env,outer.prod)
                    z
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
        }
        term[[label]] <- list(vlist=vlist,
                              iphi=iphi,nphi=nphi,phi=phi,
                              irk=irk,nrk=nrk,rk=rk)
    }
    term
}
## Make phi and rk for linear spline model terms
mkterm.linear1 <- function(mf,range)
{
    ## Obtain model terms
    mt <- attr(mf,"terms")
    xvars <- as.character(attr(mt,"variables"))[-1]
    xfacs <- attr(mt,"factors")
    term.labels <- labels(mt)
    if (attr(attr(mf,"terms"),"intercept"))
        term.labels <- c("1",term.labels)
    ## Create the phi and rk functions
    term <- list(labels=term.labels)
    iphi.wk <- irk.wk <- 1
    for (label in term.labels) {
        iphi <- irk <- phi <- rk <- NULL
        if (label=="1") {
            ## the constant term
            iphi <- iphi.wk
            iphi.wk <- iphi.wk + 1
            term[[label]] <- list(iphi=iphi,nphi=1,nrk=0)
            next
        }
        vlist <- xvars[as.logical(xfacs[,label])]
        x <- mf[,vlist]
        lmt <- range[,vlist]
        dm <- length(vlist)
        if (dm==1) {
            if (!is.factor(x)) {
                ## numeric variable
                mx <- max(lmt)
                mn <- min(lmt)
                ## phi
                nphi <- 0
                ## rk
                rk.env <- mkrk.linear(c(mn,mx))
                rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                    env$fun(x,y,env$env,outer.prod)
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
            else {
                ## factor variable
                if (!is.ordered(x)) fun.env <- mkrk.nominal(levels(x))
                else fun.env <- mkrk.ordinal(levels(x))
                if (nlevels(x)>2) {
                    ## phi
                    nphi <- 0
                    ## rk
                    rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                        env$fun(x,y,env$env,outer.prod)
                    }
                    nrk <- 1
                    irk <- irk.wk
                    irk.wk <- irk.wk + nrk
                    rk <- list(fun=rk.fun,env=fun.env)
                }
                else {
                    ## phi
                    phi.fun <- function(x,nu=1,env) {
                        wk <- as.factor(names(env$env$code)[1])
                        env$fun(x,wk,env$env)
                    }
                    nphi <- 1
                    iphi <- iphi.wk
                    iphi.wk <- iphi.wk + nphi
                    phi <- list(fun=phi.fun,env=fun.env)
                    ## rk
                    nrk <- 0
                }
            }
        }
        else {
            bin.fac <- rk.list <- NULL
            for (i in 1:dm) {
                if (!is.factor(x[[i]])) {
                    ## numeric variable
                    mx <- max(lmt[[i]])
                    mn <- min(lmt[[i]])
                    rk.wk <- mkrk.linear(c(mn,mx))
                    bin.fac <- c(bin.fac,0)
                }
                else {
                    ## factor variable
                    if (!is.ordered(x[[i]]))
                        rk.wk <- mkrk.nominal(levels(x[[i]]))
                    else rk.wk <- mkrk.ordinal(levels(x[[i]]))
                    bin.fac <- c(bin.fac,!(nlevels(x[[i]])>2))
                }
                rk.list <- c(rk.list,list(rk.wk))
            }
            rk.env <- list(dim=dm,rk=rk.list)
            if (sum(bin.fac)==dm) {
                ## phi
                phi.fun <- function(x,nu=1,env) {
                    z <- 1
                    for (i in 1:env$dim) {
                        wk <- as.factor(names(env$rk[[i]]$env$code)[1])
                        z <- z * env$rk[[i]]$fun(x[[i]],wk,env$rk[[i]]$env)
                    }
                    z
                }
                nphi <- 1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=rk.env)
                ## rk
                nrk <- 0
            }
            else {             
                ## phi
                nphi <- 0
                ## rk
                rk.fun <- function(x,y,nu,env,outer.prod=FALSE) {
                    z <- 1
                    for (i in 1:env$dim)
                        z <- z * env$rk[[i]]$fun(x[[i]],y[[i]],
                                                 env$rk[[i]]$env,outer.prod)
                    z
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
        }
        term[[label]] <- list(vlist=vlist,
                              iphi=iphi,nphi=nphi,phi=phi,
                              irk=irk,nrk=nrk,rk=rk)
    }
    term
}
## Make phi and rk for thin-plate spline model terms
mkterm.tp <- function(mf,order,mesh,weight)
{
    order <- max(order,1)
    ## Obtain model terms
    mt <- attr(mf,"terms")
    xvars <- as.character(attr(mt,"variables"))[-1]
    xfacs <- attr(mt,"factors")
    term.labels <- labels(mt)
    if (attr(attr(mf,"terms"),"intercept"))
        term.labels <- c("1",term.labels)
    ## Create the phi and rk functions
    term <- list(labels=term.labels)
    iphi.wk <- 1
    irk.wk <- 1
    for (label in term.labels) {
        iphi <- irk <- phi <- rk <- NULL
        if (label=="1") {
            ## the constant term
            iphi <- iphi.wk
            iphi.wk <- iphi.wk + 1
            term[[label]] <- list(iphi=iphi,nphi=1,nrk=0)
            next
        }
        vlist <- xvars[as.logical(xfacs[,label])]
        x <- mf[,vlist]
        xmesh <- mesh[,vlist]
        dm <- length(vlist)
        if (dm==1) {
            if (!is.factor(x)) {
                ## numeric variable
                if (is.vector(x)) xdim <- 1
                else xdim <- dim(x)[2]
                ## phi
                phi.env <- mkphi.tp(xdim,order,xmesh,weight)
                phi.fun <- function(x,nu,env) {
                    env$fun(x,nu+1,env$env)
                }
                nphi <- choose(xdim+order-1,xdim)-1
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=phi.env)
                ## rk
                rk.env <- mkrk.tp(xdim,order,xmesh,weight)
                rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                    env$fun(x,y,env$env,outer.prod)
                }
                nrk <- 1
                irk <- irk.wk
                irk.wk <- irk.wk + nrk
                rk <- list(fun=rk.fun,env=rk.env)
            }
            else {
                ## factor variable
                if (!is.ordered(x)) fun.env <- mkrk.nominal(levels(x))
                else fun.env <- mkrk.ordinal(levels(x))
                if (nlevels(x)>2) {
                    ## phi
                    nphi <- 0
                    ## rk
                    rk.fun <- function(x,y,nu=1,env,outer.prod=FALSE) {
                        env$fun(x,y,env$env,outer.prod)
                    }
                    nrk <- 1
                    irk <- irk.wk
                    irk.wk <- irk.wk + nrk
                    rk <- list(fun=rk.fun,env=fun.env)
                }
                else {
                    ## phi
                    phi.fun <- function(x,nu=1,env) {
                        wk <- as.factor(names(env$env$code)[1])
                        env$fun(x,wk,env$env)
                    }
                    nphi <- 1
                    iphi <- iphi.wk
                    iphi.wk <- iphi.wk + nphi
                    phi <- list(fun=phi.fun,env=fun.env)
                    ## rk
                    nrk <- 0
                }
            }
        }
        else {
            bin.fac <- xdim <- phi.list <- rk.list <- NULL
            for (i in 1:dm) {
                if (!is.factor(x[[i]])) {
                    ## numeric variable
                    if (is.vector(x[[i]])) xdim <- c(xdim,1)
                    else xdim <- c(xdim,dim(x[[i]])[2])
                    phi.wk <- mkphi.tp(xdim[i],order,xmesh[[i]],weight)
                    rk.wk <- mkrk.tp(xdim[i],order,xmesh[[i]],weight)
                    bin.fac <- c(bin.fac,0)
                }
                else {
                    ## factor variable
                    xdim <- c(xdim,0)
                    if (!is.ordered(x[[i]]))
                        rk.wk <- mkrk.nominal(levels(x[[i]]))
                    else rk.wk <- mkrk.ordinal(levels(x[[i]]))
                    phi.wk <- rk.wk
                    bin.fac <- c(bin.fac,!(nlevels(x[[i]])>2))
                }
                phi.list <- c(phi.list,list(phi.wk))
                rk.list <- c(rk.list,list(rk.wk))
            }
            n.phi <- choose(xdim+order-1,xdim)-1
            ## phi
            if (!all(n.phi+bin.fac)) nphi <- 0
            else {
                phi.env <- list(dim=dm,phi=phi.list,n.phi=n.phi,bin.fac=bin.fac)
                phi.fun <- function(x,nu,env) {
                    ind <- nu - 1
                    z <- 1
                    for (i in 1:env$dim) {
                        if (env$bin.fac[i]) {
                            wk <- as.factor(names(env$phi[[i]]$env$code)[1])
                            z <- z * env$phi[[i]]$fun(x[[i]],wk,env$phi[[i]]$env)
                        }
                        else {
                            code <- ind%%env$n.phi[i] + 1
                            ind <- ind%/%env$n.phi[i]
                            z <- z * env$phi[[i]]$fun(x[[i]],code+1,env$phi[[i]]$env)
                        }
                    }
                    z
                }
                nphi <- prod(n.phi+bin.fac)
                iphi <- iphi.wk
                iphi.wk <- iphi.wk + nphi
                phi <- list(fun=phi.fun,env=phi.env)
            }
            ## rk
            rk.env <- list(dim=dm,n.phi=n.phi,nphi=nphi,
                           phi=phi.list,rk=rk.list)
            rk.fun <- function(x,y,nu,env,outer.prod=FALSE) {
                n.rk <- ifelse(env$n.phi,2,1)
                ind <- nu - !env$nphi
                z <- 1
                for (i in 1:env$dim) {
                    code <- ind%%n.rk[i] + 1
                    ind <- ind%/%n.rk[i]
                    if (code==n.rk[i]) {
                        z <- z * env$rk[[i]]$fun(x[[i]],y[[i]],
                                                 env$rk[[i]]$env,outer.prod)
                    }
                    else {
                        z.wk <- 0
                        for (j in 1:env$n.phi[i]) {
                            phix <- env$phi[[i]]$fun(x[[i]],j+1,env$phi[[i]]$env)
                            phiy <- env$phi[[i]]$fun(y[[i]],j+1,env$phi[[i]]$env)
                            if (outer.prod) z.wk <- z.wk + outer(phix,phiy)
                            else z.wk <- z.wk + phix * phiy
                        }
                        z <- z * z.wk
                    }
                }
                z
            }
            n.rk <- ifelse(n.phi,2,1)
            nrk <- prod(n.rk) - as.logical(nphi)
            irk <- irk.wk
            irk.wk <- irk.wk + nrk
            rk <- list(fun=rk.fun,env=rk.env)
        }
        term[[label]] <- list(vlist=vlist,
                              iphi=iphi,nphi=nphi,phi=phi,
                              irk=irk,nrk=nrk,rk=rk)
    }
    term
}
## minimization of univariate function on finite intervals
## using 3-point quadratic fit with golden-section safe-guard
nlm0 <- function(fun,range,prec=1e-7)
{
    ratio <- 2/(sqrt(5)+1)
    ll.x <- min(range)
    uu.x <- max(range)
    if (uu.x-ll.x<prec) {
        sol <- (ll.x+uu.x)/2
        fit <- fun(sol)
        return(list(estimate=sol,minimum=fit,evaluations=1))
    }
    ml.x <- uu.x - ratio*(uu.x-ll.x)
    mu.x <- ll.x + ratio*(uu.x-ll.x)
    ## Initialization
    uu.fit <- fun(uu.x)
    mu.fit <- fun(mu.x)
    ml.fit <- fun(ml.x)
    ll.fit <- fun(ll.x)
    neval <- 4
    ## Iteration
    repeat {
        ## Fit a parabola to the 3 best points and find its minimum
        if (ll.fit<uu.fit) {
            delta.l <- ml.x-ll.x
            sigma.l <- ml.x+ll.x
            d.l <- (ml.fit-ll.fit)/delta.l
            delta.u <- mu.x-ml.x
            d.u <- (mu.fit-ml.fit)/delta.u
        }
        else {
            delta.l <- mu.x-ml.x
            sigma.l <- mu.x+ml.x
            d.l <- (mu.fit-ml.fit)/delta.l
            delta.u <- uu.x-mu.x
            d.u <- (uu.fit-mu.fit)/delta.u
        }
        a <- (d.u-d.l)/(delta.l+delta.u)
        b <- d.l-a*sigma.l
        if (a<=0) nn.x <- max(range)
        else nn.x <- -b/2/a
        ## New bracket
        if (ml.fit<mu.fit) {
            uu.x <- mu.x
            uu.fit <- mu.fit
            mm.x <- ml.x
            mm.fit <- ml.fit
        }
        else {
            ll.x <- ml.x
            ll.fit <- ml.fit
            mm.x <- mu.x
            mm.fit <- mu.fit
        }
        range.l <- mm.x-ll.x
        range.u <- uu.x-mm.x
        delta <- min(abs(nn.x-c(ll.x,mm.x,uu.x)))
        ## Safeguard
        if ((nn.x<ll.x)|(nn.x>uu.x)|(delta<prec)) {
            if (range.u>range.l) nn.x <- uu.x - ratio*range.u
            else nn.x <- ll.x + ratio*range.l
        }
        ## Update middle points
        nn.fit <- fun(nn.x)
        neval <- neval + 1
        if (nn.x<mm.x) {
            ml.x <- nn.x
            ml.fit <- nn.fit
            mu.x <- mm.x
            mu.fit <- mm.fit
        }
        else {
            ml.x <- mm.x
            ml.fit <- mm.fit
            mu.x <- nn.x
            mu.fit <- nn.fit
        }
        ## Return results
        if ((range.l+range.u<.5)&(abs(mm.x-nn.x)<sqrt(prec))) {
            if (nn.fit<mm.fit) {
                solution <- nn.x
                fit <- nn.fit
            }
            else {
                solution <- mm.x
                fit <- mm.fit
            }
            break
        }
    }
    list(estimate=solution,minimum=fit,evaluations=neval)
}
## Calculate prediction and Bayesian SE from ssanova objects
predict.ssanova <- function(object,newdata,se.fit=FALSE,
                            include=object$terms$labels,...)
{
    nnew <- dim(newdata)[1]
    nobs <- length(object$c)
    ## Extract included terms
    term <- object$terms
    philist <- rklist <- NULL
    s <- q <- NULL
    nq <- 0
    for (label in include) {
        if (label=="1") {
            philist <- c(philist,term[[label]]$iphi)
            s <- cbind(s,rep(1,len=nnew))
            next
        }
        if (label=="partial") next
        if (label=="offset") next
        xnew <- newdata[,term[[label]]$vlist]
        x <- object$mf[,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            iphi <- term[[label]]$iphi
            phi <- term[[label]]$phi
            for (i in 1:nphi) {
                philist <- c(philist,iphi+(i-1))
                s <- cbind(s,phi$fun(xnew,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            irk <- term[[label]]$irk
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                rklist <- c(rklist,irk+(i-1))
                nq <- nq+1
                q <- array(c(q,rk$fun(xnew,x,nu=i,env=rk$env,out=TRUE)),c(nnew,nobs,nq))
            }
        }
    }
    if (any(include=="partial")) {
        nphi <- term$partial$nphi
        iphi <- term$partial$iphi
        for (i in 1:nphi) philist <- c(philist,iphi+(i-1))
        s <- cbind(s,newdata$partial)
    }
    qq <- matrix(0,nnew,nobs)
    nq <- 0
    for (i in rklist) {
        nq <- nq + 1
        qq <- qq + 10^object$theta[i]*q[,,nq]
    }
    if (!is.null(object$w)) w <- object$w
    else w <- model.weights(object$mf)
    if (!is.null(w)) qq <- t(sqrt(w)*t(qq))
    ## Compute posterior mean
    nphi <- length(philist)
    pmean <- as.vector(qq%*%object$c)
    if (nphi) pmean <- pmean + as.vector(s%*%object$d[philist])
    if (any(include=="offset")) {
        if (is.null(model.offset(object$mf)))
            stop("gss error: no offset in the fit")
        offset <- newdata$offset
        if (is.null(offset)) offset <- newdata$"(offset)"
        if (is.null(offset)) stop("gss error: missing offset")
        pmean <- pmean + offset
    }
    if (se.fit) {
        b <- object$varht/10^object$nlambda
        ## Get cr, dr, and sms
        crdr <- getcrdr(object,t(qq))
        cr <- crdr$cr
        dr <- crdr$dr[philist,,drop=FALSE]
        sms <- getsms(object)[philist,philist]
        ## Compute posterior variance
        r <- 0
        for (label in include) {
            if (label=="1") next
            xnew <- newdata[,term[[label]]$vlist]
            nrk <- term[[label]]$nrk
            if (nrk) {
                irk <- term[[label]]$irk
                rk <- term[[label]]$rk
                for (i in 1:nrk) {
                    ind <- irk+(i-1)
                    r <- r + 10^object$theta[ind]*rk$fun(xnew,xnew,nu=i,env=rk$env)
                }
            }
        }
        fn2 <- function(x,n) x[1:n]%*%x[n+(1:n)]
        pvar <- r - apply(rbind(t(qq),cr),2,fn2,nobs)
        if (nphi) {
            fn1 <- function(x,sms) t(x)%*%sms%*%x
            pvar <- pvar + apply(s,1,fn1,sms)
            pvar <- pvar - 2*apply(rbind(t(s),dr),2,fn2,nphi)
        }
        pse <- as.numeric(sqrt(b*pvar))
        list(fit=pmean,se.fit=pse)
    }
    else pmean
}
## Calculate prediction and Bayesian SE from ssanova objects
predict.ssanova1 <- function(object,newdata,se.fit=FALSE,
                             include=object$terms$labels,...)
{
    nnew <- nrow(newdata)
    nbasis <- length(object$id.basis)
    nnull <- length(object$d)
    nz <- length(object$b)
    nn <- nbasis + nnull + nz
    ## Extract included terms
    term <- object$terms
    philist <- rklist <- NULL
    s <- r <- NULL
    nq <- 0
    for (label in include) {
        if (label=="1") {
            philist <- c(philist,term[[label]]$iphi)
            s <- cbind(s,rep(1,len=nnew))
            next
        }
        if (label=="partial") next
        if (label=="offset") next
        xnew <- newdata[,term[[label]]$vlist]
        x <- object$mf[object$id.basis,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            iphi <- term[[label]]$iphi
            phi <- term[[label]]$phi
            for (i in 1:nphi) {
                philist <- c(philist,iphi+(i-1))
                s <- cbind(s,phi$fun(xnew,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            irk <- term[[label]]$irk
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                rklist <- c(rklist,irk+(i-1))
                nq <- nq+1
                r <- array(c(r,rk$fun(xnew,x,nu=i,env=rk$env,out=TRUE)),c(nnew,nbasis,nq))
            }
        }
    }
    if (any(include=="partial")) {
        nphi <- term$partial$nphi
        iphi <- term$partial$iphi
        for (i in 1:nphi) philist <- c(philist,iphi+(i-1))
        s <- cbind(s,newdata$partial)
    }
    r.wk <- matrix(0,nnew,nbasis)
    nq <- 0
    for (i in rklist) {
        nq <- nq + 1
        r.wk <- r.wk + 10^object$theta[i]*r[,,nq]
    }
    ## random effects
    if (nz) {
        if (is.null(newdata$random)) z.wk <- matrix(0,nnew,nz)
        else z.wk <- newdata$random
        r.wk <- cbind(r.wk,z.wk)
    }
    ## Compute posterior mean
    nphi <- length(philist)
    pmean <- as.vector(r.wk%*%c(object$c,object$b))
    if (nphi) pmean <- pmean + as.vector(s%*%object$d[philist])
    if (any(include=="offset")) {
        if (is.null(model.offset(object$mf)))
            stop("gss error: no offset in the fit")
        offset <- newdata$offset
        if (is.null(offset)) offset <- newdata$"(offset)"
        if (is.null(offset)) stop("gss error: missing offset")
        pmean <- pmean + offset
    }
    if (se.fit) {
        b <- object$varht/10^object$nlambda
        ## Get cr, dr, and sms
        z <- .Fortran("regaux",
                      as.double(object$chol), as.integer(nn),
                      as.integer(object$jpvt), as.integer(object$rkv),
                      drcr=as.double(object$se.aux%*%t(r.wk)), as.integer(nnew),
                      sms=double(nnull^2), as.integer(nnull), double(nn*nnull),
                      PACKAGE="gss")[c("drcr","sms")]
        drcr <- matrix(z$drcr,nn,nnew)
        dr <- drcr[1:nnull,,drop=FALSE][philist,,drop=FALSE]
        cr <- drcr[(nnull+1):nn,,drop=FALSE]
        sms <- 10^object$nlambda*matrix(z$sms,nnull,nnull)[philist,philist]
        ## Compute posterior variance
        rr <- r.wk%*%object$qinv
        fn2 <- function(x,n) x[1:n]%*%x[n+(1:n)]
        pvar <- apply(t(cbind(r.wk,rr)),2,fn2,nbasis+nz)
        pvar <- pvar - apply(rbind(t(r.wk),cr),2,fn2,nbasis+nz)
        if (nphi) {
            fn1 <- function(x,sms) t(x)%*%sms%*%x
            pvar <- pvar + apply(s,1,fn1,sms)
            pvar <- pvar - 2*apply(rbind(t(s),dr),2,fn2,nphi)
        }
        pse <- as.numeric(sqrt(b*pvar))
        list(fit=pmean,se.fit=pse)
    }
    else pmean
}
## Print function for ssanova objects
print.ssanova <- function(x,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    ## terms
    cat("Terms:\n")
    print.default(x$terms$labels)
    cat("\n")
    ## terms overview
    cat("Number of unpenalized and penalized terms:\n\n")
    print.default(x$desc)
    cat("\n")
    if (x$method=="v") Method <- "GCV.\n"
    if (x$method=="m") Method <- "GML.\n"
    if (x$method=="u") Method <- "Mallows CL.\n"
    cat("Smoothing parameters are selected by",Method)
    cat("\n")
    ## the rest are suppressed
    invisible()
}

## Print function for summary.ssanova objects
print.summary.ssanova <- function (x,digits=6,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n",sep="")
    cat("\nEstimate of error standard deviation:",x$sigma,"\n")
    ## residuals
    res <- x$res
    cat("\nResiduals:\n")
    nam <- c("Min", "1Q", "Median", "3Q", "Max")
    rq <- structure(quantile(res), names = nam)
    print(rq,digits=digits)
    cat("Residual sum of squares:",x$rss)
    cat("\nR square:",x$r.squared)
    ## selected summaries
    cat("\n\nPenalty associated with the fit:",x$pen)
    cat("\n\n")
    invisible()
}

## Print function for summary.gssanova objects
print.summary.gssanova <- function (x,digits=6,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n",sep="")
    if (x$method=="u")
        cat("\n(Dispersion parameter for ",x$family,
            " family taken to be ",format(x$dispersion),")\n\n",sep="")
    if (x$method=="v")
        cat("\n(Dispersion parameter for ",x$family,
            " family estimated to be ",format(x$dispersion),")\n\n",sep="")
    ## residuals
    res <- x$res
    cat("Working residuals (weighted):\n")
    nam <- c("Min", "1Q", "Median", "3Q", "Max")
    rq <- structure(quantile(res), names = nam)
    print(rq,digits=digits)
    cat("Residual sum of squares:",x$rss,"\n")
    ## deviance residuals
    res <- x$dev.res
    cat("\nDeviance residuals:\n")
    nam <- c("Min", "1Q", "Median", "3Q", "Max")
    rq <- structure(quantile(res), names = nam)
    print(rq,digits=digits)
    cat("Deviance:",x$deviance)
    cat("\nNull deviance:",x$dev.null)
    ## selected summaries
    cat("\n\nPenalty associated with the fit:",x$pen)
    cat("\n\nNumber of performance-oriented iterations:",x$iter)
    cat("\n\n")
    invisible()
}

## Print function for ssden objects
print.ssden <- function(x,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    ## terms
    cat("Terms:\n")
    print.default(x$terms$labels)
    cat("\n")
    ## terms overview
    cat("Number of unpenalized and penalized terms:\n\n")
    print.default(x$desc)
    cat("\n")
    cat("Smoothing parameters are selected by CV with alpha=",x$alpha,".",sep="")
    cat("\n")
    ## the rest are suppressed
    invisible()
}

## Print function for ssden objects
print.sshzd <- function(x,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    ## terms
    cat("Terms:\n")
    print.default(x$terms$labels)
    cat("\n")
    ## terms overview
    cat("Number of unpenalized and penalized terms:\n\n")
    print.default(x$desc)
    cat("\n")
    cat("Smoothing parameters are selected by CV with alpha=",x$alpha,".",sep="")
    cat("\n")
    ## the rest are suppressed
    invisible()
}

## Print function for ssanova1 objects
print.ssanova1 <- function(x,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    ## terms
    cat("Terms:\n")
    print.default(x$terms$labels)
    cat("\n")
    ## terms overview
    cat("Number of unpenalized and penalized terms:\n\n")
    print.default(x$desc)
    cat("\n")
    if (x$method=="v") Method <- "GCV "
    if (x$method=="m") Method <- "GML.\n"
    if (x$method=="u") Method <- "Mallows CL "
    if (x$method=="m") cat("Smoothing parameters are selected by",Method)
    else cat("Smoothing parameters are selected by ",Method,"with alpha=",x$alpha,".",sep="")
    cat("\n")
    ## the rest are suppressed
    invisible()
}

## Print function for gssanova1 objects
print.gssanova1 <- function(x,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    ## terms
    cat("Terms:\n")
    print.default(x$terms$labels)
    cat("\n")
    ## terms overview
    cat("Number of unpenalized and penalized terms:\n\n")
    print.default(x$desc)
    cat("\n")
    cat("Smoothing parameters are selected by CV with alpha=",x$alpha,".",sep="")
    cat("\n")
    ## the rest are suppressed
    invisible()
}

## Print function for summary.gssanova objects
print.summary.gssanova1 <- function (x,digits=6,...)
{
    ## call
    cat("\nCall:\n",deparse(x$call),"\n",sep="")
    if (x$family%in%c("Gamma","inverse.gaussian")) {
        cat("\n(Dispersion parameter for ",x$family,
            " family estimated to be ",format(x$dispersion),")\n\n",sep="")
    }
    else {
        cat("\n(Dispersion parameter for ",x$family,
            " family taken to be ",format(x$dispersion),")\n\n",sep="")
    }
    ## residuals
    res <- x$res
    cat("Working residuals (weighted):\n")
    nam <- c("Min", "1Q", "Median", "3Q", "Max")
    rq <- structure(quantile(res), names = nam)
    print(rq,digits=digits)
    cat("Residual sum of squares:",x$rss,"\n")
    ## deviance residuals
    res <- x$dev.res
    cat("\nDeviance residuals:\n")
    nam <- c("Min", "1Q", "Median", "3Q", "Max")
    rq <- structure(quantile(res), names = nam)
    print(rq,digits=digits)
    cat("Deviance:",x$deviance)
    cat("\nNull deviance:",x$dev.null)
    ## selected summaries
    cat("\n\nPenalty associated with the fit:",x$pen)
    cat("\n\n")
    invisible()
}
## Calculate Kullback-Leibler projection from gssanova1 objects
project.gssanova1 <- function(object,include,...)
{
    nobs <- nrow(object$mf)
    nxi <- length(object$id.basis)
    ## evaluate full model
    family <- object$family
    eta <- object$eta
    y <- model.response(object$mf,"numeric")
    wt <- model.weights(object$mf)
    if(is.null(wt)) wt <- rep(1,nobs)
    offset <- model.offset(object$mf)
    if (!is.null(object$random)) {
        if (is.null(offset)) offset <- 0
        offset <- offset + object$random$z%*%object$b
    }
    nu <- object$nu
    y0 <- switch(family,
                 binomial=y0.binomial(y,eta,wt),
                 poisson=y0.poisson(eta),
                 Gamma=y0.Gamma(eta),
                 nbinomial=y0.nbinomial(y,eta,nu),
                 weibull=y0.weibull(y,eta,nu),
                 lognorm=y0.lognorm(y,eta,nu),
                 loglogis=y0.loglogis(y,eta,nu))
    # calculate constant fit
    cfit <- switch(family,
                   binomial=cfit.binomial(y,wt,offset),
                   poisson=cfit.poisson(y,wt,offset),
                   Gamma=cfit.Gamma(y,wt,offset),
                   nbinomial=cfit.nbinomial(y,eta,wt,nu),
                   weibull=cfit.weibull(y,wt,offset,nu),
                   lognorm=cfit.lognorm(y,wt,offset,nu),
                   loglogis=cfit.loglogis(y,wt,offset,nu))
    # calculate total entropy
    kl0 <- switch(family,
                  binomial=kl.binomial(eta,cfit,y0$wt),
                  poisson=kl.poisson(eta,cfit,wt),
                  Gamma=kl.Gamma(eta,cfit,wt),
                  nbinomial=kl.nbinomial(eta,cfit,wt,y0$nu),
                  weibull=kl.weibull(eta,cfit,wt,nu,y0$int),
                  lognorm=kl.lognorm(eta,cfit,wt,nu,y0),
                  loglogis=kl.loglogis(eta,cfit,wt,nu,y0))
    ## extract terms in subspace
    s <- matrix(1,nobs,1)
    philist <- object$term[["1"]]$iphi
    r <- NULL
    theta <- NULL
    nq.wk <- nq <- 0
    for (label in object$terms$labels) {
        if (label=="1") next
        x <- object$mf[,object$term[[label]]$vlist]
        x.basis <- object$mf[object$id.basis,object$term[[label]]$vlist]
        nphi <- object$term[[label]]$nphi
        nrk <- object$term[[label]]$nrk
        if (nphi) {
            phi <- object$term[[label]]$phi
            for (i in 1:nphi) {
                if (!any(label==include)) next
                philist <- c(philist,object$term[[label]]$iphi+(i-1))
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            rk <- object$term[[label]]$rk
            for (i in 1:nrk) {
                nq.wk <- nq.wk + 1
                if (!any(label==include)) next
                nq <- nq + 1
                theta <- c(theta,object$theta[nq.wk])
                r <- array(c(r,rk$fun(x,x.basis,nu=i,env=rk$env,out=TRUE)),
                           c(nobs,nxi,nq))
            }
        }
    }
    if (any(include=="partial")) {
        nphi <- object$term$partial$nphi
        for (i in 1:nphi)
            philist <- c(philist,object$term$partial$iphi+(i-1))
        s <- cbind(s,object$mf$partial)
    }
    ## calculate projection
    my.wls <- function(theta1=NULL) {
        if (!nq) {
            q <- matrix(0)
            sr <- cbind(s,0)
            z <- ngreg.proj(dc,family,sr,q,y0,wt,offset,nu)
        }
        else {
            theta.wk <- 1:nq
            theta.wk[fix] <- theta[fix]
            if (nq-1) theta.wk[-fix] <- theta1
            sr <- 0
            for (i in 1:nq) sr <- sr + 10^theta.wk[i]*r[,,i]
            q <- sr[object$id.basis,]
            sr <- cbind(s,sr)
            z <- ngreg.proj(dc,family,sr,q,y0,wt,offset,nu)
        }
        assign("dc",z$dc,inherit=TRUE)
        assign("eta1",z$eta,inherit=TRUE)
        z$kl
    }
    cv.wk <- function(theta) cv.scale*my.wls(theta)+cv.shift
    ## initialization
    r.wk <- 0
    for (i in 1:nq) r.wk <- r.wk + 10^theta[i]*r[,,i]
    if (is.null(s)) theta.wk <- 0
    else theta.wk <- log10(sum(s^2)/ncol(s)/sum(r.wk^2)*nxi) / 2
    theta <- theta + theta.wk
    tmp <- NULL
    for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(r[cbind(object$id.basis,1:nxi,i)]))
    fix <- rev(order(tmp))[1]
    ## projection
    if (nq) dc <- c(object$d[philist],10^(-theta.wk)*object$c)
    else dc <- c(object$d[philist],0)
    eta1 <- NULL
    if (nq>1) {
        ## scale and shift cv
        tmp <- abs(my.wls(theta[-fix]))
        cv.scale <- 1
        cv.shift <- 0
        if (tmp<1&tmp>10^(-4)) {
            cv.scale <- 10/tmp
            cv.shift <- 0
        }
        if (tmp<10^(-4)) {
            cv.scale <- 10^2
            cv.shift <- 10
        }
        zz <- nlm(cv.wk,theta[-fix],stepmax=1,ndigit=7)
        if (zz$code>3)
            warning("gss warning in project.gssanova1: theta iteration fails to converge")
        kl <- my.wls(zz$est)
    }
    else kl <- my.wls()
    ## check
    kl1 <- switch(family,
                  binomial=kl.binomial(eta1,cfit,y0$wt),
                  poisson=kl.poisson(eta1,cfit,wt),
                  Gamma=kl.Gamma(eta1,cfit,wt),
                  nbinomial=kl.nbinomial(eta1,cfit,wt,y0$nu),
                  weibull=kl.weibull(eta1,cfit,wt,nu,y0$int),
                  lognorm=kl.lognorm(eta1,cfit,wt,nu,y0),
                  loglogis=kl.loglogis(eta1,cfit,wt,nu,y0))
    list(ratio=kl/kl0,kl=kl,check=(kl+kl1)/kl0)
}

## KL projection with Non-Gaussian regression
ngreg.proj <- function(dc,family,sr,q,y0,wt,offset,nu)
{
    ## initialization
    q <- 10^(-5)*q
    eta <- sr%*%dc
    nobs <- length(eta)
    nn <- ncol(as.matrix(sr))
    nxi <- ncol(q)
    nnull <- nn-nxi
    if (!is.null(offset)) eta <- eta + offset
    iter <- 0
    flag <- 0
    adj <- 0
    fit1 <- switch(family,
                   binomial=proj0.binomial(y0,eta,offset),
                   poisson=proj0.poisson(y0,eta,wt,offset),
                   Gamma=proj0.Gamma(y0,eta,wt,offset),
                   nbinomial=proj0.nbinomial(y0,eta,wt,offset),
                   weibull=proj0.weibull(y0,eta,wt,offset,nu),
                   lognorm=proj0.lognorm(y0,eta,wt,offset,nu),
                   loglogis=proj0.loglogis(y0,eta,wt,offset,nu))
    kl <- fit1$kl
    ## Newton iteration
    repeat {
        if (!adj) iter <- iter+1
        ## weighted least squares fit
        if (!is.finite(sum(fit1$wt,fit1$ywk))) {
            if (flag) stop("gss error in project.gssanova1: Newton iteration diverges")
            eta <- rep(0,nobs)
            fit1 <- switch(family,
                           binomial=proj0.binomial(y0,eta,offset),
                           poisson=proj0.poisson(y0,eta,wt,offset),
                           Gamma=proj0.Gamma(y0,eta,wt,offset),
                           nbinomial=proj0.nbinomial(y0,eta,wt,offset),
                           weibull=proj0.weibull(y0,eta,wt,offset,nu),
                           lognorm=proj0.lognorm(y0,eta,wt,offset,nu),
                           loglogis=proj0.loglogis(y0,eta,wt,offset,nu))
            kl <- fit1$kl
            iter <- 0
            flag <- 1
            next
        }
        mumax <- max(abs(t(sr)%*%fit1$u+c(rep(0,nnull),q%*%dc[-(1:nnull)])))
        w <- sqrt(as.vector(fit1$wt))
        z <- .Fortran("reg",
                      as.double(w*sr), as.integer(nobs), as.integer(nnull),
                      as.double(q), as.integer(nxi), as.double(w*fit1$ywk),
                      as.integer(4),
                      double(1), double(1), double(1), dc=double(nn),
                      as.double(.Machine$double.eps),
                      double(nn*nn), double(nn), as.integer(rep(0,nn)),
                      double(max(nobs,nn)), integer(1), integer(1),
                      PACKAGE="gss")["dc"]
        dc.diff <- z$dc-dc
        adj <- 0
        repeat {
            dc.new <- dc + dc.diff
            eta.new <- sr%*%dc.new
            if (!is.null(offset)) eta.new <- eta.new + offset
            fit1 <- switch(family,
                           binomial=proj0.binomial(y0,eta.new,offset),
                           poisson=proj0.poisson(y0,eta.new,wt,offset),
                           Gamma=proj0.Gamma(y0,eta.new,wt,offset),
                           nbinomial=proj0.nbinomial(y0,eta.new,wt,offset),
                           weibull=proj0.weibull(y0,eta.new,wt,offset,nu),
                           lognorm=proj0.lognorm(y0,eta.new,wt,offset,nu),
                           loglogis=proj0.loglogis(y0,eta.new,wt,offset,nu))
            kl.new <- fit1$kl
            if (!is.finite(kl.new)) kl.new <- Inf
            if (kl.new-kl<(1e-4+abs(kl))*1e-1) break
            adj <- 1
            dc.diff <- dc.diff/2
        }
        disc0 <- max((mumax/(1+kl))^2,abs(kl.new-kl)/(1+kl))
        disc <- sum(fit1$wt*((eta-eta.new)/(1+abs(eta)))^2)/sum(fit1$wt)
        if (is.nan(disc)) {
            if (flag) stop("gss error in project.gssanova1: Newton iteration diverges")
            eta <- rep(0,nobs)
            fit1 <- switch(family,
                           binomial=proj0.binomial(y0,eta,offset),
                           poisson=proj0.poisson(y0,eta,wt,offset),
                           Gamma=proj0.Gamma(y0,eta,wt,offset),
                           nbinomial=proj0.nbinomial(y0,eta,wt,offset),
                           weibull=proj0.weibull(y0,eta,wt,offset,nu),
                           lognorm=proj0.lognorm(y0,eta,wt,offset,nu),
                           loglogis=proj0.loglogis(y0,eta,wt,offset,nu))
            kl <- fit1$kl
            iter <- 0
            flag <- 1
            next
        }
        dc <- dc.new
        eta <- eta.new
        kl <- kl.new
        if (adj) next
        if (disc0<1e-5) break
        if (disc<1e-5) break
        if (iter<=30) next
        warning("gss warning in project.gssanova1: Newton iteration fails to converge")
        break
    }
    fit1 <- switch(family,
                   binomial=proj0.binomial(y0,eta,offset),
                   poisson=proj0.poisson(y0,eta,wt,offset),
                   Gamma=proj0.Gamma(y0,eta,wt,offset),
                   nbinomial=proj0.nbinomial(y0,eta,wt,offset),
                   weibull=proj0.weibull(y0,eta,wt,offset,nu),
                   lognorm=proj0.lognorm(y0,eta,wt,offset,nu),
                   loglogis=proj0.loglogis(y0,eta,wt,offset,nu))
    kl <- fit1$kl
    list(dc=dc,eta=eta,kl=kl)
}
## Calculate Kullback-Leibler projection from ssanova1 objects
project.ssanova1 <- function(object,include,...)
{
    nobs <- nrow(object$mf)
    nxi <- length(object$id.basis)
    ## evaluate full model
    mf <- object$mf
    yy <- predict(object,mf)
    wt <- model.weights(object$mf)
    offset <- model.offset(object$mf)
    if (!is.null(offset)) yy <- yy - offset
    ## extract terms in subspace
    s <- matrix(1,nobs,1)
    r <- NULL
    theta <- NULL
    nq.wk <- nq <- 0
    for (label in object$terms$labels) {
        if (label=="1") next
        x <- object$mf[,object$term[[label]]$vlist]
        x.basis <- object$mf[object$id.basis,object$term[[label]]$vlist]
        nphi <- object$term[[label]]$nphi
        nrk <- object$term[[label]]$nrk
        if (nphi) {
            phi <- object$term[[label]]$phi
            for (i in 1:nphi) {
                if (!any(label==include)) next
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            rk <- object$term[[label]]$rk
            for (i in 1:nrk) {
                nq.wk <- nq.wk + 1
                if (!any(label==include)) next
                nq <- nq + 1
                theta <- c(theta,object$theta[nq.wk])
                r <- array(c(r,rk$fun(x,x.basis,nu=i,env=rk$env,out=TRUE)),
                           c(nobs,nxi,nq))
            }
        }
    }
    if (any(include=="partial")) s <- cbind(s,object$mf$partial)
    ## calculate projection
    my.ls <- function(theta1=NULL) {
        if (!nq) {
            q <- matrix(0)
            sr <- cbind(s,0)
        }
        else {
            theta.wk <- 1:nq
            theta.wk[fix] <- theta[fix]
            if (nq-1) theta.wk[-fix] <- theta1
            sr <- 0
            for (i in 1:nq) sr <- sr + 10^theta.wk[i]*r[,,i]
            q <- 10^(-5)*sr[object$id.basis,]
            sr <- cbind(s,sr)
        }
        nn <- ncol(as.matrix(sr))
        nnull <- nn-nxi
        if (!is.null(wt)) {
            wt <- sqrt(wt)
            sr <- wt*sr
            yy <- wt*yy
        }
        z <- .Fortran("reg",
                      as.double(sr), as.integer(nobs), as.integer(nnull),
                      as.double(q), as.integer(nxi), as.double(yy),
                      as.integer(4),
                      double(1), double(1), double(1), dc=double(nn),
                      as.double(.Machine$double.eps),
                      double(nn*nn), double(nn), as.integer(rep(0,nn)),
                      double(max(nobs,nn)), integer(1), integer(1),
                      PACKAGE="gss")["dc"]
        assign("yhat",sr%*%z$dc,inherit=TRUE)
        mean((yy-yhat)^2)
    }
    cv.wk <- function(theta) cv.scale*my.ls(theta)+cv.shift
    ## initialization
    r.wk <- 0
    for (i in 1:nq) r.wk <- r.wk + 10^theta[i]*r[,,i]
    if (is.null(s)) theta.wk <- 0
    else theta.wk <- log10(sum(s^2)/ncol(s)/sum(r.wk^2)*nxi) / 2
    theta <- theta + theta.wk
    tmp <- NULL
    for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(r[cbind(object$id.basis,1:nxi,i)]))
    fix <- rev(order(tmp))[1]
    ## projection    
    yhat <- NULL
    if (nq>1) {
        ## scale and shift cv
        tmp <- abs(my.ls(theta[-fix]))
        cv.scale <- 1
        cv.shift <- 0
        if (tmp<1&tmp>10^(-4)) {
            cv.scale <- 10/tmp
            cv.shift <- 0
        }
        if (tmp<10^(-4)) {
            cv.scale <- 10^2
            cv.shift <- 10
        }
        zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7)
        if (zz$code>3)
            warning("gss warning in project.ssanova1: theta iteration fails to converge")
        kl <- my.ls(zz$est)
    }
    else kl <- my.ls()
    kl0 <- mean((yy-mean(yy))^2)
    kl <- mean((yy-yhat)^2)
    kl1 <- mean((mean(yy)-yhat)^2)
    list(ratio=kl/kl0,kl=kl,check=(kl+kl1)/kl0)
}
## Calculate Kullback-Leibler projection from ssden objects
project.ssden <- function(object,include,mesh=FALSE,...)
{
    qd.pt <- object$quad$pt
    qd.wt <- object$quad$wt
    ## evaluate full model
    mesh0 <- dssden(object,qd.pt)
    ## extract terms in subspace
    nobs <- dim(object$mf)[1]
    nqd <- length(qd.wt)
    nxi <- length(object$id.basis)
    qd.s <- qd.r <- q <- NULL
    theta <- d <- NULL
    n0.wk <- nq.wk <- nq <- 0
    for (label in object$terms$labels) {
        x.basis <- object$mf[object$id.basis,object$term[[label]]$vlist]
        qd.x <- qd.pt[,object$term[[label]]$vlist]
        nphi <- object$term[[label]]$nphi
        nrk <- object$term[[label]]$nrk
        if (nphi) {
            phi <- object$term[[label]]$phi
            for (i in 1:nphi) {
                n0.wk <- n0.wk + 1
                if (!any(label==include)) next
                d <- c(d,object$d[n0.wk])
                qd.s <- cbind(qd.s,phi$fun(qd.x,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            rk <- object$term[[label]]$rk
            for (i in 1:nrk) {
                nq.wk <- nq.wk + 1
                if (!any(label==include)) next
                nq <- nq + 1
                theta <- c(theta,object$theta[nq.wk])
                qd.r <- array(c(qd.r,rk$fun(x.basis,qd.x,nu=i,env=rk$env,out=TRUE)),
                              c(nxi,nqd,nq))
                q <- cbind(q,rk$fun(x.basis,x.basis,nu=i,env=rk$env,out=FALSE))
            }
        }
    }
    if (!is.null(qd.s)) {
        nn <- nxi + ncol(qd.s)
        qd.s <- t(qd.s)
    }
    else nn <- nxi
    ## calculate projection
    rkl <- function(theta1=NULL) {
        theta.wk <- 1:nq
        theta.wk[fix] <- theta[fix]
        if (nq-1) theta.wk[-fix] <- theta1
        qd.rs <- 0
        for (i in 1:nq) qd.rs <- qd.rs + 10^theta.wk[i]*qd.r[,,i]
        qd.rs <- rbind(qd.rs,qd.s)
        z <- .Fortran("drkl",
                      cd=as.double(cd), as.integer(nn),
                      as.double(t(qd.rs)), as.integer(nqd), as.double(qd.wt),
                      mesh=as.double(mesh0), as.double(.Machine$double.eps),
                      double(nqd), double(nqd), double(nn), double(nn*nn),
                      integer(nn), double(nn), double(nn), double(nqd),
                      as.double(1e-6), as.integer(30),
                      info=integer(1), PACKAGE="gss")
        if (z$info==1)
            stop("gss error in project.ssden: Newton iteration diverges")
        if (z$info==2)
            warning("gss warning in project.ssden: Newton iteration fails to converge")
        assign("cd",z$cd,inherit=TRUE)
        assign("mesh1",z$mesh,inherit=TRUE)
        sum(qd.wt*log(mesh0/mesh1)*mesh0)
    }
    cv.wk <- function(theta) cv.scale*rkl(theta)+cv.shift
    ## initialization
    qd.r.wk <- 0
    for (i in 1:nq) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i]
    mu.r <- apply(qd.wt*t(qd.r.wk),2,sum)/sum(qd.wt)
    v.r <- apply(qd.wt*t(qd.r.wk^2),2,sum)/sum(qd.wt)
    mu.s <- apply(qd.wt*t(qd.s),2,sum)/sum(qd.wt)
    v.s <- apply(qd.wt*t(qd.s^2),2,sum)/sum(qd.wt)
    if (is.null(qd.s)) theta.wk <- 0
    else theta.wk <- log10(sum(v.s-mu.s^2)/(nn-nxi)/sum(v.r-mu.r^2)*nxi) / 2
    theta <- theta + theta.wk
    tmp <- NULL
    for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(q[,i]))
    fix <- rev(order(tmp))[1]
    ## projection
    cd <- c(10^(-theta.wk)*object$c,d)
    mesh1 <- NULL
    if (nq-1) {
        if (nq-2) {
            ## scale and shift cv
            tmp <- abs(rkl(theta[-fix]))
            cv.scale <- 1
            cv.shift <- 0
            if (tmp<1&tmp>10^(-4)) {
                cv.scale <- 10/tmp
                cv.shift <- 0
            }
            if (tmp<10^(-4)) {
                cv.scale <- 10^2
                cv.shift <- 10
            }
            zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7)
            if (zz$code>3)
                warning("gss warning in project.ssden: theta iteration fails to converge")
        }
        else {
            the.wk <- theta[-fix]
            repeat {
                mn <- the.wk-1
                mx <- the.wk+1
                zz <- nlm0(rkl,c(mn,mx))
                if (min(zz$est-mn,mx-zz$est)>=1e-3) break
                else the.wk <- zz$est
            }
        }
        kl <- rkl(zz$est)
    }
    else kl <- rkl()
    kl0 <- sum(qd.wt*log(mesh0)*mesh0) + log(sum(qd.wt))
    obj <- list(ratio=kl/kl0,kl=kl)
    if (mesh) obj$mesh <- mesh1
    obj
}
## Calculate Kullback-Leibler projection from sshzd objects
project.sshzd <- function(object,include,mesh=FALSE,...)
{
    if (!(object$tname%in%include))
        stop("gss error in project.sshzd: time main effect missing in included terms")
    quad.pt <- object$quad$pt
    quad.wt <- object$quad$wt
    nx <- dim(object$qd.wt)[2]
    nbasis <- length(object$id.basis)
    mesh0 <- object$mesh0
    ## extract terms in subspace
    nqd <- length(quad.pt)
    nxi <- length(object$id.basis)
    d <- qd.s <- q <- theta <- NULL
    qd.r <- as.list(NULL)
    n0.wk <- nu <- nq.wk <- nq <- 0
    for (label in object$terms$labels) {
        vlist <- object$terms[[label]]$vlist
        x.list <- object$xnames[object$xnames%in%vlist]
        xy.basis <- object$mf[object$id.basis,vlist]
        qd.xy <- data.frame(matrix(0,nqd,length(vlist)))
        names(qd.xy) <- vlist
        if (object$tname%in%vlist) qd.xy[,object$tname] <- quad.pt
        if (length(x.list)) xx <- object$x.pt[,x.list,drop=FALSE]
        else xx <- NULL
        nphi <- object$terms[[label]]$nphi
        nrk <- object$terms[[label]]$nrk
        if (nphi) {
            phi <- object$terms[[label]]$phi
            for (i in 1:nphi) {
                n0.wk <- n0.wk + 1
                if (label=="1") {
                    d <- object$d[n0.wk]
                    nu <- nu + 1
                    qd.wk <- matrix(1,nqd,nx)
                    qd.s <- array(c(qd.s,qd.wk),c(nqd,nx,nu))
                    next
                }
                if (!any(label==include)) next
                d <- c(d,object$d[n0.wk])
                nu <- nu + 1
                if (is.null(xx))
                    qd.wk <- matrix(phi$fun(qd.xy[,,drop=TRUE],nu=i,env=phi$env),nqd,nx)
                else {
                    qd.wk <- NULL
                    for (j in 1:nx) {
                        qd.xy[,x.list] <- xx[rep(j,nqd),]
                        for (k in x.list)
                            if (is.factor(xx[,k])) qd.xy[,k] <- as.factor(qd.xy[,k])
                        qd.wk <- cbind(qd.wk,phi$fun(qd.xy[,,drop=TRUE],i,phi$env))
                    }
                }
                qd.s <- array(c(qd.s,qd.wk),c(nqd,nx,nu))
            }
        }
        if (nrk) {
            rk <- object$terms[[label]]$rk
            for (i in 1:nrk) {
                nq.wk <- nq.wk + 1
                if (!any(label==include)) next
                nq <- nq + 1
                theta <- c(theta,object$theta[nq.wk])
                q <- cbind(q,rk$fun(xy.basis,xy.basis,i,rk$env,out=FALSE))
                if (is.null(xx))
                    qd.r[[nq]] <- rk$fun(qd.xy[,,drop=TRUE],xy.basis,i,rk$env,out=TRUE)
                else {
                    qd.wk <- NULL
                    for (j in 1:nx) {
                        qd.xy[,x.list] <- xx[rep(j,nqd),]
                        for (k in x.list)
                            if (is.factor(xx[,k])) qd.xy[,k] <- as.factor(qd.xy[,k])
                        qd.wk <- array(c(qd.wk,rk$fun(qd.xy[,,drop=TRUE],xy.basis,i,rk$env,TRUE)),
                                       c(nqd,nbasis,j))
                    }
                    qd.r[[nq]] <- qd.wk
                }
            }
        }
    }
    if (!is.null(qd.s)) nnull <- dim(qd.s)[3]
    else nnull <- 0
    nn <- nxi + nnull
    ## calculate projection
    rkl <- function(theta1=NULL) {
        theta.wk <- 1:nq
        theta.wk[fix] <- theta[fix]
        if (nq-1) theta.wk[-fix] <- theta1
        qd.r.wk <- array(0,c(nqd,nxi,nx))
        for (i in 1:nq) {
            if (length(dim(qd.r[[i]]))==3) qd.r.wk <- qd.r.wk + 10^theta.wk[i]*qd.r[[i]]
            else qd.r.wk <- qd.r.wk + as.vector(10^theta.wk[i]*qd.r[[i]])
        }
        qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
        qd.r.wk <- array(c(qd.r.wk,qd.s),c(nqd,nx,nn))
        qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
        z <- .Fortran("hrkl",
                      cd=as.double(cd), as.integer(nn),
                      as.double(qd.r.wk), as.integer(nqd), as.integer(nx),
                      as.double(object$qd.wt), mesh=as.double(object$qd.wt*mesh0),
                      as.double(.Machine$double.eps), double(nqd*nx),
                      double(nn), double(nn), double(nn*nn), integer(nn), double(nn),
                      double(nn), double(nqd*nx), as.double(1e-6), as.integer(30),
                      info=integer(1), PACKAGE="gss")
        if (z$info==1)
            stop("gss error in project.sshzd: Newton iteration diverges")
        if (z$info==2)
            warning("gss warning in project.sshzd: Newton iteration fails to converge")
        assign("cd",z$cd,inherit=TRUE)
        assign("mesh1",z$mesh,inherit=TRUE)
        sum(object$qd.wt*(log(mesh0/mesh1)*mesh0-mesh0+mesh1))
    }
    cv.wk <- function(theta) cv.scale*rkl(theta)+cv.shift
    ## initialization
    qd.r.wk <- array(0,c(nqd,nxi,nx))
    for (i in 1:nq) {
        if (length(dim(qd.r[[i]]))==3) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[[i]]
        else qd.r.wk <- qd.r.wk + as.vector(10^theta[i]*qd.r[[i]])
    }
    v.s <- v.r <- 0
    for (i in 1:nx) {
        if (nnull) v.s <- v.s + apply(object$qd.wt[,i]*qd.s[,i,,drop=FALSE]^2,2,sum)
        v.r <- v.r + apply(object$qd.wt[,i]*qd.r.wk[,,i,drop=FALSE]^2,2,sum)
    }
    if (nnull) theta.wk <- log10(sum(v.s)/nnull/sum(v.r)*nxi) / 2
    else theta.wk <- 0
    theta <- theta + theta.wk
    tmp <- NULL
    for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(q[,i]))
    fix <- rev(order(tmp))[1]
    ## projection
    cd <- c(10^(-theta.wk)*object$c,d)
    mesh1 <- NULL
    if (nq-1) {
        if (nq-2) {
            ## scale and shift cv
            tmp <- abs(rkl(theta[-fix]))
            cv.scale <- 1
            cv.shift <- 0
            if (tmp<1&tmp>10^(-4)) {
                cv.scale <- 10/tmp
                cv.shift <- 0
            }
            if (tmp<10^(-4)) {
                cv.scale <- 10^2
                cv.shift <- 10
            }
            zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7)
            if (zz$code>3)
                warning("gss warning in project.sshzd: theta iteration fails to converge")
        }
        else {
            the.wk <- theta[-fix]
            repeat {
                mn <- the.wk-1
                mx <- the.wk+1
                zz <- nlm0(rkl,c(mn,mx))
                if (min(zz$est-mn,mx-zz$est)>=1e-3) break
                else the.wk <- zz$est
            }
        }
        kl <- rkl(zz$est)
    }
    else kl <- rkl()
    kl0 <- sum(object$qd.wt*(log(mesh0/object$cfit)*mesh0-mesh0+object$cfit))
    obj <- list(ratio=kl/kl0,kl=kl)
    if (mesh) obj$mesh <- mesh1
    obj
}
## Fit Single Smoothing Parameter REGression
sspreg <- function(s,q,y,method="v",varht=1)
{
    ## Check inputs
    if (is.vector(s)) s <- as.matrix(s)
    if (!(is.matrix(s)&is.matrix(q)&is.vector(y)&is.character(method))) {
        stop("gss error in sspreg: inputs are of wrong types")
    }
    nobs <- length(y)
    nnull <- dim(s)[2]
    if (!((dim(s)[1]==nobs)&(dim(q)[1]==nobs)&(dim(q)[2]==nobs)
          &(nobs>=nnull)&(nnull>0))) {
        stop("gss error in sspreg: inputs have wrong dimensions")
    }
    ## Set method for smoothing parameter selection
    code <- (1:3)[c("v","m","u")==method]
    if (!length(code)) {
        stop("gss error: unsupported method for smoothing parameter selection")
    }
    ## Call RKPACK driver DSIDR
    z <- .Fortran("dsidr0",
                  as.integer(code),
                  swk=as.double(s), as.integer(nobs),
                  as.integer(nobs), as.integer(nnull),
                  as.double(y),
                  qwk=as.double(q), as.integer(nobs),
                  as.double(0), as.integer(0), double(2),
                  nlambda=double(1), score=double(1), varht=as.double(varht),
                  c=double(nobs), d=double(nnull),
                  qraux=double(nnull), jpvt=integer(nnull),
                  double(3*nobs),
                  info=integer(1),PACKAGE="gss")
    ## Check info for error
    if (info<-z$info) {               
        if (info>0)
            stop("gss error in sspreg: matrix s is rank deficient")
        if (info==-2)
            stop("gss error in sspreg: matrix q is indefinite")
        if (info==-1)
            stop("gss error in sspreg: input data have wrong dimensions")
        if (info==-3)
            stop("gss error in sspreg: unknown method for smoothing parameter selection.")
    }
    ## Return the fit
    c(list(method=method,theta=0),
      z[c("c","d","nlambda","score","varht","swk","qraux","jpvt","qwk")])
}

## Fit Multiple Smoothing Parameter REGression
mspreg <- function(s,q,y,method="v",varht=1,prec=1e-7,maxiter=30)
{
    ## Check inputs
    if (is.vector(s)) s <- as.matrix(s)
    if (!(is.matrix(s)&is.array(q)&(length(dim(q))==3)
          &is.vector(y)&is.character(method))) {
        stop("gss error in mspreg: inputs are of wrong types")
    }
    nobs <- length(y)
    nnull <- dim(s)[2]
    nq <- dim(q)[3]
    if (!((dim(s)[1]==nobs)&(dim(q)[1]==nobs)&(dim(q)[2]==nobs)
          &(nobs>=nnull)&(nnull>0)&(nq>1))) {
        stop("gss error in mspreg: inputs have wrong dimensions")
    }
    ## Set method for smoothing parameter selection
    code <- (1:3)[c("v","m","u")==method]
    if (!length(code)) {
        stop("gss error: unsupported method for smoothing parameter selection")
    }
    ## Call RKPACK driver DMUDR
    z <- .Fortran("dmudr0",
                  as.integer(code),
                  as.double(s),         # s
                  as.integer(nobs), as.integer(nobs), as.integer(nnull),
                  as.double(q),         # q
                  as.integer(nobs), as.integer(nobs), as.integer(nq),
                  as.double(y),         # y
                  as.double(0), as.integer(0),
                  as.double(prec), as.integer(maxiter),
                  theta=double(nq), nlambda=double(1),
                  score=double(1), varht=as.double(varht),
                  c=double(nobs), d=double(nnull),
                  double(nobs*nobs*(nq+2)),
                  info=integer(1),PACKAGE="gss")[c("theta","info")]
    ## Check info for error
    if (info<-z$info) {               
        if (info>0)
            stop("gss error in mspreg: matrix s is rank deficient")
        if (info==-2)
            stop("gss error in mspreg: matrix q is indefinite")
        if (info==-1)
            stop("gss error in mspreg: input data have wrong dimensions")
        if (info==-3)
            stop("gss error in mspreg: unknown method for smoothing parameter selection.")
        if (info==-4)
            stop("gss error in mspreg: iteration fails to converge, try to increase maxiter")
        if (info==-5)
            stop("gss error in mspreg: iteration fails to find a reasonable descent direction")
    }
    qwk <- 10^z$theta[1]*q[,,1]
    for (i in 2:nq) qwk <- qwk + 10^z$theta[i]*q[,,i]
    ## Call RKPACK driver DSIDR
    zz <- .Fortran("dsidr0",
                   as.integer(code),
                   swk=as.double(s), as.integer(nobs),
                   as.integer(nobs), as.integer(nnull),
                   as.double(y),
                   qwk=as.double(qwk), as.integer(nobs),
                   as.double(0), as.integer(0), double(2),
                   nlambda=double(1), score=double(1), varht=as.double(varht),
                   c=double(nobs), d=double(nnull),
                   qraux=double(nnull), jpvt=integer(nnull),
                   double(3*nobs),
                   info=integer(1),PACKAGE="gss")
    ## Return the fit
    c(list(method=method,theta=z$theta),
      zz[c("c","d","nlambda","score","varht","swk","qraux","jpvt","qwk")])
}

## Obtain c & d for new y's
getcrdr <- function(obj,r)
{
    ## Check inputs
    if (is.vector(r)) r <- as.matrix(r)
    if (!(any(class(obj)=="ssanova")&is.matrix(r))) {
        stop("gss error in getcrdr: inputs are of wrong types")
    }
    nobs <- length(obj$c)
    nnull <- length(obj$d)
    nr <- dim(r)[2]
    if (!((dim(r)[1]==nobs)&(nr>0))) {
        stop("gss error in getcrdr: inputs have wrong dimensions")
    }
    ## Call RKPACK ulitity DCRDR
    z <- .Fortran("dcrdr",
                  as.double(obj$swk), as.integer(nobs),
                  as.integer(nobs), as.integer(nnull),
                  as.double(obj$qraux), as.integer(obj$jpvt),
                  as.double(obj$qwk), as.integer(nobs),
                  as.double(obj$nlambda),
                  as.double(r), as.integer(nobs), as.integer(nr),
                  cr=double(nobs*nr), as.integer(nobs),
                  dr=double(nnull*nr), as.integer(nnull),
                  double(2*nobs), integer(1),PACKAGE="gss")[c("cr","dr")]
    ## Return cr and dr
    z$cr <- matrix(z$cr,nobs,nr)
    z$dr <- matrix(z$dr,nnull,nr)
    z
}

## Obtain var-cov matrix for fixed effects
getsms <- function(obj)
{
    ## Check input
    if (!any(class(obj)=="ssanova")) {
        stop("gss error in getsms: inputs are of wrong types")
    }
    nobs <- length(obj$c)
    nnull <- length(obj$d)
    ## Call RKPACK ulitity DSMS
    z <- .Fortran("dsms",
                  as.double(obj$swk), as.integer(nobs),
                  as.integer(nobs), as.integer(nnull),
                  as.integer(obj$jpvt),
                  as.double(obj$qwk), as.integer(nobs),
                  as.double(obj$nlambda),
                  sms=double(nnull*nnull), as.integer(nnull),
                  double(2*nobs), integer(1),PACKAGE="gss")["sms"]
    ## Return the nnull-by-nnull matrix
    matrix(z$sms,nnull,nnull)
}
## Fit Single Smoothing Parameter (Gaussian) REGression
sspreg1 <- function(s,r,q,y,method,alpha,varht,random)
{
    qr.trace <- FALSE
    if ((alpha<0)&(method%in%c("u","v"))) qr.trace <- TRUE
    alpha <- abs(alpha)
    ## get dimensions
    nobs <- nrow(r)
    nxi <- ncol(r)
    if (!is.null(s)) {
        if (is.vector(s)) nnull <- 1
        else nnull <- ncol(s)
    }
    else nnull <- 0
    if (!is.null(random)) nz <- ncol(as.matrix(random$z))
    else nz <- 0
    nxiz <- nxi + nz
    nn <- nxiz + nnull
    ## cv function
    cv <- function(lambda) {
        if (is.null(random)) q.wk <- 10^(lambda+theta)*q
        else {
            q.wk <- matrix(0,nxiz,nxiz)
            q.wk[1:nxi,1:nxi] <- 10^(lambda[1]+theta)*q
            q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
                10^(2*ran.scal)*random$sigma$fun(lambda[-1],random$sigma$env)
        }
        if (qr.trace) {
            qq.wk <- chol(q.wk,pivot=TRUE)
            sr <- cbind(s,10^theta*r[,attr(qq.wk,"pivot")])
            sr <- rbind(sr,cbind(matrix(0,nxiz,nnull),qq.wk))
            sr <- qr(sr,tol=0)
            rss <- mean(qr.resid(sr,c(y,rep(0,nxiz)))[1:nobs]^2)
            trc <- sum(qr.Q(sr)[1:nobs,]^2)/nobs
            if (method=="u") score <- rss + alpha*2*varht*trc
            if (method=="v") score <- rss/(1-alpha*trc)^2
            alpha.wk <- max(0,log.la0-lambda[1]-5)*(3-alpha) + alpha
            alpha.wk <- min(alpha.wk,3)
            if (alpha.wk>alpha) {
                if (method=="u") score <- score + (alpha.wk-alpha)*2*varht*trc
                if (method=="v") score <- rss/(1-alpha.wk*trc)^2
            }
            if (return.fit) {
                z <- .Fortran("reg",
                          as.double(cbind(s,10^theta*r)), as.integer(nobs), as.integer(nnull),
                          as.double(q.wk), as.integer(nxiz), as.double(y),
                          as.integer(switch(method,"u"=1,"v"=2,"m"=3)),
                          as.double(alpha), varht=as.double(varht),
                          score=double(1), dc=double(nn),
                          as.double(.Machine$double.eps),
                          chol=double(nn*nn), double(nn),
                          jpvt=as.integer(c(rep(1,nnull),rep(0,nxiz))),
                          wk=double(nobs+nnull+nz), rkv=integer(1), info=integer(1),
                          PACKAGE="gss")[c("score","varht","dc","chol","jpvt","wk","rkv","info")]
                z$score <- score
                assign("fit",z[c(1:5,7)],inherit=TRUE)
            }
        }
        else {
            z <- .Fortran("reg",
                          as.double(cbind(s,10^theta*r)), as.integer(nobs), as.integer(nnull),
                          as.double(q.wk), as.integer(nxiz), as.double(y),
                          as.integer(switch(method,"u"=1,"v"=2,"m"=3)),
                          as.double(alpha), varht=as.double(varht),
                          score=double(1), dc=double(nn),
                          as.double(.Machine$double.eps),
                          chol=double(nn*nn), double(nn),
                          jpvt=as.integer(c(rep(1,nnull),rep(0,nxiz))),
                          wk=double(nobs+nnull+nz), rkv=integer(1), info=integer(1),
                          PACKAGE="gss")[c("score","varht","dc","chol","jpvt","wk","rkv","info")]
            if (z$info) stop("gss error in ssanova: evaluation of GML score fails")
            assign("fit",z[c(1:5,7)],inherit=TRUE)
            score <- z$score
            alpha.wk <- max(0,log.la0-lambda[1]-5)*(3-alpha) + alpha
            alpha.wk <- min(alpha.wk,3)
            if (alpha.wk>alpha) {
                if (method=="u") score <- score + (alpha.wk-alpha)*2*varht*z$wk[2]
                if (method=="v") score <- z$wk[1]/(1-alpha.wk*z$wk[2])^2
            }
        }
        score
    }
    cv.wk <- function(lambda) cv.scale*cv(lambda)+cv.shift
    ## initialization
    tmp <- sum(r^2)
    if (is.null(s)) theta <- 0
    else theta <- log10(sum(s^2)/nnull/tmp*nxi) / 2
    log.la0 <- log10(tmp/sum(diag(q))) + theta
        if (!is.null(random)) {
        ran.scal <- theta - log10(sum(random$z^2)/nz/tmp*nxi) / 2
        r <- cbind(r,10^(ran.scal-theta)*random$z)
    }
    else ran.scal <- NULL
    ## lambda search
    return.fit <- FALSE
    fit <- NULL
    if (is.null(random)) la <- log.la0
    else la <- c(log.la0,random$init)
    if (length(la)-1) {
        counter <- 0
        ## scale and shift cv
        tmp <- abs(cv(la))
        cv.scale <- 1
        cv.shift <- 0
        if (tmp<1&tmp>10^(-4)) {
            cv.scale <- 10/tmp
            cv.shift <- 0
        }
        if (tmp<10^(-4)) {
            cv.scale <- 10^2
            cv.shift <- 10
        }
        repeat {
            zz <- nlm(cv.wk,la,stepmax=1,ndigit=7)
            if (zz$code<=3) break
            la <- zz$est
            counter <- counter + 1
            if (counter>=5) {
                warning("gss warning in ssanova1: iteration for model selection fails to converge")
                break
            }
        }
    }
    else {
        repeat {
            mn <- la-1
            mx <- la+1
            zz <- nlm0(cv,c(mn,mx))
            if (min(zz$est-mn,mx-zz$est)>=1e-3) break
            else la <- zz$est
        }
    }
    ## return
    return.fit <- TRUE
    jk1 <- cv(zz$est)
    if (is.null(random)) q.wk <- 10^theta*q
    else {
        q.wk <- matrix(0,nxiz,nxiz)
        q.wk[1:nxi,1:nxi] <- 10^theta*q
        q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
            10^(2*ran.scal-zz$est[1])*random$sigma$fun(zz$est[-1],random$sigma$env)
    }
    zzz <- eigen(q.wk,TRUE)
    rkq <- min(fit$rkv-nnull,sum(zzz$val/zzz$val[1]>sqrt(.Machine$double.eps)))
    val <- zzz$val[1:rkq]
    vec <- zzz$vec[,1:rkq,drop=FALSE]
    qinv <- vec%*%diag(1/val,rkq)%*%t(vec)
    se.aux <- t(cbind(s,10^theta*r))%*%(10^theta*r)%*%qinv
    c <- fit$dc[nnull+(1:nxi)]
    if (nnull) d <- fit$dc[1:nnull]
    else d <- NULL
    if (nz) b <- 10^(ran.scal)*fit$dc[nnull+nxi+(1:nz)]
    else b <- NULL
    c(list(method=method,theta=theta,ran.scal=ran.scal,c=c,d=d,b=b,
           nlambda=zz$est[1],zeta=zz$est[-1]),
      fit[-3],list(qinv=qinv,se.aux=se.aux))
}

## Fit Multiple Smoothing Parameter (Gaussian) REGression
mspreg1 <- function(s,r,q,y,method,alpha,varht,random)
{
    qr.trace <- FALSE
    if ((alpha<0)&(method%in%c("u","v"))) qr.trace <- TRUE
    alpha <- abs(alpha)
    ## get dimensions
    nobs <- nrow(r)
    nxi <- ncol(r)
    if (!is.null(s)) {
        if (is.vector(s)) nnull <- 1
        else nnull <- ncol(s)
    }
    else nnull <- 0
    if (!is.null(random)) nz <-ncol(as.matrix(random$z))
    else nz <- 0
    nxiz <- nxi + nz
    nn <- nxiz + nnull
    nq <- dim(q)[3]
    ## cv function
    cv <- function(theta) {
        r.wk <- qq.wk <- 0
        for (i in 1:nq) {
            r.wk <- r.wk + 10^theta[i]*r[,,i]
            qq.wk <- qq.wk + 10^theta[i]*q[,,i]
        }
        if (is.null(random)) q.wk <- 10^nlambda*qq.wk
        else {
            r.wk <- cbind(r.wk,10^(ran.scal)*random$z)
            q.wk <- matrix(0,nxiz,nxiz)
            q.wk[1:nxi,1:nxi] <- 10^nlambda*qq.wk
            q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
                10^(2*ran.scal)*random$sigma$fun(theta[-(1:nq)],random$sigma$env)
        }
        if (qr.trace) {
            qq.wk <- chol(q.wk,pivot=TRUE)
            sr <- cbind(s,r.wk[,attr(qq.wk,"pivot")])
            sr <- rbind(sr,cbind(matrix(0,nxiz,nnull),qq.wk))
            sr <- qr(sr,tol=0)
            rss <- mean(qr.resid(sr,c(y,rep(0,nxiz)))[1:nobs]^2)
            trc <- sum(qr.Q(sr)[1:nobs,]^2)/nobs
            if (method=="u") score <- rss + alpha*2*varht*trc
            if (method=="v") score <- rss/(1-alpha*trc)^2
            alpha.wk <- max(0,theta[1:nq]-log.th0-5)*(3-alpha) + alpha
            alpha.wk <- min(alpha.wk,3)
            if (alpha.wk>alpha) {
                if (method=="u") score <- score + (alpha.wk-alpha)*2*varht*trc
                if (method=="v") score <- rss/(1-alpha.wk*trc)^2
            }
            if (return.fit) {
                z <- .Fortran("reg",
                          as.double(cbind(s,r.wk)), as.integer(nobs), as.integer(nnull),
                          as.double(q.wk), as.integer(nxiz), as.double(y),
                          as.integer(switch(method,"u"=1,"v"=2,"m"=3)),
                          as.double(alpha), varht=as.double(varht),
                          score=double(1), dc=double(nn),
                          as.double(.Machine$double.eps),
                          chol=double(nn*nn), double(nn),
                          jpvt=as.integer(c(rep(1,nnull),rep(0,nxiz))),
                          wk=double(nobs+nnull+nz), rkv=integer(1), info=integer(1),
                          PACKAGE="gss")[c("score","varht","dc","chol","jpvt","wk","rkv","info")]
                z$score <- score
                assign("fit",z[c(1:5,7)],inherit=TRUE)
            }
        }
        else {
            z <- .Fortran("reg",
                          as.double(cbind(s,r.wk)), as.integer(nobs), as.integer(nnull),
                          as.double(q.wk), as.integer(nxiz), as.double(y),
                          as.integer(switch(method,"u"=1,"v"=2,"m"=3)),
                          as.double(alpha), varht=as.double(varht),
                          score=double(1), dc=double(nn),
                          as.double(.Machine$double.eps),
                          chol=double(nn*nn), double(nn),
                          jpvt=as.integer(c(rep(1,nnull),rep(0,nxiz))),
                          wk=double(nobs+nnull+nz), rkv=integer(1), info=integer(1),
                          PACKAGE="gss")[c("score","varht","dc","chol","jpvt","wk","rkv","info")]
            if (z$info) stop("gss error in ssanova: evaluation of GML score fails")
            assign("fit",z[c(1:5,7)],inherit=TRUE)
            score <- z$score
            alpha.wk <- max(0,theta[1:nq]-log.th0-5)*(3-alpha) + alpha
            alpha.wk <- min(alpha.wk,3)
            if (alpha.wk>alpha) {
                if (method=="u") score <- score + (alpha.wk-alpha)*2*varht*z$wk[2]
                if (method=="v") score <- z$wk[1]/(1-alpha.wk*z$wk[2])^2
            }
        }
        score
    }
    cv.wk <- function(theta) cv.scale*cv(theta)+cv.shift
    ## initialization
    theta <- -log10(apply(q,3,function(x)sum(diag(x))))
    r.wk <- q.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
    }
    ## theta adjustment
    return.fit <- FALSE
    z <- sspreg1(s,r.wk,q.wk,y,method,alpha,varht,random)
    theta <- theta + z$theta
    r.wk <- q.wk <- 0
    for (i in 1:nq) {
        theta[i] <- 2*theta[i] + log10(t(z$c)%*%q[,,i]%*%z$c)
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
    }
    log.la0 <- log10(sum(r.wk^2)/sum(diag(q.wk)))
    log.th0 <- theta-log.la0
    ## lambda search
    z <- sspreg1(s,r.wk,q.wk,y,method,alpha,varht,random)
    nlambda <- z$nlambda
    log.th0 <- log.th0 + z$lambda
    theta <- theta + z$theta
    if (!is.null(random)) ran.scal <- z$ran.scal
    ## theta search
    fit <- NULL
    if (!is.null(random)) theta <- c(theta,z$zeta)
    counter <- 0
    ## scale and shift cv
    tmp <- abs(cv(theta))
    cv.scale <- 1
    cv.shift <- 0
    if (tmp<1&tmp>10^(-4)) {
        cv.scale <- 10/tmp
        cv.shift <- 0
    }
    if (tmp<10^(-4)) {
        cv.scale <- 10^2
        cv.shift <- 10
    }
    repeat {
        zz <- nlm(cv.wk,theta,stepmax=1,ndigit=7)
        if (zz$code<=3)  break
        theta <- zz$est        
        counter <- counter + 1
        if (counter>=5) {
            warning("gss warning in ssanova1: iteration for model selection fails to converge")
            break
        }
    }
    ## return
    return.fit <- TRUE
    jk1 <- cv(zz$est)
    r.wk <- qq.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^zz$est[i]*r[,,i]
        qq.wk <- qq.wk + 10^zz$est[i]*q[,,i]
    }
    if (is.null(random)) q.wk <- qq.wk
    else {
        r.wk <- cbind(r.wk,10^(ran.scal)*random$z)
        q.wk <- matrix(0,nxiz,nxiz)
        q.wk[1:nxi,1:nxi] <- qq.wk
        q.wk[(nxi+1):nxiz,(nxi+1):nxiz] <-
            10^(2*ran.scal-nlambda)*random$sigma$fun(zz$est[-(1:nq)],random$sigma$env)
    }
    zzz <- eigen(q.wk,TRUE)
    rkq <- min(fit$rkv-nnull,sum(zzz$val/zzz$val[1]>sqrt(.Machine$double.eps)))
    val <- zzz$val[1:rkq]
    vec <- zzz$vec[,1:rkq,drop=FALSE]
    qinv <- vec%*%diag(1/val,rkq)%*%t(vec)
    se.aux <- t(cbind(s,r.wk))%*%r.wk%*%qinv
    c <- fit$dc[nnull+(1:nxi)]
    if (nnull) d <- fit$dc[1:nnull]
    else d <- NULL
    if (nz) b <- 10^(ran.scal)*fit$dc[nnull+nxi+(1:nz)]
    else b <- NULL
    c(list(method=method,theta=zz$est[1:nq],c=c,d=d,b=b,nlambda=nlambda,zeta=zz$est[-(1:nq)]),
      fit[-3],list(qinv=qinv,se.aux=se.aux))
}
smolyak.quad <- ## Generate delayed Smolyak cubature
function(d, k) {
    size <- .C("size_smolyak",
               as.integer(d),
               as.integer(d+k),
               size=integer(1),
               PACKAGE="gss")$size
    z <- .C("quad_smolyak",
            as.integer(d),
            as.integer(d+k),
            pt=double(d*size),
            wt=as.double(1:size),
            PACKAGE="gss")
    list(pt=t(matrix(z$pt,d,size)),wt=z$wt)
}

smolyak.size <- ## Get the size of delayed Smolyak cubature
function(d, k) {
    .C("size_smolyak",
       as.integer(d),
       as.integer(d+k),
       size=integer(1),
       PACKAGE="gss")$size
}
## Fit ssanova model
ssanova <- function(formula,type="cubic",data=list(),
                    weights,subset,offset,na.action=na.omit,
                    partial=NULL,method="v",varht=1,
                    prec=1e-7,maxiter=30,ext=.05,order=2)
{
    ## Obtain model frame and model terms
    mf <- match.call()
    mf$type <- mf$method <- mf$varht <- mf$partial <- NULL
    mf$prec <- mf$maxiter <- mf$ext <- mf$order <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf,sys.frame(sys.parent()))
    if (type=="cubic") term <- mkterm.cubic(mf,ext)
    if (type=="linear") term <- mkterm.linear(mf,ext)
    if (type=="tp") term <- mkterm.tp(mf,order,mf,1)
    if (is.null(term)) stop("gss error in ssanova: unknown type")
    ## Generate s, q, and y
    nobs <- dim(mf)[1]
    s <- q <- NULL
    nq <- 0
    for (label in term$labels) {
        if (label=="1") {
            s <- cbind(s,rep(1,len=nobs))
            next
        }
        x <- mf[,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            phi <- term[[label]]$phi
            for (i in 1:nphi)
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
        }
        if (nrk) {
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq+1
                q <- array(c(q,rk$fun(x,x,nu=i,env=rk$env,out=TRUE)),c(nobs,nobs,nq))
            }
        }
    }
    ## Add the partial term
    if (!is.null(partial)) {
        if (is.vector(partial)) partial <- as.matrix(partial)
        if (dim(partial)[1]!=dim(mf)[1])
            stop("gss error in ssanova: partial data are of wrong size")
        term$labels <- c(term$labels,"partial")
        term$partial <- list(nphi=dim(partial)[2],nrk=0,
                             iphi=ifelse(is.null(s),0,dim(s)[2])+1)
        s <- cbind(s,partial)
        mf$partial <- partial
    }
    ## Prepare the data
    y <- model.response(mf,"numeric")
    w <- model.weights(mf)
    offset <- model.offset(mf)
    if (!is.null(offset)) {
        term$labels <- c(term$labels,"offset")
        term$offset <- list(nphi=0,nrk=0)
        y <- y - offset
    }
    if (!is.null(w)) {
        w <- sqrt(w)
        y <- w*y
        s <- w*s
        for (i in 1:nq) q[,,i] <- w*t(w*q[,,i])
    }
    if (qr(s)$rank<dim(s)[2])
        stop("gss error in ssanova: fixed effects are linearly dependent")
    if (!nq) stop("gss error in ssanova: use lm for models with only fixed effects")
    ## Fit the model
    if (nq==1) {
        q <- q[,,1]
        z <- sspreg(s,q,y,method,varht)
    }
    else z <- mspreg(s,q,y,method,varht,prec,maxiter)
    ## Brief description of model terms
    desc <- NULL
    for (label in term$labels)
        desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
    desc <- rbind(desc,apply(desc,2,sum))
    rownames(desc) <- c(term$labels,"total")
    colnames(desc) <- c("Unpenalized","Penalized")
    ## Return the results
    obj <- c(list(call=match.call(),mf=mf,terms=term,desc=desc),z)
    class(obj) <- c("ssanova")
    obj
}
## Fit ssanova model
ssanova1 <- function(formula,type="cubic",data=list(),
                       weights,subset,offset,na.action=na.omit,
                       partial=NULL,method="v",alpha=1.4,varht=1,
                       id.basis=NULL,nbasis=NULL,seed=NULL,random=NULL,
                       ext=.05,order=2)
{
    ## Obtain model frame and model terms
    mf <- match.call()
    mf$type <- mf$method <- mf$varht <- mf$partial <- NULL
    mf$alpha <- mf$id.basis <- mf$nbasis <- mf$seed <- NULL
    mf$random <- mf$ext <- mf$order <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf,sys.frame(sys.parent()))
    wt <- model.weights(mf)
    ## Generate sub-basis
    nobs <- dim(mf)[1]
    if (is.null(id.basis)) {
        if (is.null(nbasis))  nbasis <- max(30,ceiling(10*nobs^(2/9)))
        if (nbasis>=nobs)  nbasis <- nobs
        if (!is.null(seed))  set.seed(seed)
        id.basis <- sample(nobs,nbasis,prob=wt)
    }
    else {
        if (max(id.basis)>nobs|min(id.basis)<1)
            stop("gss error in ssanova1: id.basis out of range")
        nbasis <- length(id.basis)
    }
    ## Generate terms
    if (type=="cubic") term <- mkterm.cubic(mf,ext)
    if (type=="linear") term <- mkterm.linear(mf,ext)
    if (type=="tp") term <- mkterm.tp(mf,order,mf[id.basis,],1)
    if (is.null(term)) stop("gss error in ssanova1: unknown type")
    ## Generate random
    if (!is.null(random)) {
        if (class(random)=="formula") random <- mkran(random,data)
    }
    ## Generate s, r, q, and y
    s <- r <- NULL
    nq <- 0
    for (label in term$labels) {
        if (label=="1") {
            s <- cbind(s,rep(1,len=nobs))
            next
        }
        x <- mf[,term[[label]]$vlist]
        x.basis <- mf[id.basis,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            phi <- term[[label]]$phi
            for (i in 1:nphi)
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
        }
        if (nrk) {
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq+1
                r <- array(c(r,rk$fun(x,x.basis,nu=i,env=rk$env,out=TRUE)),c(nobs,nbasis,nq))
            }
        }
    }
    if (is.null(r))
        stop("gss error in ssanova1: use lm for models with only fixed effects")
    else q <- r[id.basis,,,drop=FALSE]
    ## Add the partial term
    if (!is.null(partial)) {
        if (is.vector(partial)) partial <- as.matrix(partial)
        if (dim(partial)[1]!=dim(mf)[1])
            stop("gss error in ssanova1: partial data are of wrong size")
        term$labels <- c(term$labels,"partial")
        term$partial <- list(nphi=dim(partial)[2],nrk=0,
                             iphi=ifelse(is.null(s),0,dim(s)[2])+1)
        s <- cbind(s,partial)
        mf$partial <- partial
    }
    if (qr(s)$rank<dim(s)[2])
        stop("gss error in ssanova1: fixed effects are linearly dependent")
    ## Prepare the data
    y <- model.response(mf,"numeric")
    offset <- model.offset(mf)
    if (!is.null(offset)) {
        term$labels <- c(term$labels,"offset")
        term$offset <- list(nphi=0,nrk=0)
        y <- y - offset
    }
    if (!is.null(wt)) {
        wt <- sqrt(wt)
        y <- wt*y
        s <- wt*s
        r <- wt*r
        if (!is.null(random)) random$z <- wt*random$z
    }
    if (qr(s)$rank<dim(s)[2])
        stop("gss error in ssanova1: fixed effects are linearly dependent")
    ## Fit the model
    if (nq==1) {
        r <- r[,,1]
        q <- q[,,1]
        z <- sspreg1(s,r,q,y,method,alpha,varht,random)
    }
    else z <- mspreg1(s,r,q,y,method,alpha,varht,random)
    ## Brief description of model terms
    desc <- NULL
    for (label in term$labels)
        desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
    desc <- rbind(desc,apply(desc,2,sum))
    rownames(desc) <- c(term$labels,"total")
    colnames(desc) <- c("Unpenalized","Penalized")
    ## Return the results
    obj <- c(list(call=match.call(),mf=mf,terms=term,desc=desc,
                  alpha=alpha,id.basis=id.basis,random=random),z)
    class(obj) <- c("ssanova1","ssanova")
    obj
}
## Fit density model
ssden <- function(formula,type="cubic",data=list(),alpha=1.4,
                  weights=NULL,subset,na.action=na.omit,
                  id.basis=NULL,nbasis=NULL,seed=NULL,
                  domain=as.list(NULL),quadrature=NULL,ext=.05,order=2,
                  prec=1e-7,maxiter=30)
{
    ## Obtain model frame and model terms
    mf <- match.call()
    mf$type <- mf$alpha <- NULL
    mf$id.basis <- mf$nbasis <- mf$seed <- NULL
    mf$domain <- mf$quadrature <- mf$ext  <- NULL
    mf$prec <- mf$maxiter <- mf$order <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf,sys.frame(sys.parent()))
    cnt <- model.weights(mf)
    mf$"(weights)" <- NULL
    ## set domain
    for (i in names(mf)) {
        if (is.factor(mf[[i]])) domain[[i]] <- levels(mf[[i]])[1:2]
        else {
            if (is.null(domain[[i]])) {
                mn <- min(mf[[i]])
                mx <- max(mf[[i]])
                range <- mx-mn
                mn <- mn - ext*range
                mx <- mx + ext*range
                domain[[i]] <- c(mn,mx)
            }
        }
    }
    domain <- as.data.frame(domain)
    ## Generate sub-basis
    nobs <- dim(mf)[1]
    if (is.null(id.basis)) {
        if (is.null(nbasis))  nbasis <- max(30,ceiling(10*nobs^(2/9)))
        if (nbasis>=nobs)  nbasis <- nobs
        if (!is.null(seed))  set.seed(seed)
        id.basis <- sample(nobs,nbasis,prob=cnt)
    }
    else {
        if (max(id.basis)>nobs|min(id.basis)<1)
            stop("gss error in ssden: id.basis out of range")
        nbasis <- length(id.basis)
    }
    ## Generate terms
    if (type=="cubic") term <- mkterm.cubic1(mf,domain)
    if (type=="linear") term <- mkterm.linear1(mf,domain)
    if (type=="tp") {
        if (is.null(quadrature))
            stop("gss error in ssden: quadrature needed for type tp")
        term <- mkterm.tp(mf,order,mf[id.basis,],1)
    }
    if (is.null(term)) stop("gss error in ssden: unknown type")
    term$labels <- term$labels[term$labels!="1"]
    ## Generate default quadrature
    if (is.null(quadrature)) {
        ## TO DO: HANDLING OF FACTORS
        domain <- domain[,colnames(mf),drop=FALSE]
        mn <- apply(domain,2,min)
        mx <- apply(domain,2,max)
        if (ncol(mf)==1) {
            ## Gauss-Legendre quadrature
            quad <- gauss.quad(200,c(mn,mx))
            quad$pt <- data.frame(quad$pt)
            colnames(quad$pt) <- colnames(mf)
        }
        else {
            ## Smolyak cubature
            if (ncol(mf)>4)
                stop("gss error in ssden: dimension higher than 4 unsupported")
            code <- c(15,14,13)
            quad <- smolyak.quad(ncol(mf),code[ncol(mf)-1])
            for (i in 1:ncol(mf)) {
                wk <- mf[,i]
                jk <- ssden(~wk,domain=data.frame(wk=domain[,i]),alpha=2,
                            id.basis=id.basis)
                quad$pt[,i] <- qssden(jk,quad$pt[,i])
                quad$wt <- quad$wt/dssden(jk,quad$pt[,i])
            }
            jk <- wk <- NULL
            quad$pt <- data.frame(quad$pt)
            colnames(quad$pt) <- colnames(mf)
        }
        quadrature <- list(pt=quad$pt,wt=quad$wt)
    }
    ## Generate s, r, and q
    qd.pt <- quadrature$pt
    qd.wt <- quadrature$wt
    nmesh <- length(qd.wt)
    s <- qd.s <- r <- qd.r <- q <- NULL
    nq <- 0
    for (label in term$labels) {
        x <- mf[,term[[label]]$vlist]
        x.basis <- mf[id.basis,term[[label]]$vlist]
        qd.x <- qd.pt[,term[[label]]$vlist]
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            phi <- term[[label]]$phi
            for (i in 1:nphi) {
                s <- cbind(s,phi$fun(x,nu=i,env=phi$env))
                qd.s <- cbind(qd.s,phi$fun(qd.x,nu=i,env=phi$env))
            }
        }
        if (nrk) {
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq+1
                r <- array(c(r,rk$fun(x.basis,x,nu=i,env=rk$env,out=TRUE)),c(nbasis,nobs,nq))
                qd.r <- array(c(qd.r,rk$fun(x.basis,qd.x,nu=i,env=rk$env,out=TRUE)),
                              c(nbasis,nmesh,nq))
                q <- array(c(q,rk$fun(x.basis,x.basis,nu=i,env=rk$env,out=TRUE)),
                           c(nbasis,nbasis,nq))
            }
        }
    }
    if (!is.null(s)) {
        nnull <- dim(s)[2]
        ## Check s rank
        if (qr(s)$rank<nnull)
            stop("gss error in ssden: fixed effect MLE is not unique")
        s <- t(s)
        qd.s <- t(qd.s)
    }
    ## Fit the model
    if (nq==1) {
        r <- r[,,1]
        qd.r <- qd.r[,,1]
        q <- q[,,1]
        z <- sspdsty(s,r,q,cnt,qd.s,qd.r,qd.wt,prec,maxiter,alpha)
    }
    else z <- mspdsty(s,r,q,cnt,qd.s,qd.r,qd.wt,prec,maxiter,alpha)
    ## Brief description of model terms
    desc <- NULL
    for (label in term$labels)
        desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
    desc <- rbind(desc,apply(desc,2,sum))
    rownames(desc) <- c(term$labels,"total")
    colnames(desc) <- c("Unpenalized","Penalized")
    ## Return the results
    obj <- c(list(call=match.call(),mf=mf,cnt=cnt,terms=term,desc=desc,alpha=alpha,
                  domain=domain,quad=quadrature,id.basis=id.basis),z)
    class(obj) <- "ssden"
    obj
}
## Fit hazard model
sshzd <- function(formula,type="cubic",data=list(),alpha=1.4,
                  weights=NULL,subset,na.action=na.omit,
                  id.basis=NULL,nbasis=NULL,seed=NULL,
                  ext=.05,order=2,prec=1e-7,maxiter=30)
{
    ## Local functions handling formula
    Surv <- function(time,status,start=0) {
        tname <- as.character(as.list(match.call())$time)
        if (!is.numeric(time)|!is.vector(time))
            stop("gss error in sshzd: time should be a numerical vector")
        if ((nobs <- length(time))-length(status))
            stop("gss error in sshzd: time and status mismatch in size")
        if ((length(start)-nobs)&(length(start)-1))
            stop("gss error in sshzd: time and start mismatch in size")
        if (any(start>time))
            stop("gss error in sshzd: start after follow-up time")
        if (min(start)<0)
            warning("gss warning in sshzd: start before time 0")
        time <- cbind(start,time)
        list(tname=tname,start=time[,1],end=time[,2],status=as.logical(status))
    }
    ## Obtain model frame and model terms
    mf <- match.call()
    mf$type <- mf$alpha <- NULL
    mf$id.basis <- mf$nbasis <- mf$seed <- NULL
    mf$ext <- mf$order <- mf$prec <- mf$maxiter <- NULL
    term.wk <- terms.formula(mf$formula)
    ## response
    resp <- attr(term.wk,"variable")[[2]]
    ind.wk <- length(strsplit(deparse(resp),'')[[1]])
    if ((substr(deparse(resp),1,5)!='Surv(')
        |(substr(deparse(resp),ind.wk,ind.wk)!=')'))
        stop("gss error in sshzd: response should be Surv(...)")
    attach(data)
    yy <- eval(resp)
    tname <- yy$tname
    ## model frame
    term.labels <- attr(term.wk,"term.labels")
    if (!(tname%in%term.labels))
        stop("gss error in sshzd: time main effect missing in model")
    mf[[1]] <- as.name("model.frame")
    mf[[2]] <- eval(parse(text=paste("~",paste(term.labels,collapse="+"))))
    mf <- eval(mf,sys.frame(sys.parent()))
    ## set domain
    xnames <- names(mf)
    xnames <- xnames[!xnames%in%tname]
    domain <- as.list(NULL)
    mn <- min(yy$start)
    mx <- max(yy$end)
    domain[[tname]] <- c(mn,mx)
    for (i in xnames) {
        if (is.factor(mf[[i]])) domain[[i]] <- levels(mf[[i]])[1:2]
        else {
            mn <- min(mf[[i]])
            mx <- max(mf[[i]])
            range <- mx-mn
            mn <- mn - ext*range
            mx <- mx + ext*range
            domain[[i]] <- c(mn,mx)
        }
    }
    domain <- as.data.frame(domain)
    ## Generate sub-basis
    cnt <- model.weights(mf)
    nobs <- nrow(mf)
    if (is.null(id.basis)) {
        if (is.null(nbasis)) nbasis <- max(30,ceiling(10*nobs^(2/9)))
        if (nbasis>sum(yy$status)) nbasis <- sum(yy$status)
        if (!is.null(seed)) set.seed(seed)
        id.basis <- sample((1:nobs)[yy$status],nbasis,prob=cnt[yy$status])
    }
    else {
        if (!all(id.basis%in%(1:nobs)[yy$status]))
            stop("gss error in sshzd: id.basis not all at failure cases")
        nbasis <- length(id.basis)
    }
    ## Generate terms    
    if (type=="cubic") term <- mkterm.cubic1(mf,domain)
    if (type=="linear") term <- mkterm.linear1(mf,domain)
    if (type=="tp") term <- mkterm.tp(mf,order,mf[id.basis,],1)
    if (is.null(term)) stop("gss error in sshzd: unknown type")
    ## Generate Gauss-Legendre quadrature
    nmesh <- 200
    quad <- gauss.quad(nmesh,domain[[tname]])
    ## obtain unique covariate observations
    if (length(xnames)) {
        xx <- mf[,xnames,drop=FALSE]
        x.pt <- unique(xx)
        nx <- dim(x.pt)[1]
        x.dup.ind <- duplicated(xx)
        x.dup <- xx[x.dup.ind,,drop=FALSE]
        ## xx[i,]==x.pt[x.ind[i],]
        x.ind <- 1:nobs
        x.ind[!x.dup.ind] <- 1:nx
        if (nobs-nx) {
            x.ind.wk <- 1:(nobs-nx)
            for (i in 1:nx) {
                for (j in 1:(nobs-nx)) {
                    if (sum(duplicated(rbind(x.pt[i,],x.dup[j,]))))
                        x.ind.wk[j] <- i
                }
            }
        }
        if (nobs-nx) x.ind[x.dup.ind] <- x.ind.wk
    }
    else {
        nx <- 1
        x.ind <- rep(1,nobs)
        x.pt <- NULL
    }
    ## integration weights at x.pt[i,]
    qd.wt <- matrix(0,nmesh,nx)
    for (i in 1:nobs) {
        wk <- (quad$pt<=yy$end[i])&(quad$pt>yy$start[i])
        if (is.null(cnt)) qd.wt[,x.ind[i]] <- qd.wt[,x.ind[i]]+wk
        else qd.wt[,x.ind[i]] <- qd.wt[,x.ind[i]]+cnt[i]*wk
    }
    if (is.null(cnt)) qd.wt <- quad$wt*qd.wt/nobs
    else qd.wt <- quad$wt*qd.wt/sum(cnt)
    ## Generate s, r, and q
    s <- r <- q <- qd.s <- NULL
    qd.r <- as.list(NULL)
    nT <- sum(yy$status)
    nq <- nu <- 0
    for (label in term$labels) {
        if (label=="1") {
            nu <- nu+1
            s <- cbind(s,rep(1,len=nT))
            qd.wk <- matrix(1,nmesh,nx)
            qd.s <- array(c(qd.s,qd.wk),c(nmesh,nx,nu))
            next
        }
        vlist <- term[[label]]$vlist
        x.list <- xnames[xnames%in%vlist]
        xy <- mf[yy$status,vlist]
        xy.basis <- mf[id.basis,vlist]
        qd.xy <- data.frame(matrix(0,nmesh,length(vlist)))
        names(qd.xy) <- vlist
        if (tname%in%vlist) qd.xy[,tname] <- quad$pt
        if (length(x.list)) xx <- x.pt[,x.list,drop=FALSE]
        else xx <- NULL
        nphi <- term[[label]]$nphi
        nrk <- term[[label]]$nrk
        if (nphi) {
            phi <- term[[label]]$phi
            for (i in 1:nphi) {
                nu <- nu+1
                s <- cbind(s,phi$fun(xy,nu=i,env=phi$env))
                if (is.null(xx))
                    qd.wk <- matrix(phi$fun(qd.xy[,,drop=TRUE],nu=i,env=phi$env),nmesh,nx)
                else {
                    qd.wk <- NULL
                    for (j in 1:nx) {
                        qd.xy[,x.list] <- xx[rep(j,nmesh),]
                        qd.wk <- cbind(qd.wk,phi$fun(qd.xy[,,drop=TRUE],i,phi$env))
                    }
                }
                qd.s <- array(c(qd.s,qd.wk),c(nmesh,nx,nu))
            }
        }
        if (nrk) {
            rk <- term[[label]]$rk
            for (i in 1:nrk) {
                nq <- nq+1
                r <- array(c(r,rk$fun(xy,xy.basis,nu=i,env=rk$env,out=TRUE)),c(nT,nbasis,nq))
                q <- array(c(q,rk$fun(xy.basis,xy.basis,nu=i,env=rk$env,out=TRUE)),
                           c(nbasis,nbasis,nq))
                if (is.null(xx))
                    qd.r[[nq]] <- rk$fun(qd.xy[,,drop=TRUE],xy.basis,i,rk$env,out=TRUE)
                else {
                    qd.wk <- NULL
                    for (j in 1:nx) {
                        qd.xy[,x.list] <- xx[rep(j,nmesh),]
                        qd.wk <- array(c(qd.wk,rk$fun(qd.xy[,,drop=TRUE],xy.basis,i,rk$env,TRUE)),
                                       c(nmesh,nbasis,j))
                    }
                    qd.r[[nq]] <- qd.wk
                }
            }
        }
    }
    if (!is.null(s)) {
        nnull <- dim(s)[2]
        ## Check s rank
        if (qr(s)$rank<nnull)
            stop("gss error in sscden: fixed effect MLE is not unique")
    }
    ## Fit the model
    Nobs <- ifelse(is.null(cnt),nobs,sum(cnt))
    if (!is.null(cnt)) cntt <- cnt[yy$status]
    else cntt <- NULL
    z <- msphzd(s,r,q,Nobs,cntt,qd.s,qd.r,qd.wt,prec,maxiter,alpha)
    cfit <- sum(yy$status)/Nobs/sum(qd.wt)
    ## Brief description of model terms
    desc <- NULL
    for (label in term$labels)
        desc <- rbind(desc,as.numeric(c(term[[label]][c("nphi","nrk")])))
    desc <- rbind(desc,apply(desc,2,sum))
    rownames(desc) <- c(term$labels,"total")
    colnames(desc) <- c("Unpenalized","Penalized")
    ## Return the results
    obj <- c(list(call=match.call(),mf=mf,tname=tname,xnames=xnames,
                  terms=term,desc=desc,alpha=alpha,domain=domain,cfit=cfit,
                  quad=quad,x.pt=x.pt,qd.wt=qd.wt,id.basis=id.basis),z)
    if (is.null(cnt)) obj$se.aux$v <- sqrt(nobs)*obj$se.aux$v
    else obj$se.aux$v <- sqrt(sum(cnt))*obj$se.aux$v
    class(obj) <- c("sshzd")
    obj
}

## Fit (multiple smoothing parameter) hazard function
msphzd <- function(s,r,q,Nobs,cnt,qd.s,qd.r,qd.wt,prec,maxiter,alpha)
{
    nT <- dim(r)[1]
    nxi <- dim(r)[2]
    nqd <- dim(qd.wt)[1]
    nx <- dim(qd.wt)[2]
    if (!is.null(s)) nnull <- dim(s)[2]
    else nnull <- 0
    nxis <- nxi+nnull
    if (is.null(cnt)) cnt <- 0
    ## cv functions
    cv.s <- function(lambda) {
        fit <- .Fortran("hzdnewton",
                        cd=as.double(cd), as.integer(nxis),
                        as.double(10^lambda*q.wk), as.integer(nxi),
                        as.double(t(cbind(r.wk,s))), as.integer(nT),
                        as.integer(Nobs), as.integer(sum(cnt)), as.integer(cnt),
                        as.double(qd.r.wk), as.integer(nqd),
                        as.double(qd.wt), as.integer(nx),
                        as.double(prec), as.integer(maxiter),
                        as.double(.Machine$double.eps),
                        wk=double(2*(nqd*nx+nT)+nxis*(2*nxis+5)+max(nxis,2)),
                        info=integer(1),PACKAGE="gss")
        if (fit$info==1) stop("gss error in sshzd: Newton iteration diverges")
        if (fit$info==2) warning("gss warning in sshzd: Newton iteration fails to converge")
        assign("cd",fit$cd,inherit=TRUE)
        assign("mesh0",matrix(fit$wk[max(nxis,2)+(1:(nqd*nx))],nqd,nx),inherit=TRUE)
        cv <- alpha*fit$wk[2]-fit$wk[1]
        alpha.wk <- max(0,log.la0-lambda-5)*(3-alpha) + alpha
        alpha.wk <- min(alpha.wk,3)
        adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*fit$wk[2],0)
        cv+adj
    }
    cv.m <- function(theta) {
        r.wk <- q.wk <- 0
        qd.r.wk <- array(0,c(nqd,nxi,nx))
        for (i in 1:nq) {
            r.wk <- r.wk + 10^theta[i]*r[,,i]
            q.wk <- q.wk + 10^theta[i]*q[,,i]
            if (length(dim(qd.r[[i]]))==3) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[[i]]
            else qd.r.wk <- qd.r.wk + as.vector(10^theta[i]*qd.r[[i]])
        }
        qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
        qd.r.wk <- array(c(qd.r.wk,qd.s),c(nqd,nx,nxis))
        qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
        fit <- .Fortran("hzdnewton",
                        cd=as.double(cd), as.integer(nxis),
                        as.double(10^lambda*q.wk), as.integer(nxi),
                        as.double(t(cbind(r.wk,s))), as.integer(nT),
                        as.integer(Nobs), as.integer(sum(cnt)), as.integer(cnt),
                        as.double(qd.r.wk), as.integer(nqd),
                        as.double(qd.wt), as.integer(nx),
                        as.double(prec), as.integer(maxiter),
                        as.double(.Machine$double.eps),
                        wk=double(2*(nqd*nx+nT)+nxis*(2*nxis+5)+max(nxis,2)),
                        info=integer(1),PACKAGE="gss")
        if (fit$info==1) stop("gss error in sshzd: Newton iteration diverges")
        if (fit$info==2) warning("gss warning in sshzd: Newton iteration fails to converge")
        assign("cd",fit$cd,inherit=TRUE)
        assign("mesh0",matrix(fit$wk[max(nxis,2)+(1:(nqd*nx))],nqd,nx),inherit=TRUE)
        cv <- alpha*fit$wk[2]-fit$wk[1]
        alpha.wk <- max(0,theta-log.th0-5)*(3-alpha) + alpha
        alpha.wk <- min(alpha.wk,3)
        adj <- ifelse (alpha.wk>alpha,(alpha.wk-alpha)*fit$wk[2],0)
        cv+adj
    }
    cv.wk <- function(theta) cv.scale*cv.m(theta)+cv.shift
    ## Initialization
    theta <- -log10(apply(q,3,function(x)sum(diag(x))))
    nq <- length(theta)
    qd.r.wk <- array(0,c(nqd,nxi,nx))
    for (i in 1:nq) {
        if (length(dim(qd.r[[i]]))==3) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[[i]]
        else qd.r.wk <- qd.r.wk + as.vector(10^theta[i]*qd.r[[i]])
    }
    v.s <- v.r <- 0
    for (i in 1:nx) {
        if (nnull) v.s <- v.s + apply(qd.wt[,i]*qd.s[,i,,drop=FALSE]^2,2,sum)
        v.r <- v.r + apply(qd.wt[,i]*qd.r.wk[,,i,drop=FALSE]^2,2,sum)
    }
    if (nnull) theta.wk <- log10(sum(v.s)/nnull/sum(v.r)*nxi) / 2
    else theta.wk <- 0
    theta <- theta + theta.wk
    qd.r.wk <- aperm(10^theta.wk*qd.r.wk,c(1,3,2))
    qd.r.wk <- array(c(qd.r.wk,qd.s),c(nqd,nx,nxis))
    qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
    r.wk <- q.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
    }
    log.la0 <- log10(sum(v.r)/sum(diag(q.wk))) + 2*theta.wk
    ## fixed theta iteration
    mesh0 <- NULL
    cd <- rep(0,nxi+nnull)
    la <- log.la0
    repeat {
        mn <- la-1
        mx <- la+1
        zz <- nlm0(cv.s,c(mn,mx))
        if (min(zz$est-mn,mx-zz$est)>=1e-3) break
        else la <- zz$est
    }
    if (nq==1) {
        lambda <- zz$est
        se.aux <- .Fortran("hzdaux1",
                           as.double(cd), as.integer(nxis),
                           as.double(10^lambda*q.wk), as.integer(nxi),
                           as.double(qd.r.wk), as.integer(nqd),
                           as.double(qd.wt), as.integer(nx),
                           as.double(.Machine$double.eps), double(nqd*nx),
                           v=double(nxis*nxis), double(nxis*nxis),
                           jpvt=integer(nxis), PACKAGE="gss")[c("v","jpvt")]
        c <- cd[1:nxi]
        if (nnull) d <- cd[nxi+(1:nnull)]
        return(list(lambda=zz$est,theta=theta,c=c,d=d,cv=zz$min,mesh0=mesh0,se.aux=se.aux))
    }
    ## theta adjustment
    qd.r.wk <- array(0,c(nqd,nxi,nx))
    for (i in 1:nq) {
        theta[i] <- 2*theta[i] + log10(t(cd[1:nxi])%*%q[,,i]%*%cd[1:nxi])
        if (length(dim(qd.r[[i]]))==3) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[[i]]
        else qd.r.wk <- qd.r.wk + as.vector(10^theta[i]*qd.r[[i]])
    }
    v.s <- v.r <- 0
    for (i in 1:nx) {
        if (nnull) v.s <- v.s + apply(qd.wt[,i]*qd.s[,i,,drop=FALSE]^2,2,sum)
        v.r <- v.r + apply(qd.wt[,i]*qd.r.wk[,,i,drop=FALSE]^2,2,sum)
    }
    if (nnull) theta.wk <- log10(sum(v.s)/nnull/sum(v.r)*nxi) / 2
    else theta.wk <- 0
    theta <- theta + theta.wk
    qd.r.wk <- aperm(10^theta.wk*qd.r.wk,c(1,3,2))
    qd.r.wk <- array(c(qd.r.wk,qd.s),c(nqd,nx,nxis))
    qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
    r.wk <- q.wk <- 0
    for (i in 1:nq) {
        r.wk <- r.wk + 10^theta[i]*r[,,i]
        q.wk <- q.wk + 10^theta[i]*q[,,i]
    }
    log.la0 <- log10(sum(v.r)/sum(diag(q.wk))) + 2*theta.wk
    log.th0 <- theta-log.la0
    ## fixed theta iteration
    cd <- rep(0,nxi+nnull)
    la <- log.la0
    repeat {
        mn <- la-1
        mx <- la+1
        zz <- nlm0(cv.s,c(mn,mx))
        if (min(zz$est-mn,mx-zz$est)>=1e-3) break
        else la <- zz$est
    }
    ## theta search
    lambda <- zz$est
    counter <- 0
    tmp <- abs(cv.m(theta))
    cv.scale <- 1
    cv.shift <- 0
    if (tmp<1&tmp>10^(-4)) {
        cv.scale <- 10/tmp
        cv.shift <- 0
    }
    if (tmp<10^(-4)) {
        cv.scale <- 10^2
        cv.shift <- 10
    }
    repeat {
        zz <- nlm(cv.wk,theta,stepmax=1,ndigit=7)
        if (zz$code<=3)  break
        theta <- zz$est
        counter <- counter + 1
        if (counter>=5) {
            warning("gss warning in sscden: CV iteration fails to converge")
            break
        }
    }
    ## return
    theta <- zz$est
    cv <- (zz$min-cv.shift)/cv.scale
    q.wk <- 0
    qd.r.wk <- array(0,c(nqd,nxi,nx))
    for (i in 1:nq) {
        q.wk <- q.wk + 10^theta[i]*q[,,i]
        if (length(dim(qd.r[[i]]))==3) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[[i]]
        else qd.r.wk <- qd.r.wk + as.vector(10^theta[i]*qd.r[[i]])
    }
    qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
    qd.r.wk <- array(c(qd.r.wk,qd.s),c(nqd,nx,nxis))
    qd.r.wk <- aperm(qd.r.wk,c(1,3,2))
    se.aux <- .Fortran("hzdaux1",
                       as.double(cd), as.integer(nxis),
                       as.double(10^lambda*q.wk), as.integer(nxi),
                       as.double(qd.r.wk), as.integer(nqd),
                       as.double(qd.wt), as.integer(nx),
                       as.double(.Machine$double.eps), double(nqd*nx),
                       v=double(nxis*nxis), double(nxis*nxis),
                       jpvt=integer(nxis), PACKAGE="gss")[c("v","jpvt")]
    c <- cd[1:nxi]
    if (nnull) d <- cd[nxi+(1:nnull)]
    list(lambda=lambda,theta=theta,c=c,d=d,cv=cv,mesh0=mesh0,se.aux=se.aux)
}
## Summarize gssanova objects
summary.gssanova <- function(object,diagnostics=FALSE,...)
{
    y <- model.response(object$mf,"numeric")
    wt <- model.weights(object$mf)
    offset <- model.offset(object$mf)
    if ((object$family=="nbinomial")&(!is.null(object$nu))) y <- cbind(y,object$nu)
    dev.resid <- switch(object$family,
                        binomial=dev.resid.binomial(y,object$eta,wt),
                        nbinomial=dev.resid.nbinomial(y,object$eta,wt),
                        poisson=dev.resid.poisson(y,object$eta,wt),
                        inverse.gaussian=dev.resid.inverse.gaussian(y,object$eta,wt),
                        Gamma=dev.resid.Gamma(y,object$eta,wt),
                        weibull=dev.resid.weibull(y,object$eta,wt,object$nu),
                        lognorm=dev.resid.lognorm(y,object$eta,wt,object$nu),
                        loglogis=dev.resid.loglogis(y,object$eta,wt,object$nu))
    dev.null <- switch(object$family,
                       binomial=dev.null.binomial(y,wt,offset),
                       nbinomial=dev.null.nbinomial(y,wt,offset),
                       poisson=dev.null.poisson(y,wt,offset),
                       inverse.gaussian=dev.null.inverse.gaussian(y,wt,offset),
                       Gamma=dev.null.Gamma(y,wt,offset),
                       weibull=dev.null.weibull(y,wt,offset,object$nu),
                       lognorm=dev.null.lognorm(y,wt,offset,object$nu),
                       loglogis=dev.null.loglogis(y,wt,offset,object$nu))
    w <- object$w
    if (is.null(offset)) offset <- rep(0,length(object$eta))
    ## Residuals
    res <- 10^object$nlambda*object$c 
    ## Fitted values
    fitted <- object$eta
    fitted.off <- fitted-offset
    ## dispersion
    sigma2 <- object$varht
    ## RSS, deviance
    rss <- sum(res^2)
    dev <- sum(dev.resid)
    ## Penalty associated with the fit
    penalty <- sum(object$c*fitted.off*sqrt(w))
    penalty <- as.vector(10^object$nlambda*penalty)
    ## Calculate the diagnostics
    if (diagnostics) {
        ## Obtain retrospective linear model
        comp <- NULL
        for (label in object$terms$labels) {
            if (label=="1") next
            if (label=="offset") next
            comp <- cbind(comp,predict(object,object$mf,inc=label))
        }
        comp <- cbind(comp,yhat=fitted.off,y=fitted.off+res/sqrt(w),e=res/sqrt(w))
        term.label <- object$terms$labels[object$terms$labels!="1"]
        term.label <- term.label[term.label!="offset"]
        colnames(comp) <- c(term.label,"yhat","y","e")
        ## Sweep out constant
        comp <- sqrt(w)*comp - outer(sqrt(w),apply(w*comp,2,sum))/sum(w)
        ## Obtain pi
        comp1 <- comp[,c(term.label,"yhat")]
        decom <- t(comp1) %*% comp1[,"yhat"]
        names(decom) <- c(term.label,"yhat")
        decom <- decom[term.label]/decom["yhat"]
        ## Obtain kappa, norm, and cosines        
        corr <- t(comp)%*%comp
        corr <- t(corr/sqrt(diag(corr)))/sqrt(diag(corr))
        norm <- apply(comp,2,function(x){sqrt(sum(x^2))})
        cosines <- rbind(corr[c("y","e"),],norm)
        rownames(cosines) <- c("cos.y","cos.e","norm")
        corr <- corr[term.label,term.label,drop=FALSE]
        if (qr(corr)$rank<dim(corr)[2]) kappa <- rep(Inf,len=dim(corr)[2])
        else kappa <- as.numeric(sqrt(diag(solve(corr))))
        ## Obtain decomposition of penalty
        rough <- as.vector(10^object$nlambda*t(comp[,term.label])%*%object$c/penalty)
        names(kappa) <- names(rough) <- term.label
    }
    else decom <- kappa <- cosines <- rough <- NULL
    ## Return the summaries
    z <- list(call=object$call,family=object$family,method=object$method,iter=object$iter,
              fitted=fitted,dispersion=sigma2,residuals=res/sqrt(w),rss=rss,
              deviance=dev,dev.resid=sqrt(dev.resid)*sign(res),nu=object$nu,
              dev.null=dev.null,penalty=penalty,
              pi=decom,kappa=kappa,cosines=cosines,roughness=rough)
    class(z) <- "summary.gssanova"
    z
}
## Summarize gssanova1 objects
summary.gssanova1 <- function(object,diagnostics=FALSE,...)
{
    y <- model.response(object$mf,"numeric")
    wt <- model.weights(object$mf)
    offset <- model.offset(object$mf)
    if ((object$family=="nbinomial")&(!is.null(object$nu))) y <- cbind(y,object$nu)
    dev.null <- switch(object$family,
                       binomial=dev.null.binomial(y,wt,offset),
                       nbinomial=dev.null.nbinomial(y,wt,offset),
                       poisson=dev.null.poisson(y,wt,offset),
                       Gamma=dev.null.Gamma(y,wt,offset),
                       weibull=dev.null.weibull(y,wt,offset,object$nu),
                       lognorm=dev.null.lognorm(y,wt,offset,object$nu),
                       loglogis=dev.null.loglogis(y,wt,offset,object$nu))
    w <- object$w
    if (is.null(offset)) offset <- rep(0,length(object$eta))
    ## Residuals
    res <- residuals(object)*sqrt(w)
    dev.resid <- residuals(object,"deviance")
    ## Fitted values
    fitted <- fitted(object)
    ## dispersion
    sigma2 <- object$varht
    ## RSS, deviance
    rss <- sum(res^2)
    dev <- sum(dev.resid^2)
    ## Penalty associated with the fit
    obj.wk <- object
    obj.wk$d[] <- 0
    penalty <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,]))
    penalty <- as.vector(10^object$nlambda*penalty)
    if (!is.null(object$random)) {
        p.ran <- t(object$b)%*%object$random$sigma$fun(object$zeta,object$random$sigma$env)%*%object$b
        penalty <- penalty + p.ran
    }
    ## Calculate the diagnostics
    if (diagnostics) {
        ## Obtain retrospective linear model
        comp <- NULL
        p.dec <- NULL
        for (label in object$terms$labels) {
            if (label=="1") next
            if (label=="offset") next
            comp <- cbind(comp,predict(object,object$mf,inc=label))
            jk <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,],inc=label))
            p.dec <- c(p.dec,10^object$nlambda*jk)
        }
        term.label <- object$terms$labels[object$terms$labels!="1"]
        term.label <- term.label[term.label!="offset"]
        if (!is.null(object$random)) {
            mf <- object$mf
            mf$random <- object$random$z
            comp <- cbind(comp,predict(object,mf,inc=NULL))
            p.dec <- c(p.dec,p.ran)
            term.label <- c(term.label,"random")
        }
        fitted.off <- fitted-offset
        comp <- cbind(comp,yhat=fitted.off,y=fitted.off+res/sqrt(w),e=res/sqrt(w))
        if (any(outer(term.label,c("yhat","y","e"),"==")))
            warning("gss warning in summary.gssanova1: avoid using yhat, y, or e as variable names")
        colnames(comp) <- c(term.label,"yhat","y","e")
        ## Sweep out constant
        comp <- sqrt(w)*comp - outer(sqrt(w),apply(w*comp,2,sum))/sum(w)
        ## Obtain pi
        comp1 <- comp[,c(term.label,"yhat")]
        decom <- t(comp1) %*% comp1[,"yhat"]
        names(decom) <- c(term.label,"yhat")
        decom <- decom[term.label]/decom["yhat"]
        ## Obtain kappa, norm, and cosines        
        corr <- t(comp)%*%comp
        corr <- t(corr/sqrt(diag(corr)))/sqrt(diag(corr))
        norm <- apply(comp,2,function(x){sqrt(sum(x^2))})
        cosines <- rbind(corr[c("y","e"),],norm)
        rownames(cosines) <- c("cos.y","cos.e","norm")
        corr <- corr[term.label,term.label,drop=FALSE]
        if (qr(corr)$rank<dim(corr)[2])
            kappa <- rep(Inf,len=dim(corr)[2])
        else kappa <- as.numeric(sqrt(diag(solve(corr))))
        ## Obtain decomposition of penalty
        rough <- p.dec / penalty
        names(kappa) <- names(rough) <- term.label
    }
    else decom <- kappa <- cosines <- rough <- NULL
    ## Return the summaries
    z <- list(call=object$call,family=object$family,alpha=object$alpha,
              fitted=fitted,dispersion=sigma2,residuals=res/sqrt(w),rss=rss,
              deviance=dev,dev.resid=dev.resid,nu=object$nu,
              dev.null=dev.null,penalty=penalty,
              pi=decom,kappa=kappa,cosines=cosines,roughness=rough)
    class(z) <- "summary.gssanova1"
    z
}
## Summarize ssanova objects
summary.ssanova <- function(object,diagnostics=FALSE,...)
{
    y <- model.response(object$mf,"numeric")
    w <- model.weights(object$mf)
    offset <- model.offset(object$mf)
    if (is.null(offset)) offset <- rep(0,length(object$c))
    ## Residuals
    res <- 10^object$nlambda*object$c         
    if (!is.null(w)) res <- res/sqrt(w)
    ## Fitted values
    fitted <- as.numeric(y-res)
    fitted.off <- fitted-offset
    ## (estimated) sigma
    sigma <- sqrt(object$varht)
    ## R^2
    if (!is.null(w)) {
        r.squared <- sum(w*(fitted-sum(w*fitted)/sum(w))^2)
        r.squared <- r.squared/sum(w*(y-sum(w*y)/sum(w))^2)
    }
    else r.squared <- var(fitted)/var(y)       
    ## Residual sum of squares
    if (is.null(w)) rss <- sum(res^2)
    else rss <- sum(w*res^2)
    ## Penalty associated with the fit
    if (is.null(w)) 
        penalty <- sum(object$c*fitted.off)
    else penalty <- sum(object$c*fitted.off*sqrt(w))
    penalty <- as.vector(10^object$nlambda*penalty)
    ## Calculate the diagnostics
    if (diagnostics) {
        ## Obtain retrospective linear model
        comp <- NULL
        for (label in object$terms$labels) {
            if (label=="1") next
            if (label=="offset") next
            comp <- cbind(comp,predict(object,object$mf,inc=label))
        }
        comp <- cbind(comp,yhat=fitted.off,y=fitted.off+res,e=res)
        term.label <- object$terms$labels[object$terms$labels!="1"]
        term.label <- term.label[term.label!="offset"]
        if (any(outer(term.label,c("yhat","y","e"),"==")))
            warning("gss warning in summary.ssanova: avoid using yhat, y, or e as variable names")
        colnames(comp) <- c(term.label,"yhat","y","e")
        ## Sweep out constant
        if (!is.null(w))
            comp <- sqrt(w)*comp - outer(sqrt(w),apply(w*comp,2,sum))/sum(w)
        else comp <- sweep(comp,2,apply(comp,2,mean))
        ## Obtain pi
        comp1 <- comp[,c(term.label,"yhat")]
        decom <- t(comp1) %*% comp1[,"yhat"]
        names(decom) <- c(term.label,"yhat")
        decom <- decom[term.label]/decom["yhat"]
        ## Obtain kappa, norm, and cosines
        corr <- t(comp)%*%comp
        corr <- t(corr/sqrt(diag(corr)))/sqrt(diag(corr))
        norm <- apply(comp,2,function(x){sqrt(sum(x^2))})
        cosines <- rbind(corr[c("y","e"),],norm)
        rownames(cosines) <- c("cos.y","cos.e","norm")
        corr <- corr[term.label,term.label,drop=FALSE]
        if (qr(corr)$rank<dim(corr)[2])
            kappa <- rep(Inf,len=dim(corr)[2])
        else kappa <- as.numeric(sqrt(diag(solve(corr))))
        ## Obtain decomposition of penalty
        rough <- as.vector(10^object$nlambda*t(comp[,term.label])%*%object$c/penalty)
        names(kappa) <- names(rough) <- term.label
    }
    else decom <- kappa <- cosines <- rough <- NULL
    ## Return the summaries
    z <- list(call=object$call,method=object$method,fitted=fitted,residuals=res,
              sigma=sigma,r.squared=r.squared,rss=rss,penalty=penalty,
              pi=decom,kappa=kappa,cosines=cosines,roughness=rough)
    class(z) <- "summary.ssanova"
    z
}
## Summarize ssanova objects
summary.ssanova1 <- function(object,diagnostics=FALSE,...)
{
    y <- model.response(object$mf,"numeric")
    w <- model.weights(object$mf)
    offset <- model.offset(object$mf)
    if (is.null(offset)) offset <- rep(0,length(y))
    ## Residuals
    mf <- object$mf
    if (!is.null(object$random)) mf$random <- object$random$z
    res <- y - predict(object,mf)
    ## Fitted values
    fitted <- as.numeric(y-res)
    ## (estimated) sigma
    sigma <- sqrt(object$varht)
    ## R^2
    if (!is.null(w)) {
        r.squared <- sum(w*(fitted-sum(w*fitted)/sum(w))^2)
        r.squared <- r.squared/sum(w*(y-sum(w*y)/sum(w))^2)
    }
    else r.squared <- var(fitted)/var(y)       
    ## Residual sum of squares
    if (is.null(w)) rss <- sum(res^2)
    else rss <- sum(w*res^2)
    ## Penalty associated with the fit
    obj.wk <- object
    obj.wk$d[] <- 0
    penalty <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,]))
    penalty <- as.vector(10^object$nlambda*penalty)
    if (!is.null(object$random)) {
        p.ran <- t(object$b)%*%object$random$sigma$fun(object$zeta,object$random$sigma$env)%*%object$b
        penalty <- penalty + p.ran
    }
    ## Calculate the diagnostics
    if (diagnostics) {
        ## Obtain retrospective linear model
        comp <- NULL
        p.dec <- NULL
        for (label in object$terms$labels) {
            if (label=="1") next
            if (label=="offset") next
            comp <- cbind(comp,predict(object,object$mf,inc=label))
            jk <- sum(obj.wk$c*predict(obj.wk,obj.wk$mf[object$id.basis,],inc=label))
            p.dec <- c(p.dec,10^object$nlambda*jk)
        }
        term.label <- object$terms$labels[object$terms$labels!="1"]
        term.label <- term.label[term.label!="offset"]
        if (!is.null(object$random)) {
            comp <- cbind(comp,predict(object,mf,inc=NULL))
            p.dec <- c(p.dec,p.ran)
            term.label <- c(term.label,"random")
        }
        fitted.off <- fitted-offset
        comp <- cbind(comp,yhat=fitted.off,y=fitted.off+res,e=res)
        if (any(outer(term.label,c("yhat","y","e"),"==")))
            warning("gss warning in summary.ssanova1: avoid using yhat, y, or e as variable names")
        colnames(comp) <- c(term.label,"yhat","y","e")
        ## Sweep out constant
        if (!is.null(w))
            comp <- sqrt(w)*comp - outer(sqrt(w),apply(w*comp,2,sum))/sum(w)
        else comp <- sweep(comp,2,apply(comp,2,mean))
        ## Obtain pi
        comp1 <- comp[,c(term.label,"yhat")]
        decom <- t(comp1) %*% comp1[,"yhat"]
        names(decom) <- c(term.label,"yhat")
        decom <- decom[term.label]/decom["yhat"]
        ## Obtain kappa, norm, and cosines
        corr <- t(comp)%*%comp
        corr <- t(corr/sqrt(diag(corr)))/sqrt(diag(corr))
        norm <- apply(comp,2,function(x){sqrt(sum(x^2))})
        cosines <- rbind(corr[c("y","e"),],norm)
        rownames(cosines) <- c("cos.y","cos.e","norm")
        corr <- corr[term.label,term.label,drop=FALSE]
        if (qr(corr)$rank<dim(corr)[2])
            kappa <- rep(Inf,len=dim(corr)[2])
        else kappa <- as.numeric(sqrt(diag(solve(corr))))
        ## Obtain decomposition of penalty
        rough <- p.dec / penalty
        names(kappa) <- names(rough) <- term.label
    }
    else decom <- kappa <- cosines <- rough <- NULL
    ## Return the summaries
    z <- list(call=object$call,method=object$method,fitted=fitted,residuals=res,
              sigma=sigma,r.squared=r.squared,rss=rss,penalty=penalty,
              pi=decom,kappa=kappa,cosines=cosines,roughness=rough)
    class(z) <- "summary.ssanova"
    z
}
.First.lib <- function(lib, pkg)
{
    library.dynam("gss", pkg, lib)
}
project <- function (object,...)
{
    UseMethod("project")
}
