The R Project SVN R

#  File src/library/stats/R/dendrogram.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 1995-2017 The R Core Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

as.dendrogram <- function(object, ...) UseMethod("as.dendrogram")

as.dendrogram.dendrogram <- function(object, ...) object

as.dendrogram.hclust <- function (object, hang = -1, check = TRUE, ...)
## hang = 0.1  is default for plot.hclust
{
    nolabels <- is.null(object$labels)
    merge <- object$merge
    if(check && !isTRUE(msg <- .validity.hclust(object, merge, order=nolabels)))
    stop(msg)
    if(nolabels)
    object$labels <- seq_along(object$order)
    z <- list()
    nMerge <- length(oHgt <- object$height)
    hMax <- oHgt[nMerge]
    for (k in 1L:nMerge) {
    x <- merge[k, ]# no sort() anymore!
    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(-x)
        attr(zk, "members") <- 2L
        attr(zk, "midpoint") <- 0.5 # mean( c(0,1) )
        objlabels <- object$labels[-x]
        attr(zk[[1L]], "label") <- objlabels[1L]
        attr(zk[[2L]], "label") <- objlabels[2L]
        attr(zk[[1L]], "members") <- attr(zk[[2L]], "members") <- 1L
        attr(zk[[1L]], "height")  <- attr(zk[[2L]], "height") <- h0
        attr(zk[[1L]], "leaf")    <- attr(zk[[2L]], "leaf") <- TRUE
    }
    else if (any(neg)) {        # one leaf, one node
        X <- as.character(x)
        ## Originally had "x <- sort(..) above => leaf always left, x[1L];
        ## don't want to assume this
        isL <- x[1L] < 0 ## is leaf left?
        zk <-
        if(isL) list(-x[1L], z[[X[2L]]])
        else    list(z[[X[1L]]], -x[2L])
        attr(zk, "members") <- attr(z[[X[1 + isL]]], "members") + 1L
        attr(zk, "midpoint") <-
                (.memberDend(zk[[1L]]) + attr(z[[X[1 + isL]]], "midpoint"))/2
        attr(zk[[2 - isL]], "members") <- 1L
        attr(zk[[2 - isL]], "height") <- h0
        attr(zk[[2 - isL]], "label") <- object$labels[-x[2 - isL]]
        attr(zk[[2 - isL]], "leaf") <- TRUE
        z[[X[1 + isL]]] <- NULL
    }
    else {              # two non-leaf nodes
        x <- as.character(x)
            ## "merge" the two ('earlier') branches:
        zk <- list(z[[x[1L]]], z[[x[2L]]])
        attr(zk, "members") <- attr(z[[x[1L]]], "members") +
        attr(z[[x[2L]]], "members")
        attr(zk, "midpoint") <- (attr(z[[x[1L]]], "members") +
                     attr(z[[x[1L]]], "midpoint") +
                     attr(z[[x[2L]]], "midpoint"))/2
        z[[x[1L]]] <- z[[x[2L]]] <- NULL
    }
    attr(zk, "height") <- oHgt[k]
    z[[as.character(k)]] <- zk
    }
    structure(z[[as.character(k)]], class = "dendrogram")
}

## Count the number of leaves in a dendrogram.
nleaves <- function (node) {
    if (is.leaf(node))
    return(1L)

    todo <- NULL # Non-leaf nodes to traverse after this one.
    count <- 0L
    repeat {
    ## For each child: count iff a leaf, add to todo list otherwise.
    while (length(node)) {
        child <- node[[1L]]
        node <- node[-1L]
        if (is.leaf(child)) {
        count <- count + 1L
        } else {
        todo <- list(node=child, todo=todo)
        }
    }
    ## Advance to next node, terminating when no nodes left to count.
    if (is.null(todo)) {
        break
    } else {
        node <- todo$node
        todo <- todo$todo
    }
    }
    count
}

