multinom <- function(formula, data=sys.parent(), weights, subset, na.action,
	    contrasts=NULL, Hess=FALSE, summ=0, ...)
{
  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))
    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$... <- NULL
  m[[1]] <- as.name("model.frame")
  m <- eval(m, sys.parent())
  Terms <- attr(m, "terms")
  if(length(formals(model.matrix)) == 2) {##- For older version of R (<= 0.60):
   if(!is.null(contrasts))
     warning("`contrasts' are silently disregarded, currently")
   X <- model.matrix(Terms, m)
  } else { #-- new model.matrix --
    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]
    }
    Y <- if(length(lev) == 2) as.vector(codes(Y)) - 1 else 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,]
    Y <- if(is.matrix(Y)) Y[oX,] 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)) {
    p <- ncol(Y)
    sY <- Y %*% rep(1, p)
    if(any(sY==0)) stop("some case has no observations")
    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, 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, 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=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, ...)
      }
    } 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=TRUE, softmax=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 <- 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)
  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) {
     if(length(fit$lev)) bf <- fit$lev else bf <- fit$lab
     cf <- outer(bf[-1], cf, function(x,y)paste(x,y,sep=":"))
    }
    dimnames(fit$Hessian) <- list(cf, cf)
  }
  fit
}

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

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

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

drop1.multinom <- function(fit, scope, sorted=TRUE, trace=FALSE)
{
  if(!inherits(fit, "multinom")) stop("Not a multinom fit")
  if(missing(scope)) scope <- drop.scope(fit)
    else {
      if(!is.character(scope))
	scope <- attr(terms(update.formula(fit, scope)), "term.labels")
      if(!all(match(scope, attr(fit$terms, "term.labels"), FALSE)))
	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(fit$edf, fit$AIC)
  i <- 2
  for(tt in scope) {
    cat("trying -", tt,"\n")
    nfit <- update(fit, paste("~ . -", tt), trace=trace)
    if(nfit$edf == fit$edf) nfit$AIC <- NA
    ans[i,] <- c(nfit$edf, nfit$AIC)
    i <- i+1
  }
  if(sorted) ans <- ans[sort.list(ans[,2]),]
  ans
}

add1.multinom <- function(fit, scope, sorted=TRUE, trace=FALSE)
{
  if(!inherits(fit, "multinom")) stop("Not a multinom fit")
  if(!is.character(scope))
    scope <- add.scope(fit, update.formula(fit, scope,
					   evaluate = FALSE))
  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(fit$edf, fit$AIC)
  i <- 2
  for(tt in scope) {
    cat("trying +", tt,"\n")
    nfit <- update(fit, paste("~ . +", tt), trace=trace)
    if(nfit$edf == fit$edf) nfit$AIC <- NA
    ans[i,] <- c(nfit$edf, nfit$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(fit)
{
  if(is.null(fit$Hessian)) {
    cat("\nRe-fitting to get Hessian\n\n")
    fit <- update(fit, Hess=TRUE, trace=FALSE)
  }
  solve(fit$Hessian)
}

summary.multinom <- function(fit, correlation=TRUE, digits=NULL, ...)
{
  if(!is.null(cl <- fit$call)) {
    cat("Call:\n")
    dput(cl)
  }
  if(is.null(digits)) digits <- options()$digits
  vc <- vcov(fit)
  r <- length(fit$vcoefnames)
  cat("\nCoefficients:\n")
  if(length(fit$lev) == 2) {
    coef <- cbind(fit$wts[1+(1:r)], sqrt(diag(vc)))
    dimnames(coef) <- list(fit$vcoefnames, c("Values", "Std. Error"))
    print(coef, digits=digits)
  } else {
    coef <- matrix(fit$wts, nrow = fit$n[3], byrow=TRUE)[, 1+(1:r), drop=FALSE]
    if(length(fit$lev)) dimnames(coef) <- list(fit$lev, fit$vcoefnames)
    if(length(fit$lab)) dimnames(coef) <- list(fit$lab, fit$vcoefnames)
    coef <- coef[-1, , drop=FALSE]
    print(coef, digits=digits)
    cat("\nStd. Errors:\n")
    coef <- matrix(sqrt(diag(vc)), nrow = fit$n[3]-1, byrow=TRUE)
    if(length(fit$lev)) dimnames(coef) <- list(fit$lev[-1], fit$vcoefnames)
    if(length(fit$lab)) dimnames(coef) <- list(fit$lab[-1], fit$vcoefnames)
    print(coef, digits=digits)
  }
  cat("\nResidual Deviance:", format(fit$deviance), "\n")
  cat("AIC:", format(fit$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 = FALSE], quote = FALSE, ...)
    }
  }
  invisible()
}
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)
{
  m <- match.call(expand.dots = FALSE)
  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")
  attr(Terms, "intercept") <- 0   # Bias term is used for intercepts

  if(length(formals(model.matrix)) == 2) {##- For older version of R (<= 0.60):
   if(!is.null(contrasts))
     warning("`contrasts' are silently disregarded, currently")
   x <- model.matrix(Terms, m)
  } else { #-- new model.matrix --
    x <- model.matrix(Terms, m, contrasts)
  }
  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=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 <- 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=FALSE, entropy=FALSE, softmax=FALSE, skip=FALSE,
	 rang=0.7, decay=0, maxit=100, Hess=FALSE, trace=TRUE)
{
  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 <- TRUE; entropy <- FALSE}
  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(length(net$conn)==0) stop("No weights to fit")
  nsunits <- net$nunits
  if(linout) nsunits <- net$nunits - net$n[3]
  net$nsunits <- nsunits
  net$decay <- decay
  net$entropy <- entropy
  net$softmax <- softmax
  if(missing(Wts))
    if(rang > 0) wts <- runif(length(net$conn), -rang, rang)
    else wts <- rep(0, length(net$conn))
  else wts <- Wts
  if(length(mask) != length(wts)) stop("incorrect length of mask")
  if(length(wts) != length(net$conn))
    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)
     )
  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)
	    )
  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
  net$call <- match.call()
  if(Hess) net$Hessian <- nnet.Hess(net, x, y, weights)
  structure(net, class="nnet")
}

