Skip to content

Theory

Hifuku draws a navigational map of a gene's phylogenetic likelihood landscape over a shared two-dimensional chart of tree space. It is a survey, not an optimizer and not a posterior estimator: the aim is to show where the reasonably-likely trees lie and how they connect, with elevation read from the directly-computed alignment log-likelihood.

The map joins two fields whose vocabularies overlap in a way that can mislead. Phylogenetics models how sequences evolve along a tree. The search that fills the map is a quality-diversity algorithm, MAP-Elites, which grew out of evolutionary computation and therefore borrows biology's own words: it has a population, it applies mutation and selection, it keeps the fittest, and it speaks of genomes, parents, offspring, and lineages. In Hifuku these are terms for a search over candidate trees, not claims about the biological evolution those trees describe. The likelihood the search maximizes is a statistical quantity, not an organism's reproductive success; a "mutation" in the search is an edit to a candidate tree, not a modeled substitution in a sequence. The danger here is the reverse of the usual cross-field problem: the two fields do not use different words for the same idea, they use the same words for different ideas. This page keeps the biological meaning as the default and marks the algorithmic meaning wherever a word could be read both ways.

The same words, two meanings

Quality-diversity search descends from evolutionary computation, so it names its parts after biology. In this page the biological meaning is the default; the algorithmic meaning is used only for the mechanics of the search, and the table below is the key.

Word The science (default) The search (MAP-Elites)
evolution the substitution process along a tree the mutate-and-select loop over candidate trees
mutation a modeled sequence substitution a tree move: NNI, SPR, branch slide, branch scale
fitness, performance reproductive success (not modeled here) the alignment log-likelihood of a candidate tree
population the taxa and sequences under study the archive of elite trees
parent, offspring ancestor and descendant organisms a candidate tree and the varied copy made from it
lineage a line of descent among organisms the chain of successive elite variants in the search

One rule carries the convention: if a quantity is the measured log-likelihood or part of the substitution model, it is the science; if it only steers the search, it is the algorithm.

The rest of this page develops the design from the top down. Section 1 states the search stance (illumination, not optimization). Section 2 gives the survey algorithm that fills the map. Sections 3 and 4 define the elevation: the alignment log-likelihood and the substitution model it rests on. Sections 5 through 8 build the horizontal axes: the tree metric, the three-anchor chart, the out-of-plane residual, and the rendered surface. Section 9 explains how many genes share one frame.

1. Illumination, not optimization

Most phylogenetic inference is optimization: it returns the single highest-likelihood tree, or a posterior sample concentrated near it, and discards the rest of the search. An illumination algorithm does something different. Given a low-dimensional space of features chosen by the user, which Mouret and Clune call a feature space, it returns the highest-performing solution in each region of that space, so the output is a map of how performance varies across the space rather than a single point (Mouret and Clune 2015). Mouret and Clune call the result a phenotype-fitness map: for every combination of the chosen features, it reports the best solution that has those features and how good it is.

Hifuku applies this stance to phylogenetics. The feature space is a two-dimensional chart of tree space (Section 6). The performance is the alignment log-likelihood (Section 3). The survey keeps the best tree in each region of the chart, so the result is a map of the likelihood landscape: not just the peak, but the elevation everywhere the survey reached and the shape of the terrain between peaks. This is the quality-diversity view of search: keep a diverse archive of high performers rather than one winner.

The distinction matters because the likelihood landscape over tree space has several high-likelihood topologies separated by lower-likelihood regions. An optimizer reports one of them. A survey reports all of them that it reaches, together with the elevation of the ground between, which is the quantity a reader needs to judge how well-resolved a phylogeny is and where the support is divided.

Sources

  • Mouret, J.-B., & Clune, J. (2015). Illuminating search spaces by mapping elites. arXiv:1504.04909. Section 2 (optimization versus illumination); Figure 1 (the phenotype-fitness map).
  • DOI
  • Mouret & Clune 2015

2. The elite-archive survey

The map is filled by an elite-archive survey, the MAP-Elites algorithm of Mouret and Clune (2015). An elite, in their vocabulary, is the best solution found so far in one niche, their term for a cell of the feature space. Here the elite is the best tree found so far in one niche, and the archive is the map from niche to elite. A niche is a square cell of the barycentric chart; the grid is unbounded, so every chart point has a niche, and the filled niches form the map's footprint.

