# file nnet/multinom.q copyright (C) 1994-9 W. N. Venables and B. D. Ripley
#

multinom <-
    function(formula, data=sys.frame(sys.parent()), weights, subset, na.action,
             contrasts=NULL, Hess=FALSE, summ=0, censored=FALSE, ...)
{
    class.ind <- function(cl)
    {
        n <- length(cl)
        x <- matrix(0, n, length(levels(cl)))
        x[(1:n) + n * (as.vector(unclass(cl)) - 1)] <- 1
        dimnames(x) <- list(names(cl), levels(cl))
        x
    }
    summ2 <- function(X, Y)
    {
        X <- as.matrix(X)
        Y <- as.matrix(Y)
        n <- dim(X)[1]
        p <- dim(X)[2]
        q <- dim(Y)[2]
        Z <- t(cbind(X, Y))
        storage.mode(Z) <- "double"
        z <- .C("VR_summ2",
                as.integer(n),
                as.integer(p),
                as.integer(q),
                Z = Z,
                na = integer(1))
        Za <- t(z$Z[, 1:z$na, drop = FALSE])
        list(X = Za[, 1:p, drop = FALSE], Y = Za[, p + 1:q])
    }

    call <- match.call()
    m <- match.call(expand = FALSE)
    m$summ <- m$Hess <- m$contrasts <- m$censored <- m$... <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    X <- model.matrix(Terms, m, contrasts)
    Xr <- qr(X)$rank
    Y <- model.extract(m, response)
    if(!is.matrix(Y)) Y <- as.factor(Y)
    w <- model.extract(m, weights)
    if(length(w) == 0)
        if(is.matrix(Y)) w <- rep(1, dim(Y)[1])
        else w <- rep(1, length(Y))
    lev <- levels(Y)
    if(is.factor(Y)) {
        counts <- table(Y)
        if(any(counts == 0)) {
            warning(paste("group(s)", paste(lev[counts == 0], collapse=" "),
                          "are empty"))
            Y <- factor(Y, levels=lev[counts > 0])
            lev <- lev[counts > 0]
        }
        if(length(lev) == 2) Y <- as.vector(unclass(Y)) - 1
        else Y <- class.ind(Y)
    }
    if(summ==1) {
        Z <- cbind(X, Y)
        assign("z1", cumprod(apply(Z, 2, max)+1))
        Z1 <- apply(Z, 1, function(x) sum(z1*x))
        oZ <- order(Z1)
        Z2 <- !duplicated(Z1[oZ])
        oX <- (seq(along=Z1)[oZ])[Z2]
        X <- X[oX, , drop=FALSE]
        Y <- if(is.matrix(Y)) Y[oX, , drop=FALSE] else Y[oX]
        w <- diff(c(0,cumsum(w))[c(Z2,TRUE)])
        print(dim(X))
    }
    if(summ==2) {
        Z <- summ2(cbind(X, Y), w)
        X <- Z$X[, 1:ncol(X)]
        Y <- Z$X[, ncol(X) + 1:ncol(Y), drop = FALSE]
        w <- Z$Y
        print(dim(X))
    }
    if(summ==3) {
        Z <- summ2(X, Y*w)
        X <- Z$X
        Y <- Z$Y[, 1:ncol(Y), drop = FALSE]
        w <- rep(1, nrow(X))
        print(dim(X))
    }
    offset <- model.extract(m, offset)
    r <- ncol(X)
    if(is.matrix(Y)) {
        # 3 or more response levels or direct matrix spec.
        p <- ncol(Y)
        sY <- Y %*% rep(1, p)
        if(any(sY==0)) stop("some case has no observations")
        if(!censored) {
            Y <- Y / matrix(sY, nrow(Y), p)
            w <- w*sY
        }
        if(length(offset) > 1) {
            if(ncol(offset) !=  p) stop("ncol(offset) is wrong")
            mask <- c(rep(0, r+1+p), rep(c(0, rep(1, r), rep(0, p)), p-1) )
            X <- cbind(X, offset)
            Wts <- as.vector(rbind(matrix(0, r+1, p), diag(p)))
            fit <- nnet.default(X, Y, w, Wts=Wts, mask=mask, size=0, skip=TRUE,
                                softmax=TRUE, censored=censored, rang=0, ...)
        } else {
            mask <- c(rep(0, r+1), rep(c(0, rep(1, r)), p-1) )
            fit <- nnet.default(X, Y, w, mask=mask, size=0, skip=TRUE,
                                softmax=TRUE, censored=censored, rang=0, ...)
        }
    } else {
        # 2 response levels
        if(length(offset) <= 1) {
            mask <- c(0, rep(1, r))
            fit <- nnet.default(X, Y, w, mask=mask, size=0, skip=TRUE,
                                entropy=TRUE, rang=0, ...)
        } else {
            mask <- c(0, rep(1, r), 0)
            Wts <- c(rep(0, r+1), 1)
            X <- cbind(X, offset)
            fit <- nnet.default(X, Y, w, Wts=Wts, mask=mask, size=0, skip=TRUE,
                                entropy=TRUE, rang=0, ...)
        }
    }
    fit$formula <- as.vector(attr(Terms, "formula"))
    fit$terms <- Terms
    fit$call <- call
    fit$weights <- w
    fit$lev <- lev
    fit$deviance <- 2 * fit$value
    fit$rank <- Xr
    edf <- ifelse(length(lev) == 2, 1, length(lev)-1)*Xr
    if(is.matrix(Y)) {
        edf <- (ncol(Y)-1)*Xr
        if(length(dn <- colnames(Y)) > 0) fit$lab <- dn
        else fit$lab <- 1:ncol(Y)
    }
    fit$coefnames <- colnames(X)
    fit$vcoefnames <- fit$coefnames[1:r]# remove offset cols
    fit$edf <- edf
    fit$AIC <- fit$deviance + 2 * edf
    class(fit) <- c("multinom", "nnet")
    if(Hess) {
        mask <- as.logical(mask)
        fit$Hessian <- nnet.Hess(fit, X, Y, w)[mask, mask]
        cf <- fit$vcoefnames
        if(length(fit$lev) != 2) {
            bf <- if(length(fit$lev)) fit$lev else fit$lab
            cf <- t(outer(bf[-1], cf, function(x,y) paste(x,y,sep=":")))
        }
        dimnames(fit$Hessian) <- list(cf, cf)
    }
    fit
}

