#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 950
## file: app.R
# twist-explorer: markings and twists on a genus-2 surface ------------------
#
# The genus-2 analogue of the torus "two markings of the same lattice" app.
#
#   * A marking is a generating set (A1, B1, A2, B2) of the surface group with
#     [A1, B1][A2, B2] = 1. Its polygon has vertices p, W1 p, W2 p, ..., where
#     W_k is the product of the first k letters of the relator
#     A1 B1 A1^-1 B1^-1 A2 B2 A2^-1 B2^-1 and p is a fixed basepoint.
#     For the regular octagon's own generators this is the regular octagon.
#   * "Deform the surface" is a continuous Fenchel-Nielsen twist along one of
#     the curves a1, b1, a2, b2, or the separating curve c = [a1, b1]. The group
#     changes; at t = 1 it is the original group again, with the marking
#     changed by a Dehn twist (Buser 3.3, and iota^-1 tau^(a+1) = D tau^a).
#   * "Change the marking" applies Dehn twists to the generators. The group
#     does not change: the new polygon (orange) has its corners on the same
#     orbit points as the old one (green).
#   * The tiling of the new marking (light orange) is the set of images g(P')
#     of the new polygon P' under the same group: the tiles of the old marking
#     cut up and reassembled along different curves. Grey (old tiles) and light
#     orange (new tiles) share their vertices: both are the orbit of p.
#
# Run with  shiny::runApp("path/to/hyperbolic/twist-explorer")
# Requires hyperbolic.R and teichmuller.R in the parent folder.
# All group elements here are 2x2 complex matrices acting on the Poincare disk.
# ---------------------------------------------------------------------------

library(shiny)
library(tidyverse)
source("hyperbolic.R")
source("teichmuller.R")

col_ref   <- "grey85"   # tiling of the original (regular) surface
col_tile  <- "grey60"   # tiling of the current surface
col_mark  <- "#1D9E75"  # polygon of the current marking
col_new   <- "#D85A30"  # polygon after the change of marking
col_newtile <- "#EDA584" # tiling by the new polygon
col_traj  <- "#C0392B"  # vertex trajectories
col_curve <- c(a1 = "#185FA5", b1 = "#993556", a2 = "#534AB7", b2 = "#854F0B", c = "#2C2C2A")

# --- Group helpers (disk-model matrices) ----------------------------------------

det2 <- function(M) M[1, 1] * M[2, 2] - M[1, 2] * M[2, 1]

# Normalize to det 1 and nonnegative real trace (one representative in PSL).
psl <- function(M) {
  M <- M / sqrt(det2(M))
  if (Re(M[1, 1] + M[2, 2]) < 0) -M else M
}

minv <- function(M) psl(solve(M))

# M^t for a hyperbolic M: translation along the axis of M by t times its length.
mpow <- function(M, t) {
  if (t == 0) return(diag(2) + 0i)
  if (t == round(t)) {
    R <- diag(2) + 0i
    B <- if (t > 0) M else minv(M)
    for (k in seq_len(abs(t))) R <- R %*% B
    return(psl(R))
  }
  e <- eigen(psl(M))
  psl(e$vectors %*% diag(e$values^t) %*% solve(e$vectors))
}

comm <- function(X, Y) psl(X %*% Y %*% minv(X) %*% minv(Y))

trace_len <- function(M) {
  tr <- abs(Re(M[1, 1] + M[2, 2]))
  if (tr <= 2) 0 else 2 * acosh(tr / 2)
}

# --- The regular octagon's marking --------------------------------------------------

# Loop generators at the vertex v1 of the regular octagon. h_k is the unique
# group element with h_k(v1) = v_k, and the relator letters are h_k^-1 h_{k+1}.
regular_marking <- function() {
  P <- surface_polygon(2, "disk")
  v <- P$vertices
  g <- map(P$gens, \(x) iso_matrix(x, "disk"))
  words <- list(diag(2) + 0i)
  frontier <- words
  for (d in 1:4) {
    frontier <- unlist(lapply(frontier, \(M) lapply(g, \(G) M %*% G)), recursive = FALSE)
    words <- c(words, frontier)
  }
  h <- map(1:8, \(k) {
    words[[detect_index(words, \(M) Mod(mobius(M, v[1]) - v[k]) < 1e-9)]]
  })
  L <- map(1:8, \(k) psl(solve(h[[k]]) %*% h[[k %% 8 + 1]]))
  list(p = v[1], gens = list(A1 = L[[1]], B1 = L[[2]], A2 = L[[5]], B2 = L[[6]]))
}

relator_letters <- function(G) {
  list(G$A1, G$B1, minv(G$A1), minv(G$B1), G$A2, G$B2, minv(G$A2), minv(G$B2))
}

# Prefix products W_0 = I, W_1, ..., W_7.
prefixes <- function(G) {
  L <- relator_letters(G)
  accumulate(L[-8], \(M, X) M %*% X, .init = diag(2) + 0i)
}

marking_polygon <- function(G, p) map_vec(prefixes(G), \(W) mobius(W, p))

# Side pairings of the marking polygon. Side k runs from x_{k-1} to x_k and is
# W_{k-1}[p, L_k p]. If L_j = L_i^-1, then g = W_{j-1} L_i^-1 W_{i-1}^-1 maps
# side i onto side j (reversed).
polygon_pairings <- function(G) {
  W <- prefixes(G)
  L <- relator_letters(G)
  pairs <- list(c(1, 3), c(2, 4), c(5, 7), c(6, 8))
  gs <- map(pairs, \(ij) psl(W[[ij[2]]] %*% minv(L[[ij[1]]]) %*% minv(W[[ij[1]]])))
  c(gs, map(gs, minv))   # 8 pairings: g1..g4 and their inverses g5..g8
}

# Group elements reached by at most `depth` side crossings (no backtracking).
tile_elements <- function(G, depth) {
  gs <- polygon_pairings(G)
  inv_of <- c(5:8, 1:4)
  out <- list(list(M = diag(2) + 0i, last = 0L))
  frontier <- out
  for (d in seq_len(depth)) {
    frontier <- unlist(lapply(frontier, \(e) {
      lapply(seq_along(gs), \(j) {
        if (e$last != 0L && inv_of[j] == e$last) NULL
        else list(M = e$M %*% gs[[j]], last = j)
      })
    }), recursive = FALSE) |> compact()
    out <- c(out, frontier)
  }
  map(out, "M")
}

# Points along the disk geodesics z1[i] -> z2[i], vectorized over edges.
#
# A disk geodesic is an arc of a circle meeting the unit circle at right angles
# (or a diameter). That circle passes through z1, z2 and the inverse point
# z* = z / |z|^2 of either endpoint, so its center is the circumcenter of three
# known points. We sample the arc EVENLY IN ANGLE about that center. That gives
# a smooth curve whatever the edge's length or position. (Sampling a straight
# segment in a moved copy of the disk, as an earlier version did, leaves long
# straight chords on long edges near the boundary circle.)
# Edges through the origin, or nearly so, are drawn as straight segments.
geodesic_points <- function(z1, z2, n = 24) {
  s <- seq(0, 1, length.out = n)
  # invert the endpoint farther from 0, for numerical stability
  far <- ifelse(Mod(z1) >= Mod(z2), z1, z2)
  z3  <- far / Mod(far)^2
  # circumcenter of z1, z2, z3:  center = z1 + (|u|^2 v - |v|^2 u) / (conj(u) v - u conj(v))
  u <- z2 - z1; v <- z3 - z1
  den <- Conj(u) * v - u * Conj(v)                     # = 2i Im(conj(u) v)
  straight <- Mod(far) < 1e-9 | Mod(den) < 1e-10 * Mod(u) * Mod(v)
  ctr <- z1 + (Mod(u)^2 * v - Mod(v)^2 * u) / ifelse(straight, 1, den)
  th1 <- Arg(z1 - ctr); th2 <- Arg(z2 - ctr)
  dth <- (th2 - th1 + pi) %% (2 * pi) - pi            # the short way round
  r   <- Mod(z1 - ctr)
  arc <- ctr + r * exp(1i * (th1 + outer(dth, s)))     # rows = edges, columns = points
  line <- z1 + outer(u, s)
  arc[straight, ] <- line[straight, ]
  arc[, 1] <- z1; arc[, n] <- z2                       # exact endpoints
  arc
}

# Boundary paths of many polygons at once (a list of vertex vectors).
polygons_df <- function(polys, n = 24) {
  k  <- rep(seq_along(polys), lengths(polys))
  z1 <- unlist(polys)
  z2 <- unlist(map(polys, \(v) c(v[-1], v[1])))
  z  <- as.vector(t(geodesic_points(z1, z2, n)))   # edges in order, n points each
  tibble(x = Re(z), y = Im(z), tile = rep(k, each = n))
}

# The tiling by the polygon of the marking G: its images under the group,
# reached through the side pairings of G's own polygon. Applied to a new
# marking H of the same group, this is the "new tiling": the same group,
# so the same orbit of vertices, but different tiles.
tiling_paths <- function(G, p, depth, n = 24) {
  poly <- marking_polygon(G, p)
  polygons_df(map(tile_elements(G, depth), \(M) mobius(M, poly)), n)
}

orbit_points <- function(G, p, depth) {
  poly <- marking_polygon(G, p)
  z <- unlist(map(tile_elements(G, depth), \(M) mobius(M, poly)))
  z <- z[Mod(z) < 0.995]
  z[!duplicated(round(z, 7))]
}

# --- Labelling the curves of a marking ------------------------------------------------
#
# Side k of the marking polygon runs from x_{k-1} to x_k and is a lift of the
# based loop for the k-th relator letter: a1, b1, a1^-1, b1^-1, a2, b2, a2^-1,
# b2^-1. The diagonal from x_0 = p to x_4 = [A1,B1] p is a lift of the based
# loop c, which cuts the surface into two one-holed tori.