The survey follows the MAP-Elites loop. Each of its three roles has a direct phylogenetic reading:

  • Feature descriptor b(T). A tree is located on the chart by its barycentric coordinate (Section 6), which selects its niche. The chart projection is many-to-one: distinct trees can share a niche, just as distinct genomes share a cell in Mouret and Clune's formulation.
  • Performance f(T). The alignment log-likelihood, computed exactly by the Felsenstein recursion (Section 3). Larger is better.
  • Variation. One tree move applied to a copy of the selected elite: a nearest-neighbor interchange, a subtree prune-and-regraft, a branch slide, or a branch scale (the same move mix used elsewhere in Hifuku).

Starting from the anchor trees and any supplied start trees, the loop repeats: select a random occupied niche, vary its elite, evaluate the varied tree's niche and log-likelihood, and place it when its niche is empty or when it beats the incumbent. There is no bias potential, no temperature, and no acceptance ratio. Placement is a pure keep-if-better rule, which is what makes the recorded elevation the measured log-likelihood and nothing else.

In Mouret and Clune's formulation a genome is decoded into a phenotype before the phenotype is scored. Hifuku has no such indirection: the object varied and the object scored are the same tree, a direct encoding. Mouret and Clune observe that most elites are found by varying an elite that already sits near them on the chart, yet the chains of successive variants traverse large regions of the feature space. The diverse archive therefore acts as a reservoir of stepping stones: a high-likelihood tree in one region seeds the discovery of high-likelihood trees in regions a single local search would not reach.

2.1 The log-likelihood gate

A survey of a likelihood landscape has no interest in the vast low-likelihood interior of tree space. Hifuku bounds the surveyed domain with a log-likelihood gate. A variation creates a new niche only when its tree scores within a margin of the best tree found so far (a relative gate that follows the running best), or at or above a fixed floor (an absolute gate, useful for holding a common domain across genes). The gate bounds exploration by likelihood rather than by the grid, so the filled footprint is the reasonably-likely region and not the whole chart.

2.2 Discovery-saturation halt

Mouret and Clune measure an illumination run by its coverage, the fraction of fillable cells filled, and its precision, the performance in filled cells relative to the best possible there. Hifuku turns the same two quantities into a stopping rule. After each batch of variations the survey forms two trailing-window rates, both made scale-invariant by dividing by the filled count:

  • a coverage rate, the number of new niches per batch as a fraction of the filled set, and
  • a precision rate, the mean per-niche elite gain per batch as a fraction of the elite log-likelihood range.

The survey halts on discovery saturation, when both rates fall to or below their thresholds: new niches have stopped appearing and the elites have stopped improving. Coverage is the crisp signal; precision uses a looser threshold because within-niche refinement never fully stops. A benchmarking override runs the full budget regardless.

Sources

3. Phylogenetic likelihood: the elevation

The elevation of the map is the likelihood of a tree: the probability of the aligned sequence data given that tree, its branch lengths, and a substitution model. Computing it naively would require summing over all ancestral character states at every internal node at once, a cost that grows exponentially with tree size. Felsenstein (1981) showed that a dynamic-programming pass reduces this to linear time: because the tree is acyclic, partial likelihoods can be computed from tips to root, each internal node depending only on its children.

3.1 The pruning recurrence

The conditional (or partial) likelihood \( \mathcal{L}_k(x) \) at internal node \( k \) is the probability of the data in the subtree rooted at \( k \), given that the state at \( k \) is \( x \). It satisfies:

\[\mathcal{L}_k(x) = \prod_{c \in \mathrm{children}(k)} \left[\sum_y P(x \to y,\, t_c)\, \mathcal{L}_c(y)\right]\]

where \( P(x \to y,\, t_c) \) is the transition probability over branch length \( t_c \) (Section 4). For a tip, \( \mathcal{L}_\mathrm{tip}(x) = 1 \) when \( x \) is the observed character and 0 otherwise (a uniform vector for ambiguity codes). Weighting the root partials by the stationary frequencies and summing gives the full likelihood:

\[L = \sum_{x} \pi_{x}\,\mathcal{L}_{\mathrm{root}}(x)\]

Sources

3.2 Per-node rescaling

Float32 arithmetic underflows on long branches: partial likelihoods near the root can fall below the smallest representable float32 value before the final summation. Hifuku rescales at each internal node. After a node's children are combined, each per-pattern partial is divided by its maximum over states, and the logarithm of that maximum is added to a per-pattern log accumulator:

