biplot <- function(x, ...) UseMethod("biplot")

biplot.default <-
    function(x, y, var.axes = TRUE, col, cex = rep(par("cex"), 2),
	     xlabs = NULL, ylabs = NULL, expand=1, xlim = NULL, ylim = NULL,
	     arrow.len = 0.1,
             main = NULL, sub = NULL, xlab = NULL, ylab = NULL, ...)
{
    n <- nrow(x)
    p <- nrow(y)
    if(missing(xlabs)) {
	xlabs <- dimnames(x)[[1]]
	if(is.null(xlabs)) xlabs <- 1:n
    }
    xlabs <- as.character(xlabs)
    dimnames(x) <- list(xlabs, dimnames(x)[[2]])
    if(missing(ylabs)) {
	ylabs <- dimnames(y)[[1]]
	if(is.null(ylabs)) ylabs <- paste("Var", 1:p)
    }
    ylabs <- as.character(ylabs)
    dimnames(y) <- list(ylabs, dimnames(y)[[2]])

    if(length(cex) == 1) cex <- c(cex, cex)
    if(missing(col)) {
	col <- par("col")
	if (!is.numeric(col)) col <- match(col, palette(), nomatch=1)
	col <- c(col, col + 1)
    }
    else if(length(col) == 1) col <- c(col, col)

    unsigned.range <- function(x) c(-abs(min(x)), abs(max(x)))
    rangx1 <- unsigned.range(x[, 1])
    rangx2 <- unsigned.range(x[, 2])
    rangy1 <- unsigned.range(y[, 1])
    rangy2 <- unsigned.range(y[, 2])

    if(missing(xlim) && missing(ylim))
	xlim <- ylim <- rangx1 <- rangx2 <- range(rangx1, rangx2)
    else if(missing(xlim)) xlim <- rangx1
    else if(missing(ylim)) ylim <- rangx2
    ratio <- max(rangy1/rangx1, rangy2/rangx2)/expand
    on.exit(par(op))
    op <- par(pty = "s")
    if(!is.null(main))
        op <- c(op, par(mar = par("mar")+c(0,0,1,0)))
    plot(x, type = "n", xlim = xlim, ylim = ylim, col = col[1],
         xlab = xlab, ylab = ylab, sub = sub, main = main, ...)
    text(x, xlabs, cex = cex[1], col = col[1], ...)
    par(new = TRUE)
    plot(y, axes = FALSE, type = "n", xlim = xlim*ratio, ylim = ylim*ratio,
	 xlab = "", ylab = "", col = col[1], ...)
    axis(3, col = col[2])
    axis(4, col = col[2])
    box(col = col[1])
    text(y, labels=ylabs, cex = cex[2], col = col[2], ...)
    if(var.axes)
	arrows(0, 0, y[,1] * 0.8, y[,2] * 0.8, col = col[2], length=arrow.len)
    invisible()
}

biplot.princomp <- function(x, choices = 1:2, scale = 1, pc.biplot=FALSE, ...)
{
    if(length(choices) != 2) stop("length of choices must be 2")
    if(!length(scores <- x$scores))
	stop(paste("object", deparse(substitute(x)), "has no scores"))
    lam <- x$sdev[choices]
    if(is.null(n <- x$n.obs)) n <- 1
    lam <- lam * sqrt(n)
    if(scale < 0 || scale > 1) warning("scale is outside [0, 1]")
    if(scale != 0) lam <- lam^scale else lam <- 1
    if(pc.biplot) lam <- lam / sqrt(n)
    biplot.default(t(t(scores[, choices]) / lam),
		   t(t(x$loadings[, choices]) * lam), ...)
    invisible()
}

biplot.prcomp <- function(x, choices = 1:2, scale = 1, pc.biplot=FALSE, ...)
{
    if(length(choices) != 2) stop("length of choices must be 2")
    if(!length(scores <- x$x))
	stop(paste("object", deparse(substitute(x)), "has no scores"))
    if(is.complex(scores))
        stop("biplots are not defined for complex PCA")
    lam <- x$sdev[choices]
    n <- NROW(scores)
    lam <- lam * sqrt(n)
    if(scale < 0 || scale > 1) warning("scale is outside [0, 1]")
    if(scale != 0) lam <- lam^scale else lam <- 1
    if(pc.biplot) lam <- lam / sqrt(n)
    biplot.default(t(t(scores[, choices]) / lam),
		   t(t(x$loadings[, choices]) * lam), ...)
    invisible()
}
## Seber pages 506-507, after a Golub original

cancor <- function(x, y, xcenter=TRUE, ycenter=TRUE)
{
    x <- as.matrix(x)
    y <- as.matrix(y)
    if((nr <- nrow(x)) != nrow(y)) stop("unequal number of rows in cancor")
    ncx <- ncol(x)
    ncy <- ncol(y)
    if(!nr || !ncx || !ncy) stop("dimension 0 in x or y")
    if(is.logical(xcenter)) {
	if(xcenter) {
	    xcenter <- colMeans(x,)
	    x <- x - rep(xcenter, rep(nr, ncx))
	}
	else xcenter <- rep(0, ncx)
    }
    else {
	xcenter <- rep(xcenter, length = ncx)
	x <- x - rep(xcenter, rep(nr, ncx))
    }
    if(is.logical(ycenter)) {
	if(ycenter) {
	    ycenter <- colMeans(y)
	    y <- y - rep(ycenter, rep(nr, ncy))
	}
	else ycenter <- rep(0, ncy)
    }
    else {
	ycenter <- rep(ycenter, length = ncy)
	y <- y - rep(ycenter, rep(nr,ncy))
    }
    qx <- qr(x)
    qy <- qr(y)
    dx <- qx$rank;	if(!dx) stop("`x' has rank 0")
    dy <- qy$rank;	if(!dy) stop("`y' has rank 0")
    ## compute svd(Qx'Qy)
    z <- La.svd(qr.qty(qx, qr.qy(qy, diag(1, nr, dy)))[1:dx,, drop = FALSE],
                dx, dy)
    xcoef <- backsolve((qx$qr)[1:dx, 1:dx, drop = FALSE], z$u)
    rownames(xcoef) <- colnames(x)[qx$pivot][1:dx]
    ycoef <-  backsolve((qy$qr)[1:dy, 1:dy, drop = FALSE], t(z$vt))
    rownames(ycoef) <- colnames(y)[qy$pivot][1:dy]
    list(cor = z$d, xcoef = xcoef, ycoef = ycoef, xcenter = xcenter,
	 ycenter = ycenter)
}
cmdscale <- function (d, k = 2, eig = FALSE, add = FALSE, x.ret = FALSE) {
    if (any(is.na(d)))
	stop("NA values not allowed in d")
    if (is.null(n <- attr(d, "Size"))) {
        if(add) d <- as.matrix(d)
	x <- as.matrix(d^2)
	if ((n <- nrow(x)) != ncol(x))
	    stop("Distances must be result of dist or a square matrix")
    }
    else {
	x <- matrix(0, n, n)
        if(add) d0 <- x
	x[row(x) > col(x)] <- d^2
	x <- x + t(x)
        if(add) {
            d0[row(x) > col(x)] <- d
            d <- d0 + t(d0)
        }
    }
    if((k <- as.integer(k)) > n - 1 || k < 1)
        stop("`k' must be in {1, 2, ..  n - 1}")
    storage.mode(x) <- "double"
    ## doubly center x in-place
    .C("dblcen", x, as.integer(n), DUP = FALSE, PACKAGE="mva")

    if(add) { ## solve the additive constant problem
        ## it is c* = largest eigenvalue of 2 x 2 (n x n) block matrix Z:
        i2 <- n + (i <- 1:n)
        Z <- matrix(0, 2*n, 2*n)
        Z[cbind(i2,i)] <- -1
        Z[ i, i2] <- -x
        Z[i2, i2] <- .C("dblcen", x = 2*d, as.integer(n), PACKAGE="mva")$x
        e <- La.eigen(Z,symmetric = FALSE, only.val = TRUE)$values
        add.c <- max(Re(e))
	x <- matrix(double(n*n), n, n)
        non.diag <- row(d) != col(d)
        x[non.diag] <- (d[non.diag] + add.c)^2
    }
    e <- La.eigen(-x/2, symmetric = TRUE)
    ev <- e$values[1:k]
    if(any(ev < 0))
        warning("some of the first ", k, " eigenvalues are < 0")
    points <- e$vectors[, 1:k, drop = FALSE] %*% diag(sqrt(ev), k)
    rn <- if(is.matrix(d)) rownames(d) else names(d)
    dimnames(points) <- list(rn, NULL)
    if (eig || x.ret || add) {
        evalus <- e$values[-n]
        list(points = points, eig = if(eig) ev, x = if(x.ret) x,
             ac = if(add) add.c else 0,
             GOF = sum(ev)/c(sum(abs(evalus)),
                             sum(evalus[evalus > 0])))
    } else points
}
cutree <- function(tree, k=NULL, h=NULL)
{
    if(is.null(n1 <- nrow(tree$merge)) || n1 < 1)
        stop("invalid tree (merge component)")
    n <- n1 + 1
    if(is.null(k) && is.null(h))
        stop("Either k or h must be specified")
    if(is.null(k)) {
        ## h |--> k
        k <- integer(length(h))
        ## S+6 help(cutree) says k(h) = k(h+), but does k(h-) [continuity]
        ## h < min() should give k = n;
        k <- n+1 - apply(outer(c(tree$height,Inf), h, ">"), 2, which.max)
        if(getOption("verbose")) cat("cutree(): k(h) = ",k,"\n")
    }
    else {
        k <- as.integer(k)
        if(min(k) < 1 || max(k) > n)
            stop(paste("Elements of k must be between 1 and", n))
    }

    ans <- .Call("R_cutree", tree$merge, k, PACKAGE = "mva")

    if(length(k) == 1) {
        ans <- as.vector(ans)
        names(ans) <- tree$labels
    }
    else{
        colnames(ans) <- if(!is.null(h)) h else k
        rownames(ans) <- tree$labels
    }
    return(ans)
}
as.dendrogram <- function(object, ...) UseMethod("as.dendrogram")