side_names <- c("a1", "b1", "a1", "b1", "a2", "b2", "a2", "b2")
side_labels <- c("a1", "b1", "a1\u207b\u00b9", "b1\u207b\u00b9",
                 "a2", "b2", "a2\u207b\u00b9", "b2\u207b\u00b9")

# Hyperbolic midpoint of the segment z1 z2 (vectorized): move z1 to 0, take the
# point halfway out in hyperbolic distance, move back.
hmidpoint <- function(z1, z2) {
  w <- (z2 - z1) / (1 - Conj(z1) * z2)
  u <- w / Mod(w) * tanh(atanh(Mod(w)) / 2)
  (u + z1) / (1 + Conj(z1) * u)
}

marking_sides_df <- function(G, p, n = 30) {
  x <- marking_polygon(G, p)
  z1 <- x; z2 <- c(x[-1], x[1])
  pts <- geodesic_points(z1, z2, n)
  centre <- mean(x)
  mid <- hmidpoint(z1, z2)
  lab <- mid + 0.13 * (mid - centre) / Mod(mid - centre)   # nudge outward
  list(
    paths = tibble(x = Re(as.vector(t(pts))), y = Im(as.vector(t(pts))),
                   side = rep(1:8, each = n), curve = rep(side_names, each = n)),
    labels = tibble(x = Re(lab), y = Im(lab), label = side_labels, curve = side_names),
    c_path = { z <- geodesic_points(x[1], x[5], n)[1, ]; tibble(x = Re(z), y = Im(z)) },
    c_label = { z <- hmidpoint(x[1], hmidpoint(x[1], x[5]))   # a quarter of the way down
                tibble(x = Re(z) + 0.05, y = Im(z), label = "c") }
  )
}

# Fixed points of a hyperbolic disk isometry: roots of c z^2 + (d - a) z - b = 0.
fixed_points <- function(M) {
  a <- M[1, 1]; b <- M[1, 2]; c <- M[2, 1]; d <- M[2, 2]
  disc <- sqrt((d - a)^2 + 4 * b * c)
  (-(d - a) + c(-1, 1) * disc) / (2 * c)
}

# Points (x, y) inside a polygon given by its boundary path (ray casting).
inside_polygon <- function(x, y, px, py) {
  n <- length(px); j <- c(n, seq_len(n - 1)); inside <- logical(length(x))
  for (i in seq_len(n)) {
    cross <- ((py[i] > y) != (py[j[i]] > y)) &
      (x < (px[j[i]] - px[i]) * (y - py[i]) / (py[j[i]] - py[i]) + px[i])
    inside <- xor(inside, cross)
  }
  inside
}

# The closed geodesics a1, b1, a2, b2, c of the marking, drawn on its polygon D.
# Each lifts to the axes of the conjugates g X g^-1 (g in the group); we keep the
# parts of those axes that lie inside D. On the surface these pieces join up into
# the closed geodesic, so their lengths add up to its length.
closed_geodesics_df <- function(G, p, depth = 3, n = 400) {
  D <- ngon_path(marking_polygon(G, p), "disk", n = 40)
  rmax <- max(Mod(complex(real = D$x, imaginary = D$y)))
  els <- tile_elements(G, depth)
  gens <- list(a1 = G$A1, b1 = G$B1, a2 = G$A2, b2 = G$B2, c = comm(G$A1, G$B1))
  imap_dfr(gens, \(X, nm) {
    e <- fixed_points(psl(X))
    ends <- map(els, \(g) mobius(g, e))
    key <- map_chr(ends, \(z) paste(sort(round(c(Arg(z[1]), Arg(z[2])), 6)), collapse = ","))
    ends <- ends[!duplicated(key)]
    # Skip axes that never come within Euclidean radius rmax of 0: for ideal
    # endpoints an angle phi apart, the geodesic's closest approach to 0 is
    # (1 - sin(phi/2)) / cos(phi/2).
    phi <- map_dbl(ends, \(z) abs(.wrap(Arg(z[2]) - Arg(z[1]))))
    near <- (1 - sin(phi / 2)) / cos(phi / 2) < rmax
    ends <- ends[near]
    map_dfr(seq_along(ends), \(k) {
      pts <- .disk_geodesic_1(ends[[k]][1], ends[[k]][2], TRUE, n)
      keep <- inside_polygon(Re(pts), Im(pts), D$x, D$y) & Mod(pts) < 0.9999
      if (!any(keep)) return(NULL)
      run <- cumsum(c(TRUE, diff(keep) != 0))
      tibble(z = pts, run = run)[keep, ] |>
        mutate(x = Re(z), y = Im(z), piece = paste(nm, k, run), curve = nm) |>
        select(-z, -run)
    })
  })
}

geodesic_piece_lengths <- function(cg) {
  cg |>
    group_by(curve, piece) |>
    summarise(len = sum(hdist(complex(real = head(x, -1), imaginary = head(y, -1)),
                              complex(real = x[-1], imaginary = y[-1]), "disk")),
              .groups = "drop") |>
    group_by(curve) |>
    summarise(length_on_D = sum(len), pieces = n(), .groups = "drop")
}

geodesic_labels_df <- function(cg) {
  cg |>
    group_by(curve, piece) |>
    mutate(npts = n()) |>
    ungroup() |>
    group_by(curve) |>
    filter(npts == max(npts)) |>
    slice(ceiling(n() / 2)) |>
    ungroup() |>
    mutate(label = curve)
}

# --- Twists -------------------------------------------------------------------------

curve_choices <- c(
  "a1  (b1 → b1·a1ᵗ)" = "a1",
  "b1  (a1 → a1·b1ᵗ)" = "b1",
  "a2  (b2 → b2·a2ᵗ)" = "a2",
  "b2  (a2 → a2·b2ᵗ)" = "b2",
  "c = [a1,b1], separating" = "c"
)

# Twist of power t about a curve. Integer t is a Dehn twist (change of marking,
# same group); fractional t is a Fenchel-Nielsen deformation (new group).
# Each substitution preserves [A1,B1][A2,B2] = 1.
twist <- function(G, curve, t) {
  if (t == 0) return(G)
  switch(curve,
    a1 = { G$B1 <- psl(G$B1 %*% mpow(G$A1, t)) },
    b1 = { G$A1 <- psl(G$A1 %*% mpow(G$B1, t)) },
    a2 = { G$B2 <- psl(G$B2 %*% mpow(G$A2, t)) },
    b2 = { G$A2 <- psl(G$A2 %*% mpow(G$B2, t)) },
    c  = {
      C  <- comm(G$A1, G$B1)
      Ct <- mpow(C, t)
      G$A2 <- psl(Ct %*% G$A2 %*% minv(Ct))
      G$B2 <- psl(Ct %*% G$B2 %*% minv(Ct))
    },
    stop("unknown curve ", curve)
  )
  G
}

# Parse "a1 b1^-1 c^2" into a list of (curve, integer power), applied left to right.
parse_mapping_class <- function(text) {
  toks <- str_split(str_trim(text), "[\\s,]+")[[1]]
  toks <- toks[toks != ""]
  map(toks, \(tok) {
    m <- str_match(tok, "^(a1|b1|a2|b2|c)(\\^(-?\\d+))?$")
    if (is.na(m[1, 1])) stop("can't read '", tok, "': use a1, b1, a2, b2 or c, optionally with ^n")
    list(curve = m[1, 2], power = if (is.na(m[1, 4])) 1L else as.integer(m[1, 4]))
  })
}

apply_mapping_class <- function(G, steps) {
  reduce(steps, \(H, s) twist(H, s$curve, s$power), .init = G)
}

relator_error <- function(G) {
  R <- comm(G$A1, G$B1) %*% comm(G$A2, G$B2)
  max(Mod(psl(R) - diag(2)))
}

length_table <- function(G) {
  w <- list(
    a1 = G$A1, b1 = G$B1, a2 = G$A2, b2 = G$B2,
    `a1 b1` = G$A1 %*% G$B1, `a2 b2` = G$A2 %*% G$B2, `a1 a2` = G$A1 %*% G$A2
  )
  tibble(curve = names(w), length = map_dbl(w, trace_len))
}

polygon_status <- function(z) {
  convex <- tryCatch(is_convex_polygon(z, "disk"), error = \(e) FALSE)
  area <- if (convex) ngon_area(z, "disk") else NA_real_
  list(convex = convex, area = area)
}


# --- Vertex trajectories and lunes ------------------------------------------------
#
# Under a twist of amount t along one of a1, b1, a2, b2, the marking polygon's
# vertices x_k = W_k p either stay put or move as
#
#   x_k(t) = U X^t q            (U, q fixed; X the generator of the twisted curve).
#
# X^t is translation by t * l(X) along the axis of X, so X^t q runs along the
# EQUIDISTANT CURVE through q at distance d = dist(q, axis X), and U carries that
# curve to an equidistant curve of the lift U(axis X). In Fermi coordinates
# ds^2 = drho^2 + l^2 cosh^2(rho) dt^2, so over t in [0, 1] the vertex travels
#
#   l(X) * cosh(d).
#
# This is the earthquake picture: the vertices on one side of the lifts of the
# curve stay put, the others slide along lifts, farther the farther they are.
# (Which side stays put is a normalization; distances, lengths and area are not.)
#
# slide_spec() lists (U, X, q) for each moving vertex, read off from the relator
# A1 B1 A1^-1 B1^-1 A2 B2 A2^-1 B2^-1 with the substitution made by twist(). For
# example, for a1 (B1 -> B1 A1^t): x2 = A1 B1(t) p = (A1 B1) A1^t p and
# x3 = A1 B1(t) A1^-1 p = (A1 B1) A1^t (A1^-1 p), while x4 = [A1, B1 A1^t] p
# = [A1, B1] p does not move. The separating twist c moves x5, x6, x7 by two
# slides each (conjugation by C^t), so it has no single-slide prediction.

