tropical_grassmannian/Figure 2.r at main · bhattsamir/tropical_grassmannian · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
library(ape)
library(phylo2vec)
library(phangorn)
library(distory)

# Normalize vector so that max(D) becomes h and negatives are truncated at 0
normalize <- function(D, h = 1) {
  a <- max(D) - h
  pmax(D - a, 0)
}

# Draw a random coalescent tree with n tips
draw_random_tree <- function(n) {
  Tr2 <- ape::rcoal(n)
  return(Tr2)
}

# Map a tree to its "tropical" representation (normalized cophenetic distances)
to_tropical <- function(Tr2) {
  D2 <- cophenetic.phylo(Tr2)
  D2 <- as.matrix(D2)
  e2 <- D2[lower.tri(D2)]
  e2 <- e2 / max(e2)
  e2 <- normalize(e2)
  return(e2)
}

# Reconstruct full symmetric matrix from lower-triangular vector (excluding diagonal)
lower_to_full <- function(v) {
  n <- (1 + sqrt(1 + 8 * length(v))) / 2
  m <- matrix(0, n, n)
  m[lower.tri(m)] <- v
  m <- m + t(m)
  return(m)
}

## Parameters --------------------------------------------------------------

set.seed(12345)
n    <- 50
num  <- 100

tree_samples <- vector("list", num)

## Initial trees and tropical representations ------------------------------

e_init1 <- to_tropical(draw_random_tree(n))
tree_samples_end <- ladderize(upgma(as.dist(lower_to_full(e_init1))))

e_init2 <- to_tropical(draw_random_tree(n))
tree_samples_start <- ladderize(upgma(as.dist(lower_to_full(e_init2))))

lambda <- e_init2 - e_init1

l <- seq(min(lambda), max(lambda), length.out = num)

## Quick checks for endpoints ----------------------------------------------

l_min <- min(lambda)
Dlmin <- pmax(l_min + e_init1, e_init2)
Dlmin <- normalize(Dlmin)

l_max <- max(lambda)
Dlmax <- pmax(l_max + e_init1, e_init2)
Dlmax <- normalize(Dlmax)

par(mfrow = c(1, 2))
plot(Dlmin, e_init2)
plot(Dlmax, e_init1)
all.equal(Dlmin, e_init2)
all.equal(Dlmax, e_init1)

par(mfrow = c(2, 2))
plot(ladderize(upgma(as.dist(lower_to_full(Dlmin)))),  main = "Start tree tropical")
plot(ladderize(upgma(as.dist(lower_to_full(e_init2)))), main = "Start tree")
plot(ladderize(upgma(as.dist(lower_to_full(Dlmax)))),  main = "End tree tropical")
plot(ladderize(upgma(as.dist(lower_to_full(e_init1)))), main = "End tree")

## Path in tropical space and tree samples --------------------------------

tropical_distance <- c()

for (i in 1:num) {
  Dls <- pmax(l[i] + e_init1, e_init2)
  Dls <- normalize(Dls)
  tropical_distance[i] <- max(e_init2 - Dls) - min(e_init2 - Dls)
  tree_samples[[i]] <- ladderize(upgma(as.dist(lower_to_full(Dls))))
}

class(tree_samples) <- "multiPhylo"
tree_samples <- c(tree_samples_start, tree_samples, tree_samples_end)

## Distances: BHV, SPR, RF, tropical --------------------------------------

d   <- dist.multiPhylo(tree_samples)
d   <- as.matrix(d)
BHV <- d[, 1]

SPR <- c()
for (i in 1:length(tree_samples)) {
  SPR[i] <- SPR.dist(tree_samples[[1]], tree_samples[[i]])
}

RF <- c()
for (i in 1:length(tree_samples)) {
  RF[i] <- RF.dist(tree_samples[[1]], tree_samples[[i]], normalize = TRUE)
}

par(mfrow = c(2, 3))
plot(ladderize(upgma(as.dist(lower_to_full(e_init2)))), main = "Start tree",show.tip.label = FALSE)
plot(ladderize(upgma(as.dist(lower_to_full(e_init1)))), main = "End tree",show.tip.label = FALSE)
plot(BHV, type = "b", pch = 16, main = "BHV distance",      ylab = "distance", xlab = "tree")
plot(SPR, type = "b", pch = 16, main = "SPR distance",      ylab = "distance", xlab = "tree")
plot(RF,  type = "b", pch = 16, main = "RF distance",       ylab = "distance", xlab = "tree")
plot(tropical_distance,
     type = "b", pch = 16, main = "Tropical distance", ylab = "distance", xlab = "tree")