# file nnet/multinom.q copyright (C) 1994-8 W. N. Venables and B. D. Ripley
#
multinom <- function(formula, data=sys.parent(), weights, subset, na.action,
	    contrasts=NULL, Hess=F, summ=0, censored=F, ...)
{
  class.ind <- function(cl)
  {
    n <- length(cl)
    x <- matrix(0, n, length(levels(cl)))
    x[(1:n) + n * (as.vector(codes(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 = F])
    list(X = Za[, 1:p, drop = F], Y = Za[, p + 1:q])
  }

  call <- match.call()
  m <- match.call(expand = F)
  m$summ <- m$Hess <- m$contrasts <- m$censored <- m$... <- NULL
  m[[1]] <- as.name("model.frame")
  m <- eval(m, 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(codes(Y)) - 1
    else Y <- class.ind(Y)
  }
  if(summ==1) {
    Z <- cbind(X,Y)
    assign("z1", cumprod(apply(Z, 2, max)+1), frame=1)
    Z1 <- apply(Z, 1, function(x) sum(z1*x))
    oZ <- sort.list(Z1)
    Z2 <- !duplicated(Z1[oZ])
    oX <- ((1:length(Z1))[oZ])[Z2]
    X <- X[oX,]
    if(is.matrix(Y)) Y <- Y[oX,] else Y <- Y[oX]
    w <- diff(c(0,cumsum(w))[c(Z2,T)])
    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 = F]
    w <- Z$Y
    print(dim(X))
  }
  if(summ==3) {
    Z <- summ2(X, Y*w)
    X <- Z$X
    Y <- Z$Y[,1:ncol(Y), drop = F]
    w <- rep(1, nrow(X))
    print(dim(X))
  }
  offset <- model.extract(m, offset)
  r <- ncol(X)
  if(is.matrix(Y)) {
    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=T, softmax=T, 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=T, softmax=T, censored=censored, rang=0, ...)
    }
  } else {
    if((p <- length(lev)) == 2) {
      if(length(offset) == 1) {
	mask <- c(0, rep(1, r))
	fit <- nnet.default(X, Y, w, mask=mask, size=0, skip=T, entropy=T, 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=T, entropy=T, rang=0, ...)
      }
    } else {
      mask <- c(rep(0, ncol(X)+1), rep(1, (p-1)*(ncol(X)+1)))
      fit <- nnet.default(X, Y, w, mask=mask, size=0, skip=T, softmax=T, 
			  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 <- dimnames(Y)[[2]]) > 0) fit$lab <- dn
      else fit$lab <- 1:ncol(Y)
  }
  fit$coefnames <- dimnames(X)[[2]]
  fit$vcoefnames <- fit$coefnames[1:r] # remove offset cols
  fit$edf <- edf
  fit$AIC <- 2 * (fit$value + edf)
  attr(fit, "class") <- 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) {
     if(length(fit$lev)) bf <- fit$lev else bf <- 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.values
  else {
    form <- delete.response(terms(object$form))
    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)
  },)
  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=T)[, 1+(1:r), drop=F]
    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=F]
  }
  coef
}

drop1.multinom <- function(object, scope, sorted=T, trace=F)
{
  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"), F)))
	stop("scope is not a subset of term labels")
    }
  n <- length(scope)
  ans <- matrix(nrow=n+1, ncol=2)
  dimnames(ans) <- 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[sort.list(ans[,2]),]
  ans
}

add1.multinom <- function(object, scope, sorted=T, trace=F)
{
  if(!inherits(object, "multinom")) stop("Not a multinom fit")
  if(!is.character(scope))
    scope <- add.scope(object, update.formula(object, scope, 
					   evaluate = F))
  if(!length(scope))
    stop("no terms in scope for adding to object")	
  n <- length(scope)
  ans <- matrix(nrow=n+1, ncol=2)
  dimnames(ans) <- 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[sort.list(ans[,2]),]
  ans
}
extractAIC.multinom <- function(fit, ...) c(fit$edf, c(fit$AIC))

vcov.multinom <- function(object)
{
  if(is.null(object$Hessian)) {
    cat("\nRe-fitting to get Hessian\n\n")
    object <- update(object, Hess=T, trace=F)
  }
  solve(object$Hessian)
}