slide_spec <- function(G, curve) {
  p <- P0
  C <- comm(G$A1, G$B1)
  switch(curve,
    a1 = list(`2` = list(U = G$A1 %*% G$B1, X = G$A1, q = p),
              `3` = list(U = G$A1 %*% G$B1, X = G$A1, q = mobius(minv(G$A1), p))),
    b1 = list(`1` = list(U = G$A1, X = G$B1, q = p),
              `2` = list(U = G$A1, X = G$B1, q = mobius(G$B1, p))),
    a2 = list(`6` = list(U = C %*% G$A2 %*% G$B2, X = G$A2, q = p),
              `7` = list(U = C %*% G$A2 %*% G$B2, X = G$A2, q = mobius(minv(G$A2), p))),
    b2 = list(`5` = list(U = C %*% G$A2, X = G$B2, q = p),
              `6` = list(U = C %*% G$A2, X = G$B2, q = mobius(G$B2, p))),
    c  = list()
  )
}

# The generator whose axis is the twisted curve.
twist_generator <- function(G, curve)
  switch(curve, a1 = G$A1, b1 = G$B1, a2 = G$A2, b2 = G$B2, c = comm(G$A1, G$B1))

# Hyperbolic distance from a disk point z to the axis of the hyperbolic M.
# Move z to 0; the axis then has ideal endpoints an angle phi apart (as seen
# from 0), and its closest approach to 0 is at Euclidean radius
# (1 - sin(phi/2)) / cos(phi/2).
dist_to_axis <- function(M, z) {
  e <- fixed_points(psl(M))
  w <- (e - z) / (1 - Conj(z) * e)
  phi <- abs(.wrap(Arg(w[2]) - Arg(w[1])))
  2 * atanh((1 - sin(phi / 2)) / cos(phi / 2))
}

# Positions of all 8 vertices for t in [0, 1] (rows = t, columns = x0..x7).
# n is odd so that every other sample is itself an even sampling (used for the
# length extrapolation below).
vertex_paths <- function(G, curve, n = 1201) {
  ts <- seq(0, 1, length.out = n)
  V <- t(vapply(ts, \(t) marking_polygon(twist(G, curve, t), P0), complex(8)))
  # a vertex "moves" if it ever gets more than 1e-6 from where it started
  # (rounding noise in the fractional matrix powers is far smaller)
  moved <- apply(Mod(V - V[rep(1, n), , drop = FALSE]), 2, max)
  list(t = ts, V = V, moving = which(moved > 1e-6))
}

trajectories_df <- function(vp, every = 5) {
  keep <- unique(c(seq(1, length(vp$t), by = every), length(vp$t)))
  map_dfr(vp$moving, \(j) tibble(x = Re(vp$V[keep, j]), y = Im(vp$V[keep, j]), vertex = j - 1))
}

# One row per moving vertex: measured path length over t in [0, 1], the
# single-slide prediction l(X) cosh d, d, and where the vertex lands at t = 1.
trajectory_table <- function(G, curve, vp, depth = 2) {
  spec <- slide_spec(G, curve)
  poly0 <- marking_polygon(G, P0)
  els <- tile_elements(G, depth)
  depths <- c(0, rep(1:depth, 8 * 7^(0:(depth - 1))))[seq_along(els)]
  # Path length from the sampled points. Each step is measured by the geodesic
  # chord, which undercounts a curved path by O(step^2), so we extrapolate from
  # the full sampling and every other sample (Richardson): L ~ L_h + (L_h - L_2h)/3.
  # This matters for the separating twist c, whose vertices run far out toward
  # the boundary circle.
  chord_len <- \(z) sum(hdist(head(z, -1), z[-1], "disk"))
  map_dfr(vp$moving, \(j) {
    z <- vp$V[, j]
    L_h <- chord_len(z); L_2h <- chord_len(z[seq(1, length(z), by = 2)])
    meas <- L_h + (L_h - L_2h) / 3
    sp <- spec[[as.character(j - 1)]]
    if (!is.null(sp)) {
      d <- dist_to_axis(sp$X, sp$q)
      pred <- trace_len(sp$X) * cosh(d)
    } else { d <- NA_real_; pred <- NA_real_ }
    # landing point at t = 1: which vertex of which tile of the original tiling?
    end <- z[length(z)]
    hit <- NULL
    for (i in seq_along(els)) {
      k <- which(Mod(mobius(els[[i]], poly0) - end) < 1e-6)
      if (length(k)) { hit <- c(i, k[1]); break }
    }
    tibble(
      vertex = sprintf("x%d", j - 1),
      # single slides: accurate to ~1e-4. Two slides (c): the path runs out
      # toward the boundary circle, and the estimate is good to a few percent.
      `path length` = if (is.null(sp)) sprintf("≈ %.0f", meas) else sprintf("%.4f", meas),
      `ℓ·cosh d` = if (is.na(pred)) "two slides" else sprintf("%.4f", pred),
      d = if (is.na(d)) "" else sprintf("%.4f", d),
      `lands on (t = 1)` = if (is.null(hit)) sprintf("beyond %d crossings", depth) else
        if (hit[1] == 1) sprintf("x%d (same tile)", hit[2] - 1) else
        sprintf("x%d of a tile %d crossing%s away", hit[2] - 1, depths[hit[1]],
                if (depths[hit[1]] == 1) "" else "s")
    )
  })
}

# The lune {z : dist(z, axis Y) < d} around the axis of a hyperbolic Y.
# Its two boundary curves are the equidistant curves {Y^s q_+} and {Y^s q_-},
# where q_+- are the points at distance d on either side of the axis, along
# the perpendicular from the axis's closest point m to the origin.
lune_poly <- function(Y, d, n = 200) {
  Y <- psl(Y)
  e <- fixed_points(Y)
  phi <- abs(.wrap(Arg(e[2]) - Arg(e[1])))
  if (abs(phi - pi) < 1e-9) {                 # axis is a diameter: m = 0
    u <- 1i * e[1]; rho_m <- 0
  } else {
    u <- (e[1] + e[2]) / Mod(e[1] + e[2])     # direction of the closest point m
    rho_m <- 2 * atanh((1 - sin(phi / 2)) / cos(phi / 2))
  }
  q_plus  <- u * tanh((rho_m + d) / 2)
  q_minus <- u * tanh((rho_m - d) / 2)
  S <- 9 / trace_len(Y)                       # run about 9 units each way along the axis
  ss <- seq(-S, S, length.out = n)
  side <- \(q) vapply(ss, \(s) mobius(mpow(Y, s), q), complex(1))
  axis <- vapply(ss, \(s) mobius(mpow(Y, s), u * tanh(rho_m / 2)), complex(1))
  list(poly = tibble(z = c(side(q_plus), rev(side(q_minus)))),
       axis = tibble(z = axis))
}

# The lifts to draw lunes around: the axis of the twisted generator X, and the
# lifts U(axis X) along which the moving vertices slide.
lune_axes <- function(G, curve) {
  X <- twist_generator(G, curve)
  ys <- c(list(X), map(slide_spec(G, curve), \(sp) psl(sp$U %*% sp$X %*% minv(sp$U))))
  keys <- map_chr(ys, \(Y) paste(round(sort(Arg(fixed_points(psl(Y)))), 6), collapse = ","))
  ys[!duplicated(keys)]
}

# --- Precomputed reference --------------------------------------------------------

REG <- regular_marking()
P0 <- REG$p

# --- UI -------------------------------------------------------------------------------

ui <- fluidPage(
  titlePanel("Markings and twists on a genus-2 surface"),
  sidebarLayout(
    sidebarPanel(
      width = 3,
      h4("1. Deform the surface"),
      selectInput("curve", "Twist along", choices = curve_choices, selected = "a1"),
      sliderInput("t", "Twist amount t", min = 0, max = 1, value = 0, step = 0.02,
                  animate = animationOptions(interval = 150)),
      helpText("Fractional t gives a different surface. At t = 1 the group is",
               "the original one again, with the marking changed by a Dehn twist."),
      checkboxInput("show_traj", "Trace the vertex trajectories for t from 0 to 1", FALSE),
      checkboxInput("show_lune", "Show lunes around the twisted geodesic's lifts", FALSE),
      conditionalPanel("input.show_lune",
        sliderInput("lune_d", "Lune half-width d", min = 0.05, max = 3, value = 1, step = 0.05)),
      hr(),
      h4("2. Change the marking"),
      textInput("mc", "Dehn twists (applied left to right)", value = "a1"),
      helpText("Examples: a1   b1^-1   a1 b1   c^2. Each token twists the",
               "current generators; the group, and so the orbit, is unchanged."),
      uiOutput("mc_msg"),
      hr(),
      h4("Display"),
      sliderInput("depth", "Tiling depth", min = 1, max = 3, value = 2, step = 1),
      checkboxInput("show_ref", "Original regular tiling (light)", TRUE),
      checkboxInput("show_tiles", "Tiling by the green polygon (grey)", TRUE),
      checkboxInput("show_newtiles", "Tiling by the orange polygon (light orange)", TRUE),
      checkboxInput("show_orbit", "Orbit of the basepoint", FALSE),
      checkboxInput("show_mark", "Polygon of the original marking (green)", TRUE),
      checkboxInput("show_new", "Polygon of the new marking (orange)", TRUE),
      checkboxInput("show_labels", "Label the curves a1, b1, a2, b2, c of the green marking", TRUE),
      checkboxInput("show_axes", "Closed geodesics a1, b1, a2, b2, c, drawn on the green polygon", FALSE)
    ),
    mainPanel(
      width = 9,
      fluidRow(
        column(8, plotOutput("disk", height = "640px")),
        column(4,
               h4("Polygons"),
               tableOutput("status"),
               h4("Marked lengths"),
               tableOutput("lengths"),
               conditionalPanel("input.show_traj",
                 h4("Vertex trajectories (t from 0 to 1)"),
                 tableOutput("traj"),
                 helpText("A vertex that slides along one lift of the twisted",
                          "geodesic travels ℓ·cosh d, where d is its distance from",
                          "that lift (Fermi coordinates). At t = 1 it lands on a",
                          "vertex of the original tiling: the Dehn twist acting on",
                          "the orbit.")),
               helpText("Green = current marking, orange = after the Dehn twists.",
                        "The lengths of individual curves get shuffled, but both",
                        "markings describe the same surface.",
                        "Grey lines tile the disk by copies of the green polygon,",
                        "light orange lines by copies of the orange one. Both",
                        "tilings use the same group, so their corners are the same",
                        "orbit points. If the orange polygon is not convex, its",
                        "copies can overlap and the picture is no longer a tiling."))
      )
    )
  )
)