\[\mathcal{L}_k^*(x,\,p) = \frac{\mathcal{L}_k(x,\,p)}{M_p}, \qquad M_p = \max_x \mathcal{L}_k(x,\,p), \qquad \ell_p \mathrel{+}= \log M_p\]

The site log-likelihood at the root is then \( \ell_p^{\mathrm{site}} = \log\!\left(\sum_x \pi_x \mathcal{L}_r^*(x,\,p)\right) + \ell_p \). Hifuku holds the partials in float32 (see CONSTRAINTS.md) and accumulates \( \ell_p \) in float64, so the dynamic range of the final sum is preserved without float64 partials.

Sources

3.3 Total log-likelihood

Summing the site terms over all \( n_\mathrm{pat} \) alignment patterns, weighted by their multiplicity \( w_p \):

\[\log P(D \mid T) = \sum_{p=1}^{n_\mathrm{pat}} w_p\, \ell_p^{\mathrm{site}}\]

Pattern compression groups identical alignment columns so each is evaluated once. Site-repeat detection extends the saving to patterns that produce identical partials at internal nodes, allowing those intermediate results to be shared.

Sources

4. Substitution models

The transition probability \( P(t) \) in the pruning recurrence comes from a substitution model: the stochastic process by which characters change along a branch. Hifuku uses time-reversible models, for which the probability flux from state \( i \) to state \( j \) in the stationary distribution equals the reverse flux. Reversibility is well-supported empirically and it yields the Pulley Principle: the likelihood does not depend on where the root is placed on an unrooted tree.

4.1 The reversible rate matrix

A time-reversible continuous-time Markov chain has a rate matrix \( Q \) of the form:

\[Q_{ij} = r_{ij}\,\pi_j \quad (i \neq j), \qquad Q_{ii} = -\sum_{j \neq i} Q_{ij}\]

where \( r_{ij} = r_{ji} \) are exchangeability parameters (six free values for four DNA states) and \( \pi \) is the stationary distribution. By convention \( Q \) is normalized to one expected substitution per unit branch length. The general four-state reversible model, GTR, leaves all six exchangeabilities free; JC69 sets them equal and fixes \( \pi_i = 1/4 \).

4.2 Eigendecomposition and the transition matrix

Computing \( P(t) = e^{Qt} \) directly at every evaluation would be expensive. Reversibility gives a faster route: the symmetrized matrix \( S = D_{\sqrt{\pi}}\, Q\, D_{\sqrt{\pi}}^{-1} \) is symmetric, so it has real eigenvalues and orthonormal eigenvectors, \( S = U \Lambda U^T \), computed once per model. The transition matrix follows by rescaling:

\[P(t) = D_{\sqrt{\pi}}^{-1}\; U\, \mathrm{diag}(e^{\lambda_k t})\; U^T\; D_{\sqrt{\pi}}\]

Rows sum to 1 and \( P(0) = I \). All eigenvalues of a valid \( Q \) are \( \leq 0 \), so entries decay toward the stationary frequencies as \( t \) grows. For protein data, Hifuku uses empirically derived rate matrices (JTT, WAG, LG), fixed from large protein databases and substituted at the eigendecomposition step; the rest of the computation is unchanged.

Sources

  • Tavaré, S. (1986). Some probabilistic and statistical problems in the analysis of DNA sequences. Lectures on Mathematics in the Life Sciences, 17, 57-86. Matrix exponential Eq. (2.7); reversibility Eqs. (2.11)-(2.13).
  • WorldCat
  • Felsenstein, J. (1981). Journal of Molecular Evolution, 17(6), 368-376. Eqs. (6)-(7) (JC69 closed form).
  • DOI
  • Tavaré 1986 · JTT · WAG · LG
  • EigenDecomp · gtr · jc69

5. The tree metric

The horizontal axes of the map measure distance between trees. Hifuku uses a normalized branch-score metric that combines branch-length difference with topological difference under one weight:

\[d^2 = (1 - w)\,\frac{\mathrm{CBS}^2}{C_\mathrm{ref}} + w\,\frac{\mathrm{RF}}{R_\mathrm{ref}}\]