## Reversing   as.dendrogram.hclust() above (as much as possible)
## is only possible for dendrograms with *binary* splits
as.hclust.dendrogram <- function(x, ...)
{
    stopifnot(is.list(x), length(x) == 2L)
    n <- nleaves(x)
    stopifnot(n == attr(x, "members"))

    ## Ord and labels for each leaf node (in preorder).
    ord <- integer(n)
    labsu <- character(n)

    ## Height and (parent,index) for each internal node (in preorder).
    n.h <- n - 1L
    height <- numeric(n.h)
    myIdx <- matrix(NA_integer_, 2L, n.h)

    ## Record merges initially in preorder traversal
    ## We will resort into merge order at end.
    merge <- matrix(NA_integer_, 2L, n.h)

    ## Starting at root, traverse dendrogram recording
    ## information above about leaves and nodes encountered
    position <- 0L  # position within current node
    stack <- NULL   # parents of current node plus saved state
    leafCount <- 0L # number of leaves seen
    nodeCount <- 0L # number of nodes seen
    repeat {
        ## Pre-order traversal of the current node.
        ## Will descend into non-leaf children pushing parents onto stack.
    while (length(x)) {
            ## Record height and index list on first visit to each internal node.
        if (position == 0L) {
        nodeCount <- nodeCount + 1L
                myNodeIndex <- nodeCount
                if (nodeCount != 1L) {
                    myIdx[,nodeCount] <- c(stack$position, stack$myNodeIndex)
                }
        height[nodeCount] <- attr(x, "height")
        }
        position <- position + 1L
        child <- x[[1L]]
        x <- x[-1L]
        if (is.leaf(child)) {
                ## Record information about leaf nodes.
                leafCount <- leafCount + 1L
                labsu[leafCount] <- attr(child,'label')
                ord[leafCount] <- as.integer(child)
                merge[position,myNodeIndex] <- - ord[leafCount]
            } else {
                stopifnot (length(child)==2L)
                ## Descend into non-leaf nodes, saving state on stack.
        stack <- list(node=x, position=position,
                              myNodeIndex=myNodeIndex, stack=stack)
        x <- child
        position <- 0L
        }
    }
        ## All children of current node have been traversed.

        ## Terminate if current node was the root node.
    if (is.null(stack)) {
        break
    }

        ## Otherwise, pop parent node and state.
    position <- stack$position   # Restore position in parent node.
    x <- stack$node
        myNodeIndex <- stack$myNodeIndex
    stack <- stack$stack
    }

    iOrd <- sort.list(ord)
    if(!identical(ord[iOrd], seq_len(n)))
    stop(gettextf(
        "dendrogram entries must be 1,2,..,%d (in any order), to be coercible to \"hclust\"",
        n), domain=NA)

    ## ties: break ties "compatibly" with above preorder traversal -- relies on stable sort here:
    ii <- sort.list(height, decreasing=TRUE)[n.h:1L]
    stopifnot(ii[n.h] == 1L)

    ## Record internal merges
    k <- seq_len(n.h-1L)
    merge[t(myIdx[,ii[k]])] <- + k

    if (getOption("as.hclust.dendr", FALSE)) { # be verbose
    for(k in seq_len(n.h)) {
        cat(sprintf("ii[k=%2d]=%2d ", k, ii[k]))
        cat("-> s=merge[[,ii[k]]]=")
        str(merge[,ii[k]])
    }
    }

    structure(list(merge = t(merge[,ii]),  # Resort into merge order
           height = height[ii], # Resort into merge order
           order = ord,
           labels = labsu[iOrd],
           call = match.call(),
           method = NA_character_,
           dist.method = NA_character_),
          class = "hclust")
}

