-## SDM.R (2011-10-11)
+## SDM.R (2012-04-02)
## Construction of Consensus Distance Matrix With SDM
-## Copyright 2011 Andrei-Alin Popescu
+## Copyright 2011-2012 Andrei-Alin Popescu
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
{
st <- list(...) # first half contains matrices, second half s_p
k <- length(st)/2
- labels <- attr(as.dist(st[[1]]), "Labels")
- tot <- length(rownames(st[[1]]))
- for (i in 2:k) {
- labels <- union(labels, attr(as.dist(st[[i]]), "Labels"))
- tot <- tot + length(rownames(st[[i]]))
- }
- sp <- mat.or.vec(1,k)
- for (i in c(k + 1:k))
- sp[i - k] <- st[[i]]
+ ONEtoK <- seq_len(k)
+
+ ## make sure we have only matrices:
+ for (i in ONEtoK) st[[i]] <- as.matrix(st[[i]])
+
+ ## store the rownames of each matrix in a list because they are often called:
+ ROWNAMES <- lapply(st[ONEtoK], rownames)
+
+ ## the number of rows of each matrix:
+ NROWS <- sapply(ROWNAMES, length)
+ tot <- sum(NROWS)
+
+ labels <- unique(unlist(ROWNAMES))
+ sp <- unlist(st[k + ONEtoK])
- astart <- mat.or.vec(1, tot) # start of aip, astart[p] is start of aip
+ astart <- numeric(tot) # start of aip, astart[p] is start of aip
astart[1] <- k
for (i in 2:k)
- astart[i] <- astart[i - 1] + length(rownames(st[[i - 1]]))
+ astart[i] <- astart[i - 1] + NROWS[i - 1]
- a <- mat.or.vec(1, k + tot + k + length(labels))
- ## first k are alphas, subseqeunt ones aip
+ ## apparently erased by the operation below so no need to initialize:
+ ## a <- mat.or.vec(1, k + tot + k + length(labels))
+
+ ## first k are alphas, subsequent ones aip
## each matrix p starting at astart[p], next are
## Lagrange multipliers, miu, niu, lambda in that order
n <- length(labels)
niustart <- miustart + n
lambstart <- niustart + k - 1
- X <- mat.or.vec(n, n)
- V <- mat.or.vec(n, n)
- w <- mat.or.vec(n, n)
-
- col <- mat.or.vec(k + tot + k + length(labels), 1) # free terms of system
+ X <- matrix(0, n, n, dimnames = list(labels, labels))
+ V <- w <- X
- dimnames(X) <- list(labels, labels)
- dimnames(V) <- list(labels, labels)
- dimnames(W) <- list(labels, labels)
+ tmp <- 2 * k + tot + n
+ col <- numeric(tmp) # free terms of system
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
- for (p in 1:k) {
- d <- st[[p]]
- if (is.element(rownames(X)[i], rownames(d)) & is.element(colnames(X)[j], colnames(d))) {
+ for (p in ONEtoK) {
+ ## d <- st[[p]] # not needed anymore
+ if (is.element(labels[i], ROWNAMES[[p]]) && is.element(labels[j], ROWNAMES[[p]])) {
w[i, j] <- w[j, i] <- w[i, j] + sp[p]
}
}
}
}
- Q <- mat.or.vec(length(labels) + k + k + tot, length(labels) + k + k + tot)
+ ONEtoN <- seq_len(n)
+
+ Q <- matrix(0, tmp, tmp)
## first decompose first sum in paper
- for (p in 1:k) {
+ for (p in ONEtoK) {
d_p <- st[[p]]
- for (l in 1:k) { # first compute coeficients of alphas
+ for (l in ONEtoK) { # first compute coefficients of alphas
+ d <- st[[l]]
sum <- 0
dijp <- -1
if (l == p) { # calculate alpha_p
- for (i in 1:n) {
- for (j in 1:n) { # check if {i,j}\subset L_l
- d <- st[[l]]
- if (i != j & is.element(labels[i],rownames(d)) & is.element(labels[j],colnames(d))) {
- dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
- sum <- sum + ((dij*dij) - sp[l]*dij*dij/w[i,j])
- ipos <- which(rownames(d) == labels[i])
- jpos <- which(rownames(d) == labels[j])
- Q[p, astart[l] + ipos] <- Q[p, astart[l] + ipos] + (dij - (sp[l]*dij/w[i,j]))
- Q[p, astart[l] + jpos] <- Q[p, astart[l] + jpos] + (dij - (sp[l]*dij/w[i,j]))
+ for (i in ONEtoN) {
+ for (j in ONEtoN) { # check if {i,j}\subset L_l
+ if (i == j) next # make sure i != j
+ ## d <- st[[l]] # <- moved-up
+ pos <- match(labels[c(i, j)], ROWNAMES[[l]]) # <- returns NA if not in this matrix
+ if (all(!is.na(pos))) {
+ ipos <- pos[1L]
+ jpos <- pos[2L]
+ dij <- d[ipos, jpos]
+ sum <- sum + dij * dij - sp[l] * dij * dij / w[i,j]
+ tmp2 <- dij - sp[l] * dij / w[i,j]
+ Q[p, astart[l] + ipos] <- Q[p, astart[l] + ipos] + tmp2
+ Q[p, astart[l] + jpos] <- Q[p, astart[l] + jpos] + tmp2
}
}
}
} else {
- for (i in 1:n) {
- for (j in 1:n) { # check if {i,j}\subset L_l
- d <- st[[l]]
- 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))) {
- dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
- dijp <- d_p[rownames(d_p) == labels[i], colnames(d_p) == labels[j]]
- sum <- sum - sp[l]*dij*dijp/w[i, j]
- ipos <- which(rownames(d) == labels[i])
- jpos <- which(rownames(d) == labels[j])
- Q[p,astart[l] + ipos] <- Q[p, astart[l] + ipos] - sp[l]*dijp/w[i, j]
- Q[p,astart[l] + jpos] <- Q[p, astart[l] + jpos] - sp[l]*dijp/w[i, j]
+ for (i in ONEtoN) {
+ for (j in ONEtoN) { # check if {i,j}\subset L_l
+ if (i == j) next
+ ## d <- st[[l]] # <- moved-up
+ pos <- match(labels[c(i, j)], ROWNAMES[[l]])
+ posp <- match(labels[c(i, j)], ROWNAMES[[p]])
+ if (all(!is.na(pos)) && all(!is.na(posp))) {
+ ipos <- pos[1L]
+ jpos <- pos[2L]
+ dij <- d[ipos, jpos]
+ dijp <- d_p[posp[1L], posp[2L]]
+ sum <- sum - sp[l] * dij * dijp / w[i, j]
+ tmp2 <- sp[l] * dijp / w[i, j]
+ Q[p,astart[l] + ipos] <- Q[p, astart[l] + ipos] - tmp2
+ Q[p,astart[l] + jpos] <- Q[p, astart[l] + jpos] - tmp2
}
}
}
}
Q[p, lambstart + 1] <- 1
}
+
r <- k
- for (p in 1:k) {
+
+ for (p in ONEtoK) {
dp <- st[[p]]
- for (i in 1:n) {
- if (is.element(labels[i], rownames(dp))) {
+ for (i in ONEtoN) {
+ if (is.element(labels[i], ROWNAMES[[p]])) {
r <- r + 1
- for (l in 1:k) {
+ for (l in ONEtoK) {
d <- st[[l]]
if (l == p) {
- for (j in 1:n) {
- if (i != j & is.element(labels[j], rownames(dp))) {
- dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
- Q[r, l] <- Q[r, l] + (dij - sp[l]*dij/w[i, j])
- ipos <- which(rownames(d) == labels[i])
- jpos <- which(rownames(d) == labels[j])
- Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] + (1 - sp[l]/w[i, j])
- Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] + (1 - sp[l]/w[i, j])
+ ipos <- match(labels[i], ROWNAMES[[p]])
+ for (j in ONEtoN) {
+ if (i == j) next
+ jpos <- match(labels[j], ROWNAMES[[p]])
+ if (!is.na(jpos)) {
+ dij <- d[ipos, jpos]
+ Q[r, l] <- Q[r, l] + dij - sp[l] * dij / w[i, j]
+ tmp2 <- 1 - sp[l] / w[i, j]
+ Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] + tmp2
+ Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] + tmp2
}
}
} else {
- for (j in 1:n) {
- if (i != j & is.element(labels[j], rownames(dp)) & is.element(labels[i], rownames(d)) & is.element(labels[j], colnames(d))) {
- dij <- d[rownames(d) == labels[i], colnames(d) == labels[j]]
- Q[r,l] <- Q[r,l] - sp[l]*dij/w[i, j]
- ipos <- which(rownames(d) == labels[i])
- jpos <- which(rownames(d) == labels[j])
- Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] - sp[l]/w[i, j]
- Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] - sp[l]/w[i, j]
+ for (j in ONEtoN) {
+ if (i == j) next
+ if (!is.element(labels[j], rownames(dp))) next
+ pos <- match(labels[c(i, j)], ROWNAMES[[l]])
+ if (all(!is.na(pos))) {
+ ipos <- pos[1L]
+ jpos <- pos[2L]
+ dij <- d[ipos, jpos]
+ Q[r, l] <- Q[r, l] - sp[l] * dij / w[i, j]
+ tmp2 <- sp[l]/w[i, j]
+ Q[r, astart[l] + ipos] <- Q[r, astart[l] + ipos] - tmp2
+ Q[r, astart[l] + jpos] <- Q[r, astart[l] + jpos] - tmp2
}
}
}
}
}
}
+
r <- r + 1
col[r] <- k
- for (i in 1:k) Q[r,i] <- 1
+ Q[r, ONEtoK] <- 1
+ ## for (i in 1:k) Q[r, i] <- 1
- for (i in 1:n) {
+ for (i in ONEtoN) {
r <- r + 1
- for (p in 1:k) {
- d <- st[[p]]
- if (is.element(labels[i], rownames(d))) {
- ipos <- which(rownames(d) == labels[i])
- Q[r, astart[p] + ipos] <- 1
- }
+ for (p in ONEtoK) {
+ ## d <- st[[p]] # not needed
+ ipos <- match(labels[i], ROWNAMES[[p]])
+ if (!is.na(ipos)) Q[r, astart[p] + ipos] <- 1
}
}
+
for (p in 1:(k - 1)) {
r <- r + 1
- for (i in 1:n) {
- d <- st[[p]]
- if (is.element(labels[i], rownames(d))) {
- ipos <- which(rownames(d) == labels[i])
- Q[r, astart[p] + ipos] <- 1
- }
+ for (i in ONEtoN) {
+ ## d <- st[[p]]
+ ipos <- match(labels[i], ROWNAMES[[p]])
+ if (!is.na(ipos)) Q[r, astart[p] + ipos] <- 1
}
}
a <- solve(Q, col, 1e-19)
- for (i in 1:n) {
- for (j in 1:n) {
+ for (i in ONEtoN) {
+ for (j in ONEtoN) {
+ if (i == j) {
+ X[i, j] <- V[i, j] <- 0
+ next
+ }
sum <- 0
sumv <- 0
- for (p in 1:k) {
+ for (p in ONEtoK) {
d <- st[[p]]
- if (is.element(labels[i], rownames(d)) & is.element(labels[j], rownames(d))) {
- ipos <- which(rownames(d) == labels[i])
- jpos <- which(rownames(d) == labels[j])
- sum <- sum + sp[p]*(a[p]*d[ipos, jpos] + a[astart[p] + ipos] + a[astart[p] + jpos])
- sumv <- sumv + sp[p]*(a[p]*d[ipos, jpos])*(a[p]*d[ipos, jpos])
+ pos <- match(labels[c(i, j)], ROWNAMES[[p]])
+ if (all(!is.na(pos))) {
+ ipos <- pos[1L]
+ jpos <- pos[2L]
+ dij <- d[ipos, jpos]
+ sum <- sum + sp[p] * (a[p] * dij + a[astart[p] + ipos] + a[astart[p] + jpos])
+ sumv <- sumv + sp[p] * (a[p] * dij)^2
}
}
- X[i, j] <- sum/w[i, j]
- V[i, j] <- sumv/(w[i, j] * w[i, j])
- if (i == j)
- X[i, j] <- V[i, j] <- 0
+ X[i, j] <- sum / w[i, j]
+ V[i, j] <- sumv / (w[i, j])^2
}
}
list(X, V)
projects / ape.git / blobdiff