hyperbolic tree

R code of the ggplot2 example

library(ggplot2)
library(ape)
library(hydra)

## The 137 bird families of Sibley & Ahlquist in the Poincaré disk. The
## layout is hydraPlus() (Keller-Ressel & Nargang 2020): a hyperbolic
## embedding of the path lengths along the tree, which keeps the branches
## apart as they crowd toward the rim. A Möbius transformation then puts the
## focus on the root, then on the passerines, which open up while the rest
## of the tree is pushed to the edge. Links are hyperbolic geodesics: arcs
## of circles that meet the rim at right angles.
data(bird.families)
tree <- reorder(bird.families, "cladewise")
n_tip <- Ntip(tree)
tree$edge.length <- rep(1, nrow(tree$edge))
## Path lengths halved: at full length the embedding overflows (hydra warns).
embedding <- hydraPlus(dist.nodes(tree) * 0.5, dim = 2, curvature = 1,
                       control = list(trace = 0))
z <- complex(modulus = embedding$r, argument = embedding$theta)

passerines <- getMRCA(tree, c("Tyrannidae", "Fringillidae"))
## Cladewise order: a parent comes before its children.
in_passerines <- seq_along(z) == passerines
for (i in seq_len(nrow(tree$edge))) {
  in_passerines[tree$edge[i, 2]] <- in_passerines[tree$edge[i, 2]] || in_passerines[tree$edge[i, 1]]
}

## Geodesic between p and q: a straight segment through the centre, otherwise
## the circle through p, q and the inverse of p in the unit circle.
geodesic <- function(p, q, n = 30) {
  if (abs(Im(Conj(p) * q)) < 1e-9) {
    return(seq(0, 1, length.out = n) * (q - p) + p)
  }
  r <- p / Mod(p)^2
  d <- 2 * (Re(p) * (Im(q) - Im(r)) + Re(q) * (Im(r) - Im(p)) + Re(r) * (Im(p) - Im(q)))
  cx <- (Mod(p)^2 * (Im(q) - Im(r)) + Mod(q)^2 * (Im(r) - Im(p)) + Mod(r)^2 * (Im(p) - Im(q))) / d
  cy <- (Mod(p)^2 * (Re(r) - Re(q)) + Mod(q)^2 * (Re(p) - Re(r)) + Mod(r)^2 * (Re(q) - Re(p))) / d
  centre <- complex(real = cx, imaginary = cy)
  start <- Arg(p - centre)
  turn <- Arg((q - centre) / (p - centre))
  centre + Mod(p - centre) * exp(1i * (start + seq(0, turn, length.out = n)))
}

layout <- function(z, focus) {
  edges <- do.call(rbind, lapply(seq_len(nrow(tree$edge)), function(i) {
    w <- geodesic(z[tree$edge[i, 1]], z[tree$edge[i, 2]])
    data.frame(edge = i, x = Re(w), y = Im(w), passerine = in_passerines[tree$edge[i, 2]])
  }))
  nodes <- data.frame(x = Re(z), y = Im(z), passerine = in_passerines)[seq_len(n_tip), ]
  edges$focus <- focus
  nodes$focus <- focus
  list(edges = edges, nodes = nodes)
}
## The Möbius transformation that brings node `a` to the centre.
focus <- function(a) {
  (z - z[a]) / (1 - Conj(z[a]) * z)
}
views <- list(layout(focus(n_tip + 1), "Focus on the root"),
              layout(focus(passerines), "Focus on the passerines"))
edges <- do.call(rbind, lapply(views, `[[`, "edges"))
nodes <- do.call(rbind, lapply(views, `[[`, "nodes"))
edges$focus <- factor(edges$focus, levels = c("Focus on the root", "Focus on the passerines"))
nodes$focus <- factor(nodes$focus, levels = levels(edges$focus))
rim <- data.frame(x = cos(seq(0, 2 * pi, length.out = 200)), y = sin(seq(0, 2 * pi, length.out = 200)))

p <- ggplot() +
  geom_path(data = rim, aes(x, y), colour = "#d9d9d9", linewidth = 0.4) +
  geom_path(data = edges, aes(x, y, group = edge, colour = passerine), linewidth = 0.35) +
  geom_point(data = nodes, aes(x, y, colour = passerine), size = 0.7) +
  scale_colour_manual(values = c(`TRUE` = "#9b2367", `FALSE` = "#2e7891")) +
  facet_wrap(~focus) +
  coord_equal(xlim = c(-1.02, 1.02), ylim = c(-1.02, 1.02), expand = FALSE) +
  labs(title = "Moving the focus spreads out the 42 passerine families",
       subtitle = paste("The 137 bird families in the Poincaré disk;",
                        "<span style='color:#9b2367'>**passerines**</span> and",
                        "<span style='color:#2e7891'>**other orders**</span>"),
       caption = "Data: Sibley & Ahlquist 1990, via the ape R package") +
  theme_lex_viz(grid = "", arrows = "") +
  theme(legend.position = "none", axis.text = element_blank(), axis.title = element_blank(),
        axis.line = element_blank(), axis.ticks = element_blank())