# --- Server --------------------------------------------------------------------------

server <- function(input, output, session) {

  ref_tiles <- reactive({
    tiling_paths(REG$gens, P0, input$depth)
  })

  G_now <- reactive({
    twist(REG$gens, input$curve, input$t)
  })

  steps <- reactive({
    tryCatch(parse_mapping_class(input$mc), error = \(e) e)
  })

  output$mc_msg <- renderUI({
    s <- steps()
    if (inherits(s, "error")) tags$p(style = "color:#A32D2D;", conditionMessage(s))
    else if (length(s) == 0) tags$p(style = "color:#5F5E5A;", "No twists: the orange polygon is the green one.")
    else if (relator_error(G_new()) > 1e-3)
      tags$p(style = "color:#854F0B;", "The new generators are very long, so their matrices are",
             "huge and floating-point precision is running out (see 'relator error').",
             "The picture near the boundary may be unreliable.")
    else NULL
  })

  G_new <- reactive({
    s <- steps()
    if (inherits(s, "error")) G_now() else apply_mapping_class(G_now(), s)
  })

  output$disk <- renderPlot({
    G <- G_now(); H <- G_new()
    poly_now <- marking_polygon(G, P0)
    poly_new <- marking_polygon(H, P0)

    p <- hyp_canvas("disk")
    if (isTRUE(input$show_ref)) {
      p <- p + geom_path(data = ref_tiles(), aes(x, y, group = tile),
                         colour = col_ref, linewidth = 0.35)
    }
    if (isTRUE(input$show_tiles)) {
      p <- p + geom_path(data = tiling_paths(G, P0, input$depth), aes(x, y, group = tile),
                         colour = col_tile, linewidth = 0.35)
    }
    # The new tiling: the images of the orange polygon under the (same) group.
    if (isTRUE(input$show_newtiles)) {
      p <- p + geom_path(data = tiling_paths(H, P0, input$depth), aes(x, y, group = tile),
                         colour = col_newtile, linewidth = 0.45)
    }
    if (isTRUE(input$show_lune)) {
      col_l <- col_curve[[input$curve]]
      for (Y in lune_axes(REG$gens, input$curve)) {
        lp <- lune_poly(Y, input$lune_d)
        p <- p + geom_polygon(data = lp$poly, aes(Re(z), Im(z)), fill = col_l, alpha = 0.13,
                              colour = col_l, linewidth = 0.3, inherit.aes = FALSE) +
          geom_path(data = lp$axis, aes(Re(z), Im(z)), colour = col_l, linewidth = 0.7,
                    linetype = "22", inherit.aes = FALSE)
      }
    }
    if (isTRUE(input$show_orbit)) {
      p <- p + geom_hpoints(orbit_points(G, P0, input$depth), "disk",
                            size = 0.9, colour = "grey30")
    }
    if (isTRUE(input$show_mark)) {
      p <- p + geom_hngon(poly_now, "disk", fill = col_mark, alpha = 0.25,
                          colour = col_mark, linewidth = 0.9)
    }
    if (isTRUE(input$show_new)) {
      p <- p + geom_hngon(poly_new, "disk", fill = col_new, alpha = 0.2,
                          colour = col_new, linewidth = 0.9) +
        geom_hpoints(poly_new, "disk", colour = col_new, size = 1.8)
    }
    if (isTRUE(input$show_axes)) {
      cg <- closed_geodesics_df(G, P0)
      p <- p + geom_path(data = cg, aes(x, y, group = piece, colour = curve),
                         linewidth = 0.8, linetype = "22", inherit.aes = FALSE) +
        geom_label(data = geodesic_labels_df(cg), aes(x, y, label = label, colour = curve),
                   size = 3.6, fill = "white", alpha = 0.85, inherit.aes = FALSE)
    }
    if (isTRUE(input$show_labels)) {
      ms <- marking_sides_df(G, P0)
      p <- p + geom_path(data = ms$paths, aes(x, y, group = side, colour = curve),
                         linewidth = 1.3, inherit.aes = FALSE) +
        geom_path(data = ms$c_path, aes(x, y), colour = col_curve[["c"]],
                  linewidth = 0.9, linetype = "42", inherit.aes = FALSE) +
        geom_text(data = ms$labels, aes(x, y, label = label, colour = curve),
                  size = 4.6, fontface = "bold", inherit.aes = FALSE) +
        geom_text(data = ms$c_label, aes(x, y, label = label),
                  colour = col_curve[["c"]], size = 4.6, fontface = "bold",
                  inherit.aes = FALSE)
    }
    if (isTRUE(input$show_mark)) {
      p <- p + geom_hpoints(poly_now, "disk", colour = col_mark, size = 1.8)
    }
    if (isTRUE(input$show_traj)) {
      tr <- trajectories_df(vpaths())
      p <- p + geom_path(data = tr, aes(x, y, group = vertex), colour = col_traj,
                         linewidth = 1, arrow = arrow(length = unit(0.08, "inches")),
                         inherit.aes = FALSE)
    }
    p + geom_hpoints(P0, "disk", colour = "black", size = 2.6) +
      scale_colour_manual(values = col_curve, guide = "none")
  }, res = 96)

  # Vertex paths for t in [0, 1], for the chosen curve (independent of the t slider).
  vpaths <- reactive(vertex_paths(REG$gens, input$curve))

  output$traj <- renderTable(trajectory_table(REG$gens, input$curve, vpaths()),
                             striped = TRUE, spacing = "s")

  output$status <- renderTable({
    G <- G_now(); H <- G_new()
    s1 <- polygon_status(marking_polygon(G, P0))
    s2 <- polygon_status(marking_polygon(H, P0))
    fmt_area <- \(s) if (s$convex) sprintf("%.4f", s$area) else "n/a"
    tibble(
      ` ` = c("convex", "area (4π = 12.566)", "relator error"),
      green = c(if (s1$convex) "yes" else "no", fmt_area(s1), sprintf("%.1e", relator_error(G))),
      orange = c(if (s2$convex) "yes" else "no", fmt_area(s2), sprintf("%.1e", relator_error(H)))
    )
  }, striped = TRUE, spacing = "s")

  output$lengths <- renderTable({
    left_join(length_table(G_now()) |> rename(green = length),
              length_table(G_new()) |> rename(orange = length), by = "curve")
  }, digits = 4, striped = TRUE, spacing = "s")
}

shinyApp(ui, server)

## file: hyperbolic.R
# hyperbolic.R --------------------------------------------------------------
#
# Tools for computing and plotting in the hyperbolic plane, in two models:
#
#   model = "disk"  Poincare disk      { w : |w| < 1 }
#   model = "uhp"   upper half-plane   { z : Im z > 0 }
#
# Conventions
#   * Points are R complex numbers: hpt(0.3, 0.2) or 0.3 + 0.2i.
#   * Ideal (boundary) points: in the disk, complex numbers of modulus 1;
#     in the UHP, real numbers (imaginary part 0) or HINF (the point at infinity).
#   * Every function takes a `model` argument; points are given and returned in
#     that model's coordinates.
#   * Isometries are stored internally as SL(2,R) matrices acting on the UHP and
#     can be applied in either model.
#   * Curve-producing functions return tibbles with columns x, y (and an id or
#     edge column), ready for ggplot2.
#
# Usage:  library(tidyverse); source("hyperbolic.R")
# ---------------------------------------------------------------------------

suppressPackageStartupMessages(library(tidyverse))

# --- Basics ----------------------------------------------------------------

HINF <- complex(real = Inf, imaginary = 0) # the ideal point at infinity (UHP)

hpt <- function(x, y = 0) complex(real = x, imaginary = y)

is_hinf <- function(z) is.infinite(Re(z)) | is.infinite(Im(z))

as_xy <- function(z, ..., .before = NULL) {
  out <- tibble(x = Re(z), y = Im(z), ...)
  if (!is.null(.before)) out <- relocate(out, -c(x, y))   # extra columns first
  out
}

.model <- function(model) match.arg(model, c("disk", "uhp"))

.mat <- function(a, b, c, d) matrix(as.complex(c(a, c, b, d)), 2, 2)

.normalize <- function(M) M / sqrt(M[1, 1] * M[2, 2] - M[1, 2] * M[2, 1])

.wrap <- function(t) (t + pi) %% (2 * pi) - pi # angle into [-pi, pi)

