R/mst.R

   1 ## mst.R (2006-11-08)
   2
   3 ##   Minimum Spanning Tree
   4
   5 ## Copyright 2002-2006 Yvonnick Noel, Julien Claude, and Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 mst <- function(X)
  11 {
  12     if (class(X) == "dist") X <- as.matrix(X)
  13     n <- dim(X)[1]
  14     N <- matrix(0, n, n)
  15     tree <- NULL
  16     large.value <- max(X) + 1
  17     diag(X) <- large.value
  18     index.i <- 1
  19
  20     for (i in 1:(n - 1)) {
  21         tree <- c(tree, index.i)
  22         m <- apply(as.matrix(X[, tree]), 2, min)  #calcul les minimum par colonne
  23         a <- sortIndex(X[, tree])[1, ]
  24         b <- sortIndex(m)[1]
  25         index.j <- tree[b]
  26         index.i <- a[b]
  27
  28         N[index.i, index.j] <- 1
  29         N[index.j, index.i] <- 1
  30
  31         for (j in tree) {
  32             X[index.i, j] <- large.value
  33             X[j, index.i] <- large.value
  34         }
  35     }
  36     dimnames(N) <- dimnames(X)
  37     class(N) <- "mst"
  38     return(N)
  39 }
  40
  41 ### Function returning an index matrix for an increasing sort
  42 sortIndex <- function(X)
  43 {
  44     if(length(X) == 1) return(1)                  # sorting a scalar?
  45     if(!is.matrix(X)) X <- as.matrix(X)           # force vector into matrix
  46     ## n <- nrow(X)
  47     apply(X, 2, function(v) order(rank(v)))       # find the permutation
  48 }
  49
  50 plot.mst <- function(x, graph = "circle", x1 = NULL, x2 = NULL, ...)
  51 {
  52     n <- nrow(x)
  53     if (is.null(x1) || is.null(x2)) {
  54         if (graph == "circle") {
  55             ang <- seq(0, 2 * pi, length = n + 1)
  56             x1 <- cos(ang)
  57             x2 <- sin(ang)
  58             plot(x1, x2, type = "n", xlab = "", ylab = "",
  59                  xaxt = "n", yaxt = "n", bty = "n", ...)
  60         }
  61         if (graph == "nsca") {
  62             XY <- nsca(x)
  63             x1 <- XY[, 1]
  64             x2 <- XY[, 2]
  65             plot(XY, type = "n", xlab = "\"nsca\" -- axis 1",
  66                  ylab = "\"nsca\" -- axis 2", ...)
  67         }
  68     } else plot(x1, x2, type = "n", xlab = deparse(substitute(x1)),
  69                 ylab = deparse(substitute(x2)), ...)
  70     for (i in 1:n) {
  71         w1 <- which(x[i, ] == 1)
  72         segments(x1[i], x2[i], x1[w1], x2[w1])
  73     }
  74     points(x1, x2, pch = 21, col = "black", bg = "white", cex = 3)
  75     text(x1, x2, 1:n, cex = 0.8)
  76 }
  77
  78 nsca <- function(A)
  79 {
  80     Dr <- apply(A, 1, sum)
  81     Dc <- apply(A, 2, sum)
  82
  83     eig.res <- eigen(diag(1 / sqrt(Dr)) %*% A %*% diag(1 / sqrt(Dc)))
  84     r <- diag(1 / Dr) %*% (eig.res$vectors)[, 2:4]
  85     ## The next line has been changed by EP (20-02-2003), since
  86     ## it does not work if 'r' has no dimnames already defined
  87     ## dimnames(r)[[1]] <- dimnames(A)[[1]]
  88     rownames(r) <- rownames(A)
  89     r
  90 }