predict.nnet.formula <- function(net, newdata, type=c("raw", "class"))
{
  if(!inherits(net, "nnet.formula"))
    stop("object not of class nnet.formula")
  type <- match.arg(type)
  if(missing(newdata))
    switch(type,
	   raw = net$fitted.values,
	   class = {
	     z <- net$fitted.values
	     net$lev[if(ncol(z) > 1) max.col(z) else 1 + (z > 0.5)]
	   })
  else
    predict.nnet(net=net,
                 x=model.matrix(delete.response(net$terms), newdata),
                 type=type)
}

predict.nnet <- function(net, x, type=c("raw","class"))
{
  if(!inherits(net, "nnet")) stop("object not of class nnet")
  type <- match.arg(type)
  z <- if(missing(x)) net$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 <- net$n[3]
    ##? if(length(net)==10) net$softmax <- FALSE
    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(net$decay),
       as.integer(net$nsunits),
       as.integer(0),
       as.integer(net$softmax),
       )
    matrix(.C("VR_nntest",
              as.integer(ntr),
              as.double(x),
              tclass = double(ntr*nout),
              as.double(net$wts)
              )$tclass,
           ntr, nout,
           dimnames= list(dimnames(x)[[1]], dimnames(net$fitted)[[2]]))
  }
  switch(type,
         raw = z,
	 class = net$lev[if(ncol(z) > 1) max.col(z) else 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]
    na <- length(cadd <-c(if(nconn[i+1] == ns) 0,from))
    con <- c(if(ns > 1) conn[1:ns] else NULL,
             cadd,
             if(length(conn) > ns) conn[(ns+1):length(conn)])
    for(j in (i+1):net$nunits)
      nconn[j+1] <- nconn[j+1]+na
    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)
     )
  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)
	  )
  matrix(.C("VR_nnHessian",as.double(net$wts),H = double(nw*nw))$H,  nw, nw)
}

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

##--- These two functions  print.nnet & print.nnet.formula
##--- should really be `one'

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

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

coef.nnet <- function(net, ...)
{
  wts <- format(round(net$wts,2))#- this is ugly..
  wm <- c("b", paste("i", seq(length=net$n[1]), sep=""))
  if(net$n[2] > 0)
    wm <- c(wm, paste("h", seq(length=net$n[2]), sep=""))
  wm <- c(wm, paste("o", if(net$n[3] > 1) seq(length=net$n[3]), sep=""))
  names(wts) <- apply(cbind(wm[1+net$conn], wm[1+rep(1:net$nunits -1, tconn)]),
		      1, function(x)  paste(x, collapse = "->"))
  wts
}

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

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

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),
     n,
     ncol(m),
     rmax = integer(n))$rmax
}
which.na <- function(x) seq(along = x)[is.na(x)]

library.dynam("nnet")