### 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 <- 1L
    }
    r
}

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

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

    type <- match.arg(type) ## currently only "hclust"
    stopifnot( inherits(x, "dendrogram") )
    verbose <- getOption("verbose", 0) >= 2 # non-public
    setmid <- function(d, type) {
    depth <- 0L
    kk <- integer()
    jj <- integer()
    dd <- list()
    repeat {
        if(!is.leaf(d)) {# no "midpoint" for leaf
        k <- length(d)
        if(k < 1)
            stop("dendrogram node with non-positive #{branches}")
        depth <- depth + 1L
        if(verbose) cat(sprintf(" depth(+)=%4d, k=%d\n", depth, k))
        kk[depth] <- k
        if(storage.mode(jj) != storage.mode(kk)) # (long vectors)
            storage.mode(jj) <- storage.mode(kk)
        dd[[depth]] <- d
        d <- d[[jj[depth] <- 1L]]
        next
        }
        while(depth) {
        k <- kk[depth]
        j <- jj[depth]
        r <- dd[[depth]] # incl. attributes!
        r[[j]] <- unclass(d)
        if(j < k) break
        depth <- depth - 1L
        if(verbose) cat(sprintf(" depth(-)=%4d, k=%d\n", depth, k))
        midS <- sum(vapply(r, .midDend, 0))
        if(!quiet && type == "hclust" && k != 2)
            warning("midcache() of non-binary dendrograms only partly implemented")
        ## compatible to as.dendrogram.hclust() {MM: doubtful if k > 2}
        attr(r, "midpoint") <- (.memberDend(r[[1L]]) + midS) / 2
        d <- r
        }
        if(!depth) break
        dd[[depth]] <- r
        d <- r[[jj[depth] <- j + 1L]]
    }
    d
    }
    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(is.leaf(x))
    cat("leaf '", format(attr(x, "label"), digits = digits),"'", sep = "")
    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 = NA, digits.d = 3L, give.attr = FALSE,
          wid = getOption("width"), nest.lev = 0L, indent.str = "",
          last.str = getOption("str.dendrogram.last"), stem = "--", ...)
{
## TO DO: when object is part of a larger structure which is str()ed
##    with default max.level= NA, it should not be str()ed to all levels,
##   but only to e.g. level 2
## Implement via smarter default for 'max.level' (?)

    pasteLis <- function(lis, dropNam, sep = " = ") {
    ## drop uninteresting "attributes" here
    lis <- lis[!(names(lis) %in% dropNam)]
    fl <- sapply(lis, format, digits = digits.d)
    paste(paste(names(fl), fl, sep = sep), collapse = ", ")
    }

    todo <- NULL # Nodes to process after this one
    repeat {
        ## when  indent.str  ends in a blank, i.e. "last" (see below)
    istr <- sub(" $", last.str, indent.str)
    cat(istr, stem, sep = "")

    at <- attributes(object)
    memb <- at[["members"]]
    hgt  <- at[["height"]]
    if(!is.leaf(object)) {
        le <- length(object)
        if(give.attr) {
        if(nzchar(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, "]",
        if(!is.na(max.level) && nest.lev == max.level)" ..", "\n", sep = "")
        if (is.na(max.level) || nest.lev < max.level) {
        ## Push children onto todo list in reverse order.
        ## Assumes at least one child.
        nest.lev <- nest.lev + 1L
        todo <- list(object=object[[le]], nest.lev = nest.lev,
                 indent.str = paste(indent.str, "  "), todo = todo)
        indent.str <- paste (indent.str, " |")
        while ((le <- le - 1L) > 0L) {
            todo <- list(object=object[[le]], nest.lev = nest.lev,
                 indent.str = indent.str, todo = todo)
        }
        }
    } 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 || nzchar(at)) cat(if(!any.at)"(", at, ")")
        cat("\n")
    }
        ## Advance to next node, if any.
    if (is.null(todo)) {
        break
    } else {
        object <- todo$object
        nest.lev <- todo$nest.lev
        indent.str <- todo$indent.str
        todo <- todo$todo
    }
    }
    invisible()
}


## The ``generic'' method for "[["  (analogous to e.g., "[[.POSIXct"):
## --> subbranches (including leafs!) are dendrograms as well!
`[[.dendrogram` <- function(x, ..., drop = TRUE) {
    if(!is.null(r <- NextMethod("[[")))
        structure(r, class = "dendrogram")
}

nobs.dendrogram <- function(object, ...) attr(object, "members")