as.dendrogram.hclust <- function (object, hang = -1, ...)
## hang = 0.1  is default for plot.hclust
{
    if (is.null(object$labels))
	object$labels <- (1:length(object$order))
    z <- list()
    nMerge <- length(oHgt <- object$height)
    if (nMerge != nrow(object$merge))
	stop("`merge' and `height' do not fit!")
    hMax <- oHgt[nMerge]
    one <- 1:1;	   two <- 2:2 # integer!
    for (k in 1:nMerge) {
	x <- sort(object$merge[k, ])
	if (any(neg <- x < 0))
	    h0 <- if (hang < 0) 0 else max(0, oHgt[k] - hang * hMax)
	if (all(neg)) {			# two leaves
	    zk <- as.list(rev(-x))
	    attr(zk, "members") <- two
	    attr(zk, "midpoint") <- 0.5
	    objlabels <- rev(object$labels[-x])
	    attr(zk[[1]], "label") <- objlabels[1]
	    attr(zk[[2]], "label") <- objlabels[2]
	    attr(zk[[1]], "members") <- attr(zk[[2]], "members") <- one
	    attr(zk[[1]], "height") <- attr(zk[[2]], "height") <- h0
	    attr(zk[[1]], "leaf") <- attr(zk[[2]], "leaf") <- TRUE
	}
	else if (any(neg)) {		# one leaf, one node
	    X <- as.character(x)
	    ## For original "hclust" node ordering, the leaf is always left,
	    ## i.e. x[1]; don't want to assume this
	    isL <- x[1] < 0 ##is leaf left?
	    zk <-
		if(isL) list(-x[1], z[[X[2]]])
		else	list(z[[X[1]]], object$labels[-x[2]])
	    attr(zk, "members") <- attr(z[[X[1 + isL]]], "members") + one
	    attr(zk, "midpoint") <- (1 + attr(z[[X[1 + isL]]], "midpoint"))/2
	    attr(zk[[2 - isL]], "members") <- one
	    attr(zk[[2 - isL]], "height") <- h0
	    attr(zk[[2 - isL]], "label") <- object$labels[-x[1]]
	    attr(zk[[2 - isL]], "leaf") <- TRUE
	}
	else {				# two nodes
	    x <- as.character(x)
	    zk <- list(z[[x[1]]], z[[x[2]]])
	    attr(zk, "members") <- attr(z[[x[1]]], "members") +
		attr(z[[x[2]]], "members")
	    attr(zk, "midpoint") <- (attr(z[[x[1]]], "members") +
				     attr(z[[x[1]]], "midpoint") +
				     attr(z[[x[2]]], "midpoint"))/2
	}
	attr(zk, "height") <- oHgt[k]
	z[[k <- as.character(k)]] <- zk
    }
    z <- z[[k]]
    class(z) <- "dendrogram"
    z
}

### MM: `FIXME'	 (2002-05-14):
###	 =====
## We currently (mis)use a node's "members" attribute for two things:
## 1) #{sub nodes}
## 2) information about horizontal layout of the given node
## Because of that, cut.dend..() cannot correctly set "members" as it should!

## ==> start using "x.member" and the following function :

.memberDend <- function(x) {
    r <- attr(x,"x.member")
    if(is.null(r)) {
	r <- attr(x,"members")
	if(is.null(r)) r <- 1:1
    }
    r
}

.midDend <- function(x)
    if(is.null(mp <- attr(x, "midpoint"))) 0 else mp

midcache.dendrogram <- function (x, type = "hclust")
{
    ## Recompute the  "midpoint" attributes of a dendrogram, e.g. after reorder().

    type <- match.arg(type) ## currently only "hclust"
    if( !inherits(x, "dendrogram") )
        stop("we require a dendrogram")
    setmid <- function(x, type) {
        if(isLeaf(x))# no "midpoint"
            return(x)
        k <- length(x)
        if(k < 1)
            stop("dendrogram node with non-positive #{branches}")
        if(type == "hclust" && k != 2)
            warning("midcache()ing of non-binary dendrograms may be buggy")
        r <- x # incl. attributes!

        ## for(j in 1:k) {  r[[j]] <- remid(x[[j]]) .... }
        r[[1]] <- setmid(x[[1]], type)
        r[[2]] <- setmid(x[[2]], type)

        attr(r, "midpoint") <- (.memberDend(x[[1]]) +
                                .midDend(x[[1]]) + .midDend(x[[2]])) / 2
        r
    }
    setmid(x, type=type)
}


### Define a very concise print() method for dendrograms:
##  Martin Maechler, 15 May 2002
print.dendrogram <- function(x, digits = getOption("digits"), ...)
{
    cat("`dendrogram' ")
    if(isLeaf(x))
	cat("leaf", format(attr(x, "label"), digits = digits))
    else
	cat("with", length(x), "branches and",
	    attr(x,"members"), "members total")

    cat(", at height", format(attr(x,"height"), digits = digits), "\n")
    invisible(x)
}

str.dendrogram <-
function (object, max.level = 0, digits.d = 3, give.attr = FALSE,
          wid = getOption("width"), nest.lev = 0, indent.str = "", ...)
{
    ## FIXME: `wid' argument is currently disregarded
    pasteLis <- function(lis, dropNam, sep = " = ") {
	## drop uninteresting "attributes" here
	lis <- lis[!(names(lis) %in% dropNam)]
        fl <- sapply(lis, format, digits=digits.d, wid=wid)
        paste(paste(names(fl), fl, sep=sep), collapse = ", ")
    }

    nind <- nchar(istr <- indent.str)
    if(substr(istr, nind,nind) == " ")
       substr(istr, nind,nind) <- "`"
    cat(istr, "--", sep="")

    nAt <- names(at <- attributes(object))
    memb <- at[["members"]]
    hgt  <- at[["height"]]
    if(!isLeaf(object)) {
	le <- length(object)
        if(give.attr) {
            if(nchar(at <- pasteLis(at, c("class", "height", "members"))))
                at <- paste(",", at)
        }
	cat("[dendrogram w/ ", le, " branches and ", memb, " members at h = ",
            format(hgt, digits=digits.d), if(give.attr) at, "]\n", sep="")
	if (max.level==0 || nest.lev < max.level) {
	    for(i in 1:le) {
		##cat(indent.str, nam.ob[i], ":", sep="")
		str(object[[i]], nest.lev = nest.lev + 1,
		    indent.str= paste(indent.str, if(i < le) " |" else "  "),
		    max.level=max.level, digits.d=digits.d,
		    give.attr= give.attr, wid=wid)
	    }
	}
    } else { ## leaf
	cat("leaf",
	    if(is.character(at$label)) paste("",at$label,"",sep='"') else
	    format(object, digits=digits.d),"")
	any.at <- hgt != 0
	if(any.at) cat("(h=",format(hgt, digits=digits.d))
	if(memb != 1) #MM: when can this happen?
	    cat(if(any.at)", " else {any.at <- TRUE; "("}, "memb= ",memb,sep="")
        at <- pasteLis(at, c("class", "height", "members", "leaf", "label"))
        if(any.at || nchar(at)) cat(if(!any.at)"(", at, ")")
        cat("\n")
    }
}

## The ``generic'' method for "[["  (identical to e.g., "[[.POSIXct"):
## --> subbranches are dendrograms as well!
"[[.dendrogram" <- function(x, ..., drop = TRUE)
{
    cl <- class(x)
    class(x) <- NULL
    val <- NextMethod("[[")
    class(val) <- cl
    val
}