# Apply the Mobius map z -> (a z + b) / (c z + d), handling the point at infinity.
mobius <- function(M, z) {
  a <- M[1, 1]
  b <- M[1, 2]
  c <- M[2, 1]
  d <- M[2, 2]
  inf <- is_hinf(z)
  den <- c * z + d
  out <- (a * z + b) / den
  out[!inf & Mod(den) < 1e-14] <- HINF
  if (any(inf)) {
    out[inf] <- if (Mod(c) < 1e-14) HINF else a / c
  }
  out
}

# --- Moving between the models (Cayley transform) ---------------------------

# UHP -> disk:  w = (z - i) / (z + i)
to_disk <- function(z) {
  w <- (z - 1i) / (z + 1i)
  w[is_hinf(z)] <- 1 + 0i
  w
}

# disk -> UHP:  z = i (1 + w) / (1 - w)
to_uhp <- function(w) {
  z <- 1i * (1 + w) / (1 - w)
  ideal <- abs(Mod(w) - 1) < 1e-12
  z[ideal] <- complex(real = Re(z[ideal]), imaginary = 0)
  z[Mod(w - 1) < 1e-12] <- HINF
  z
}

.to_uhp_from <- function(z, model) if (.model(model) == "disk") to_uhp(z) else z
.from_uhp_to <- function(z, model) {
  if (.model(model) == "disk") to_disk(z) else z
}

# Default "center" of each model.
horigin <- function(model = c("disk", "uhp")) {
  if (.model(model) == "disk") 0 + 0i else 1i
}

# --- Distance ----------------------------------------------------------------

# Hyperbolic distance (curvature -1). Vectorized over z1, z2.
hdist <- function(z1, z2, model = c("disk", "uhp")) {
  if (.model(model) == "uhp") {
    acosh(1 + Mod(z1 - z2)^2 / (2 * Im(z1) * Im(z2)))
  } else {
    acosh(1 + 2 * Mod(z1 - z2)^2 / ((1 - Mod(z1)^2) * (1 - Mod(z2)^2)))
  }
}

# --- Geodesics ---------------------------------------------------------------

# Ideal endpoints of the geodesic through z1 and z2, ordered so that moving
# from z1 toward z2 heads toward `forward`.
hgeodesic_endpoints <- function(z1, z2, model = c("disk", "uhp")) {
  model <- .model(model)
  p <- .to_uhp_from(z1, model)
  q <- .to_uhp_from(z2, model)
  if (is_hinf(p) || is_hinf(q) || abs(Re(p) - Re(q)) < 1e-12 * (1 + Mod(p))) {
    x0 <- if (is_hinf(p)) Re(q) else Re(p)
    up <- is_hinf(q) || (!is_hinf(p) && Im(q) > Im(p))
    ends <- if (up) c(hpt(x0), HINF) else c(HINF, hpt(x0))
  } else {
    c0 <- (Mod(q)^2 - Mod(p)^2) / (2 * (Re(q) - Re(p)))
    R <- Mod(p - c0)
    left <- abs(Arg(q - c0)) > abs(Arg(p - c0)) # angle increasing -> toward c0 - R
    ends <- if (left) hpt(c(c0 + R, c0 - R)) else hpt(c(c0 - R, c0 + R))
  }
  ends <- .from_uhp_to(ends, model)
  c(back = ends[1], forward = ends[2])
}

.uhp_geodesic_1 <- function(p, q, full, n, ymax) {
  vertical <- is_hinf(p) ||
    is_hinf(q) ||
    abs(Re(p) - Re(q)) < 1e-12 * (1 + Mod(p))
  if (vertical) {
    x0 <- if (is_hinf(p)) Re(q) else Re(p)
    yp <- if (is_hinf(p)) ymax else Im(p)
    yq <- if (is_hinf(q)) ymax else Im(q)
    y <- if (full) seq(0, ymax, length.out = n) else seq(yp, yq, length.out = n)
    return(hpt(rep(x0, n), y))
  }
  c0 <- (Mod(q)^2 - Mod(p)^2) / (2 * (Re(q) - Re(p)))
  R <- Mod(p - c0)
  th <- if (full) {
    seq(0, pi, length.out = n)
  } else {
    seq(abs(Arg(p - c0)), abs(Arg(q - c0)), length.out = n)
  }
  c0 + R * exp(1i * th)
}

.disk_geodesic_1 <- function(a, b, full, n) {
  det <- Re(a) * Im(b) - Im(a) * Re(b)
  if (abs(det) < 1e-12) {
    # a, b, 0 collinear: a diameter
    if (!full) {
      return(a + seq(0, 1, length.out = n) * (b - a))
    }
    u <- if (Mod(a) >= Mod(b)) a / Mod(a) else b / Mod(b)
    return(seq(-1, 1, length.out = n) * u)
  }
  # circle orthogonal to the unit circle: 2 Re(c conj(p)) = 1 + |p|^2 for p = a, b
  A <- matrix(c(Re(a), Re(b), Im(a), Im(b)), 2)
  cc <- solve(A, c(1 + Mod(a)^2, 1 + Mod(b)^2) / 2)
  c0 <- hpt(cc[1], cc[2])
  R <- sqrt(Mod(c0)^2 - 1)
  if (full) {
    beta <- acos(1 / Mod(c0))
    a <- exp(1i * (Arg(c0) - beta))
    b <- exp(1i * (Arg(c0) + beta))
  }
  t1 <- Arg(a - c0)
  th <- t1 + seq(0, 1, length.out = n) * .wrap(Arg(b - c0) - t1)
  c0 + R * exp(1i * th)
}

# Points along geodesics, as a tibble (id, x, y). Vectorized over pairs (z1, z2).
#   full = FALSE: the segment from z1 to z2.  full = TRUE: the whole geodesic.
#   ymax: where vertical UHP geodesics are cut off.
hgeodesic <- function(
  z1,
  z2,
  model = c("disk", "uhp"),
  full = FALSE,
  n = 200,
  ymax = 100
) {
  model <- .model(model)
  stopifnot(length(z1) == length(z2))
  map_dfr(seq_along(z1), function(k) {
    pts <- if (model == "uhp") {
      .uhp_geodesic_1(z1[k], z2[k], full, n, ymax)
    } else {
      .disk_geodesic_1(z1[k], z2[k], full, n)
    }
    as_xy(pts, id = k, .before = 1)
  })
}

# --- Hyperbolic circles -------------------------------------------------------

# Circle of hyperbolic radius r about `center`. Vectorized; returns (id, x, y).
hcircle <- function(center, radius, model = c("disk", "uhp"), n = 200) {
  model <- .model(model)
  k <- vctrs::vec_recycle_common(center, radius)
  th <- seq(0, 2 * pi, length.out = n)
  map_dfr(seq_along(k[[1]]), function(j) {
    p <- k[[1]][j]
    r <- k[[2]][j]
    pts <- if (model == "uhp") {
      hpt(Re(p), Im(p) * cosh(r)) + Im(p) * sinh(r) * exp(1i * th)
    } else {
      w <- tanh(r / 2) * exp(1i * th)
      (w + p) / (1 + Conj(p) * w)
    }
    as_xy(pts, id = j, .before = 1)
  })
}

# --- Isometries ---------------------------------------------------------------
#
# An isometry is a list of class "hiso" holding M, a 2x2 matrix with det 1
# acting on the UHP by z -> (a z + b) / (c z + d).

.hiso <- function(M, label) {
  structure(list(M = .normalize(M), label = label), class = "hiso")
}

iso_identity <- function() .hiso(diag(2) + 0i, "identity")

# Elliptic: rotation by `angle` (radians, counterclockwise) about `center`.
iso_rotation <- function(center, angle, model = c("disk", "uhp")) {
  p <- .to_uhp_from(center, model)
  x <- Re(p)
  y <- Im(p)
  A <- .mat(sqrt(y), x / sqrt(y), 0, 1 / sqrt(y)) # i -> p
  R <- .mat(cos(angle / 2), sin(angle / 2), -sin(angle / 2), cos(angle / 2))
  .hiso(A %*% R %*% solve(A), sprintf("rotation by %.4g about a point", angle))
}

# Parabolic: fixes the ideal point `fixed`. Conjugate to z -> z + t with `fixed`
# sent to infinity; t controls how far points move, and its sign the direction.
iso_parabolic <- function(fixed, t, model = c("disk", "uhp")) {
  xi <- .to_uhp_from(fixed, model)
  T <- .mat(1, t, 0, 1)
  if (is_hinf(xi)) {
    return(.hiso(T, "parabolic fixing infinity"))
  }
  B <- .mat(Re(xi), -1, 1, 0) # infinity -> xi
  .hiso(B %*% T %*% solve(B), "parabolic")
}

# Hyperbolic: translation by distance d along the geodesic with ideal endpoints
# `from` and `to`, moving points toward `to`.
iso_hyperbolic <- function(from, to, d, model = c("disk", "uhp")) {
  a <- .to_uhp_from(from, model)
  b <- .to_uhp_from(to, model)
  C <- if (is_hinf(b)) {
    # 0 -> a, inf -> b
    .mat(1, Re(a), 0, 1)
  } else if (is_hinf(a)) {
    .mat(Re(b), -1, 1, 0)
  } else {
    .mat(Re(b), Re(a), 1, 1)
  }
  T <- .mat(exp(d / 2), 0, 0, exp(-d / 2)) # z -> e^d z
  .hiso(C %*% T %*% solve(C), sprintf("translation by %.4g", d))
}

# Hyperbolic translation taking point p to point q (along the geodesic pq).
iso_translate <- function(p, q, model = c("disk", "uhp")) {
  model <- .model(model)
  if (Mod(p - q) < 1e-14) {
    return(iso_identity())
  }
  e <- hgeodesic_endpoints(p, q, model)
  iso_hyperbolic(e[["back"]], e[["forward"]], hdist(p, q, model), model)
}

