"""
Quartet distance module.
This module provides functions for computing distance matrices from quartet profiles
using quartet distance metrics.
"""
import itertools
from typing import TYPE_CHECKING
import numpy as np
from ...utils.exceptions import PhyloZooValueError
from ..distance import DistanceMatrix
if TYPE_CHECKING:
from .qprofile import QuartetProfile
from .qprofileset import QuartetProfileSet
[docs]
def quartet_distance(
profileset: "QuartetProfileSet",
rho: tuple[float, float, float, float],
) -> DistanceMatrix:
"""
Compute a distance matrix between taxa based on quartet profiles.
This function computes pairwise distances between taxa by aggregating contributions
from all quartet profiles. The distance is based on the quartet topology and uses
a rho vector to weight different quartet types.
Parameters
----------
profileset : QuartetProfileSet
The quartet profile set to compute distances from. Must be dense.
Each profile must contain exactly 1 or 2 resolved quartets.
rho : tuple[float, float, float, float]
Rho vector (rho_c, rho_s, rho_a, rho_o) specifying distance contributions
for different quartet topologies. rho_c: contribution when two leaves are
on the same side of a split (used for profiles with 1 quartet). rho_s:
contribution when two leaves are on different sides of a split (used for
profiles with 1 quartet). rho_a: contribution when two leaves are adjacent
in a circular ordering (used for profiles with 2 quartets). rho_o:
contribution when two leaves are opposite in a circular ordering (used for
profiles with 2 quartets). For each quartet profile containing a pair of
leaves, 2*rho_xy is added to the distance between leaves x and y, where
rho_xy is the appropriate rho value based on the profile type and leaf
positions.
Returns
-------
DistanceMatrix
A distance matrix with pairwise distances between taxa.
The matrix is symmetric with zero diagonal.
Raises
------
ValueError
If rho vector has invalid length or values.
If the profile set is not dense.
If any profile contains more than 2 quartets or contains unresolved quartets.
Examples
--------
>>> from phylozoo.core.split.base import Split
>>> from phylozoo.core.quartet.base import Quartet
>>> from phylozoo.core.quartet.qprofileset import QuartetProfileSet
>>> # Create profile set and compute distance
>>> q1 = Quartet(Split({'A', 'B'}, {'C', 'D'}))
>>> q2 = Quartet(Split({'A', 'C'}, {'B', 'D'}))
>>> profileset = QuartetProfileSet(profiles=[q1, q2])
>>> dist_matrix = quartet_distance(profileset, rho=(0.5, 1.0, 0.5, 1.0))
>>> # Access distances
>>> dist_matrix.get_distance('A', 'B')
5.0
Notes
-----
See the manual for more details.
"""
# Validate rho vector
if len(rho) != 4:
raise PhyloZooValueError(f"Rho vector must have 4 elements, got {len(rho)}")
rho_c, rho_s, rho_a, rho_o = rho
# Validate rho constraints
if rho_a > rho_o or rho_c > rho_s:
raise PhyloZooValueError("Rho vector must satisfy: rho_a <= rho_o and rho_c <= rho_s")
# Check if profile set is dense
if not profileset.is_dense:
raise PhyloZooValueError(
"Profile set must be dense (have a profile for every 4-taxon combination)"
)
# Get all taxa and initialize distance matrix as numpy array
taxa = sorted(profileset.taxa)
n = len(taxa)
D = np.zeros((n, n), dtype=np.float64)
# Create mapping from taxon to index
taxon_to_index = {taxon: i for i, taxon in enumerate(taxa)}
rho_c, rho_s, rho_a, rho_o = (float(value) for value in rho)
# The set is dense, so its profiles are exactly the C(n, 4) taxon sets; walk them
# directly. Each profile's six pair contributions follow from its split(s)
# alone: with one quartet a pair is on the same side or split; with two, the
# profile fixes one circular ordering, in which a pair is adjacent exactly when
# the two splits disagree about it and opposite when both splits separate it.
for four_taxa_set, (profile, _profile_weight) in profileset.profiles.items():
quartets = list(profile.quartets)
num_quartets = len(quartets)
if num_quartets == 0:
raise PhyloZooValueError(f"Profile for {four_taxa_set} has no quartets")
if num_quartets > 2:
raise PhyloZooValueError(
f"Profile for {four_taxa_set} has {num_quartets} quartets. "
"Each profile must contain exactly 1 or 2 quartets."
)
sides: list[set[str]] = []
for quartet in quartets:
if not quartet.is_resolved():
raise PhyloZooValueError(
f"Profile for {four_taxa_set} contains unresolved quartet (star tree). "
"All quartets must be resolved."
)
sides.append(quartet.split.set1) # type: ignore[union-attr]
if num_quartets == 1:
side = sides[0]
for leaf1, leaf2 in itertools.combinations(four_taxa_set, 2):
same_side = (leaf1 in side) == (leaf2 in side)
delta = 2.0 * (rho_c if same_side else rho_s)
i, j = taxon_to_index[leaf1], taxon_to_index[leaf2]
D[i, j] += delta
D[j, i] += delta
else:
first, second = sides
for leaf1, leaf2 in itertools.combinations(four_taxa_set, 2):
together_first = (leaf1 in first) == (leaf2 in first)
together_second = (leaf1 in second) == (leaf2 in second)
# Both splits separating the pair puts it across the four-cycle.
opposite = not together_first and not together_second
delta = 2.0 * (rho_o if opposite else rho_a)
i, j = taxon_to_index[leaf1], taxon_to_index[leaf2]
D[i, j] += delta
D[j, i] += delta
# Add constant 2*n - 4 to all off-diagonal entries
constant = 2 * n - 4
D += constant
np.fill_diagonal(D, 0) # Keep diagonal at 0
# Create and return DistanceMatrix
return DistanceMatrix(D, labels=taxa)
def _rho_distance(
profile: "QuartetProfile",
leaf1: str,
leaf2: str,
rho: tuple[float, float, float, float],
) -> float:
"""
Compute rho-distance between two leaves in a quartet profile.
If the profile has 1 quartet (split quarnet), uses split-based logic:
- Same side of split: rho_c
- Different sides: rho_s
If the profile has 2 quartets (four-cycle), uses circular ordering logic:
- Adjacent in ordering: rho_a
- Opposite in ordering: rho_o
Parameters
----------
profile : QuartetProfile
The quartet profile. Must have 1 or 2 resolved quartets.
leaf1 : str
First leaf.
leaf2 : str
Second leaf.
rho : tuple[float, float, float, float]
Rho vector (rho_c, rho_s, rho_a, rho_o).
Returns
-------
float
The rho-distance.
"""
rho_c, rho_s, rho_a, rho_o = rho
num_quartets = len(profile.quartets)
if num_quartets == 1:
# Single quartet: use split-based logic
quartet = next(iter(profile.quartets))
split = quartet.split
if split is None:
raise PhyloZooValueError("Quartet must be resolved (have a split)")
# Check if both leaves are on the same side of the split
same_side = (leaf1 in split.set1 and leaf2 in split.set1) or (
leaf1 in split.set2 and leaf2 in split.set2
)
if same_side:
return float(rho_c)
else:
return float(rho_s)
elif num_quartets == 2:
# Two quartets: use circular ordering logic
circular_orderings = profile.circular_orderings
if circular_orderings is None or len(circular_orderings) == 0:
raise PhyloZooValueError("Profile with 2 quartets must have a circular ordering")
# Use the first circular ordering (they should all be equivalent)
ordering = next(iter(circular_orderings))
# Check if leaves are adjacent or opposite in the circular ordering
if ordering.are_neighbors(leaf1, leaf2):
# Adjacent: use rho_a
return float(rho_a)
else:
# Opposite: use rho_o
return float(rho_o)
else:
raise PhyloZooValueError(f"Profile must have 1 or 2 quartets, got {num_quartets}")