## Also: need larger par("mar")[1] or [4] for longish labels !
## [probably don't change, just print a warning or so !!
plot.dendrogram <-
    function (x, type = c("rectangle", "triangle"), center = FALSE,
	      edge.root = !is.null(attr(x, "edgetext")), nodePar = NULL,
	      xaxt="n", yaxt="s",
	      edgePar = list(), leaflab= c("perpendicular", "textlike", "none"),
	      xlab = "", ylab = "", horiz = FALSE, frame.plot = FALSE, ...)
{
    type <- match.arg(type)
    leaflab <- match.arg(leaflab)
    hgt <- attr(x, "height")
    if (edge.root && is.logical(edge.root))
	edge.root <- 0.0625 * hgt
    mem.x <- .memberDend(x)
    yTop <- hgt + edge.root
    if(center) { x1 <- 0.5 ; x2 <- mem.x + 0.5 }
    else       { x1 <- 1   ; x2 <- mem.x }
    xlim <- c(x1 - 1/2, x2 + 1/2)
    ylim <- c(0, yTop)
    if (horiz) {## swap and reverse direction on `x':
	xl <- xlim; xlim <- rev(ylim); ylim <- xl
	tmp <- xaxt; xaxt <- yaxt; yaxt <- tmp
    }
    plot(0, xlim = xlim, ylim = ylim, type = "n", xlab = xlab, ylab = ylab,
	 xaxt = xaxt, yaxt = yaxt, frame.plot = frame.plot, ...)
    if (edge.root) {
	x0 <- plotNodeLimit(x1, x2, x, center)$x
	if (horiz)
	    segments(hgt, x0, yTop, x0)
	else segments(x0, hgt, x0, yTop)
	if (!is.null(et <- attr(x, "edgetext"))) {
	    my <- mean(hgt, yTop)
	    if (horiz)
		text(my, x0, et)
	    else text(x0, my, et)
	}
    }
    plotNode(x1, x2, x, type = type, center = center, leaflab = leaflab,
	     nodePar = nodePar, edgePar = edgePar, horiz = horiz)
}

### the work horse: plot node (if pch) and lines to all children
plotNode <-
    function(x1, x2, subtree, type, center, leaflab,
             nodePar, edgePar, horiz = FALSE)
{
    inner <- !isLeaf(subtree) && x1 != x2
    yTop <- attr(subtree, "height")
    bx <- plotNodeLimit(x1, x2, subtree, center)
    xTop <- bx$x

    ## handle node specific parameters in "nodePar":
    hasP <- !is.null(nPar <- attr(subtree, "nodePar"))
    if(!hasP) nPar <- nodePar

    if(getOption("verbose")) {
	cat(if(inner)"inner node" else "leaf", ":")
	if(!is.null(nPar)) { cat(" with node pars\n"); str(nPar) }
	cat(if(inner)paste(" height", formatC(yTop),"; "),
	    "(x1,x2)= (",formatC(x1,wid=4),",",formatC(x2,wid=4),")",
	    "--> xTop=", formatC(xTop, wid=8),"\n", sep="")
    }

    Xtract <- function(nam, L, default, indx)
	rep(if(any(nam == names(L))) L[[nam]] else default, length = indx)[indx]

    if(!is.null(nPar)) { ## draw this node
	i <- if(inner || hasP) 1 else 2 # only 1 node specific par
	pch <- Xtract("pch", nPar, default = 1:2,	 i)
	cex <- Xtract("cex", nPar, default = c(1,1),	 i)
	col <- Xtract("col", nPar, default = par("col"), i)
	bg <- Xtract("bg", nPar, default = par("bg"), i)
	points(if (horiz) cbind(yTop, xTop) else cbind(xTop, yTop),
	       pch = pch, bg = bg, col = col, cex = cex)
    }
    if (isLeaf(subtree)) {
	## label leaf
	if (leaflab == "perpendicular") { # somewhat like plot.hclust
	    if(horiz) {
	       X <- yTop; Y <- xTop; srt <- 0; adj <- c(0, 0.5)
	    }
	    else{
	       X <- xTop; Y <- yTop; srt <- 90; adj <- 1
	    }
	    text(X, Y, attr(subtree,"label"), xpd = TRUE, srt = srt, adj = adj)
	}
    }
    else if (inner) {
	for (k in 1:length(subtree)) {
	    child <- subtree[[k]]
	    ## draw lines to the children and draw themselves recursively
	    yBot <- attr(child, "height")
	    if (getOption("verbose"))
		cat("ch.", k, "@ h=", yBot, "; ")
	    if (is.null(yBot))
		yBot <- 0
	    xBot <-
                if (center) mean(bx$limit[k:(k + 1)])
                else bx$limit[k] + .midDend(child)

	    hasE <- !is.null(ePar <- attr(child, "edgePar"))
	    if (!hasE)
		ePar <- edgePar
	    i <- if (!isLeaf(child) || hasE) 1 else 2
	    col <- Xtract("col", ePar, default = par("col"), i)
	    lty <- Xtract("lty", ePar, default = par("lty"), i)
	    lwd <- Xtract("lwd", ePar, default = par("lwd"), i)
	    segmentsHV <- function(x0, y0, x1, y1) {
		if (horiz)
		     segments(y0, x0, y1, x1, col = col, lty = lty, lwd = lwd)
		else segments(x0, y0, x1, y1, col = col, lty = lty, lwd = lwd)
	    }
	    if (type == "triangle") {
		segmentsHV(xTop, yTop, xBot, yBot)
	    }
	    else { # rectangle
		segmentsHV(xTop,yTop, xBot,yTop)# h
		segmentsHV(xBot,yTop, xBot,yBot)# v
	    }
	    vln <- NULL
	    if (isLeaf(child) && leaflab == "textlike") {
		nodeText <- as.character(attr(child,"label"))
		if(getOption("verbose")) cat('-- with "label"',nodeText)
		hln <- 0.6 * strwidth(nodeText)/2
		vln <- 1.5 * strheight(nodeText)/2
		rect(xBot - hln, yBot,
		     xBot + hln, yBot + 2 * vln, col = "white")
		text(xBot, yBot + vln, nodeText)
	    }
	    if (!is.null(attr(child, "edgetext"))) {
		edgeText <- as.character(attr(child, "edgetext"))
		if(getOption("verbose")) cat('-- with "edgetext"',edgeText)
		hlm <- 1.2 * strwidth(edgeText)/2
		vlm <- 1.5 * strheight(edgeText)/2
		if (!is.null(vln)) {
		  mx <- (xTop + xBot + ((xTop - xBot)/(yTop - yBot)) * vln)/2
		  my <- (yTop + yBot + 2 * vln)/2
		}
		else {
		  mx <- (xTop + xBot)/2
		  my <- (yTop + yBot)/2
	      }
		if (type == "triangle") {
		    polygon(c(mx - hlm, mx - 1.3*hlm, mx - hlm,
			      mx + hlm, mx + 1.3*hlm, mx + hlm),
			    c(my - vlm, my, my + vlm, my + vlm, my, my - vlm),
			    col = "white")
		  text(mx, my, edgeText)
		}
		else {
		  polygon(c(xBot - hlm, xBot - 1.3*hlm, xBot - hlm,
			    xBot + hlm, xBot + 1.3*hlm, xBot + hlm),
			  c(my - vlm, my, my + vlm, my + vlm, my, my - vlm),
			  col = "white")
		  text(xBot, my, edgeText)
		}
	    }
	    plotNode(bx$limit[k], bx$limit[k + 1], subtree = child,
                     type, center, leaflab, nodePar, edgePar, horiz)
	}
    }
}

plotNodeLimit <- function(x1, x2, subtree, center)
{
    ## get the left borders limit[k] of all children k=1..K, and
    ## the handle point `x' for the edge connecting to the parent.
    inner <- !isLeaf(subtree) && x1 != x2
    if(inner) {
	K <- length(subtree)
	mTop <- .memberDend(subtree)
	limit <- integer(K)
	xx1 <- x1
	for(k in 1:K) {
	    m <- .memberDend(subtree[[k]])
	    ##if(is.null(m)) m <- 1
	    xx1 <- xx1 + (if(center) (x2-x1) * m/mTop else m)
	    limit[k] <- xx1
	}
	limit <- c(x1, limit)
    } else { ## leaf
	limit <- c(x1, x2)
    }
    mid <- attr(subtree, "midpoint")
    center <- center || (inner && !is.numeric(mid))
    x <- if(center) mean(c(x1,x2)) else x1 + (if(inner) mid else 0)
    list(x = x, limit = limit)
}

