R/SDM.R

   1 ## SDM.R (2011-10-11)
   2
   3 ## Construction of Consensus Distance Matrix With SDM
   4
   5 ## Copyright 2011 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     labels <- attr(as.dist(st[[1]]), "Labels")
  15     tot <- length(rownames(st[[1]]))
  16     for (i in 2:k) {
  17         labels <- union(labels, attr(as.dist(st[[i]]), "Labels"))
  18         tot <- tot + length(rownames(st[[i]]))
  19     }
  20     sp <- mat.or.vec(1,k)
  21     for (i in c(k+1:k))
  22         sp[i - k] <- st[[i]]
  23
  24     astart <- mat.or.vec(1, tot) # start of aip, astart[p] is start of aip
  25     astart[1] <- k
  26     for (i in 2:k)
  27         astart[i] <- astart[i - 1] + length(rownames(st[[i - 1]]))
  28
  29     a <- mat.or.vec(1, k + tot + k + length(labels))
  30     ## first k are alphas, subseqeunt ones aip
  31     ## each matrix p starting at astart[p], next are
  32     ## Lagrange multipliers, miu, niu, lambda in that order
  33     n <- length(labels)
  34     miustart <- k + tot
  35     niustart <- miustart + n
  36     lambstart <- niustart + k - 1
  37
  38     X <- mat.or.vec(n, n)
  39     V <- mat.or.vec(n, n)
  40     w <- mat.or.vec(n, n)
  41
  42     col <- mat.or.vec(k + tot + k + length(labels), 1) # free terms of system
  43
  44     dimnames(X) <- list(labels, labels)
  45     dimnames(V) <- list(labels, labels)
  46     dimnames(W) <- list(labels, labels)
  47
  48     for (i in 1:(n - 1)) {
  49         for (j in (i + 1):n) {
  50             for (p in 1:k) {
  51                 d <- st[[p]]
  52                 if (is.element(rownames(X)[i], rownames(d)) & is.element(colnames(X)[j], colnames(d))) {
  53                     w[i, j] <- w[j, i] <- w[i, j] + sp[p]
  54                 }
  55             }
  56         }
  57     }
  58
  59     Q <- mat.or.vec(length(labels) + k + k + tot, length(labels) + k + k + tot)
  60     ## first decompose first sum in paper
  61     for (p in 1:k) {
  62         d_p <- st[[p]]
  63         for (l in 1:k) { # first compute coeficients of alphas
  64             sum <- 0
  65             dijp <- -1
  66             if (l == p) { # calculate alpha_p
  67                 for (i in 1:n) {
  68                     for (j in 1:n) { #check if {i,j}\subset L_l
  69                         d <- st[[l]]
  70                         if (i != j & is.element(labels[i],rownames(d)) & is.element(labels[j],colnames(d))) {
  71                             dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
  72                             sum <- sum + ((dij*dij) - sp[l]*dij*dij/w[i,j])
  73                             ipos <- which(rownames(d) == labels[i])
  74                             jpos <- which(rownames(d) == labels[j])
  75                             Q[p, astart[l] + ipos] <- Q[p, astart[l] + ipos] + (dij - (sp[l]*dij/w[i,j]))
  76                             Q[p, astart[l] + jpos] <- Q[p, astart[l] + jpos] + (dij - (sp[l]*dij/w[i,j]))
  77                         }
  78                     }
  79                 }
  80             } else {
  81                 for (i in 1:n) {
  82                     for (j in 1:n) { #check if {i,j}\subset L_l
  83                         d <- st[[l]]
  84                         if (i != j & is.element(labels[i], rownames(d)) & is.element(labels[j], colnames(d)) & is.element(labels[i], rownames(d_p)) & is.element(labels[j], colnames(d_p))) {
  85                             dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
  86                             dijp <- d_p[rownames(d_p) == labels[i], colnames(d_p) == labels[j]]
  87                             sum <- sum - sp[l]*dij*dijp/w[i, j]
  88                             ipos <- which(rownames(d) == labels[i])
  89                             jpos <- which(rownames(d) == labels[j])
  90                             Q[p,astart[l] + ipos] <- Q[p, astart[l] + ipos] - sp[l]*dijp/w[i, j]
  91                             Q[p,astart[l] + jpos] <- Q[p, astart[l] + jpos] - sp[l]*dijp/w[i, j]
  92                         }
  93                     }
  94                 }
  95             }
  96             Q[p, l] <- sum
  97         }
  98         Q[p, lambstart + 1] <- 1
  99     }
 100     r <- k
 101     for (p in 1:k) {
 102         dp <- st[[p]]
 103         for (i in 1:n) {
 104             if (is.element(labels[i], rownames(dp))) {
 105                 r <- r + 1
 106                 for (l in 1:k) {
 107                     d <- st[[l]]
 108                     if (l == p) {
 109                         for (j in 1:n) {
 110                             if (i != j & is.element(labels[j], rownames(dp))) {
 111                                 dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
 112                                 Q[r, l] <- Q[r, l] + (dij - sp[l]*dij/w[i, j])
 113                                 ipos <- which(rownames(d) == labels[i])
 114                                 jpos <- which(rownames(d) == labels[j])
 115                                 Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] + (1 - sp[l]/w[i, j])
 116                                 Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] + (1 - sp[l]/w[i, j])
 117                             }
 118                         }
 119                     } else {
 120                         for (j in 1:n) {
 121                             if (i != j & is.element(labels[j], rownames(dp)) & is.element(labels[i], rownames(d)) & is.element(labels[j], colnames(d))) {
 122                                 dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
 123                                 Q[r,l] <- Q[r,l] - sp[l]*dij/w[i, j]
 124                                 ipos <- which(rownames(d) == labels[i])
 125                                 jpos <- which(rownames(d) == labels[j])
 126                                 Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] - sp[l]/w[i, j]
 127                                 Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] - sp[l]/w[i, j]
 128                             }
 129                         }
 130                     }
 131                 }
 132                 if (p < k) Q[r, ] <- Q[r, ] * sp[p]
 133                 Q[r, miustart + i] <- 1
 134                 if (p < k) Q[r, niustart + p] <- 1
 135             }
 136         }
 137     }
 138     r <- r + 1
 139     col[r] <- k
 140     for (i in 1:k) Q[r,i] <- 1
 141
 142     for (i in 1:n) {
 143         r <- r + 1
 144         for (p in 1:k) {
 145             d <- st[[p]]
 146             if (is.element(labels[i], rownames(d))) {
 147                 ipos <- which(rownames(d) == labels[i])
 148                 Q[r, astart[p] + ipos] <- 1
 149             }
 150         }
 151     }
 152     for (p in 1:(k - 1)) {
 153         r <- r + 1
 154         for (i in 1:n) {
 155             d <- st[[p]]
 156             if (is.element(labels[i], rownames(d))) {
 157                 ipos <- which(rownames(d) == labels[i])
 158                 Q[r, astart[p] + ipos] <- 1
 159             }
 160         }
 161     }
 162     a <- solve(Q, col, 1e-19)
 163     for(i in 1:n) {
 164         for(j in 1:n) {
 165             sum <- 0
 166             sumv <- 0
 167             for(p in 1:k) {
 168                 d <- st[[p]]
 169                 if (is.element(labels[i], rownames(d)) & is.element(labels[j], rownames(d))) {
 170                     ipos <- which(rownames(d) == labels[i])
 171                     jpos <- which(rownames(d) == labels[j])
 172                     sum <- sum + sp[p]*(a[p]*d[ipos, jpos] + a[astart[p] + ipos] + a[astart[p] + jpos])
 173                     sumv <- sumv + sp[p]*(a[p]*d[ipos, jpos])*(a[p]*d[ipos, jpos])
 174                 }
 175             }
 176             X[i, j] <- sum/w[i, j]
 177             V[i, j] <- sumv/(w[i, j] * w[i, j])
 178             if (i == j)
 179                 X[i, j] <- V[i, j] <- 0
 180         }
 181     }
 182     list(X, V)
 183 }