predict.multinom <- function(object, newdata, type=c("class","probs"), ...)
{
    if(!inherits(object, "multinom")) stop("Not a multinom fit")
    type <- match.arg(type)
    if(missing(newdata)) Y <- object$fitted
    else {
        form <- delete.response(object$terms)
        X <- model.matrix(form, newdata)
        Y <- predict.nnet(object, X)
    }
    switch(type, class={
        if(length(object$lev) > 2)
            Y <- factor(max.col(Y), levels=seq(along=object$lev),
                        labels=object$lev)
        if(length(object$lev) == 2)
            Y <- factor(1 + (Y > 0.5), levels=1:2, labels=object$lev)
        if(length(object$lev) == 0)
            Y <- factor(max.col(Y), levels=seq(along=object$lab),
                        labels=object$lab)
    }, probs={})
    drop(Y)
}

print.multinom <- function(x, ...)
{
    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl)
    }
    cat("\nCoefficients:\n")
    print(coef(x), ...)
    cat("\nResidual Deviance:", format(x$deviance), "\n")
    cat("AIC:", format(x$AIC), "\n")
    invisible(x)
}

coef.multinom <- function(object, ...)
{
  r <- length(object$vcoefnames)
  if(length(object$lev) == 2) {
    coef <- object$wts[1+(1:r)]
    names(coef) <- object$vcoefnames
  } else {
    coef <- matrix(object$wts, nrow = object$n[3], byrow=TRUE)[, 1+(1:r), drop=FALSE]
    if(length(object$lev))
        dimnames(coef) <- list(object$lev, object$vcoefnames)
    if(length(object$lab))
        dimnames(coef) <- list(object$lab, object$vcoefnames)
    coef <- coef[-1, , drop=FALSE]
  }
  coef
}