cut.dendrogram <- function(x, h, ...)
{
    LOWER <- list()
    X <- 1

    assignNodes <- function(subtree, h) {
	if(!isLeaf(subtree)) {
	    if(!(K <- length(subtree)))
		warning("`subtree' of length 0 !!")
	    new.mem <- 0
	    for(k in 1:K) {
		if(attr(subtree[[k]], "height") <= h) {
		    ## cut it, i.e. save to LOWER[] and make a leaf
		    sub <- subtree[[k]]
		    at <- attributes(sub)
		    at$leaf <- TRUE
		    at$x.member <- at$members
		    new.mem <- new.mem + (at$members <- 1)
		    at$label <- paste("Branch", X)
		    subtree[[k]] <- paste("Branch", X)
		    attributes(subtree[[k]]) <- at
		    class(sub) <- "dendrogram"
		    LOWER[[X]] <<- sub
		    X <<- X+1
		}
		else { ## don't cut up here, possibly its children:
		    subtree[[k]] <- assignNodes(subtree[[k]], h)
		    new.mem <- new.mem + attr(subtree[[k]], "members")
		}
	    }
	    ## re-count members:
	    attr(subtree,"x.member") <- attr(subtree,"members")
	    attr(subtree,"members") <- new.mem
	}
	subtree
    }# assignNodes()

    list(upper = assignNodes(x, h), lower = LOWER)
}## cut.dendrogram()

isLeaf <- function(object) (is.logical(L <- attr(object, "leaf"))) && L


## *Not* a method (yet):
order.dendrogram <- function(x) {
    if( !inherits(x, "dendrogram") )
        stop("order.dendrogram requires a dendrogram")
    unlist(x)
}

##RG's first version -- for posterity
# order.dendrogram <- function(x) {
#    if( !inherits(x, "dendrogram") )
#       stop("order.dendrogram requires a dendrogram")
#    ord <- function(x) {
#      if( isLeaf(x) ) return(x[1])
#      return(c(ord(x[[1]]), ord(x[[2]])))
#    }
#   return(ord(x))
# }

reorder <- function(x, ...) UseMethod("reorder")

reorder.dendrogram <- function(x, wts, ...)
{
    if( !inherits(x, "dendrogram") )
        stop("we require a dendrogram")
    oV <- function(x, wts) {
        if( isLeaf(x) ) {
            attr(x, "value") <- wts[x[1]]
            return(x)
        }
###--- FIXME: This works only for binary dendrograms !!!

        if(length(x) != 2)
            stop("reorder()ing of non-binary dendrograms not yet implemented")
        ## recurse:
        left  <- oV(x[[1]], wts)
        right <- oV(x[[2]], wts)
        lV <- attr(left, "value")
        rV <- attr(right, "value")
        attr(x, "value") <- lV+rV
        if( lV > rV ) {
            x[[1]] <- right
            x[[2]] <- left
        } else {
            x[[1]] <- left
            x[[2]] <- right
        }
        x
    }
    midcache.dendrogram( oV(x, wts) )
}

rev.dendrogram <- function(x) {
    if(isLeaf(x))
	return(x)

    k <- length(x)
    if(k < 1)
	stop("dendrogram non-leaf node with non-positive #{branches}")
    r <- x # incl. attributes!
    for(j in 1:k) ## recurse
 	r[[j]] <- rev(x[[k+1-j]])
    midcache.dendrogram( r )
}

## original Andy Liaw; modified RG, MM