# The point(s) at (signed) distance d from p along the geodesic from p toward q.
hpoint_toward <- function(p, q, d, model = c("disk", "uhp")) {
  e <- hgeodesic_endpoints(p, q, model)
  map_vec(d, \(t) {
    iso_apply(p, iso_hyperbolic(e[["back"]], e[["forward"]], t, model), model)
  })
}

iso_compose <- function(f, g) .hiso(f$M %*% g$M, "composition") # f after g
iso_inverse <- function(f) .hiso(solve(f$M), "inverse")

# Matrix of the isometry in a given model's coordinates.
iso_matrix <- function(iso, model = c("disk", "uhp")) {
  if (.model(model) == "uhp") {
    return(iso$M)
  }
  K <- .mat(1, -1i, 1, 1i) # Cayley map
  .normalize(K %*% iso$M %*% solve(K))
}

iso_apply <- function(z, iso, model = c("disk", "uhp")) {
  mobius(iso_matrix(iso, model), z)
}

# Apply an isometry to the x, y columns of a tibble (e.g. a curve).
iso_apply_xy <- function(df, iso, model = c("disk", "uhp")) {
  w <- iso_apply(hpt(df$x, df$y), iso, model)
  mutate(df, x = Re(w), y = Im(w))
}

iso_type <- function(iso, tol = 1e-9) {
  M <- iso$M
  if (max(Mod(M - diag(2))) < tol || max(Mod(M + diag(2))) < tol) {
    return("identity")
  }
  tr <- abs(Re(M[1, 1] + M[2, 2]))
  if (tr < 2 - tol) {
    "elliptic"
  } else if (tr > 2 + tol) {
    "hyperbolic"
  } else {
    "parabolic"
  }
}

# Translation length (hyperbolic) or rotation angle (elliptic).
iso_amount <- function(iso) {
  tr <- abs(Re(iso$M[1, 1] + iso$M[2, 2]))
  switch(
    iso_type(iso),
    hyperbolic = c(translation_length = 2 * acosh(tr / 2)),
    elliptic = c(rotation_angle = 2 * acos(tr / 2)),
    c(none = 0)
  )
}

# Fixed points in the closure of the model (interior point for elliptic,
# one ideal point for parabolic, two ideal points for hyperbolic).
iso_fixed_points <- function(iso, model = c("disk", "uhp"), tol = 1e-9) {
  M <- iso$M
  a <- M[1, 1]
  b <- M[1, 2]
  c <- M[2, 1]
  d <- M[2, 2]
  if (Mod(c) < tol) {
    fp <- if (Mod(a - d) < tol) HINF else c(b / (d - a), HINF)
  } else {
    disc <- sqrt((a - d)^2 + 4 * b * c + 0i)
    fp <- c((a - d + disc) / (2 * c), (a - d - disc) / (2 * c))
    fp <- complex(
      real = Re(fp),
      imaginary = ifelse(abs(Im(fp)) < tol, 0, Im(fp))
    )
    fp <- if (iso_type(iso) == "elliptic") fp[Im(fp) > 0] else unique(fp)
  }
  .from_uhp_to(fp, model)
}

print.hiso <- function(x, ...) {
  cat("<hyperbolic isometry>", iso_type(x), "\n")
  amt <- iso_amount(x)
  if (names(amt) != "none") {
    cat(" ", names(amt), "=", signif(amt, 6), "\n")
  }
  fp <- iso_fixed_points(x, "uhp")
  fp <- if (all(Im(fp) == 0)) {
    format(signif(Re(fp), 6))
  } else {
    format(signif(fp, 6))
  }
  cat("  fixed points (uhp):", fp, "\n")
  cat("  matrix (uhp):\n")
  print(signif(Re(x$M), 6))
  invisible(x)
}

# --- Polygons -----------------------------------------------------------------

# Vertices of a regular hyperbolic n-gon. Give either
#   radius: hyperbolic distance from center to each vertex, or
#   angle:  interior angle (must be < (n - 2) pi / n),
# using  cosh(radius) = cot(pi / n) cot(angle / 2).
# rotation turns the polygon about its center; the default (0) puts the first
# vertex straight "up" from the center (toward i in the disk, toward
# infinity in the UHP).
regular_ngon <- function(
  n,
  radius = NULL,
  angle = NULL,
  center = NULL,
  rotation = 0,
  model = c("disk", "uhp")
) {
  model <- .model(model)
  rotation <- rotation + if (model == "disk") pi / 2 else 0
  stopifnot(n >= 3, xor(is.null(radius), is.null(angle)))
  if (!is.null(angle)) {
    if (angle >= (n - 2) * pi / n) {
      stop(
        "a regular hyperbolic ",
        n,
        "-gon needs angle < ",
        signif((n - 2) * pi / n, 5)
      )
    }
    radius <- acosh(1 / (tan(pi / n) * tan(angle / 2)))
  }
  v <- tanh(radius / 2) * exp(1i * (rotation + 2 * pi * (0:(n - 1)) / n))
  v <- .from_uhp_to(to_uhp(v), model)
  if (!is.null(center)) {
    v <- iso_apply(v, iso_translate(horigin(model), center, model), model)
  }
  v
}

# Boundary of the polygon with the given vertices (in order), as (edge, x, y).
ngon_path <- function(vertices, model = c("disk", "uhp"), n = 100) {
  nxt <- c(vertices[-1], vertices[1])
  hgeodesic(vertices, nxt, model, n = n) |> rename(edge = id)
}

# Interior angles, assuming a convex polygon with vertices in order.
ngon_angles <- function(vertices, model = c("disk", "uhp")) {
  v <- if (.model(model) == "uhp") to_disk(vertices) else vertices
  n <- length(v)
  map_dbl(seq_len(n), function(k) {
    p <- v[k]
    phi <- function(w) (w - p) / (1 - Conj(p) * w) # moves p to 0
    prev <- v[(k - 2) %% n + 1]
    nxt <- v[k %% n + 1]
    abs(.wrap(Arg(phi(nxt)) - Arg(phi(prev))))
  })
}

ngon_sides <- function(vertices, model = c("disk", "uhp")) {
  hdist(vertices, c(vertices[-1], vertices[1]), model)
}

# Area by Gauss-Bonnet: (n - 2) pi - (sum of interior angles).
ngon_area <- function(vertices, model = c("disk", "uhp")) {
  (length(vertices) - 2) * pi - sum(ngon_angles(vertices, model))
}

ngon_summary <- function(vertices, model = c("disk", "uhp")) {
  as_xy(vertices, vertex = seq_along(vertices), .before = 1) |>
    mutate(
      angle = ngon_angles(vertices, model),
      side_to_next = ngon_sides(vertices, model)
    )
}

# --- Plotting -----------------------------------------------------------------

# A blank canvas for a model: the unit circle, or the real axis.
hyp_canvas <- function(
  model = c("disk", "uhp"),
  xlim = c(-3, 3),
  ylim = c(0, 4)
) {
  if (.model(model) == "disk") {
    th <- seq(0, 2 * pi, length.out = 400)
    ggplot() +
      geom_polygon(
        data = tibble(x = cos(th), y = sin(th)),
        aes(x, y),
        fill = "grey97",
        colour = "grey35",
        linewidth = 0.6
      ) +
      coord_equal(
        xlim = c(-1.05, 1.05),
        ylim = c(-1.05, 1.05),
        expand = FALSE
      ) +
      theme_void()
  } else {
    ggplot() +
      annotate(
        "rect",
        xmin = xlim[1],
        xmax = xlim[2],
        ymin = 0,
        ymax = ylim[2],
        fill = "grey97"
      ) +
      geom_hline(yintercept = 0, colour = "grey35", linewidth = 0.6) +
      coord_equal(xlim = xlim, ylim = ylim, expand = FALSE) +
      theme_minimal() +
      theme(axis.title = element_blank(), panel.grid.minor = element_blank())
  }
}

# Layers to add to a canvas. `...` goes to the underlying ggplot2 geom.
geom_hgeodesic <- function(
  z1,
  z2,
  model = c("disk", "uhp"),
  full = FALSE,
  n = 200,
  ymax = 100,
  ...
) {
  geom_path(
    data = hgeodesic(z1, z2, model, full, n, ymax),
    aes(x, y, group = id),
    inherit.aes = FALSE,
    ...
  )
}

geom_hcircle <- function(center, radius, model = c("disk", "uhp"), ...) {
  geom_path(
    data = hcircle(center, radius, model),
    aes(x, y, group = id),
    inherit.aes = FALSE,
    ...
  )
}

geom_hngon <- function(vertices, model = c("disk", "uhp"), fill = NA, ...) {
  geom_polygon(
    data = ngon_path(vertices, model),
    aes(x, y),
    inherit.aes = FALSE,
    fill = fill,
    ...
  )
}

geom_hpoints <- function(z, model = c("disk", "uhp"), ...) {
  geom_point(data = as_xy(z[!is_hinf(z)]), aes(x, y), inherit.aes = FALSE, ...)
}

# --- Surface groups -----------------------------------------------------------
#
# A closed orientable surface of genus g >= 2 is the quotient of the regular
# 4g-gon with interior angles 2 pi / (4g) (angle sum 2 pi) by side pairings.
# Side k of a polygon runs from vertex k to vertex k + 1 (indices mod n).
#
# For each side j, gens[[j]] is the side pairing that maps the partner side
# pair[j] onto side j. It sends the polygon D to the tile across side j, and
# gens[[pair[j]]] is its inverse.
#
# Tiles are named by words: crossing sides j1, j2, ..., jk in turn, starting
# from D, reaches the tile  g_{j1} g_{j2} ... g_{jk} (D)  (composition written
# in crossing order, so the last side crossed is applied first). The neighbor of
# the tile h(D) across its side h(side j) is  h g_j (D), not  h g_j h^-1 (D).