drop1.multinom <- function(object, scope, sorted = TRUE, trace = FALSE)
{
    if(!inherits(object, "multinom")) stop("Not a multinom fit")
    if(missing(scope)) scope <- drop.scope(object)
    else {
        if(!is.character(scope))
            scope <- attr(terms(update.formula(object, scope)), "term.labels")
        if(!all(match(scope, attr(object$terms, "term.labels"), FALSE)))
            stop("scope is not a subset of term labels")
    }
    ns <- length(scope)
    ans <- matrix(nrow = ns+1, ncol = 2,
                  dimnames = list(c("<none>", paste("-",scope,sep="")),
                  c("df", "AIC")))
    ans[1, ] <- c(object$edf, object$AIC)
    i <- 2
    for(tt in scope) {
        cat("trying -", tt,"\n")
        nobject <- update(object, paste("~ . -", tt), trace = trace)
        if(nobject$edf == object$edf) nobject$AIC <- NA
        ans[i, ] <- c(nobject$edf, nobject$AIC)
        i <- i+1
    }
    if(sorted) ans <- ans[order(ans[, 2]), ]
    ans
}

add1.multinom <- function(object, scope, sorted = TRUE, trace = FALSE)
{
    if(!inherits(object, "multinom")) stop("Not a multinom fit")
    if(!is.character(scope))
        scope <- add.scope(object, update.formula(object, scope,
                                                  evaluate = FALSE))
    if(!length(scope))
        stop("no terms in scope for adding to object")
    ns <- length(scope)
    ans <- matrix(nrow = ns+1, ncol = 2,
                  dimnames = list(c("<none>",paste("+",scope,sep="")),
                  c("df", "AIC")))
    ans[1, ] <- c(object$edf, object$AIC)
    i <- 2
    for(tt in scope) {
        cat("trying +", tt,"\n")
        nobject <- update(object, paste("~ . +", tt), trace=trace)
        if(nobject$edf == object$edf) nobject$AIC <- NA
        ans[i, ] <- c(nobject$edf, nobject$AIC)
        i <- i+1
    }
    if(sorted) ans <- ans[order(ans[, 2]), ]
    ans
}

extractAIC.multinom <- function(fit, scale, k = 2, ...)
    c(fit$edf, fit$AIC + (k-2)*fit$edf)

vcov.multinom <- function(object)
{
    ginv <- function(X, tol = sqrt(.Machine$double.eps))
    {
        #
        # based on suggestions of R M Heiberger, TRUE M Hesterburg and WNV
        #
        Xsvd <- svd(X)
        Positive <- Xsvd$d > max(tol * Xsvd$d[1], 0)
        if(!any(Positive)) array(0, dim(X)[2:1])
        else Xsvd$v[, Positive] %*% ((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive]))
    }

    if(is.null(object$Hessian)) {
        cat("\nRe-fitting to get Hessian\n\n")
        object <- update(object, Hess=TRUE, trace=FALSE)
    }
    structure(ginv(object$Hessian), dimnames = dimnames(object$Hessian))
}

summary.multinom <-
function(object, correlation = TRUE, digits = options()$digits,
         Wald.ratios = FALSE, ...)
{
    vc <- vcov(object)
    r <- length(object$vcoefnames)
    se <- sqrt(diag(vc))
    if(length(object$lev) == 2) {
        coef <- object$wts[1 + (1:r)]
        stderr <- se
        names(coef) <- names(stderr) <- object$vcoefnames
    } else {
        coef <- matrix(object$wts, nrow = object$n[3],
                       byrow = TRUE)[-1, 1 + (1:r), drop = FALSE]
        stderr <- matrix(se, nrow = object$n[3] - 1, byrow = TRUE)
        if(length(l <- object$lab) || length(l <- object$lev))
            dimnames(coef) <- dimnames(stderr) <- list(l[-1], object$vcoefnames)
    }
    object$is.binomial <- (length(object$lev) == 2)
    object$digits <- digits
    object$coefficients <- coef
    object$standard.errors <- stderr
    if(Wald.ratios) object$Wald.ratios <- coef/stderr
    if(correlation) object$correlation <- vc/outer(se, se)
    class(object) <- "summary.multinom"
    object
}