summary.multinom <- function(object, correlation=T, digits=NULL, ...)
{        
  if(!is.null(cl <- object$call)) {
    cat("Call:\n")
    dput(cl)
  }
  if(is.null(digits)) digits <- options()$digits
  vc <- vcov(object)
  r <- length(object$vcoefnames)
  cat("\nCoefficients:\n")
  if(length(object$lev) == 2) {
    coef <- cbind(object$wts[1+(1:r)], sqrt(diag(vc)))
    dimnames(coef) <- list(object$vcoefnames, c("Values", "Std. Error"))
    print(coef, digits=digits)
  } else {
    coef <- matrix(object$wts, nrow = object$n[3], byrow=T)[, 1+(1:r), drop=F]
    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=F]
    print(coef, digits=digits)
    cat("\nStd. Errors:\n")
    coef <- matrix(sqrt(diag(vc)), nrow = object$n[3]-1, byrow=T)
    if(length(object$lev)) dimnames(coef) <- list(object$lev[-1], object$vcoefnames)
    if(length(object$lab)) dimnames(coef) <- list(object$lab[-1], object$vcoefnames)
    print(coef, digits=digits)    
  }
  cat("\nResidual Deviance:", format(object$deviance), "\n")
  cat("AIC:", format(object$AIC), "\n")
  if(correlation) {
    stddev<- sqrt(diag(vc))
    correl <- vc / outer(stddev, stddev)
    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, drop = F], quote = F, ...)
    }
  }
  invisible()
}
# file nnet/nnet.q copyright (C) 1994-8 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(codes(cl)) - 1)] <- 1
    dimnames(x) <- list(names(cl), levels(cl))
    x
  }
  m <- match.call(expand.dots = F)
  if(is.matrix(eval(m$data, sys.parent())))
      m$data <- as.data.frame(data)
  m$... <- m$contrasts <- NULL
  m[[1]] <- as.name("model.frame")
  m <- eval(m, sys.parent())
  Terms <- attr(m, "terms")
  x <- model.matrix(Terms, m, contrasts)
  xint <- match("(Intercept)", dimnames(x)[[2]], nomatch=0)
  if(xint > 0) x <- x[, -xint, drop=F] # 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(codes(y)) - 1
      res <- nnet.default(x, y, w, entropy=T, ...)
      res$lev <- lev
    } else {
      y <- class.ind(y)
      res <- nnet.default(x, y, w, softmax=T, ...)
      res$lev <- lev
    } 
  } else res <- nnet.default(x, y, w, ...)
  res$terms <- Terms
  res$coefnames <- dimnames(x)[[2]]
  res$call <- match.call()
  class(res) <- c("nnet.formula", "nnet")
  res
}

nnet.default <-
function(x, y, weights, size, Wts, mask=rep(1, length(wts)), 
	 linout=F, entropy=F, softmax=F, censored=F, skip=F, 
	 rang=0.7, decay=0, maxit=100, Hess=F, trace=T, MaxNWts=1000)
{
  net <- NULL
  x <- as.matrix(x)
  y <- as.matrix(y)
  if(length(which.na(x))) stop("missing values in x")
  if(length(which.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 <- T; entropy <- F}
  if(censored) {linout <- T; entropy <- F; softmax <- T}
  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(mask) != length(wts)) stop("incorrect length of mask")
  if(length(wts) != nwts) stop("weights vector of incorrect length")
  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 <- .C("VR_set_train",
	  as.integer(ntr),
	  as.double(cbind(x,y)),
	  as.double(weights)
	  )
  tmp <- .C("VR_dovm",
	    as.integer(length(wts)),
	    wts=as.double(wts),
	    val=double(1),
	    as.integer(maxit),
	    trace,
	    as.integer(mask)
	    )
  .C("VR_unset_train")
  net$value <- tmp$val
  net$wts <- tmp$wts
  tmp <- matrix(.C("VR_nntest",
		   as.integer(ntr),
		   as.double(cbind(x,y)),
		   tclass = double(ntr*nout),
		   as.double(net$wts)
		   )$tclass,  ntr, nout)
  dimnames(tmp) <- list(dimnames(x)[[1]], dimnames(y)[[2]])
  net$fitted.values <- tmp
  .C("VR_unset_net")
  if(entropy) net$lev <- c("0","1")
  if(softmax) net$lev <- dimnames(y)[[2]]
  net$call <- match.call()
  if(Hess) net$Hessian <- nnet.Hess(net, x, y, weights)
  structure(net, class="nnet")
}

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)", dimnames(x)[[2]], nomatch=0)
  if(xint > 0) x <- x[, -xint, drop=F] # 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(length(which.na(x))) stop("missing values in x")
    x <- as.matrix(x)
    ntr <- dim(x)[1]
    nout <- object$n[3]
    if(length(object)==10) object$softmax <- F
    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(dimnames(x)[[1]], 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 <- .C("VR_set_train",
	  as.integer(ntr),
	  as.double(cbind(x,y)),
	  as.double(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*(codes(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")
  if(length(x)==10) x$softmax <- F
  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 <- format(round(object$wts,2))
  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) <- "summary.nnet"
  object
}

print.summary.nnet <- function(x, ...)
{
  if(length(x)==10) x$softmax <- F
  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(x$wts,2))
  wm <- c("b", paste("i", seq(length=x$n[1]), sep=""))
  if(x$n[2] > 0) wm <- c(wm, paste("h", seq(length=x$n[2]), sep=""))
  if(x$n[3] > 1) 
    wm <- c(wm, paste("o", seq(length=x$n[3]), sep=""))
  else wm <- c(wm, "o")
  names(wts) <- apply(cbind(wm[1+x$conn], wm[1+rep(1:x$nunits - 1, tconn)]), 
		      1, function(x)  paste(x, collapse = "->"))
  lapply(split(wts,rep(1:x$nunits, tconn)), 
	 function(x) print(x, quote=F))
  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
}
which.na <- function(x) seq(along = x)[is.na(x)]

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

if (version$minor < "0.62")
  library.dynam("nnet")
