R/SDM.R

   1 ## SDM.R (2012-04-02)
   2
   3 ## Construction of Consensus Distance Matrix With SDM
   4
   5 ## Copyright 2011-2012 Andrei-Alin Popescu
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 SDM <- function(...)
  11 {
  12     st <- list(...) # first half contains matrices, second half s_p
  13     k <- length(st)/2
  14     ONEtoK <- seq_len(k)
  15
  16     ## make sure we have only matrices:
  17     for (i in ONEtoK) st[[i]] <- as.matrix(st[[i]])
  18
  19     ## store the rownames of each matrix in a list because they are often called:
  20     ROWNAMES <- lapply(st[ONEtoK], rownames)
  21
  22     ## the number of rows of each matrix:
  23     NROWS <- sapply(ROWNAMES, length)
  24     tot <- sum(NROWS)
  25
  26     labels <- unique(unlist(ROWNAMES))
  27     sp <- unlist(st[k + ONEtoK])
  28
  29     astart <- numeric(tot) # start of aip, astart[p] is start of aip
  30     astart[1] <- k
  31     for (i in 2:k)
  32         astart[i] <- astart[i - 1] + NROWS[i - 1]
  33
  34     ## apparently erased by the operation below so no need to initialize:
  35     ## a <- mat.or.vec(1, k + tot + k + length(labels))
  36
  37     ## first k are alphas, subsequent ones aip
  38     ## each matrix p starting at astart[p], next are
  39     ## Lagrange multipliers, miu, niu, lambda in that order
  40     n <- length(labels)
  41     miustart <- k + tot
  42     niustart <- miustart + n
  43     lambstart <- niustart + k - 1
  44
  45     X <- matrix(0, n, n, dimnames = list(labels, labels))
  46     V <- w <- X
  47
  48     tmp <- 2 * k + tot + n
  49     col <- numeric(tmp) # free terms of system
  50
  51     for (i in 1:(n - 1)) {
  52         for (j in (i + 1):n) {
  53             for (p in ONEtoK) {
  54                 ## d <- st[[p]] # not needed anymore
  55                 if (is.element(labels[i], ROWNAMES[[p]]) && is.element(labels[j], ROWNAMES[[p]])) {
  56                     w[i, j] <- w[j, i] <- w[i, j] + sp[p]
  57                 }
  58             }
  59         }
  60     }
  61
  62     ONEtoN <- seq_len(n)
  63
  64     Q <- matrix(0, tmp, tmp)
  65     ## first decompose first sum in paper
  66     for (p in ONEtoK) {
  67         d_p <- st[[p]]
  68         for (l in ONEtoK) { # first compute coefficients of alphas
  69             d <- st[[l]]
  70             sum <- 0
  71             dijp <- -1
  72             if (l == p) { # calculate alpha_p
  73                 for (i in ONEtoN) {
  74                     for (j in ONEtoN) { # check if {i,j}\subset L_l
  75                         if (i == j) next # make sure i != j
  76                         ## d <- st[[l]] # <- moved-up
  77                         pos <- match(labels[c(i, j)], ROWNAMES[[l]]) # <- returns NA if not in this matrix
  78                         if (all(!is.na(pos))) {
  79                             ipos <- pos[1L]
  80                             jpos <- pos[2L]
  81                             dij <- d[ipos, jpos]
  82                             sum <- sum + dij * dij - sp[l] * dij * dij / w[i,j]
  83                             tmp2 <- dij - sp[l] * dij / w[i,j]
  84                             Q[p, astart[l] + ipos] <- Q[p, astart[l] + ipos] + tmp2
  85                             Q[p, astart[l] + jpos] <- Q[p, astart[l] + jpos] + tmp2
  86                         }
  87                     }
  88                 }
  89             } else {
  90                 for (i in ONEtoN) {
  91                     for (j in ONEtoN) { # check if {i,j}\subset L_l
  92                         if (i == j) next
  93                         ## d <- st[[l]] # <- moved-up
  94                         pos <- match(labels[c(i, j)], ROWNAMES[[l]])
  95                         posp <- match(labels[c(i, j)], ROWNAMES[[p]])
  96                         if (all(!is.na(pos)) && all(!is.na(posp))) {
  97                             ipos <- pos[1L]
  98                             jpos <- pos[2L]
  99                             dij <- d[ipos, jpos]
 100                             dijp <- d_p[posp[1L], posp[2L]]
 101                             sum <- sum - sp[l] * dij * dijp / w[i, j]
 102                             tmp2 <- sp[l] * dijp / w[i, j]
 103                             Q[p,astart[l] + ipos] <- Q[p, astart[l] + ipos] - tmp2
 104                             Q[p,astart[l] + jpos] <- Q[p, astart[l] + jpos] - tmp2
 105                         }
 106                     }
 107                 }
 108             }
 109             Q[p, l] <- sum
 110         }
 111         Q[p, lambstart + 1] <- 1
 112     }
 113
 114     r <- k
 115
 116     for (p in ONEtoK) {
 117         dp <- st[[p]]
 118         for (i in ONEtoN) {
 119             if (is.element(labels[i], ROWNAMES[[p]])) {
 120                 r <- r + 1
 121                 for (l in ONEtoK) {
 122                     d <- st[[l]]
 123                     if (l == p) {
 124                         ipos <- match(labels[i], ROWNAMES[[p]])
 125                         for (j in ONEtoN) {
 126                             if (i == j) next
 127                             jpos <- match(labels[j], ROWNAMES[[p]])
 128                             if (!is.na(jpos)) {
 129                                 dij <- d[ipos, jpos]
 130                                 Q[r, l] <- Q[r, l] + dij - sp[l] * dij / w[i, j]
 131                                 tmp2 <- 1 - sp[l] / w[i, j]
 132                                 Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] + tmp2
 133                                 Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] + tmp2
 134                             }
 135                         }
 136                     } else {
 137                         for (j in ONEtoN) {
 138                             if (i == j) next
 139                             if (!is.element(labels[j], rownames(dp))) next
 140                             pos <- match(labels[c(i, j)], ROWNAMES[[l]])
 141                             if (all(!is.na(pos))) {
 142                                 ipos <- pos[1L]
 143                                 jpos <- pos[2L]
 144                                 dij <- d[ipos, jpos]
 145                                 Q[r, l] <- Q[r, l] - sp[l] * dij / w[i, j]
 146                                 tmp2 <- sp[l]/w[i, j]
 147                                 Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] - tmp2
 148                                 Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] - tmp2
 149                             }
 150                         }
 151                     }
 152                 }
 153                 if (p < k) Q[r, ] <- Q[r, ] * sp[p]
 154                 Q[r, miustart + i] <- 1
 155                 if (p < k) Q[r, niustart + p] <- 1
 156             }
 157         }
 158     }
 159
 160     r <- r + 1
 161     col[r] <- k
 162     Q[r, ONEtoK] <- 1
 163     ## for (i in 1:k) Q[r, i] <- 1
 164
 165     for (i in ONEtoN) {
 166         r <- r + 1
 167         for (p in ONEtoK) {
 168             ## d <- st[[p]] # not needed
 169             ipos <- match(labels[i], ROWNAMES[[p]])
 170             if (!is.na(ipos)) Q[r, astart[p] + ipos] <- 1
 171         }
 172     }
 173
 174     for (p in 1:(k - 1)) {
 175         r <- r + 1
 176         for (i in ONEtoN) {
 177             ## d <- st[[p]]
 178             ipos <- match(labels[i], ROWNAMES[[p]])
 179             if (!is.na(ipos)) Q[r, astart[p] + ipos] <- 1
 180         }
 181     }
 182     a <- solve(Q, col, 1e-19)
 183     for (i in ONEtoN) {
 184         for (j in ONEtoN) {
 185             if (i == j) {
 186                 X[i, j] <- V[i, j] <- 0
 187                 next
 188             }
 189             sum <- 0
 190             sumv <- 0
 191             for (p in ONEtoK) {
 192                 d <- st[[p]]
 193                 pos <- match(labels[c(i, j)], ROWNAMES[[p]])
 194                 if (all(!is.na(pos))) {
 195                     ipos <- pos[1L]
 196                     jpos <- pos[2L]
 197                     dij <- d[ipos, jpos]
 198                     sum <- sum + sp[p] * (a[p] * dij + a[astart[p] + ipos] + a[astart[p] + jpos])
 199                     sumv <- sumv + sp[p] * (a[p] * dij)^2
 200                 }
 201             }
 202             X[i, j] <- sum / w[i, j]
 203             V[i, j] <- sumv / (w[i, j])^2
 204         }
 205     }
 206     list(X, V)
 207 }