print.summary.multinom <- function(x, digits = x$digits, ...)
{
    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl)
    }
    cat("\nCoefficients:\n")
    if(x$is.binomial) {
        print(cbind(Values = x$coefficients,
                    "Std. Err." = x$standard.errors,
                    "Value/SE" = x$Wald.ratios),
              digits = digits)
    } else {
        print(x$coefficients, digits = digits)
        cat("\nStd. Errors:\n")
        print(x$standard.errors, digits = digits)
        if(!is.null(x$Wald.ratios)) {
            cat("\nValue/SE (Wald statistics):\n")
            print(x$coefficients/x$standard.errors, digits = digits)
        }
    }
    cat("\nResidual Deviance:", format(x$deviance), "\n")
    cat("AIC:", format(x$AIC), "\n")
    if(!is.null(correl <- x$correlation)) {
        p <- dim(correl)[2]
        if(p > 1) {
            cat("\nCorrelation of Coefficients:\n")
            ll <- lower.tri(correl)
            correl[ll] <- format(round(correl[ll], digits))
            correl[!ll] <- ""
            print(correl[-1, -p], quote = FALSE, ...)
        }
    }
    invisible(x)
}

anova.multinom <- function(object, ..., test = c("Chisq", "none"))
{
    test <- match.arg(test)
    dots <- list(...)
    if(length(dots) == 0)
        stop("anova is not implemented for a single multinom object")
    mlist <- list(object, ...)
    nt <- length(mlist)
    dflis <- sapply(mlist, function(x) x$edf)
    s <- order(dflis)
    dflis <- nrow(residuals(object)) * (ncol(residuals(object))-1) - dflis
    mlist <- mlist[s]
    if(any(!sapply(mlist, inherits, "multinom")))
        stop("not all objects are of class `multinom'")
    rsp <- unique(sapply(mlist, function(x) paste(formula(x)[2])))
    mds <- sapply(mlist, function(x) paste(formula(x)[3]))
    dfs <- dflis[s]
    lls <- sapply(mlist, function(x) deviance(x))
    tss <- c("", paste(1:(nt - 1), 2:nt, sep = " vs "))
    df <- c(NA, -diff(dfs))
    x2 <- c(NA, -diff(lls))
    pr <- c(NA, 1 - pchisq(x2[-1], df[-1]))
    out <- data.frame(Model = mds, Resid.df = dfs,
                      Deviance = lls, Test = tss, Df = df, LRtest = x2,
                      Prob = pr)
    names(out) <- c("Model", "Resid. df", "Resid. Dev", "Test",
                    "   Df", "LR stat.", "Pr(Chi)")
    if(test=="none") out <- out[, 1:6]
    class(out) <- c("Anova", "data.frame")
    attr(out, "heading") <-
        c("Likelihood ratio tests of Multinomial Models\n",
          paste("Response:", rsp))
    out
}

# file nnet/nnet.q copyright (C) 1994-9 W. N. Venables and B. D. Ripley
#
nnet <- function(object, ...)
{
    if(is.null(class(object))) class(object) <- data.class(object)
    UseMethod("nnet")
}

nnet.formula <- function(formula, data = NULL, weights, ...,
			subset, na.action = na.fail, contrasts=NULL)
{
    class.ind <- function(cl)
    {
        n <- length(cl)
        x <- matrix(0, n, length(levels(cl)))
        x[(1:n) + n * (as.vector(unclass(cl)) - 1)] <- 1
        dimnames(x) <- list(names(cl), levels(cl))
        x
    }
    m <- match.call(expand.dots = FALSE)
    if(is.matrix(eval(m$data, sys.frame(sys.parent()))))
        m$data <- as.data.frame(data)
    m$... <- m$contrasts <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval(m, sys.frame(sys.parent()))
    Terms <- attr(m, "terms")
    x <- model.matrix(Terms, m, contrasts)
    xint <- match("(Intercept)", colnames(x), nomatch=0)
    if(xint > 0) x <- x[, -xint, drop=FALSE]# Bias term is used for intercepts
    w <- model.extract(m, weights)
    if(length(w) == 0) w <- rep(1, nrow(x))
    y <- model.extract(m, response)
    if(is.factor(y)) {
        lev <- levels(y)
        counts <- table(y)
        if(any(counts == 0)) {
            warning(paste("group(s)", paste(lev[counts == 0], collapse=" "),
                          "are empty"))
            y <- factor(y, levels=lev[counts > 0])
        }
        if(length(lev) == 2) {
            y <- as.vector(unclass(y)) - 1
            res <- nnet.default(x, y, w, entropy=TRUE, ...)
            res$lev <- lev
        } else {
            y <- class.ind(y)
            res <- nnet.default(x, y, w, softmax=TRUE, ...)
            res$lev <- lev
        }
    } else res <- nnet.default(x, y, w, ...)
    res$terms <- Terms
    res$coefnames <- colnames(x)
    res$call <- match.call()
    class(res) <- c("nnet.formula", "nnet")
    res
}