The clade branch score (CBS) is a variant of the branch-score distance of Kuhner and Felsenstein (1994). Trees are unrooted, so a branch is a split of the taxa. For every split that appears in either tree, CBS sums the squared difference of the two branch lengths, counting an absent split as zero length. The Robinson-Foulds count (RF) is the number of splits present in one tree but not the other, the purely topological distance. The weight \( w \in [0, 1] \) sets topology against branch length: at \( w = 0 \) the metric is the pure branch score, at \( w = 1 \) it is the topological count.

The normalizers \( C_\mathrm{ref} \) and \( R_\mathrm{ref} \) are the mean \( \mathrm{CBS}^2 \) and mean \( \mathrm{RF} \) over the anchor set, so \( w \) means the same thing across datasets: it always mixes the two terms after each has been brought to a common scale. The GPU computes the RF term from a shared-split count with a root-pair deduplication, folded into the squared distance in every kernel, and matches the CPU reference.

Alongside each distance the metric returns a saturation score \( B_\mathrm{frac} \in [0, 1] \): 1 for identical trees, near 0 when the two trees share so few splits that their branch-score distance has lost resolution. This is metric saturation, a property of the tree pair under the branch-score metric, and is distinct from substitution saturation, which is a property of the sequence model. The chart rejects anchor pairs whose \( B_\mathrm{frac} \) falls below a threshold before any sampling starts (Section 6).

The chart owns one metric instance, built from its anchor set. A tree cannot be placed under a metric inconsistent with the chart it lands on: there is no per-call metric override, and the normalizers travel with the chart.

Sources

  • Kuhner, M. K., & Felsenstein, J. (1994). A simulation comparison of phylogeny algorithms. Molecular Biology and Evolution, 11(3), 459-468. (branch-score distance)
  • DOI
  • Robinson, D. F., & Foulds, L. R. (1981). Comparison of phylogenetic trees. Mathematical Biosciences, 53(1-2), 131-147. (topological distance)
  • DOI
  • Kuhner & Felsenstein 1994 · Robinson & Foulds 1981
  • The weighted form, its w-limits, portability under branch rescaling, and the RF shared-split identity are verified in docs/theory/algebra/kuhner1994-robinson1981-normalized-weight.py and docs/theory/algebra/robinson1981rf-shared-count.py.
  • CBSMetric · compute_normalizers

6. The three-anchor 2D chart

A single scalar distance from one reference tree is a fragile chart. It collapses all of tree space onto a line, so two topologies that differ in different directions can land at the same coordinate, and a saddle between two likelihood modes is invisible. Hifuku uses a two-dimensional chart fixed by three anchor trees \( r_1, r_2, r_3 \), each a neighbor-joining tree built from one gene alignment.

6.1 Embedding the anchors

The three pairwise tree distances \( d_{12}, d_{13}, d_{23} \) fix a triangle in the plane by the law of cosines (classical multidimensional scaling for three points), placing the anchors at \( P_1, P_2, P_3 \) with \( |P_i - P_j| = d_{ij} \). The embedding exists only when the three distances satisfy the triangle inequality, which for three points is the Schoenberg condition for Euclidean embeddability. The plane is a local Euclidean chart of tree space, which is globally non-Euclidean; the chart is faithful where the sampled region is approximately flat, and the out-of-plane residual of Section 7 measures the departure.

6.2 Triangle quality

The chart is trustworthy only when the anchors are well-spread. The triangle quality \( Q = \lambda_\mathrm{min}/\lambda_\mathrm{max} \) is the ratio of the eigenvalues of the anchor-position covariance: \( Q = 1 \) at an equilateral triangle, \( Q \to 0 \) as the anchors become collinear. \( Q \) is the a-priori proxy for geometric dilution of precision (GDOP): trilateration is most accurate when the query sees the three anchors about 120 degrees apart, which happens near the centroid of an equilateral triangle. Anchor sets below a quality floor are rejected, because a near-collinear triangle gives the second axis too little information to carry.

6.3 Locating a query tree

A query tree \( Q \) at distances \( s_1, s_2, s_3 \) from the anchors is placed by linear trilateration. Writing \( f_i(X) = |X - P_i|^2 - s_i^2 \) and subtracting \( f_1 \) from \( f_2 \) and \( f_3 \) cancels \( |X|^2 \) and leaves a 2x2 linear system for the in-plane position \( X \):

\[2\,(P_2 - P_1)\cdot X = (|P_2|^2 - |P_1|^2) - (s_2^2 - s_1^2), \qquad 2\,(P_3 - P_1)\cdot X = (|P_3|^2 - |P_1|^2) - (s_3^2 - s_1^2)\]

