drawNNet <- function(nnet,gray=FALSE) {

  W <- nnet$W
  v <- W[[1]] # hidden layer
  vmax <- max(abs(v))
  v <- v / vmax
  w <- W[[2]] # output layer
  wmax <- max(abs(w))
  w <- w / wmax

  ni <- nrow(v) - 1
  nh <- ncol(v)
  no <- ncol(w)
  allw <- rbind( cbind(NA, v, NA), rep(NA,nh+2), cbind(t(w),NA))
  nr <- nrow(allw)
  nc <- ncol(allw)

  ## Set up plot
  plot(NA, NA, type = "n", xlim=c(0.5,nc+1.5), ylim=c(0.5,nr+0.5), axes=FALSE,
       asp=1, xlab="",ylab="")

  ## vertical lines
  xs <- 1:nh + 1
  y1 <- 0.5
  y2 <- nr+0.5
  segments(1,0.5,1,no+0.5)
  text(1,no+1,"1")
  segments(xs,y1,xs,y2)
  ## horizontal lines
  ## hidden layer
  ys <- 1:(ni+1) + no + 1
  x1 <- 1
  x2 <- nh + 1.5
  segments(x1,ys,x2,ys)
  ## input labels
  for (i in 1:ni) 
   text(x1-0.4,ys[ni+no-1 - i],bquote(x[.(i)]))
  text(x1-0.5,ys[ni+1],"1")
  ## output layer
  ys <- 1:no
  x2 <- nh + 2.5
  segments(x1,ys,x2,ys)
  ## output labels
  for (i in 1:no) 
   text(x2+0.4,ys[no-i+1],bquote(y[.(i)]))
  ## cell bodies
  r <- 0.2
  hTri <- list(x=c(-r,r,0), y=c(r,r,-r))
  hTri$x <- hTri$x + 1
  hTri$y <- hTri$y + no + 1
  for (i in 1:nh) {
    hTri$x <- hTri$x + 1
    polygon(hTri,col="gray")
  }
  oTri <- list(x=c(-r,-r,r), y=c(-r,r,0))
  oTri$x <- oTri$x + nh + 2
  for (i in 1:no) {
    oTri$y <- oTri$y + 1
    polygon(oTri,col="gray")
  }

  ## Draw the weights
  allwc <- c(t(allw[nr:1,]))
  if (gray)
    colors <- ifelse(allwc < 0, "black", "gray")
  else
    colors <- ifelse(allwc < 0, "red", "green")

  symbols(expand.grid(1:nc,1:nr), squares=abs(allwc), bg=colors, #fg=NA,
          inches=FALSE, add=TRUE)

  ## Draw max values
  text(nh+1.5,1+no+1+ni/2,bquote(v[max]==.(round(vmax,3))),pos=4)
  text(1+(nh+1)/2,0.1,bquote(w[max]==.(round(wmax,3))),pos=3)
}

drawNNetDemo <- function() {
  
  ni <- 5
  nh <- 10
  no <- 2
  makew <- function(nr,nc) { matrix(1:(nr*nc) - nr*nc/2, nr,nc,byrow=TRUE)}
  makewr <- function(nr,nc) { matrix(runif(nr*nc)-0.5, nr,nc,byrow=TRUE)}
  net <- list(W = list(makew(ni+1,nh), makew(nh+1,no)))

#  x11(type="Xlib")
  
  p <- par(mar=c(0,0,0,0))
  drawNNet(net)
  par(p)

  print("Press enter")
  scan()

  for (i in 1:100) {
    net <- list(W = list(makewr(ni+1,nh), makewr(nh+1,no)))
    drawNNet(net)
    system("sleep 0.01")
  }
  
}

makeStandardizeF <- function(X) {
  if (missing(X)) {
    cat("Usage:
         standardize <- makeStandardizeF(X)  ## X is nSamples x nDimensions
         Xs <- standardize(X)
         X2s <- standardize(X2)\n")
    return(invisible())
  }
  ## X is nSamples x nDimensions
  mu <- colMeans(X)
  sigma <- sd(X) ##sd should be named colSds

  function(newX) {
    nr <- nrow(newX)
    nc <- ncol(newX)
    (newX - matrix(mu,nr,nc,byrow=TRUE)) / matrix(sigma,nr,nc,byrow=TRUE)
  }
}