nnet.default <-
function(x, y, weights, size, Wts, mask=rep(TRUE, length(wts)),
	 linout=FALSE, entropy=FALSE, softmax=FALSE, censored=FALSE,
         skip=FALSE, rang=0.7, decay=0, maxit=100, Hess=FALSE, trace=TRUE,
         MaxNWts=1000, abstol=1.0e-4, reltol=1.0e-8, ...)
{
    net <- NULL
    x <- as.matrix(x)
    y <- as.matrix(y)
    if(any(is.na(x))) stop("missing values in x")
    if(any(is.na(y))) stop("missing values in y")
    if(dim(x)[1] != dim(y)[1]) stop("nrows of x and y must match")
    if(linout && entropy) stop("entropy fit only for logistic units")
    if(softmax) {linout <- TRUE; entropy <- FALSE}
    if(censored) {linout <- TRUE; entropy <- FALSE; softmax <- TRUE}
    net$n <- c(dim(x)[2], size, dim(y)[2])
    net$nunits <- 1 + sum(net$n)
    net$nconn <- rep(0, net$nunits+1)
    net$conn <- numeric(0)
    net <- norm.net(net)
    if(skip) net <- add.net(net, seq(1,net$n[1]),
                            seq(1+net$n[1]+net$n[2], net$nunits-1))
    if((nwts <- length(net$conn))==0) stop("No weights to fit")
    if(nwts > MaxNWts)
        stop(paste("Too many (", nwts, ") weights", sep=""))
    nsunits <- net$nunits
    if(linout) nsunits <- net$nunits - net$n[3]
    net$nsunits <- nsunits
    net$decay <- decay
    net$entropy <- entropy
    net$softmax <- softmax
    net$censored <- censored
    if(missing(Wts))
        if(rang > 0) wts <- runif(nwts, -rang, rang)
        else wts <- rep(0, nwts)
    else wts <- Wts
    if(length(wts) != nwts) stop("weights vector of incorrect length")
    if(length(mask) != length(wts)) stop("incorrect length of mask")
    if(trace) {
        cat("                           # weights: ", length(wts))
        nw <- sum(mask != 0)
        if(nw < length(wts)) cat(" (", nw, " variable)\n",sep="")
        else cat("\n")
    }
    if(length(decay) == 1) decay <- rep(decay, length(wts))
    if(!is.loaded(symbol.C("VR_set_net")))
        stop("Compiled code has not been dynamically loaded")
    .C("VR_set_net",
       as.integer(net$n),
       as.integer(net$nconn),
       as.integer(net$conn),
       as.double(decay),
       as.integer(nsunits),
       as.integer(entropy),
       as.integer(softmax),
       as.integer(censored)
       )
    ntr <- dim(x)[1]
    nout <- dim(y)[2]
    if(missing(weights)) weights <- rep(1, ntr)
    if(length(weights) != ntr || any(weights < 0))
        stop("invalid weights vector")
    Z <- as.double(cbind(x,y))
    storage.mode(weights) <- "double"
    z <- .C("VR_set_train", as.integer(ntr), Z, weights)
    tmp <- .C("VR_dovm",
              as.integer(length(wts)),
              wts=as.double(wts),
              val=double(1),
              as.integer(maxit),
              as.logical(trace),
              as.integer(mask),
              as.double(abstol), as.double(reltol)
              )
    .C("VR_unset_train")
    net$value <- tmp$val
    net$wts <- tmp$wts
    tmp <- matrix(.C("VR_nntest",
                     as.integer(ntr), Z, tclass = double(ntr*nout),
                     as.double(net$wts))$tclass,  ntr, nout)
    dimnames(tmp) <- list(rownames(x), colnames(y))
    net$fitted.values <- tmp
    tmp <- y - tmp
    dimnames(tmp) <- list(rownames(x), colnames(y))
    net$residuals <- tmp
    .C("VR_unset_net")
    if(entropy) net$lev <- c("0","1")
    if(softmax) net$lev <- colnames(y)
    net$call <- match.call()
    if(Hess) net$Hessian <- nnet.Hess(net, x, y, weights)
    class(net) <- "nnet"
    net
}