heatmap <-
function (x, Rowv=NULL, Colv=if(symm)"Rowv" else NULL,
          distfun = dist, hclustfun = hclust, add.expr,
          symm = FALSE, revC = identical(Colv, "Rowv"),
	  scale = c("row", "column", "none"), na.rm=TRUE,
	  margins = c(5, 5), ColSideColors, RowSideColors,
          cexRow = 0.2 + 1/log10(nr), cexCol = 0.2 + 1/log10(nc),
          labRow = NULL, labCol = NULL, main = NULL, xlab = NULL, ylab = NULL,
          ...)
{
    scale <- if(symm && missing(scale)) "none" else match.arg(scale)
    if(length(di <- dim(x)) != 2 || !is.numeric(x))
	stop("`x' must be a numeric matrix")
    nr <- di[1]
    nc <- di[2]
    if(nr <= 1 || nc <= 1)
	stop("`x' must have at least 2 rows and 2 columns")
    if(!is.numeric(margins) || length(margins) != 2)
	stop("`margins' must be a numeric vector of length 2")

    ## by default order by row/col means
    if(is.null(Rowv)) Rowv <- rowMeans(x, na.rm = na.rm)
    if(is.null(Colv)) Colv <- colMeans(x, na.rm = na.rm)

    ## get the dendrograms and reordering indices

    if(inherits(Rowv, "dendrogram"))
	ddr <- Rowv
    else {
	hcr <- hclustfun(distfun(x))
	ddr <- as.dendrogram(hcr)
	if(!is.logical(Rowv) || Rowv)
	    ddr <- reorder(ddr, Rowv)
    }
    if(nr != length(rowInd <- order.dendrogram(ddr)))
	stop("row dendrogram ordering gave index of wrong length")

    if(inherits(Colv, "dendrogram"))
	ddc <- Colv
    else if(identical(Colv, "Rowv")) {
        if(nr != nc)
            stop('Colv = "Rowv" but nrow(x) != ncol(x)')
        ddc <- ddr
    }
    else {
	hcc <- hclustfun(distfun(if(symm)x else t(x)))
	ddc <- as.dendrogram(hcc)
	if(!is.logical(Colv) || Colv)
	    ddc <- reorder(ddc, Colv)
    }
    if(nc != length(colInd <- order.dendrogram(ddc)))
	stop("column dendrogram ordering gave index of wrong length")

    ## reorder x
    x <- x[rowInd, colInd]

    if(is.null(labRow))
        labRow <- if(is.null(rownames(x))) (1:nr)[rowInd] else rownames(x)
    if(is.null(labCol))
        labCol <- if(is.null(colnames(x))) (1:nc)[colInd] else colnames(x)

    if(scale == "row") {
	x <- sweep(x, 1, rowMeans(x, na.rm = na.rm))
	sx <- apply(x, 1, sd, na.rm = na.rm)
	x <- sweep(x, 1, sx, "/")
    }
    else if(scale == "column") {
	x <- sweep(x, 2, colMeans(x, na.rm = na.rm))
	sx <- apply(x, 2, sd, na.rm = na.rm)
	x <- sweep(x, 2, sx, "/")
    }

    ## Calculate the plot layout
    lmat <- rbind(c(NA, 3), 2:1)
    lwid <- lhei <- c(1, 4)
    if(!missing(ColSideColors)) { ## add middle row to layout
	if(!is.character(ColSideColors) || length(ColSideColors) != nc)
	    stop("'ColSideColors' must be a character vector of length ncol(x)")
	lmat <- rbind(lmat[1,]+1, c(NA,1), lmat[2,]+1)
	lhei <- c(lhei[1], 0.2, lhei[2])
    }
    if(!missing(RowSideColors)) { ## add middle column to layout
	if(!is.character(RowSideColors) || length(RowSideColors) != nr)
	    stop("'RowSideColors' must be a character vector of length nrow(x)")
	lmat <- cbind(lmat[,1]+1, c(rep(NA, nrow(lmat)-1), 1), lmat[,2]+1)
	lwid <- c(lwid[1], 0.2, lwid[2])
    }
    lmat[is.na(lmat)] <- 0

    ## Graphics `output' -----------------------

    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    layout(lmat, widths = lwid, heights = lhei, respect = TRUE)
    ## draw the side bars
    if(!missing(RowSideColors)) {
	par(mar = c(margins[1],0, 0,0.5))
	image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
    }
    if(!missing(ColSideColors)) {
	par(mar = c(0.5,0, 0,margins[2]))
	image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
    }
    ## draw the main carpet
    par(mar = c(margins[1], 0, 0, margins[2]))
    if(!symm || scale != "none") x <- t(x)
    if(revC) { # x columns reversed
        iy <- nr:1
        ddr <- rev(ddr)
        x <- x[,iy]
    } else iy <- 1:nr

    image(1:nc, 1:nr, x, xlim = 0.5+ c(0, nc), ylim = 0.5+ c(0, nr),
          axes = FALSE, xlab = "", ylab = "", ...)
    axis(1, 1:nc, labels= labCol, las= 2, line= -0.5, tick= 0, cex.axis= cexCol)
    if(!is.null(xlab)) mtext(xlab, side = 1, line = margins[1] - 1.25)
    axis(4, iy, labels= labRow, las= 2, line= -0.5, tick= 0, cex.axis= cexRow)
    if(!is.null(ylab)) mtext(ylab, side = 4, line = margins[2] - 1.25)

    if (!missing(add.expr))
	eval(substitute(add.expr))

    ## the two dendrograms :
    par(mar = c(margins[1], 0, 0, 0))
    plot(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")

    par(mar = c(0, 0, if(!is.null(main)) 1 else 0, margins[2]))
    plot(ddc,		    axes = FALSE, xaxs = "i", leaflab = "none")
    if(!is.null(main)) title(main, cex.main = 1.5*op[["cex.main"]])

    invisible(list(rowInd = rowInd, colInd = colInd))
}
dist <- function(x, method="euclidean", diag=FALSE, upper=FALSE)
{
    ## account for possible spellings of euclid?an
    if(!is.na(pmatch(method, "euclidian")))
	method <- "euclidean"

    METHODS <- c("euclidean", "maximum",
		 "manhattan", "canberra", "binary")
    method <- pmatch(method, METHODS)
    if(is.na(method))
	stop("invalid distance method")
    if(method == -1)
	stop("ambiguous distance method")

    N <- nrow(x <- as.matrix(x))
    d <- .C("R_distance",
	    x = as.double(x),
	    nr= N,
	    nc= ncol(x),
	    d = double(N*(N - 1)/2),
	    diag  = as.integer(FALSE),
	    method= as.integer(method),
	    DUP = FALSE, NAOK=TRUE, PACKAGE="mva")$d
    attr(d, "Size") <- N
    attr(d, "Labels") <- dimnames(x)[[1]]
    attr(d, "Diag") <- diag
    attr(d, "Upper") <- upper
    attr(d, "method") <- METHODS[method]
    attr(d, "call") <- match.call()
    class(d) <- "dist"
    return(d)
}

names.dist <- function(x) attr(x, "Labels")

"names<-.dist" <- function(x, value)
{
    if(length(value) != attr(x, "Size"))
	stop("invalid names for dist object")
    attr(x, "Labels") <- value
    x
}

## Because names(d) != length(d) for "dist"-object d, we need
format.dist <- function(x, ...) format(as.vector(x), ...)

as.matrix.dist <- function(x)
{
    size <- attr(x, "Size")
    df <- matrix(0, size, size)
    df[row(df) > col(df)] <- x
    df <- df + t(df)
    labels <- attr(x, "Labels")
    dimnames(df) <-
	if(is.null(labels)) list(1:size,1:size) else list(labels,labels)
    df
}


as.dist <- function(m, diag = FALSE, upper = FALSE)
{
    if (inherits(m,"dist"))
	ans <- m
    else { ## matrix |-> dist
	m <- as.matrix(m)
	ans <- m[row(m) > col(m)]
	attributes(ans) <- NULL
	if(!is.null(rownames(m)))
	    attr(ans,"Labels") <- rownames(m)
	else if(!is.null(colnames(m)))
	    attr(ans,"Labels") <- colnames(m)
	attr(ans,"Size") <- nrow(m)
	attr(ans, "call") <- match.call()
	class(ans) <- "dist"
    }
    if(is.null(attr(ans,"Diag")) || !missing(diag))
	attr(ans,"Diag") <- diag
    if(is.null(attr(ans,"Upper")) || !missing(upper))
	attr(ans,"Upper") <- upper
    ans
}


print.dist <- function(x, diag = NULL, upper = NULL, ...)
{
    if(is.null(diag))
	diag <-	 if(is.null(a <- attr(x, "Diag"))) FALSE else a
    if(is.null(upper))
	upper <- if(is.null(a <- attr(x,"Upper"))) FALSE else a

    size <- attr(x, "Size")
    df <- as.matrix.dist(x)
    if(!upper)
	df[row(df) < col(df)] <- NA
    if(!diag)
	df[row(df) == col(df)] <- NA
    print(if(diag || upper) df else df[-1, -size], na = "", ...)
    invisible(x)
}
## Hmm, MM thinks diag(.) needs checking { diag(vec) when length(vec)==1 !}
factanal <-
    function (x, factors, data = NULL, covmat = NULL, n.obs = NA,
              subset, na.action, start = NULL,
              scores = c("none", "regression", "Bartlett"),
              rotation = "varimax",
              control = NULL, ...)
{
    sortLoadings <- function(Lambda)
    {
        cn <- colnames(Lambda)
        Phi <- attr(Lambda, "covariance")
        ssq <- apply(Lambda, 2, function(x) -sum(x^2))
        Lambda <- Lambda[, order(ssq), drop = FALSE]
        colnames(Lambda) <- cn
        neg <- colSums(Lambda) < 0
        Lambda[, neg] <- -Lambda[, neg]
        if(!is.null(Phi)) {
            unit <- ifelse(neg, -1, 1)
            attr(Lambda, "covariance") <-
                unit %*% Phi[order(ssq), order(ssq)] %*% unit
        }
        Lambda
    }
    cl <- match.call()
    na.act <- NULL
    if (is.list(covmat)) {
        if (any(is.na(match(c("cov", "n.obs"), names(covmat)))))
            stop("covmat is not a valid covariance list")
        cv <- covmat$cov
        n.obs <- covmat$n.obs
        have.x <- FALSE
    }
    else if (is.matrix(covmat)) {
        cv <- covmat
        have.x <- FALSE
    }
    else if (is.null(covmat)) {
        if(missing(x)) stop("neither x nor covmat supplied")
        have.x <- TRUE
        if(inherits(x, "formula")) {
            mt <- terms(x, data = data)
            if(attr(mt, "response") > 0)
                stop("response not allowed in formula")
            attr(mt, "intercept") <- 0
            mf <- match.call(expand.dots = FALSE)
            names(mf)[names(mf) == "x"] <- "formula"
            mf$factors <- mf$covmat <- mf$scores <- mf$start <-
                mf$rotation <- mf$control <- mf$... <- NULL
            mf[[1]] <- as.name("model.frame")
            mf <- eval(mf, parent.frame())
            na.act <- attr(mf, "na.action")
            z <- model.matrix(mt, mf)
        } else {
            z <- as.matrix(x)
            if(!missing(subset)) z <- z[subset, , drop = FALSE]
        }
        covmat <- cov.wt(z)
        cv <- covmat$cov
        n.obs <- covmat$n.obs
    }
    else stop("covmat is of unknown type")
    scores <- match.arg(scores)
    if(scores != "none" && !have.x)
        stop("requested scores without an x matrix")
    sds <- sqrt(diag(cv))
    cv <- cv/(sds %o% sds)
    p <- ncol(cv)
    dof <- 0.5 * ((p - factors)^2 - p - factors)
    if(dof < 0)
        stop(paste(factors, "factors is too many for", p, "variables"))

    cn <- list(nstart = 1, trace = FALSE, lower = 0.005)
    cn[names(control)] <- control
    more <- list(...)[c("nstart", "trace", "lower", "opt", "rotate")]
    if(length(more)) cn[names(more)] <- more

    if(is.null(start)) {
        start <- (1 - 0.5*factors/p)/diag(solve(cv))
        if((ns <- cn$nstart) > 1)
            start <- cbind(start, matrix(runif(ns-1), p, ns-1, byrow=TRUE))
    }
    start <- as.matrix(start)
    if(nrow(start) != p) stop(paste("start must have", p, "rows"))
    nc <- ncol(start)
    if(nc < 1) stop("no starting values supplied")

    best <- Inf
    for (i in 1:nc) {
        nfit <- factanal.fit.mle(cv, factors, start[, i],
                                 max(cn$lower, 0), cn$opt)
        if(cn$trace)
            cat("start", i, "value:", format(nfit$criteria[1]),
                "uniqs:", format(as.vector(round(nfit$uniquenesses, 4))), "\n")
        if(nfit$converged && nfit$criteria[1] < best) {
            fit <- nfit
            best <- fit$criteria[1]
        }
    }
    if(best == Inf) stop("Unable to optimize from these starting value(s)")
    load <- fit$loadings
    if(rotation != "none") {
        rot <- do.call(rotation, c(list(load), cn$rotate))
        load <- if(is.list(rot)) rot$loadings else rot
    }
    fit$loadings <- sortLoadings(load)
    class(fit$loadings) <- "loadings"
    fit$na.action <- na.act
    if(have.x && scores != "none") {
        Lambda <- fit$loadings
        zz <- scale(z, TRUE, TRUE)
        switch(scores,
               regression = {
                   sc <- zz %*% solve(cv, Lambda)
                   if(!is.null(Phi <- attr(Lambda, "covariance")))
                       sc <- sc %*% Phi
               },
               Bartlett = {
                   d <- 1/fit$uniquenesses
                   tmp <- t(Lambda * d)
                   sc <- t(solve(tmp %*% Lambda, tmp %*% t(zz)))
               })
        rownames(sc) <- rownames(z)
        colnames(sc) <- colnames(Lambda)
        if(!is.null(na.act)) sc <- napredict(na.act, sc)
        fit$scores <- sc
    }
    if(!is.na(n.obs) && dof > 0) {
        fit$STATISTIC <- (n.obs - 1 - (2 * p + 5)/6 -
                     (2 * factors)/3) * fit$criteria["objective"]
        fit$PVAL <- pchisq(fit$STATISTIC, dof, lower.tail = FALSE)
    }
    fit$n.obs <- n.obs
    fit$call <- cl
    fit
}

factanal.fit.mle <-
    function(cmat, factors, start=NULL, lower = 0.005, control = NULL, ...)
{
    FAout <- function(Psi, S, q)
    {
        sc <- diag(1/sqrt(Psi))
        Sstar <- sc %*% S %*% sc
        E <- La.eigen(Sstar, symmetric = TRUE)
        L <- E$vectors[, 1:q, drop = FALSE]
        load <- L %*% diag(sqrt(pmax(E$values[1:q] - 1, 0)), q)
        diag(sqrt(Psi)) %*% load
    }
    FAfn <- function(Psi, S, q)
    {
        sc <- diag(1/sqrt(Psi))
        Sstar <- sc %*% S %*% sc
        E <- La.eigen(Sstar, symmetric = TRUE, only.values = TRUE)
        e <- E$values[-(1:q)]
        e <- sum(log(e) - e) - q + nrow(S)
##        print(round(c(Psi, -e), 5))  # for tracing
        -e
    }
    FAgr <- function(Psi, S, q)
    {
        sc <- diag(1/sqrt(Psi))
        Sstar <- sc %*% S %*% sc
        E <- La.eigen(Sstar, symmetric = TRUE)
        L <- E$vectors[, 1:q, drop = FALSE]
        load <- L %*% diag(sqrt(pmax(E$values[1:q] - 1, 0)), q)
        load <- diag(sqrt(Psi)) %*% load
        g <- load %*% t(load) + diag(Psi) - S
        diag(g)/Psi^2
    }
    p <- ncol(cmat)
    if(is.null(start))
        start <- (1 - 0.5*factors/p)/diag(solve(cmat))
    res <- optim(start, FAfn, FAgr, method = "L-BFGS-B",
                 lower = lower, upper = 1,
                 control = c(list(fnscale=1,
                 parscale = rep(0.01, length(start))), control),
                 q = factors, S = cmat)
    Lambda <- FAout(res$par, cmat, factors)
    dimnames(Lambda) <- list(dimnames(cmat)[[1]],
                             paste("Factor", 1:factors, sep = ""))
    p <- ncol(cmat)
    dof <- 0.5 * ((p - factors)^2 - p - factors)
    un <- res$par
    names(un) <- colnames(cmat)
    class(Lambda) <- "loadings"
    ans <- list(converged = res$convergence == 0,
                loadings = Lambda, uniquenesses = un,
                correlation = cmat,
                criteria = c(objective = res$value, counts = res$counts),
                factors = factors, dof = dof, method = "mle")
    class(ans) <- "factanal"
    ans
}

print.loadings <- function(x, digits = 3, cutoff = 0.1, sort = FALSE, ...)
{
    Lambda <- unclass(x)
    p <- nrow(Lambda)
    factors <- ncol(Lambda)
    if (sort) {
        mx <- max.col(abs(Lambda))
        ind <- cbind(1:p, mx)
        mx[abs(Lambda[ind]) < 0.5] <- factors + 1
        Lambda <- Lambda[order(mx, 1:p),]
    }
    cat("\nLoadings:\n")
    fx <- format(round(Lambda, digits))
    names(fx) <- NULL
    nc <- nchar(fx[1])
    fx[abs(Lambda) < cutoff] <- paste(rep(" ", nc), collapse = "")
    print(fx, quote = FALSE, ...)
    vx <- colSums(x^2)
    varex <- rbind("SS loadings" = vx)
    if(is.null(attr(x, "covariance"))) {
        varex <- rbind(varex, "Proportion Var" = vx/p)
        if(factors > 1)
            varex <- rbind(varex, "Cumulative Var" = cumsum(vx/p))
    }
    cat("\n")
    print(round(varex, digits))
    invisible(x)
}

print.factanal <- function(x, digits = 3, ...)
{
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    cat("Uniquenesses:\n")
    print(round(x$uniquenesses, digits), ...)
    print(x$loadings, digits = digits, ...)
    if(!is.null(x$STATISTIC)) {
        factors <- x$factors
        cat("\nTest of the hypothesis that", factors, if(factors == 1)
            "factor is" else "factors are", "sufficient.\n")
        cat("The chi square statistic is", round(x$STATISTIC, 2), "on", x$dof,
            if(x$dof == 1) "degree" else "degrees",
            "of freedom.\nThe p-value is", signif(x$PVAL, 3), "\n")
    } else {
        cat(paste("\nThe degrees of freedom for the model is",
                  x$dof, "and the fit was", round(x$criteria["objective"], 4),
                  "\n"))
    }
    invisible(x)
}

varimax <- function(x, normalize = TRUE, eps = 1e-5)
{
    nc <- ncol(x)
    if(nc < 2) return(x)
    if(normalize) {
        sc <- sqrt(drop(apply(x, 1, function(x) sum(x^2))))
        x <- x/sc
    }
    p <- nrow(x)
    TT <- diag(nc)
    d <- 0
    for(i in 1:1000) {
        z <- x %*% TT
        B  <- t(x) %*% (z^3 - z %*% diag(drop(rep(1, p) %*% z^2))/p)
        sB <- La.svd(B)
        TT <- sB$u %*% sB$vt
        dpast <- d
        d <- sum(sB$d)
        if(d < dpast * (1 + eps)) break
    }
    z <- x %*% TT
    if(normalize) z <- z * sc
    dimnames(z) <- dimnames(x)
    list(loadings = z, rotmat = TT)
}

promax <- function(x, m = 4)
{
    if(ncol(x) < 2) return(x)
    dn <- dimnames(x)
    xx <- varimax(x)
    x <- xx$loadings
    Q <- x * abs(x)^(m-1)
    U <- lm.fit(x, Q)$coefficients
    d <- diag(solve(t(U) %*% U))
    U <- U %*% diag(sqrt(d))
    dimnames(U) <- NULL
    z <- x %*% U
    U <- xx$rotmat %*% U
    dimnames(z) <- dn
    list(loadings = z, rotmat = U)
}
## Hierarchical clustering, on raw input data; we will use Euclidean
## distance.  A range of criteria are supported; also there is a
## storage-economic option.
##
## We use the general routine, `hc', which caters for 7 criteria,
## using a half dissimilarity matrix; (BTW, this uses the very efficient
## nearest neighbor chain algorithm, which makes this algorithm of
## O(n^2) computational time, and differentiates it from the less
## efficient -- i.e. O(n^3) -- implementations in all commercial
## statistical packages -- as far as I am aware -- except Clustan.)
##
## Clustering Methods:
##
## 1. Ward's minimum variance or error sum of squares method.
## 2. single linkage or nearest neighbor method.
## 3. complete linkage or diameter.
## 4. average linkage, group average, or UPGMA method.
## 5. McQuitty's or WPGMA method.
## 6. median, Gower's or WPGMC method.
## 7. centroid or UPGMC method (7).
##
## Original author: F. Murtagh, May 1992
## R Modifications: Ross Ihaka, Dec 1996
##		    Friedrich Leisch, Apr 1998, Jun 2000

hclust <- function(d, method="complete", members=NULL)
{
    METHODS <- c("ward", "single",
                 "complete", "average", "mcquitty",
                 "median", "centroid")
    method <-  pmatch(method, METHODS)
    if(is.na(method))
	stop("invalid clustering method")
    if(method == -1)
	stop("ambiguous clustering method")

    n <- as.integer(attr(d, "Size"))
    if(is.null(n))
	stop("invalid dissimilarities")
    if(n < 2)
        stop("Must have n >= 2 objects to cluster")
    len <- as.integer(n*(n-1)/2)
    if(length(d) != len)
        (if (length(d) < len) stop else warning
         )("dissimilarities of improper length")

    if(is.null(members))
        members <- rep(1, n)
    else if(length(members) != n)
        stop("Invalid length of members")

    hcl <- .Fortran("hclust",
		    n = n,
		    len = len,
		    method = as.integer(method),
		    ia = integer(n),
		    ib = integer(n),
		    crit = double(n),
		    members = as.double(members),
		    nn = integer(n),
		    disnn = double(n),
		    flag = logical(n),
		    diss = as.double(d), PACKAGE="mva")

    ## 2nd step: interpret the information that we now have
    ## as merge, height, and order lists.

    hcass <- .Fortran("hcass2",
		      n = as.integer(n),
		      ia = as.integer(hcl$ia),
		      ib = as.integer(hcl$ib),
  		      order = integer(n),
		      iia = integer(n),
		      iib = integer(n), PACKAGE="mva")

    tree <- list(merge = cbind(hcass$iia[1:(n-1)], hcass$iib[1:(n-1)]),
		 height= hcl$crit[1:(n-1)],
		 order = hcass$order,
		 labels=attr(d, "Labels"),
                 method=METHODS[method],
                 call = match.call(),
                 dist.method = attr(d, "method"))
    class(tree) <- "hclust"
    tree
}

plot.hclust <-
    function (x, labels = NULL, hang = 0.1,
              axes = TRUE, frame.plot = FALSE, ann = TRUE,
              main = "Cluster Dendrogram",
              sub = NULL, xlab = NULL, ylab = "Height", ...)
{
    merge <- x$merge
    if (!is.matrix(merge) || ncol(merge) != 2)
	stop("invalid dendrogram")
    n <- nrow(merge)
    height <- as.double(x$height)
    ord <- as.double(order(x$order))

    labels <-
	if(missing(labels) || is.null(labels)) {
	    if (is.null(x$labels))
		paste(1:(n+1))
	    else
		as.character(x$labels)
	} else {
	    if(is.logical(labels) && !labels)# FALSE
		character(n+1)
	    else
		as.character(labels)
	}

    plot.new()
    .Internal(dend.window(n, merge, height, ord, hang, labels, ...))
    .Internal(dend       (n, merge, height, ord, hang, labels, ...))
    if(axes)
        axis(2, at=pretty(range(height)))
    if (frame.plot)
        box(...)
    if (ann) {
        if(!is.null(cl <- x$call) && is.null(sub))
            sub <- paste(deparse(cl[[1]])," (*, \"", x$method,"\")",sep="")
        if(is.null(xlab) && !is.null(cl))
            xlab <- deparse(cl[[2]])
        title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...)
    }
    invisible()
}

