3 ## Construction of Consensus Distance Matrix With SDM
5 ## Copyright 2011 Andrei-Alin Popescu
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
12 st <- list(...) # first half contains matrices, second half s_p
14 labels <- attr(as.dist(st[[1]]), "Labels")
15 tot <- length(rownames(st[[1]]))
17 labels <- union(labels, attr(as.dist(st[[i]]), "Labels"))
18 tot <- tot + length(rownames(st[[i]]))
24 astart <- mat.or.vec(1, tot) # start of aip, astart[p] is start of aip
27 astart[i] <- astart[i - 1] + length(rownames(st[[i - 1]]))
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
35 niustart <- miustart + n
36 lambstart <- niustart + k - 1
42 col <- mat.or.vec(k + tot + k + length(labels), 1) # free terms of system
44 dimnames(X) <- list(labels, labels)
45 dimnames(V) <- list(labels, labels)
46 dimnames(W) <- list(labels, labels)
48 for (i in 1:(n - 1)) {
49 for (j in (i + 1):n) {
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]
59 Q <- mat.or.vec(length(labels) + k + k + tot, length(labels) + k + k + tot)
60 ## first decompose first sum in paper
63 for (l in 1:k) { # first compute coeficients of alphas
66 if (l == p) { # calculate alpha_p
68 for (j in 1:n) { #check if {i,j}\subset L_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]))
82 for (j in 1:n) { #check if {i,j}\subset L_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]
98 Q[p, lambstart + 1] <- 1
104 if (is.element(labels[i], rownames(dp))) {
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])
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]
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
140 for (i in 1:k) Q[r,i] <- 1
146 if (is.element(labels[i], rownames(d))) {
147 ipos <- which(rownames(d) == labels[i])
148 Q[r, astart[p] + ipos] <- 1
152 for (p in 1:(k - 1)) {
156 if (is.element(labels[i], rownames(d))) {
157 ipos <- which(rownames(d) == labels[i])
158 Q[r, astart[p] + ipos] <- 1
162 a <- solve(Q, col, 1e-19)
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])
176 X[i, j] <- sum/w[i, j]
177 V[i, j] <- sumv/(w[i, j] * w[i, j])
179 X[i, j] <- V[i, j] <- 0