predict.nnet.formula <- function(object, newdata, type=c("raw", "class"), ...)
{
    if(!inherits(object, "nnet.formula"))
        stop("object not of class nnet.formula")
    type <- match.arg(type)
    if(missing(newdata))
        switch(type,
               raw = return(object$fitted.values),
               class = {
                   if(is.null(object$lev)) stop("inappropriate fit for class")
                   z <- object$fitted.values
                   if(ncol(z) > 1) object$lev[max.col(z)]
                   else object$lev[1 + (z > 0.5)]
               })
    x <- model.matrix(delete.response(object$terms), newdata)
    xint <- match("(Intercept)", colnames(x), nomatch=0)
    if(xint > 0) x <- x[, -xint, drop=FALSE]# Bias term is used for intercepts
    predict.nnet(object=object, x=x, type=type)
}

predict.nnet <- function(object, x, type=c("raw","class"), ...)
{
    if(!inherits(object, "nnet")) stop("object not of class nnet")
    type <- match.arg(type)
    if(missing(x)) z <- object$fitted
    else {
        if(is.null(dim(x))) dim(x) <- c(1, length(x))
        if(any(is.na(x))) stop("missing values in x")
        x <- as.matrix(x)
        ntr <- dim(x)[1]
        nout <- object$n[3]
        if(!is.loaded(symbol.C("VR_set_net")))
            stop("Compiled code has not been dynamically loaded")
        .C("VR_set_net",
           as.integer(object$n),
           as.integer(object$nconn),
           as.integer(object$conn),
           as.double(object$decay),
           as.integer(object$nsunits),
           as.integer(0),
           as.integer(object$softmax),
           as.integer(object$censored)
           )
        z <- .C("VR_nntest",
                as.integer(ntr),
                as.double(x),
                tclass = double(ntr*nout),
                as.double(object$wts)
                )
        z <- matrix(z$tclass, ntr, nout)
        dimnames(z) <- list(rownames(x), dimnames(object$fitted)[[2]])
        .C("VR_unset_net")
    }
    switch(type, raw = z,
           class = {
               if(is.null(object$lev)) stop("inappropriate fit for class")
               if(ncol(z) > 1) object$lev[max.col(z)]
               else object$lev[1 + (z > 0.5)]
           })
}


eval.nn <- function(wts)
{
    z <- .C("VR_dfunc",
            as.double(wts), df = double(length(wts)), fp = as.double(1))
    fp <- z$fp
    attr(fp, "gradient") <- z$df
    fp
}

add.net <- function(net, from, to)
{
    nconn <- net$nconn
    conn <- net$conn
    for(i in to){
        ns <- nconn[i+2]
        cadd <- from
        if(nconn[i+1] == ns) cadd <- c(0,from)
        con <- NULL
        if(ns > 1) con <- conn[1:ns]
        con <- c(con, cadd)
        if(length(conn) > ns) con <- c(con, conn[(ns+1):length(conn)])
        for(j in (i+1):net$nunits) nconn[j+1] <- nconn[j+1]+length(cadd)
        conn <- con
    }
    net$nconn <- nconn
    net$conn <- con
    net
}

norm.net <- function(net)
{
    n <- net$n; n0 <- n[1]; n1 <- n0+n[2]; n2 <- n1+n[3];
    if(n[2] <= 0) return(net)
    net <- add.net(net, 1:n0,(n0+1):n1)
    add.net(net, (n0+1):n1, (n1+1):n2)
}

which.is.max <- function(x)
{
    y <- seq(length(x))[x == max(x)]
    if(length(y) > 1) sample(y,1) else y
}