The position is then expressed in barycentric coordinates \( (\lambda_1, \lambda_2, \lambda_3) \) with \( \sum_k \lambda_k = 1 \), solved from \( [\lambda_1, \lambda_2]^T = [P_1 - P_3,\, P_2 - P_3]^{-1}(X - P_3) \) and \( \lambda_3 = 1 - \lambda_1 - \lambda_2 \). The barycentric form removes the sign ambiguity of a raw distance and gives the two chart axes the map is drawn on. Trees inside the triangle have all \( \lambda_k \ge 0 \); trees outside have at least one negative coordinate. A tree's position outside the triangle carries no special meaning: the anchors are coordinate references, not landmarks, and the maximum-likelihood tree can lie anywhere on the chart.

Sources

7. The out-of-plane residual

Because the plane is a local chart of a non-Euclidean space, three distances that came from a real tree need not be exactly co-planar. After solving for the in-plane position \( X \), the common residual \( c = |X - P_1|^2 - s_1^2 \) satisfies \( c = -z^2 \), where \( z \ge 0 \) is the distance of the tree from the anchor plane:

\[z = \sqrt{\max(0,\, -c)}\]

When the three distances are planar-consistent, \( z \approx 0 \). When the query lies off the plane, \( z > 0 \), and the flat chart is a lossy picture of that tree's neighborhood. The residual is stored for every placed tree and is used offline only. No decision inside the survey loop reads it: it is a measure-only diagnostic, and the chart geometry never gates on it.

Offline, \( z \) separates two causes of spread within a chart cell. Because the projection is many-to-one, several trees can share a cell. The variance of their residuals splits, by the law of total variance, into a within-source part and a between-source part. The within-source part is genuine out-of-plane roughness of the landscape; the between-source part is aliasing, several distinct sheets stacked at one chart point. Their ratio is a planar-fidelity indicator that tells the reader where the flat chart is faithful and where it has folded distinct regions together. \( z \) is a scalar probe of a high-dimensional orthogonal complement, so it is a cheap and partial diagnostic, not a full account of the missing dimensions.

Sources

  • The residual identity \( c = -z^2 \) and the variance decomposition are verified in docs/theory/algebra/zfidelity-variance-decomposition.py.
  • barycentric · z_fidelity_map

8. Rendering the map

The archive holds one elite per niche: its log-likelihood, its chart coordinate, and its out-of-plane residual. Two renderings turn the archive into a map.

Elevation. The elevation of a chart cell is a high quantile of the log-likelihood over the samples that fell in it, rather than the mean. A high quantile reports the best tree found near that chart point, which is the quantity a survey of a likelihood landscape is meant to show. The value comes from the likelihood alone, so the drawn elevation does not depend on any survey parameter.

Interpolated surface. For a smooth surface the scattered elite log-likelihoods are interpolated onto a regular grid by Shepard's inverse-distance weighting with compact support. Each grid point is a weighted average of the elites within a support radius, the weight falling with distance, and the grid is masked beyond the radius where no elite constrains it. The result is a plain gridded field that drops into a contour, mesh, or surface plot. Hifuku computes the field; the caller draws it.

Sources

  • Shepard, D. (1968). A two-dimensional interpolation function for irregularly-spaced data. Proc. 23rd ACM National Conference, 517-524.
  • DOI
  • Shepard 1968
  • The Shepard weight is verified in docs/theory/algebra/shepard1968-weight.py.
  • elevation_map · EliteArchive.surface

9. One frame for many genes

Every alignment shares a single taxon table. Coordinates are indexed in this global namespace and masked for taxa absent from a given tree, so the patristic and branch-score comparisons are well-defined even when gene trees differ in taxon coverage, and maps from different genes are directly comparable once produced from the same table.

The neighbor-joining anchors and the chart they fix are shared across genes, so each per-gene survey lands in one common frame. The chart's metric derives its normalizers from the anchor set, which can be the three anchors of one triplet or a larger global set when several triplet charts are stitched into one frame; in the latter case every chart shares one metric, so distances mean the same thing across the whole map. A completed survey is written to a self-contained HDF5 file, with every niche's elite stored as a newick tree so that any elite is a real tree offline, and reloaded into one archive per gene on the shared chart.

Sources