3 ## Construction of Consensus Distance Matrix With SDM
5 ## Copyright 2011-2012 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
16 ## make sure we have only matrices:
17 for (i in ONEtoK) st[[i]] <- as.matrix(st[[i]])
19 ## store the rownames of each matrix in a list because they are often called:
20 ROWNAMES <- lapply(st[ONEtoK], rownames)
22 ## the number of rows of each matrix:
23 NROWS <- sapply(ROWNAMES, length)
26 labels <- unique(unlist(ROWNAMES))
27 sp <- unlist(st[k + ONEtoK])
29 astart <- numeric(tot) # start of aip, astart[p] is start of aip
32 astart[i] <- astart[i - 1] + NROWS[i - 1]
34 ## apparently erased by the operation below so no need to initialize:
35 ## a <- mat.or.vec(1, k + tot + k + length(labels))
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
42 niustart <- miustart + n
43 lambstart <- niustart + k - 1
45 X <- matrix(0, n, n, dimnames = list(labels, labels))
48 tmp <- 2 * k + tot + n
49 col <- numeric(tmp) # free terms of system
51 for (i in 1:(n - 1)) {
52 for (j in (i + 1):n) {
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]
64 Q <- matrix(0, tmp, tmp)
65 ## first decompose first sum in paper
68 for (l in ONEtoK) { # first compute coefficients of alphas
72 if (l == p) { # calculate alpha_p
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))) {
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
91 for (j in ONEtoN) { # check if {i,j}\subset L_l
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))) {
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
111 Q[p, lambstart + 1] <- 1
119 if (is.element(labels[i], ROWNAMES[[p]])) {
124 ipos <- match(labels[i], ROWNAMES[[p]])
127 jpos <- match(labels[j], ROWNAMES[[p]])
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
139 if (!is.element(labels[j], rownames(dp))) next
140 pos <- match(labels[c(i, j)], ROWNAMES[[l]])
141 if (all(!is.na(pos))) {
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
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
163 ## for (i in 1:k) Q[r, i] <- 1
168 ## d <- st[[p]] # not needed
169 ipos <- match(labels[i], ROWNAMES[[p]])
170 if (!is.na(ipos)) Q[r, astart[p] + ipos] <- 1
174 for (p in 1:(k - 1)) {
178 ipos <- match(labels[i], ROWNAMES[[p]])
179 if (!is.na(ipos)) Q[r, astart[p] + ipos] <- 1
182 a <- solve(Q, col, 1e-19)
186 X[i, j] <- V[i, j] <- 0
193 pos <- match(labels[c(i, j)], ROWNAMES[[p]])
194 if (all(!is.na(pos))) {
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
202 X[i, j] <- sum / w[i, j]
203 V[i, j] <- sumv / (w[i, j])^2