nnet.Hess <- function(net, x, y, weights)
{
    x <- as.matrix(x)
    y <- as.matrix(y)
    if(dim(x)[1] != dim(y)[1]) stop("dims of x and y must match")
    nw <- length(net$wts)
    decay <- net$decay
    if(length(decay) == 1) decay <- rep(decay, nw)
    if(!is.loaded(symbol.C("VR_set_net")))
        stop("Compiled code has not been dynamically loaded")
    .C("VR_set_net",
       as.integer(net$n),
       as.integer(net$nconn),
       as.integer(net$conn),
       as.double(decay),
       as.integer(net$nsunits),
       as.integer(net$entropy),
       as.integer(net$softmax),
       as.integer(net$censored)
       )
    ntr <- dim(x)[1]
    nout <- dim(y)[2]
    if(missing(weights)) weights <- rep(1, ntr)
    if(length(weights) != ntr || any(weights < 0))
        stop("invalid weights vector")
    Z <- as.double(cbind(x,y))
    storage.mode(weights) <- "double"
    z <- .C("VR_set_train",
            as.integer(ntr),
            Z,
            weights
            )
    z <- matrix(.C("VR_nnHessian", as.double(net$wts), H = double(nw*nw))$H,
                nw, nw)
    .C("VR_unset_train"); .C("VR_unset_net")
    z
}

class.ind <- function(cl)
{
    n <- length(cl)
    cl <- as.factor(cl)
    x <- matrix(0, n, length(levels(cl)) )
    x[(1:n) + n*(unclass(cl)-1)] <- 1
    dimnames(x) <- list(names(cl), levels(cl))
    x
}

print.nnet <- function(x, ...)
{
    if(!inherits(x, "nnet")) stop("Not legitimate a neural net fit")
    cat("a ",x$n[1],"-",x$n[2],"-",x$n[3]," network", sep="")
    cat(" with", length(x$wts),"weights\n")
    if(length(x$coefnames))  cat("inputs:", x$coefnames, "\noutput(s):",
                                 deparse(formula(x)[[2]]), "\n")
    cat("options were -")
    tconn <- diff(x$nconn)
    if(tconn[length(tconn)] > x$n[2]+1) cat(" skip-layer connections ")
    if(x$nunits > x$nsunits && !x$softmax) cat(" linear output units ")
    if(x$entropy) cat(" entropy fitting ")
    if(x$softmax) cat(" softmax modelling ")
    if(x$decay[1] > 0) cat(" decay=", x$decay[1], sep="")
    cat("\n")
    invisible(x)
}

coef.nnet <- function(object, ...)
{
    wts <- object$wts
    wm <- c("b", paste("i", seq(length=object$n[1]), sep=""))
    if(object$n[2] > 0)
        wm <- c(wm, paste("h", seq(length=object$n[2]), sep=""))
    if(object$n[3] > 1)  wm <- c(wm,
                                 paste("o", seq(length=object$n[3]), sep=""))
    else wm <- c(wm, "o")
    names(wts) <- apply(cbind(wm[1+object$conn],
                              wm[1+rep(1:object$nunits - 1, diff(object$nconn))]),
                        1, function(x)  paste(x, collapse = "->"))
    wts
}

summary.nnet <- function(object, ...)
{
    class(object) <- c("summary.nnet", class(object))
    object
}


print.summary.nnet <- function(x, ...)
{
    cat("a ",x$n[1],"-",x$n[2],"-",x$n[3]," network", sep="")
    cat(" with", length(x$wts),"weights\n")
    cat("options were -")
    tconn <- diff(x$nconn)
    if(tconn[length(tconn)] > x$n[2]+1) cat(" skip-layer connections ")
    if(x$nunits > x$nsunits && !x$softmax) cat(" linear output units ")
    if(x$entropy) cat(" entropy fitting ")
    if(x$softmax) cat(" softmax modelling ")
    if(x$decay[1] > 0) cat(" decay=", x$decay[1], sep="")
    cat("\n")
    wts <- format(round(coef.nnet(x),2))
    lapply(split(wts,rep(1:x$nunits, tconn)),
           function(x) print(x, quote=FALSE))
    invisible(x)
}

max.col <- function(m)
{
    if(!is.loaded(symbol.C("VR_max_col")))
        stop("Compiled code has not been dynamically loaded")
    m <- as.matrix(m)
    n <- nrow(m)
    .C("VR_max_col",
       as.double(m),
       as.integer(n),
       as.integer(ncol(m)),
       rmax = integer(n)
       )$rmax
}

residuals.nnet <- function(object, ...) object$residuals
.First.lib <- function(lib, pkg)
    library.dynam("nnet", pkg, lib)