# Orientation-preserving isometry sending p1 -> q1 and p2 -> q2.
# Requires d(p1, p2) = d(q1, q2).
iso_segment <- function(p1, p2, q1, q2, model = c("disk", "uhp")) {
  model <- .model(model)
  if (abs(hdist(p1, p2, model) - hdist(q1, q2, model)) > 1e-8)
    stop("iso_segment: the two segments have different lengths")
  f <- iso_translate(p1, q1, model)
  disk <- function(z) if (model == "uhp") to_disk(z) else z
  q <- disk(q1)
  at_q <- function(z) { w <- disk(z); (w - q) / (1 - Conj(q) * w) }   # q1 -> 0
  ang <- Arg(at_q(q2)) - Arg(at_q(iso_apply(p2, f, model)))
  iso_compose(iso_rotation(q1, ang, model), f)
}

polygon_side <- function(vertices, k) {
  n <- length(vertices)
  vertices[c((k - 1) %% n + 1, k %% n + 1)]
}

# Orientation-preserving pairing of side i onto side j, reversing direction
# along the boundary (as an orientable gluing requires): v_i -> v_{j+1},
# v_{i+1} -> v_j.
side_pairing <- function(vertices, i, j, model = c("disk", "uhp")) {
  s <- polygon_side(vertices, i); t <- polygon_side(vertices, j)
  iso_segment(s[1], s[2], t[2], t[1], model)
}

# The regular 4g-gon with its surface-group side pairings.
#   word = "standard": a1 b1 a1^-1 b1^-1 ... ag bg ag^-1 bg^-1
#                      (side 4m+1 <-> 4m+3, side 4m+2 <-> 4m+4)
#   word = "opposite": side k <-> side k + 2g
surface_polygon <- function(genus, model = c("disk", "uhp"), center = NULL,
                            word = c("standard", "opposite")) {
  model <- .model(model); word <- match.arg(word)
  stopifnot(genus >= 2)
  n <- 4 * genus
  v <- regular_ngon(n, angle = 2 * pi / n, center = center, model = model)
  pair <- if (word == "standard") {
    as.vector(sapply(0:(genus - 1), \(m) 4 * m + c(3, 4, 1, 2)))
  } else {
    (seq_len(n) + 2 * genus - 1) %% n + 1
  }
  gens <- map(seq_len(n), \(j) side_pairing(v, pair[j], j, model))
  structure(list(vertices = v, pair = pair, gens = gens, genus = genus,
                 word = word, model = model), class = "surface_polygon")
}

# The group element for a word of side crossings, and the tile it reaches.
tile_iso <- function(P, word) {
  reduce(P$gens[word], iso_compose, .init = iso_identity())
}
tile <- function(P, word) iso_apply(P$vertices, tile_iso(P, word), P$model)

# All tiles reachable by crossing at most `depth` sides (no immediate
# backtracking), as a tibble of boundary paths with a `word` column.
tiles_upto <- function(P, depth = 1, n = 60) {
  words <- list(integer(0))
  frontier <- words
  for (d in seq_len(depth)) {
    frontier <- unlist(lapply(frontier, \(w) lapply(seq_along(P$gens), \(j) {
      if (length(w) && P$pair[j] == w[length(w)]) NULL else c(w, j)
    })), recursive = FALSE)
    frontier <- compact(frontier)
    words <- c(words, frontier)
  }
  map_dfr(words, \(w) ngon_path(tile(P, w), P$model, n = n) |>
            mutate(word = paste(w, collapse = " "), depth = length(w)))
}

print.surface_polygon <- function(x, ...) {
  cat("<surface polygon> genus", x$genus, "|", 4 * x$genus, "-gon |", x$word,
      "pairing |", x$model, "\n")
  cat("  side j is paired with side:", x$pair, "\n")
  cat("  translation lengths:", signif(map_dbl(x$gens, iso_amount), 6), "\n")
  invisible(x)
}

## file: teichmuller.R
# teichmuller.R -------------------------------------------------------------
#
# Fundamental polygons for points of Teichmuller space other than the regular
# one, and tools for comparing the surfaces they give.
#
#   source("hyperbolic.R"); source("teichmuller.R")
#
# Background (Poincare's polygon theorem). A convex hyperbolic polygon with a
# side pairing is a fundamental domain for the group generated by the pairing
# transformations, and the quotient is a closed surface, provided
#   (1) paired sides have equal length, and
#   (2) for every cycle of identified vertices, the interior angles sum to 2 pi.
# For a 4g-gon with the standard word a1 b1 a1^-1 b1^-1 ... (or with opposite
# sides paired) all vertices form one cycle, so (2) is "angle sum = 2 pi".
#
# Counting: 4g vertices = 8g real parameters, minus 3 for isometries, minus
# 2g + 1 equations, leaves 6g - 4 = dim Teich(S_g) + 2. The extra 2 is the
# position of the polygon's vertex on the surface.
#
# Polygons here are worked with in the disk model; the resulting surface
# polygons can be converted with to_uhp() like any others.
# ---------------------------------------------------------------------------

# --- Pairing patterns ---------------------------------------------------------

pairing_pattern <- function(genus, word = c("standard", "opposite")) {
  word <- match.arg(word)
  n <- 4 * genus
  if (word == "standard") as.vector(sapply(0:(genus - 1), \(m) 4 * m + c(3, 4, 1, 2)))
  else (seq_len(n) + 2 * genus - 1) %% n + 1
}

# Cycles of vertices identified by the side pairing. Side i (v_i -> v_{i+1}) is
# glued to side j = pair[i] reversed, so v_i ~ v_{j+1} and v_{i+1} ~ v_j.
vertex_cycles <- function(pair) {
  n <- length(pair)
  nxt <- \(k) k %% n + 1
  parent <- seq_len(n)
  find <- function(k) { while (parent[k] != k) k <- parent[k]; k }
  join <- function(a, b) { ra <- find(a); rb <- find(b); if (ra != rb) parent[ra] <<- rb }
  for (i in seq_len(n)) { j <- pair[i]; join(i, nxt(j)); join(nxt(i), j) }
  roots <- vapply(seq_len(n), find, 1)
  unname(split(seq_len(n), roots))
}

# Interior angles of a polygon with vertices in counterclockwise order,
# in (0, 2 pi): angles > pi mean the polygon is not convex there.
signed_angles <- function(vertices, model = c("disk", "uhp")) {
  v <- if (.model(model) == "uhp") to_disk(vertices) else vertices
  n <- length(v)
  vapply(seq_len(n), function(k) {
    p <- v[k]
    phi <- function(w) (w - p) / (1 - Conj(p) * w)          # p -> 0
    prev <- v[(k - 2) %% n + 1]; nxt <- v[k %% n + 1]
    (Arg(phi(prev)) - Arg(phi(nxt))) %% (2 * pi)
  }, 0)
}

is_convex_polygon <- function(vertices, model = c("disk", "uhp")) {
  all(signed_angles(vertices, model) < pi)
}

# --- From any polygon to a surface polygon -----------------------------------

# Build a surface polygon (same structure as surface_polygon()) from vertices.
# Checks Poincare's conditions and stops with a clear message if they fail.
surface_from_polygon <- function(vertices, model = c("disk", "uhp"),
                                 word = c("standard", "opposite"), tol = 1e-8) {
  model <- .model(model); word <- match.arg(word)
  n <- length(vertices)
  if (n %% 4 != 0 || n < 8) stop("need a 4g-gon with g >= 2")
  genus <- n / 4
  pair <- pairing_pattern(genus, word)
  len <- ngon_sides(vertices, model)
  bad <- which(abs(len - len[pair]) > tol)
  if (length(bad)) stop("paired sides differ in length: sides ", paste(bad, collapse = ", "))
  ang <- signed_angles(vertices, model)
  if (any(ang >= pi)) stop("polygon is not convex")
  for (cyc in vertex_cycles(pair)) {
    if (abs(sum(ang[cyc]) - 2 * pi) > tol)
      stop("angles at vertex cycle {", paste(cyc, collapse = ", "), "} sum to ",
           signif(sum(ang[cyc]), 8), ", not 2 pi")
  }
  gens <- map(seq_len(n), \(j) side_pairing(vertices, pair[j], j, model))
  structure(list(vertices = vertices, pair = pair, gens = gens, genus = genus,
                 word = word, model = model), class = "surface_polygon")
}

# --- Constraint equations and projection -------------------------------------

.to_x   <- function(v) c(Re(v), Im(v))
.from_x <- function(x) { n <- length(x) / 2; complex(real = x[1:n], imaginary = x[(n + 1):(2 * n)]) }

# Residuals of Poincare's conditions (disk coordinates): one equation per pair
# of sides, one per vertex cycle.
poincare_residuals <- function(vertices, pair) {
  len <- ngon_sides(vertices, "disk")
  reps <- which(seq_along(pair) < pair)
  ang <- signed_angles(vertices, "disk")
  cyc <- vapply(vertex_cycles(pair), \(k) sum(ang[k]) - 2 * pi, 0)
  c(len[reps] - len[pair[reps]], cyc)
}

.jacobian <- function(f, x, h = 1e-7) {
  f0 <- f(x)
  vapply(seq_along(x), function(k) {
    e <- replace(numeric(length(x)), k, h)
    (f(x + e) - f(x - e)) / (2 * h)
  }, numeric(length(f0)))
}