## For S ``compatibility'': was in cluster just as
## plclust <- plot.hclust ## .Alias
plclust <- function(tree, hang = 0.1, unit = FALSE, level = FALSE, hmin = 0,
                    square = TRUE, labels = NULL, plot. = TRUE,
                    axes = TRUE, frame.plot = FALSE, ann = TRUE,
                    main = "", sub = NULL, xlab = NULL, ylab = "Height")
{
    if(!missing(level) && level)	.NotYetUsed("level", error = FALSE)
    if(!missing(hmin) && hmin != 0)	.NotYetUsed("hmin",  error = FALSE)
    if(!missing(square) && !square)	.NotYetUsed("square",error = FALSE)
    if(!missing(plot.) && !plot.)	.NotYetUsed("plot.", error = TRUE)
    if(!missing(hmin)) tree$height <- pmax(tree$height, hmin)
    if(unit) tree$height <- rank(tree$height)
    plot.hclust(x = tree, labels = labels, hang = hang,
                axes = axes, frame.plot = frame.plot, ann = ann,
                main = main, sub = sub, xlab = xlab, ylab = ylab)
}


as.hclust <- function(x, ...) UseMethod("as.hclust")
## need *.default for idempotency:
as.hclust.default <- function(x, ...) {
    if(inherits(x, "hclust")) x
    else
	stop("argument `x' cannot be coerced to class `hclust'.",
             if(!is.null(oldClass(x)))
             "\n Consider providing an as.hclust.",oldClass(x)[1],"() method")
}