## FIXME: need larger par("mar")[1L] or [4L] for longish labels !
## {probably don't change, just print a warning ..}
plot.dendrogram <-
    function (x, type = c("rectangle", "triangle"), center = FALSE,
          edge.root = is.leaf(x) || !is.null(attr(x, "edgetext")),
          nodePar = NULL, edgePar = list(),
          leaflab = c("perpendicular", "textlike", "none"), dLeaf = NULL,
          xlab = "", ylab = "", xaxt="n", yaxt="s",
          horiz = FALSE, frame.plot = FALSE, xlim, ylim, ...)
{
    type <- match.arg(type)
    leaflab <- match.arg(leaflab)
    hgt <- attr(x, "height")
    if (edge.root && is.logical(edge.root))
    edge.root <- 0.0625 * if(is.leaf(x)) 1 else 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 }
    xl. <- c(x1 - 1/2, x2 + 1/2)
    yl. <- c(0, yTop)
    if (horiz) {## swap and reverse direction on `x':
    tmp <- xl.; xl. <- rev(yl.); yl. <- tmp
    tmp <- xaxt; xaxt <- yaxt; yaxt <- tmp
    }
    if(missing(xlim) || is.null(xlim)) xlim <- xl.
    if(missing(ylim) || is.null(ylim)) ylim <- yl.
    dev.hold(); on.exit(dev.flush())
    plot(0, xlim = xlim, ylim = ylim, type = "n", xlab = xlab, ylab = ylab,
     xaxt = xaxt, yaxt = yaxt, frame.plot = frame.plot, ...)
    if(is.null(dLeaf))
        dLeaf <- .75*(if(horiz) strwidth("w") else strheight("x"))

    if (edge.root) {
### FIXME: the first edge + edgetext is drawn here, all others in plotNode()
### -----  maybe use trick with adding a single parent node to the top ?
    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,
             dLeaf = dLeaf, 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, dLeaf,
         nodePar, edgePar, horiz = FALSE)
{
  wholetree <- subtree
  depth <- 0L
  llimit <- list()
  KK <- integer()
  kk <- integer()

  repeat {
    inner <- !is.leaf(subtree) && x1 != x2
    yTop <- attr(subtree, "height")
    bx <- plotNodeLimit(x1, x2, subtree, center)
    xTop <- bx$x
    depth <- depth + 1L
    llimit[[depth]] <- bx$limit

    ## 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, width = 4), ",", formatC(x2, width = 4), ")",
        "--> xTop=", formatC(xTop, width = 8), "\n", sep = "")
    }

    Xtract <- function(nam, L, default, indx)
    rep(if(nam %in% names(L)) L[[nam]] else default,
        length.out = indx)[indx]
    asTxt <- function(x) # to allow 'plotmath' labels:
    if(is.character(x) || is.expression(x) || is.null(x)) x else as.character(x)

    i <- if(inner || hasP) 1 else 2 # only 1 node specific par

    if(!is.null(nPar)) { ## draw this node
    pch <- Xtract("pch", nPar, default = 1L: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(leaflab == "textlike")
        p.col <- Xtract("p.col", nPar, default = "white", i)
    lab.col <- Xtract("lab.col", nPar, default = par("col"), i)
    lab.cex <- Xtract("lab.cex", nPar, default = c(1,1), i)
    lab.font <- Xtract("lab.font", nPar, default = par("font"), i)
    lab.xpd <- Xtract("xpd", nPar, default = c(TRUE, TRUE), i)
    if (is.leaf(subtree)) {
    ## label leaf
    if (leaflab == "perpendicular") { # somewhat like plot.hclust
        if(horiz) {
                X <- yTop + dLeaf * lab.cex
                Y <- xTop; srt <- 0; adj <- c(0, 0.5)
        }
        else {
                Y <- yTop - dLeaf * lab.cex
                X <- xTop; srt <- 90; adj <- 1
        }
            nodeText <- asTxt(attr(subtree,"label"))
        text(X, Y, nodeText, xpd = lab.xpd, srt = srt, adj = adj,
                 cex = lab.cex, col = lab.col, font = lab.font)
    }
    }
    else if (inner) {
    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)
    }
    for (k in seq_along(subtree)) {
        child <- subtree[[k]]
        ## draw lines to the children and draw them 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 (!is.leaf(child) || hasE) 1 else 2
        ## define line attributes for segmentsHV():
        col <- Xtract("col", ePar, default = par("col"), i)
        lty <- Xtract("lty", ePar, default = par("lty"), i)
        lwd <- Xtract("lwd", ePar, default = par("lwd"), i)
        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 (is.leaf(child) && leaflab == "textlike") {
        nodeText <- asTxt(attr(child,"label"))
        if(getOption("verbose"))
            cat('-- with "label"',format(nodeText))
        hln <- 0.6 * strwidth(nodeText, cex = lab.cex)/2
        vln <- 1.5 * strheight(nodeText, cex = lab.cex)/2
        rect(xBot - hln, yBot,
             xBot + hln, yBot + 2 * vln, col = p.col)
        text(xBot, yBot + vln, nodeText, xpd = lab.xpd,
                     cex = lab.cex, col = lab.col, font = lab.font)
        }
        if (!is.null(attr(child, "edgetext"))) {
        edgeText <- asTxt(attr(child, "edgetext"))
        if(getOption("verbose"))
            cat('-- with "edgetext"',format(edgeText))
        if (!is.null(vln)) {
            mx <-
            if(type == "triangle")
                (xTop+ xBot+ ((xTop - xBot)/(yTop - yBot)) * vln)/2
            else xBot
            my <- (yTop + yBot + 2 * vln)/2
        }
        else {
            mx <- if(type == "triangle") (xTop + xBot)/2 else xBot
            my <- (yTop + yBot)/2
        }
        ## Both for "triangle" and "rectangle" : Diamond + Text

                p.col <- Xtract("p.col", ePar, default = "white", i)
                p.border <- Xtract("p.border", ePar, default = par("fg"), i)
                ## edge label pars: defaults from the segments pars
                p.lwd <- Xtract("p.lwd", ePar, default = lwd, i)
                p.lty <- Xtract("p.lty", ePar, default = lty, i)
                t.col <- Xtract("t.col", ePar, default = col, i)
                t.cex <- Xtract("t.cex", ePar, default =  1,  i)
                t.font <- Xtract("t.font", ePar, default = par("font"), i)

        vlm <- strheight(c(edgeText,"h"), cex = t.cex)/2
        hlm <- strwidth (c(edgeText,"m"), cex = t.cex)/2
        hl3 <- c(hlm[1L], hlm[1L] + hlm[2L], hlm[1L])
                if(horiz) {
                    polygon(my+ c(-hl3, hl3), mx + sum(vlm)*c(-1L:1L, 1L:-1L),
                            col = p.col, border = p.border,
                            lty = p.lty, lwd = p.lwd)
                    text(my, mx, edgeText, cex = t.cex, col = t.col,
                         font = t.font)
                } else {
                    polygon(mx+ c(-hl3, hl3), my + sum(vlm)*c(-1L:1L, 1L:-1L),
                            col = p.col, border = p.border,
                            lty = p.lty, lwd = p.lwd)
                    text(mx, my, edgeText, cex = t.cex, col = t.col,
                         font = t.font)
                }
        }
    }
    }

    if (inner && length(subtree)) {
    KK[depth] <- length(subtree)
    if (storage.mode(kk) != storage.mode(KK))
        storage.mode(kk) <- storage.mode(KK)

    ## go to first child
    kk[depth] <- 1L
    x1 <- bx$limit[1L]
    x2 <- bx$limit[2L]
    subtree <- subtree[[1L]]
    }
    else {
    repeat {
        depth <- depth - 1L
        if (!depth || kk[depth] < KK[depth]) break
    }
    if (!depth) break
    length(kk) <- depth
    kk[depth] <- k <- kk[depth] + 1L
    x1 <- llimit[[depth]][k]
    x2 <- llimit[[depth]][k + 1L]
    subtree <- wholetree[[kk]]
    }
  } ## repeat
  invisible()
}

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 <- !is.leaf(subtree) && x1 != x2
    limit <- c(x1,
           if(inner) {
           K <- length(subtree)
           mTop <- .memberDend(subtree)
           limit <- integer(K)
           xx1 <- x1
           for(k in 1L: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
           } else ## leaf
           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(!is.leaf(subtree)) {
        if(!(K <- length(subtree)))
        stop("non-leaf subtree of length 0")
        new.mem <- 0L
        for(k in 1L:K) {
                sub <- subtree[[k]]
        if(attr(sub, "height") <= h) {
            ## cut it, i.e. save to LOWER[] and make a leaf
            at <- attributes(sub)
            at$leaf <- TRUE
                    at$class <- NULL# drop it from leaf
            at$x.member <- at$members
            new.mem <- new.mem + (at$members <- 1L)
            at$label <- paste("Branch", X)
            subtree[[k]] <- X #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(sub, 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()

is.leaf <- 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")
    if(is.list(x))
    unlist(x)
    else ## leaf
    as.vector(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( is.leaf(x) ) return(x[1L])
#      return(c(ord(x[[1L]]), ord(x[[2L]])))
#    }
#   return(ord(x))
# }

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

reorder.dendrogram <- function(x, wts, agglo.FUN = sum, ...)
{
    if( !inherits(x, "dendrogram") )
    stop("'reorder.dendrogram' requires a dendrogram")
    agglo.FUN <- match.fun(agglo.FUN)
    oV <- function(x, wts) {
    depth <- 0L
    kk <- jj <- integer()
    xx <- list()
    repeat {
        if(is.leaf(x))
        attr(x, "value") <- wts[x[1L]]
        else {
        k <- length(x)
        if(k == 0L) stop("invalid (length 0) node in dendrogram")
        depth <- depth + 1L
        kk[depth] <- k
        if(storage.mode(jj) != storage.mode(kk))
            storage.mode(jj) <- storage.mode(kk)
        ## insert/compute 'wts' recursively down the branches:
        xx[[depth]] <- x
        x <- x[[jj[depth] <- 1L]]
        next
        }
        while(depth) {
        b <- x
        x <- xx[[depth]]
        j <- jj[depth]
        x[[j]] <- b
        if(j < kk[depth]) break
        depth <- depth - 1L
        vals <- vapply(x, attr, numeric(1L), which="value")
        iOrd <- sort.list(vals)
        attr(x, "value") <- agglo.FUN(vals[iOrd])
        x[] <- x[iOrd]
        }
        if(!depth) break
        xx[[depth]] <- x
        x <- x[[jj[depth] <- j + 1L]]
    }
        x
    }
    midcache.dendrogram( oV(x, wts) )
}

rev.dendrogram <- function(x) {
    if(is.leaf(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 1L:k) ## recurse
    r[[j]] <- rev(x[[k+1-j]])
    midcache.dendrogram( r )
}

labels.dendrogram <- function(object, ...) {
    if(is.list(object))
        rapply(object, function(n) attr(n,"label"))
    else # can "end" in a leaf here
        attr(object, "label")
}

merge.dendrogram <- function(x, y, ..., height,
                             adjust = c("auto", "add.max", "none"))
{
    stopifnot(inherits(x,"dendrogram"), inherits(y,"dendrogram"))
    if((adjust <- match.arg(adjust)) == "auto")
        adjust <-
            ## dendrograms as from hclust(), have entries {1,2,..,n}; "cheap" check:
            if(min(unlist(x)) == 1 && min(unlist(y)) == 1)
                "add.max"
            else # for now, can imagine more:
                "none"
    if(adjust == "add.max") {
        add.ifleaf <- function(i, add) if(is.leaf(i)) i + add else i
        add <- max(unlist(x))
        y <- dendrapply(y, add.ifleaf, add=add)
    }
    r <- list(x,y)
    if(length(xtr <- list(...))) {
    if(!all(is.d <- vapply(xtr, inherits, NA, what="dendrogram"))) {
        xpr <- substitute(c(...))
        nms <- sapply(xpr[-1][!is.d], deparse, nlines = 1L)
            ## do not simplify: xgettext needs this form
            msg <- ngettext(length(nms),
                            "extra argument %s is not of class \"%s\"",
                            "extra arguments %s are not of class \"%s\"s")
        stop(sprintf(msg, paste(nms, collapse=", "), "dendrogram"),
                 domain = NA)
    }
    if(adjust == "add.max") {
        add <- max(add, unlist(y))
        for(i in seq_along(xtr)) {
        if(i > 1L) add <- max(add, unlist(xtr[i-1L]))
        xtr[[i]] <- dendrapply(xtr[[i]], add.ifleaf, add=add)
        }
    }
    r <- c(r, xtr)
    }
    attr(r, "members") <- sum(vapply(r, attr, 0L, which="members"))
    h.max <- max(vapply(r, attr, 0., which="height"))
    if(missing(height) || is.null(height))
    height <- 1.1 * h.max
    else if(height < h.max) {
        msg <- gettextf("'height' must be at least %g, the maximal height of its components", h.max)
        stop(msg, domain = NA)
    }
    attr(r, "height") <- height
    class(r) <- "dendrogram"
    midcache.dendrogram(r, quiet=TRUE)
}

dendrapply <- function(X, FUN, ...)
{
    ## Purpose: "dendrogram" recursive apply {to each node}
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 26 Jun 2004, 22:43
    FUN <- match.fun(FUN)
    if( !inherits(X, "dendrogram") ) stop("'X' is not a dendrogram")

    ## Node apply recursively:
    Napply <- function(d) {
    r <- FUN(d, ...)
    if(!is.leaf(d)) {
        if(!is.list(r)) r <- as.list(r) # fixing unsafe FUN()s
        if(length(r) < (n <- length(d))) r[seq_len(n)] <- vector("list", n)
        ## and overwrite recursively, possibly keeping "attr"
        r[] <- lapply(d, Napply)
        }
    r
    }
    Napply(X)
}


## original Andy Liaw; modified RG, MM :
heatmap <-
function (x, Rowv=NULL, Colv=if(symm)"Rowv" else NULL,
      distfun = dist, hclustfun = hclust,
          reorderfun = function(d,w) reorder(d,w),
          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,
      keep.dendro = FALSE,
      verbose = getOption("verbose"), ...)
{
    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[1L]
    nc <- di[2L]
    if(nr <= 1 || nc <= 1)
    stop("'x' must have at least 2 rows and 2 columns")
    if(!is.numeric(margins) || length(margins) != 2L)
    stop("'margins' must be a numeric vector of length 2")

    doRdend <- !identical(Rowv,NA)
    doCdend <- !identical(Colv,NA)
    if(!doRdend && identical(Colv, "Rowv")) doCdend <- FALSE
    ## 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(doRdend) {
    if(inherits(Rowv, "dendrogram"))
        ddr <- Rowv
    else {
        hcr <- hclustfun(distfun(x))
        ddr <- as.dendrogram(hcr)
        if(!is.logical(Rowv) || Rowv)
        ddr <- reorderfun(ddr, Rowv)
    }
    if(nr != length(rowInd <- order.dendrogram(ddr)))
        stop("row dendrogram ordering gave index of wrong length")
    }
    else rowInd <- 1L:nr

    if(doCdend) {
    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 <- reorderfun(ddc, Colv)
    }
    if(nc != length(colInd <- order.dendrogram(ddc)))
        stop("column dendrogram ordering gave index of wrong length")
    }
    else colInd <- 1L:nc

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

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

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

    ## Calculate the plot layout
    lmat <- rbind(c(NA, 3), 2:1)
    lwid <- c(if(doRdend) 1 else 0.05, 4)
    lhei <- c((if(doCdend) 1 else 0.05) + if(!is.null(main)) 0.2 else 0, 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[1L], 0.2, lhei[2L])
    }
    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[1L], 0.2, lwid[2L])
    }
    lmat[is.na(lmat)] <- 0
    if(verbose) {
    cat("layout: widths = ", lwid, ", heights = ", lhei,"; lmat=\n")
    print(lmat)
    }

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

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

    image(1L:nc, 1L:nr, x, xlim = 0.5+ c(0, nc), ylim = 0.5+ c(0, nr),
      axes = FALSE, xlab = "", ylab = "", ...)
    axis(1, 1L:nc, labels = labCol, las = 2, line = -0.5, tick = 0,
         cex.axis = cexCol)
    if(!is.null(xlab)) mtext(xlab, side = 1, line = margins[1L] - 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[2L] - 1.25)

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

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

    par(mar = c(0, 0, if(!is.null(main)) 1 else 0, margins[2L]))
    if(doCdend)
    plot(ddc,       axes = FALSE, xaxs = "i", leaflab = "none")
    else if(!is.null(main)) frame()

    ## title
    if(!is.null(main)) {
        par(xpd = NA)# {we have room on the left}
        title(main, cex.main = 1.5*op[["cex.main"]])
    }

    invisible(list(rowInd = rowInd, colInd = colInd,
           Rowv = if(keep.dendro && doRdend) ddr,
           Colv = if(keep.dendro && doCdend) ddc ))
}