# Move vertices (disk) to the nearest polygon satisfying Poincare's conditions,
# by minimum-norm Newton steps.
project_to_poincare <- function(vertices, pair, tol = 1e-13, maxit = 100) {
  f <- \(x) poincare_residuals(.from_x(x), pair)
  x <- .to_x(vertices)
  for (it in seq_len(maxit)) {
    r <- f(x)
    if (max(abs(r)) < tol) break
    J <- .jacobian(f, x)
    step <- drop(t(J) %*% solve(J %*% t(J), r))
    # damp the step if it would leave the disk or make things worse
    for (k in 1:30) {
      x_new <- x - step
      r_new <- f(x_new)
      if (all(Mod(.from_x(x_new)) < 1) && all(is.finite(r_new)) &&
          max(abs(r_new)) < max(abs(r))) break
      step <- step / 2
    }
    x <- x_new
  }
  if (max(abs(f(x))) > 1e-9) stop("projection did not converge")
  .from_x(x)
}

# A new point of Teichmuller space near P: move each vertex v_k to
# (delta_k + v_k) / (1 + conj(v_k) delta_k), a point at hyperbolic distance
# about 2 |delta_k| from v_k, then project back onto Poincare's conditions.
# delta is random with size `scale`, or given as a complex vector.
deform_surface_polygon <- function(P, scale = 0.05, delta = NULL, seed = NULL) {
  v <- if (P$model == "uhp") to_disk(P$vertices) else P$vertices
  if (is.null(delta)) {
    if (!is.null(seed)) set.seed(seed)
    delta <- complex(real = rnorm(length(v)), imaginary = rnorm(length(v))) * scale
  }
  w <- project_to_poincare((delta + v) / (1 + Conj(v) * delta), P$pair)
  if (any(Mod(w) >= 1)) stop("deformation left the disk; try a smaller scale")
  if (!is_convex_polygon(w, "disk")) stop("deformed polygon is not convex; try a smaller scale")
  w <- if (P$model == "uhp") to_uhp(w) else w
  surface_from_polygon(w, P$model, P$word)
}

# --- Polygons from side lengths and angles ("turtle" construction) -----------

# Walk counterclockwise: from each vertex go `lengths[k]` along the current
# direction, then turn left by pi - angles[k + 1]. Starts at the origin of the
# disk heading along the positive real axis. Returns the vertices, with the
# closing error (0 for a closed polygon) as an attribute.
turtle_polygon <- function(lengths, angles, model = c("disk", "uhp")) {
  model <- .model(model)
  n <- length(lengths); stopifnot(length(angles) == n)
  step <- function(s) iso_hyperbolic(-1, 1, s, "disk")           # along the real axis
  turn <- function(a) iso_rotation(0, a, "disk")
  F <- iso_identity(); v <- complex(n)
  for (k in seq_len(n)) {
    v[k] <- iso_apply(0, F, "disk")
    F <- iso_compose(F, iso_compose(step(lengths[k]), turn(pi - angles[k %% n + 1])))
  }
  err <- min(max(Mod(F$M - diag(2))), max(Mod(F$M + diag(2))))
  v <- if (model == "uhp") to_uhp(v) else v
  structure(v, closing_error = err)
}

# Genus 2: the equiangular octagon (all angles pi/4) whose sides alternate
# a, b, a, b, ... Paired sides (1,3), (2,4), (5,7), (6,8) automatically have
# equal length, so this is a one-parameter family of points of Teichmuller
# space through the regular octagon. Given a, b is determined by closing up.
alternating_octagon <- function(a, model = c("disk", "uhp")) {
  model <- .model(model)
  ext <- pi - pi / 4
  X <- function(b) iso_compose(iso_compose(iso_hyperbolic(-1, 1, a, "disk"), iso_rotation(0, ext, "disk")),
                               iso_compose(iso_hyperbolic(-1, 1, b, "disk"), iso_rotation(0, ext, "disk")))
  # Four copies of X must close up, so X is a rotation by pi/2 about the
  # octagon's center: |trace| = sqrt(2). That equation has several roots; some
  # give star octagons (winding 3 times, like the Euclidean {8/3} star, which
  # also has angles pi/4). Keep the simple one: vertices wind once around the
  # center.
  g  <- \(b) abs(Re(sum(diag(X(b)$M)))) - sqrt(2)
  bs <- exp(seq(log(1e-3), log(20), length.out = 400))
  gb <- vapply(bs, g, 0)
  roots <- map_dbl(which(diff(sign(gb)) != 0),
                   \(k) uniroot(g, bs[c(k, k + 1)], tol = 1e-14)$root)
  centered <- function(b) {
    v <- turtle_polygon(rep(c(a, b), 4), rep(pi / 4, 8), "disk")
    ctr <- iso_fixed_points(X(b), "disk")
    w <- iso_apply(v, iso_translate(ctr, 0, "disk"), "disk")
    structure(w, closing_error = attr(v, "closing_error"))
  }
  simple <- function(b) {
    w <- centered(b)
    winding <- sum(.wrap(diff(Arg(c(w, w[1]))))) / (2 * pi)
    attr(w, "closing_error") < 1e-8 && abs(winding - 1) < 1e-6 &&
      all(abs(signed_angles(w, "disk") - pi / 4) < 1e-8)
  }
  ok <- keep(roots, simple)
  if (!length(ok)) stop("no convex equiangular octagon with a = ", a,
                        " (a must be larger than about 1.77, and at most about 12 numerically)")
  b <- ok[1]
  v <- centered(b)
  P <- surface_from_polygon(if (model == "uhp") to_uhp(v) else v, model)
  P$family <- c(a = a, b = b)
  P
}

# --- Comparing surfaces ---------------------------------------------------------

.tr <- function(M) abs(Re(M[1, 1] + M[2, 2]))

# Generators a_m, b_m for the standard word: a_m = g_{4m-3}, b_m = g_{4m}.
standard_generators <- function(P) {
  stopifnot(P$word == "standard")
  idx <- as.vector(rbind(4 * seq_len(P$genus) - 3, 4 * seq_len(P$genus)))
  set_names(P$gens[idx], as.vector(rbind(paste0("a", seq_len(P$genus)),
                                         paste0("b", seq_len(P$genus)))))
}

# Traces of the generators and of all their pairwise products: coordinates
# of the *marked* surface (a point of Teichmuller space). Two marked surfaces
# are the same point iff all these traces agree (for enough words).
marked_traces <- function(P) {
  G <- if (P$word == "standard") standard_generators(P)
       else set_names(P$gens[seq_len(2 * P$genus)], paste0("g", seq_len(2 * P$genus)))
  nm <- names(G)
  singles <- tibble(word = nm, trace = map_dbl(G, \(g) .tr(g$M)))
  pairs <- combn(seq_along(G), 2, simplify = FALSE) |>
    map_dfr(\(ij) tibble(word = paste(nm[ij], collapse = " "),
                         trace = .tr(G[[ij[1]]]$M %*% G[[ij[2]]]$M)))
  bind_rows(singles, pairs)
}

# Translation lengths 2 arccosh(|tr|/2) of all cyclically reduced words in the
# side pairings of length <= max_word. Marking-free: it depends only on the
# surface, so it is a function on moduli space. Returns distinct lengths with
# how many words give each (words, not geodesics: conjugates and inverses are
# counted separately) and one example word.
length_spectrum <- function(P, max_word = 4, digits = 8) {
  n <- length(P$gens); pair <- P$pair
  Ms <- map(P$gens, "M")
  out <- list(); k <- 0
  rec <- function(w, M) {
    if (length(w) && P$pair[w[length(w)]] != w[1]) {        # cyclically reduced
      t <- .tr(M)
      if (t > 2 + 1e-9) { k <<- k + 1; out[[k]] <<- list(w = w, len = 2 * acosh(t / 2)) }
    }
    if (length(w) == max_word) return(invisible())
    for (j in seq_len(n)) {
      if (length(w) && pair[j] == w[length(w)]) next
      rec(c(w, j), M %*% Ms[[j]])
    }
  }
  rec(integer(0), diag(2) + 0i)
  tibble(length = map_dbl(out, "len"),
         word = map_chr(out, \(o) paste0("g", o$w, collapse = " "))) |>
    mutate(length = round(length, digits)) |>
    group_by(length) |>
    summarise(words = n(), example = first(word), .groups = "drop") |>
    arrange(length)
}

systole <- function(P, max_word = 4) min(length_spectrum(P, max_word)$length)

# Relabel: start the side labels at vertex 1 + shift instead of vertex 1.
# Same polygon, same group, different marking (a mapping class acts).
relabel_polygon <- function(P, shift) {
  n <- length(P$vertices)
  v <- P$vertices[(seq_len(n) + shift - 1) %% n + 1]
  surface_from_polygon(v, P$model, P$word)
}

# --- The dimension of Teichmuller space, numerically -----------------------------

# Tangent space to the space of Poincare polygons at P, and the rank of the
# map  polygon -> marked traces  on it. The rank should be 6g - 6.
teichmuller_dimension <- function(P, tol = 1e-6) {
  v <- if (P$model == "uhp") to_disk(P$vertices) else P$vertices
  pair <- P$pair; word <- P$word
  f <- \(x) poincare_residuals(.from_x(x), pair)
  x0 <- .to_x(v)
  J <- .jacobian(f, x0)
  sv <- svd(J, nv = length(x0))
  r_con <- sum(sv$d > tol * max(sv$d))
  N <- sv$v[, (r_con + 1):length(x0), drop = FALSE]          # tangent directions
  trace_vec <- function(x) {
    w <- project_to_poincare(.from_x(x), pair)
    marked_traces(surface_from_polygon(w, "disk", word))$trace
  }
  h <- 1e-5
  T <- vapply(seq_len(ncol(N)), \(k) (trace_vec(x0 + h * N[, k]) - trace_vec(x0 - h * N[, k])) / (2 * h),
              numeric(length(trace_vec(x0))))
  st <- svd(T)$d
  c(parameters = length(x0),
    equations = r_con,
    polygon_space = ncol(N),
    teichmuller = sum(st > 1e-4 * max(st)),
    expected = 6 * P$genus - 6)
}