as.hclust.twins <- function(x, ...)
{
    r <- list(merge = x$merge,
	      height = sort(x$height),
	      order = x$order,
	      labels = if(!is.null(lb <- x$order.lab)) {
                  lb[sort.list(x$order)] } else rownames(x$data),# may be NULL
	      call = if(!is.null(cl <- x$call)) cl else match.call(),
	      method = if(!is.null(mt <- x$method)) mt else NA,
	      dist.method = attr(x$diss, "Metric"))
    class(r) <- "hclust"
    r
}

print.hclust <- function(x, ...)
{
    if(!is.null(x$call))
        cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
    if(!is.null(x$method))
        cat("Cluster method   :", x$method, "\n")
    if(!is.null(x$dist.method))
        cat("Distance         :", x$dist.method, "\n")
    cat("Number of objects:", length(x$height)+1, "\n")
    cat("\n")
}

cophenetic <- function(x) {
    x <- as.hclust(x)
    nobs <- length(x$order)
    ilist <- vector("list", length=nobs)
    names(ilist) <- 1:nobs # FIXME: do better when you can!
    rmat <- matrix(NA, nr=nobs, nc=nobs)
    for( i in 1:(nobs-1)) {
        inds <- x$merge[i,]
        ids1 <- if(inds[1] < 0) -inds[1] else ilist[[inds[1]]]
        ids2 <- if(inds[2] < 0) -inds[2] else ilist[[inds[2]]]
        ilist[[i]] <- c(ids1, ids2)
        for( ival1 in ids1)
            for( ival2 in ids2 ){
                if( ival1 > ival2 )
                    rmat[ival1, ival2] <- x$height[i]
                else
                    rmat[ival2, ival1] <- x$height[i]
            }
    }
    return(as.dist(rmat))
}
rect.hclust <- function(tree, k=NULL, which=NULL,
                        x=NULL, h=NULL, border=2, cluster=NULL)
{
    if(length(h)>1 | length(k)>1)
        stop("k and h must be a scalar")

    if(!is.null(h)){
        if(!is.null(k))
            stop("specify exactly one of k and h")
        k <- min(which(rev(tree$height)<h))
        k <- max(k, 2)
    }
    else
        if(is.null(k))
            stop("specify exactly one of k and h")

    if(k < 2 | k > length(tree$height))
        stop(paste("k must be between 2 and", length(tree$height)))

    if(is.null(cluster))
        cluster <- cutree(tree, k=k)
    ## cutree returns classes sorted by data, we need classes
    ## as occurring in the tree (from left to right)
    clustab <- table(cluster)[unique(cluster[tree$order])]
    m <- c(0, cumsum(clustab))

    if(!is.null(x)){
        if(!is.null(which))
            stop("specify exactly one of which and x")
        which <- x
        for(n in 1:length(x))
            which[n] <- max(which(m<x[n]))
    }
    else
        if(is.null(which))
            which <- 1:k

    if(any(which>k))
        stop(paste("all elements of which must be between 1 and", k))

    border <- rep(border, length=length(which))

    retval <- list()
    for(n in 1:length(which)){
        rect(m[which[n]]+0.66, par("usr")[3],
             m[which[n]+1]+0.33, mean(rev(tree$height)[(k-1):k]),
             border = border[n])
        retval[[n]] <- which(cluster==as.integer(names(clustab)[which[n]]))
    }
    invisible(retval)
}

identify.hclust <- function(x, FUN = NULL, N = 20, MAXCLUSTER = 20,
                            DEV.FUN = NULL, ...)
{
    cluster <- cutree(x, k = 2:MAXCLUSTER)

    retval <- list()
    oldk <- NULL
    oldx <- NULL
    DEV.x <- dev.cur()

    for(n in 1:N){

        dev.set(DEV.x)
        X <- locator(1)
        if(is.null(X))
            break

        k <- min(which(rev(x$height) < X$y), MAXCLUSTER)
        k <- max(k, 2)
        if(!is.null(oldx)){
            rect.hclust(x, k = oldk, x = oldx, cluster = cluster[, oldk-1],
                        border = "grey")
        }
        retval[[n]] <- unlist(rect.hclust(x, k = k, x = X$x,
                                          cluster = cluster[, k-1],
                                          border = "red"))
        if(!is.null(FUN)){
            if(!is.null(DEV.FUN)){
                dev.set(DEV.FUN)
            }
            retval[[n]] <- FUN(retval[[n]], ...)
        }

        oldx <- X$x
        oldk <- k
    }
    dev.set(DEV.x)
    invisible(retval)
}




kmeans <- function(x, centers, iter.max = 10)
{
    x <- as.matrix(x)
    m <- nrow(x)
    if(missing(centers))
	stop("centers must be a number or a matrix")
    if(length(centers) == 1) {
	k <- centers
	if(m < k)
	    stop("more cluster centers than data points.")
	centers <- x[sample(1:m, k), , drop=FALSE]
    } else {
	centers <- as.matrix(centers)
	k <- nrow(centers)
    }
    if(iter.max < 1) stop("iter.max must be positive")
    if(m < k)
	stop("more cluster centers than data points")
    if(ncol(x) != ncol(centers))
	stop("must have same number of columns in x and centers")
    Z <- .Fortran("kmns",
		  as.double(x),
		  as.integer(m),
		  as.integer(ncol(x)),
		  centers = as.double(centers),
		  as.integer(k),
		  c1 = integer(m),
		  integer(m),
		  nc =integer(k),
		  double(k),
		  double(k),
		  integer(k),
		  double(m),
		  integer(k),
		  integer(k),
		  as.integer(iter.max),
		  wss = double(k),
		  ifault = as.integer(0), PACKAGE="mva")
    switch(Z$ifault,
	   stop("empty cluster: try a better set of initial centers",
                call.=FALSE),
	   warning("did not converge in iter.max iterations", call.=FALSE),
	   stop("number of cluster centres must lie between 1 and nrow(x)",
                call.=FALSE)
	   )
    centers <- matrix(Z$centers, k)
    dimnames(centers) <- list(1:k, dimnames(x)[[2]])
    list(cluster = Z$c1, centers = centers, withinss = Z$wss, size = Z$nc)
}

prcomp <- function(x, retx = TRUE, center = TRUE, scale. = FALSE,
                   tol = NULL) {
    x <- as.matrix(x)
    x <- scale(x, center = center, scale = scale.)
#  as from 1.7.0 svd uses LAPACK: LINPACK was slow in this case.
#     if(dn[1] < dn[2]) {
#         s <- La.svd(x, nu = 0)
#         s$v <- t(s$vt)
#         s$vt <- NULL
#     } else
    s <- svd(x, nu = 0)
    if (!is.null(tol)) {
        rank <- sum(s$d > (s$d[1]*tol))
        if (rank < ncol(x))
            s$v <- s$v[, 1:rank, drop = FALSE]
    }
    s$d <- s$d / sqrt(max(1, nrow(x) - 1))
    dimnames(s$v) <-
        list(colnames(x), paste("PC", seq(len = ncol(s$v)), sep = ""))
    r <- list(sdev = s$d, rotation = s$v)
    if (retx) r$x <- x %*% s$v
    class(r) <- "prcomp"
    r
}

plot.prcomp <- function(x, main = deparse(substitute(x)), ...)
    screeplot(x, main = main, ...)

print.prcomp <- function(x, print.x = FALSE, ...) {
    cat("Standard deviations:\n")
    print(x$sdev, ...)
    cat("\nRotation:\n")
    print(x$rotation, ...)
    if (print.x && length(x$x)) {
        cat("\nRotated variables:\n")
        print(x$x, ...)
    }
    invisible(x)
}

summary.prcomp <- function(object, ...)
{
    vars <- object$sdev^2
    vars <- vars/sum(vars)
    importance <- rbind("Standard deviation" = object$sdev,
                        "Proportion of Variance" = round(vars, 5),
                        "Cumulative Proportion" = round(cumsum(vars), 5))
    colnames(importance) <- colnames(object$rotation)
    object$importance <- importance
    class(object) <- "summary.prcomp"
    object
}

print.summary.prcomp <-
function(x, digits = min(3, getOption("digits")-3), ...) {
    cat("Importance of components:\n")
    print(x$importance, digits = digits)
    invisible(x)
}
## copyright (C) 1998 W. N. Venables and B. D. Ripley
##
predict.princomp <- function(object, newdata, ...) {
    if (missing(newdata)) return(object$scores)
    scale(newdata, object$center, object$scale) %*% object$loadings
}

summary.princomp <- function(object, loadings = FALSE, cutoff = 0.1, ...)
{
    object$cutoff <- cutoff
    object$print.loadings <- loadings
    class(object) <- "summary.princomp"
    object
}

print.summary.princomp <-
    function(x, digits = 3, loadings = x$print.loadings, cutoff = x$cutoff,
             ...)
{
    vars <- x$sdev^2
    vars <- vars/sum(vars)
    cat("Importance of components:\n")
    print(rbind("Standard deviation" = x$sdev,
                "Proportion of Variance" = vars,
                "Cumulative Proportion" = cumsum(vars)))
    if(loadings) {
        cat("\nLoadings:\n")
        cx <- format(round(x$loadings, digits = digits))
        cx[abs(x$loadings) < cutoff] <-
            substring("       ", 1, nchar(cx[1,1]))
        print(cx, quote = FALSE, ...)
    }
    invisible(x)
}

plot.princomp <- function(x, main = deparse(substitute(x)), ...)
    screeplot(x, main = main, ...)

screeplot <-
function(x, npcs = min(10, length(x$sdev)),
         type = c("barplot", "lines"),
         main = deparse(substitute(x)), ...)
{
    main
    type <- match.arg(type)
    pcs <- x$sdev^2
    xp <- seq(length=npcs)
    if(type=="barplot")
        barplot(pcs[xp], names = names(pcs[xp]), main = main,
                ylab = "Variances", ...)
    else {
        plot(xp, pcs[xp], type = "b", axes = FALSE, main = main,
             xlab = "", ylab = "Variances", ...)
        axis(2)
        axis(1, at = xp, labels = names(pcs[xp]))
    }
    invisible()
}

loadings <- function(x) x$loadings
princomp <- function(x, ...) UseMethod("princomp")

## use formula to allow update() to be used.
princomp.formula <- function(formula, data = NULL, subset, na.action, ...)
{
    mt <- terms(formula, data = data)
    if(attr(mt, "response") > 0) stop("response not allowed in formula")
    attr(mt, "intercept") <- 0
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    mf$... <- NULL
    mf[[1]] <- as.name("model.frame")
    mf <- eval(mf, parent.frame())
    na.act <- attr(mf, "na.action")
    x <- model.matrix(mt, mf)
    res <- princomp.default(x, ...)
    ## fix up call to refer to the generic, but leave arg name as `formula'
    cl[[1]] <- as.name("princomp")
    res$call <- cl
    if(!is.null(na.act)) {
        res$na.action <- na.act
        if(!is.null(sc <- res$scores))
            res$scores <- napredict(na.act, sc)
    }
    res
}

princomp.default <-
    function(x, cor = FALSE, scores = TRUE, covmat = NULL,
             subset = rep(TRUE, nrow(as.matrix(x))), ...)
{
    cl <- match.call()
    cl[[1]] <- as.name("princomp")
    z <- if(!missing(x)) as.matrix(x)[subset, , drop = FALSE]
    if (is.list(covmat)) {
        if(any(is.na(match(c("cov", "n.obs"), names(covmat)))))
            stop("covmat is not a valid covariance list")
        cv <- covmat$cov
        n.obs <- covmat$n.obs
        cen <- covmat$center
    } else if(is.matrix(covmat)) {
        cv <- covmat
        n.obs <- NA
        cen <- NULL
    } else if(is.null(covmat)){
        dn <- dim(z)
        if(dn[1] < dn[2])
            stop("princomp can only be used with more units than variables")
        covmat <- cov.wt(z)             # returns list, cov() does not
        n.obs <- covmat$n.obs
        cv <- covmat$cov * (1 - 1/n.obs)# for S-PLUS compatibility
        cen <- covmat$center
    } else stop("covmat is of unknown type")
    if (cor) {
        sds <- sqrt(diag(cv))
        cv <- cv/(sds %o% sds)
    }
    edc <- La.eigen(cv, symmetric = TRUE)
    ev <- edc$values
    if (any(neg <- ev < 0)) { # S-PLUS sets all := 0
        ## 9 * : on Solaris found case where 5.59 was needed (MM)
        if (any(ev[neg] < - 9 * .Machine$double.eps * ev[1]))
            stop("covariance matrix is not non-negative definite")
        else
            ev[neg] <- 0
    }
    cn <- paste("Comp.", 1:ncol(cv), sep = "")
    names(ev) <- cn
    dimnames(edc$vectors) <- if(missing(x))
        list(dimnames(cv)[[2]], cn) else list(dimnames(x)[[2]], cn)
    sdev <- sqrt(ev)
    sc <- if (cor) sds else rep(1, ncol(cv))
    names(sc) <- colnames(cv)
    scr <- if (scores && !missing(x))
        scale(z, center = TRUE, scale = sc) %*% edc$vectors
    if (is.null(cen)) cen <- rep(NA, nrow(cv))
    edc <- list(sdev = sdev,
                loadings = structure(edc$vectors, class="loadings"),
                center = cen, scale = sc, n.obs = n.obs,
                scores = scr, call = cl)
    ## The Splus function also return list elements factor.sdev,
    ## correlations and coef, but these are not documented in the help.
    ## coef seems to equal load.  The Splus function also returns list
    ## element terms which is not supported here.
    class(edc) <- "princomp"
    edc
}

print.princomp <- function(x, ...)
{
    cat("Call:\n"); dput(x$call)
    cat("\nStandard deviations:\n")
    print(x$sdev, ...)
    cat("\n", length(x$scale), " variables and ", x$n.obs,
        "observations.\n")
    invisible(x)
}
.noGenerics <- TRUE

.onUnload <- function(libpath)
    library.dynam.unload("mva", libpath)
