Anchor Triangle and Chart

The chart is fixed by three anchor trees. This module builds the anchor triangle from pairwise clade-branch-score distances, validates its conditioning, and places a query tree in barycentric coordinates. See The Three-Anchor 2D Chart.

AnchorTriangle dataclass

Geometry of the three-anchor 2D chart.

The plain constructor stores already-computed geometry. To build a chart from three anchor trees (deriving the metric, measuring distances, running the validation checks, and solving the embedding), use :meth:from_anchors.

Attributes:

Name	Type	Description
P1, P2, P3	(ndarray, shape(2))	Anchor positions in the local Euclidean chart.
r1, r2, r3	HifukuTree	The three anchor trees.
dist	(ndarray, shape(3, 3))	Pairwise metric distances.
qual	(ndarray, shape(3, 3))	Pairwise B_frac quality scores.
metric	TreeMetric	The chart's metric instance, built from an anchor set (this triangle's own anchors, or a global set when stitching several triplets into one frame). Required: a chart without a metric is broken.

Source code in src/hifuku/anchor_triangle.py

@dataclass
class AnchorTriangle:
    """
    Geometry of the three-anchor 2D chart.

    The plain constructor stores already-computed geometry.  To build a chart
    from three anchor trees (deriving the metric, measuring distances, running
    the validation checks, and solving the embedding), use
    :meth:`from_anchors`.

    Attributes
    ----------
    P1, P2, P3 : np.ndarray, shape (2,)
        Anchor positions in the local Euclidean chart.
    r1, r2, r3 : HifukuTree
        The three anchor trees.
    dist : np.ndarray, shape (3, 3)
        Pairwise metric distances.
    qual : np.ndarray, shape (3, 3)
        Pairwise B_frac quality scores.
    metric : TreeMetric
        The chart's metric instance, built from an anchor set (this triangle's
        own anchors, or a global set when stitching several triplets into one
        frame).  Required: a chart without a metric is broken.
    """
    P1: np.ndarray
    P2: np.ndarray
    P3: np.ndarray
    r1: HifukuTree
    r2: HifukuTree
    r3: HifukuTree
    dist: np.ndarray
    qual: np.ndarray
    metric: object = None

    @property
    def tree_metric(self) -> TreeMetric:
        """The chart's metric instance.  Raises if the chart is missing its
        metric or its anchors, both of which a usable chart must have."""
        if self.metric is None:
            raise ValueError(
                "AnchorTriangle has no metric; build the chart with validate_anchors."
            )
        if self.r1 is None or self.r2 is None or self.r3 is None:
            raise ValueError("AnchorTriangle is missing anchor trees.")
        return self.metric

    @classmethod
    def from_anchors(
        cls,
        r1: HifukuTree,
        r2: HifukuTree,
        r3: HifukuTree,
        *,
        w: float = METRIC_W,
        stitch_anchors=None,
        b_frac_threshold: float = B_FRAC_THRESHOLD,
        anchor_q_min: float = ANCHOR_Q_MIN,
    ) -> "AnchorTriangle":
        """
        Build a validated chart from three anchor trees.

        This is the entry point that turns three anchor trees into a chart.  It
        does the work the plain constructor does not: it derives the chart's
        metric from the anchor set, measures the three pairwise distances,
        checks that the anchors form a usable chart (any failure raises), solves
        the triangle embedding, and returns a fully populated instance.  The
        plain constructor is the low-level counterpart: it takes the already
        computed fields (``P1``, ``P2``, ``P3``, ``dist``, ``qual``,
        ``metric``, ...) and only stores them.  Use it when the geometry is
        already known (loading a saved chart, or a test with fixed geometry);
        use ``from_anchors`` when starting from trees.

        Three checks run in order (any failure raises ValueError and stops):

        1. **Metric saturation** - the B_frac quality score for each of the
           three anchor pairs must be >= b_frac_threshold.  Near-zero B_frac
           means the branch-score metric has lost resolution for that pair
           (metric saturation; distinct from substitution saturation, which is
           a property of the sequence model).

        2. **Chartability** (triangle inequality / Euclidean-embeddability) -
           the three distances must form a valid planar triangle.  Equivalently,
           the double-centered squared-distance matrix must be positive
           semidefinite (Schoenberg's condition for 3 points).  Checked before
           the sqrt in ``embed_anchor_triangle``.

        3. **Triangle quality** Q >= anchor_q_min - the anchor-position
           covariance must be non-degenerate (the three anchors are not
           near-collinear).  The second axis carries too little information to
           form a useful 2D CV when Q is below this floor.

        Parameters
        ----------
        r1, r2, r3 : HifukuTree
            Anchor trees restricted to the shared global namespace, with
            ``global_leaf_indices`` set.
        w : float
            Topology weight in [0, 1] for the chart's metric (default METRIC_W).
        stitch_anchors : sequence of HifukuTree, optional
            Anchor set the metric normalizers are derived from.  Default is the
            triangle's own three anchors.  Pass the full global anchor set when
            stitching several triplet charts into one common frame, so every
            chart shares one metric.
        b_frac_threshold : float
            Minimum B_frac quality score per anchor pair (default
            B_FRAC_THRESHOLD).
        anchor_q_min : float
            Minimum triangle quality Q (default ANCHOR_Q_MIN).

        Returns
        -------
        AnchorTriangle
            The fully-built chart: embedding positions, distance/quality
            matrices, the three anchors, and the metric instance the chart owns.

        Raises
        ------
        ValueError
            On metric saturation, chartability failure, or insufficient triangle
            quality.  The error message names the failing pair/triple and the
            relevant diagnostic values.
        """
        metric = CBSMetric.for_anchors(
            list(stitch_anchors) if stitch_anchors is not None else [r1, r2, r3],
            w=w,
        )
        dist, qual = compute_anchor_distances(r1, r2, r3, metric)

        # 1. Metric saturation check.
        pair_names = [("r1", "r2"), ("r1", "r3"), ("r2", "r3")]
        pair_indices = [(0, 1), (0, 2), (1, 2)]
        for (i, j), (na, nb) in zip(pair_indices, pair_names):
            q = qual[i, j]
            if q < b_frac_threshold:
                raise ValueError(
                    f"Metric saturation: anchor pair ({na}, {nb}) has B_frac = "
                    f"{q:.4f} < b_frac_threshold = {b_frac_threshold}.  "
                    f"The branch-score metric has lost resolution for this pair "
                    f"(metric saturation; distinct from substitution saturation).  "
                    f"Distance d = {dist[i,j]:.4f}."
                )

        # 2. Chartability check (triangle inequality / Schoenberg condition).
        d12 = float(dist[0, 1])
        d13 = float(dist[0, 2])
        d23 = float(dist[1, 2])
        P1, P2, P3 = embed_anchor_triangle(d12, d13, d23)

        # 3. Triangle quality check.
        Q, lmin, lmax = triangle_quality(P1, P2, P3)
        if Q < anchor_q_min:
            raise ValueError(
                f"Triangle quality Q = {Q:.4f} < anchor_q_min = {anchor_q_min}.  "
                f"Anchor distances: d12={d12:.4g}, d13={d13:.4g}, d23={d23:.4g}.  "
                f"Covariance eigenvalues: lambda_min={lmin:.4g}, lambda_max={lmax:.4g}.  "
                f"The three anchors are near-collinear; the second CV axis carries "
                f"too little information."
            )

        return cls(
            P1=P1,
            P2=P2,
            P3=P3,
            r1=r1,
            r2=r2,
            r3=r3,
            dist=dist,
            qual=qual,
            metric=metric,
        )

tree_metric property

tree_metric: TreeMetric

The chart's metric instance. Raises if the chart is missing its metric or its anchors, both of which a usable chart must have.

from_anchors classmethod

from_anchors(r1: HifukuTree, r2: HifukuTree, r3: HifukuTree, *, w: float = METRIC_W, stitch_anchors=None, b_frac_threshold: float = B_FRAC_THRESHOLD, anchor_q_min: float = ANCHOR_Q_MIN) -> 'AnchorTriangle'

Build a validated chart from three anchor trees.

This is the entry point that turns three anchor trees into a chart. It does the work the plain constructor does not: it derives the chart's metric from the anchor set, measures the three pairwise distances, checks that the anchors form a usable chart (any failure raises), solves the triangle embedding, and returns a fully populated instance. The plain constructor is the low-level counterpart: it takes the already computed fields (P1, P2, P3, dist, qual, metric, ...) and only stores them. Use it when the geometry is already known (loading a saved chart, or a test with fixed geometry); use from_anchors when starting from trees.

Three checks run in order (any failure raises ValueError and stops):

1. Metric saturation - the B_frac quality score for each of the three anchor pairs must be >= b_frac_threshold. Near-zero B_frac means the branch-score metric has lost resolution for that pair (metric saturation; distinct from substitution saturation, which is a property of the sequence model).

2. Chartability (triangle inequality / Euclidean-embeddability) - the three distances must form a valid planar triangle. Equivalently, the double-centered squared-distance matrix must be positive semidefinite (Schoenberg's condition for 3 points). Checked before the sqrt in embed_anchor_triangle.

3. Triangle quality Q >= anchor_q_min - the anchor-position covariance must be non-degenerate (the three anchors are not near-collinear). The second axis carries too little information to form a useful 2D CV when Q is below this floor.

Parameters:

Name	Type	Description	Default
r1	HifukuTree	Anchor trees restricted to the shared global namespace, with global_leaf_indices set.	required
r2	HifukuTree	Anchor trees restricted to the shared global namespace, with global_leaf_indices set.	required
r3	HifukuTree	Anchor trees restricted to the shared global namespace, with global_leaf_indices set.	required
w	float	Topology weight in [0, 1] for the chart's metric (default METRIC_W).	METRIC_W
stitch_anchors	sequence of HifukuTree	Anchor set the metric normalizers are derived from. Default is the triangle's own three anchors. Pass the full global anchor set when stitching several triplet charts into one common frame, so every chart shares one metric.	None
b_frac_threshold	float	Minimum B_frac quality score per anchor pair (default B_FRAC_THRESHOLD).	B_FRAC_THRESHOLD
anchor_q_min	float	Minimum triangle quality Q (default ANCHOR_Q_MIN).	ANCHOR_Q_MIN

Returns:

Type	Description
AnchorTriangle	The fully-built chart: embedding positions, distance/quality matrices, the three anchors, and the metric instance the chart owns.

Raises:

Type	Description
ValueError	On metric saturation, chartability failure, or insufficient triangle quality. The error message names the failing pair/triple and the relevant diagnostic values.

Source code in src/hifuku/anchor_triangle.py

@classmethod
def from_anchors(
    cls,
    r1: HifukuTree,
    r2: HifukuTree,
    r3: HifukuTree,
    *,
    w: float = METRIC_W,
    stitch_anchors=None,
    b_frac_threshold: float = B_FRAC_THRESHOLD,
    anchor_q_min: float = ANCHOR_Q_MIN,
) -> "AnchorTriangle":
    """
    Build a validated chart from three anchor trees.

    This is the entry point that turns three anchor trees into a chart.  It
    does the work the plain constructor does not: it derives the chart's
    metric from the anchor set, measures the three pairwise distances,
    checks that the anchors form a usable chart (any failure raises), solves
    the triangle embedding, and returns a fully populated instance.  The
    plain constructor is the low-level counterpart: it takes the already
    computed fields (``P1``, ``P2``, ``P3``, ``dist``, ``qual``,
    ``metric``, ...) and only stores them.  Use it when the geometry is
    already known (loading a saved chart, or a test with fixed geometry);
    use ``from_anchors`` when starting from trees.

    Three checks run in order (any failure raises ValueError and stops):

    1. **Metric saturation** - the B_frac quality score for each of the
       three anchor pairs must be >= b_frac_threshold.  Near-zero B_frac
       means the branch-score metric has lost resolution for that pair
       (metric saturation; distinct from substitution saturation, which is
       a property of the sequence model).

    2. **Chartability** (triangle inequality / Euclidean-embeddability) -
       the three distances must form a valid planar triangle.  Equivalently,
       the double-centered squared-distance matrix must be positive
       semidefinite (Schoenberg's condition for 3 points).  Checked before
       the sqrt in ``embed_anchor_triangle``.

    3. **Triangle quality** Q >= anchor_q_min - the anchor-position
       covariance must be non-degenerate (the three anchors are not
       near-collinear).  The second axis carries too little information to
       form a useful 2D CV when Q is below this floor.

    Parameters
    ----------
    r1, r2, r3 : HifukuTree
        Anchor trees restricted to the shared global namespace, with
        ``global_leaf_indices`` set.
    w : float
        Topology weight in [0, 1] for the chart's metric (default METRIC_W).
    stitch_anchors : sequence of HifukuTree, optional
        Anchor set the metric normalizers are derived from.  Default is the
        triangle's own three anchors.  Pass the full global anchor set when
        stitching several triplet charts into one common frame, so every
        chart shares one metric.
    b_frac_threshold : float
        Minimum B_frac quality score per anchor pair (default
        B_FRAC_THRESHOLD).
    anchor_q_min : float
        Minimum triangle quality Q (default ANCHOR_Q_MIN).

    Returns
    -------
    AnchorTriangle
        The fully-built chart: embedding positions, distance/quality
        matrices, the three anchors, and the metric instance the chart owns.

    Raises
    ------
    ValueError
        On metric saturation, chartability failure, or insufficient triangle
        quality.  The error message names the failing pair/triple and the
        relevant diagnostic values.
    """
    metric = CBSMetric.for_anchors(
        list(stitch_anchors) if stitch_anchors is not None else [r1, r2, r3],
        w=w,
    )
    dist, qual = compute_anchor_distances(r1, r2, r3, metric)

    # 1. Metric saturation check.
    pair_names = [("r1", "r2"), ("r1", "r3"), ("r2", "r3")]
    pair_indices = [(0, 1), (0, 2), (1, 2)]
    for (i, j), (na, nb) in zip(pair_indices, pair_names):
        q = qual[i, j]
        if q < b_frac_threshold:
            raise ValueError(
                f"Metric saturation: anchor pair ({na}, {nb}) has B_frac = "
                f"{q:.4f} < b_frac_threshold = {b_frac_threshold}.  "
                f"The branch-score metric has lost resolution for this pair "
                f"(metric saturation; distinct from substitution saturation).  "
                f"Distance d = {dist[i,j]:.4f}."
            )

    # 2. Chartability check (triangle inequality / Schoenberg condition).
    d12 = float(dist[0, 1])
    d13 = float(dist[0, 2])
    d23 = float(dist[1, 2])
    P1, P2, P3 = embed_anchor_triangle(d12, d13, d23)

    # 3. Triangle quality check.
    Q, lmin, lmax = triangle_quality(P1, P2, P3)
    if Q < anchor_q_min:
        raise ValueError(
            f"Triangle quality Q = {Q:.4f} < anchor_q_min = {anchor_q_min}.  "
            f"Anchor distances: d12={d12:.4g}, d13={d13:.4g}, d23={d23:.4g}.  "
            f"Covariance eigenvalues: lambda_min={lmin:.4g}, lambda_max={lmax:.4g}.  "
            f"The three anchors are near-collinear; the second CV axis carries "
            f"too little information."
        )

    return cls(
        P1=P1,
        P2=P2,
        P3=P3,
        r1=r1,
        r2=r2,
        r3=r3,
        dist=dist,
        qual=qual,
        metric=metric,
    )

TreeMetric

Bases: Protocol

Protocol for tree distance metrics used by the 2D CV.

A callable that takes two HifukuTree objects and returns a (distance, quality) pair. The quality score is a reliability / confidence value interpretable as a likelihood-like weight: high quality = the metric distinguishes the two trees well; low quality = metric saturation (the metric has lost resolution; distinct from substitution saturation).

Only the CBS concrete implementation (CBSMetric below) is provided. Code outside it must not assume CBS internals.

Source code in src/hifuku/metric/branch_score.py

@runtime_checkable
class TreeMetric(Protocol):
    """
    Protocol for tree distance metrics used by the 2D CV.

    A callable that takes two HifukuTree objects and returns a (distance,
    quality) pair.  The quality score is a reliability / confidence value
    interpretable as a likelihood-like weight: high quality = the metric
    distinguishes the two trees well; low quality = metric saturation (the
    metric has lost resolution; distinct from substitution saturation).

    Only the CBS concrete implementation (CBSMetric below) is provided.  Code
    outside it must not assume CBS internals.
    """

    def __call__(
        self, t1: HifukuTree, t2: HifukuTree
    ) -> tuple[float, float]:
        """
        Parameters
        ----------
        t1, t2 : HifukuTree
            Trees with ``global_leaf_indices`` set and restricted to the same
            global taxon namespace.

        Returns
        -------
        (distance, quality) : (float, float)
            ``distance >= 0``.  ``quality`` in [0, 1], where 0 = metric
            saturation and 1 = identical trees.
        """
        ...

CBSMetric dataclass

Normalized branch-score implementation of the TreeMetric protocol.

distance = d = sqrt((1 - w) * CBS^2 / c_ref + w * RF / r_ref): the Kuhner-Felsenstein (1994) clade branch score combined with the Robinson-Foulds topological distance, weighted by w in [0, 1] (0 = pure branch score, 1 = pure topology). c_ref and r_ref normalize the two terms; the defaults 1.0 give the raw branch score at w = 0. Use :meth:for_anchors to derive the normalizers from an anchor set. quality = B_frac in [0, 1]; the length-overlap fraction. Near 0 = branch- score saturation; near 1 = similar trees.

kuhner1994branch, robinson1981rf.

Source code in src/hifuku/metric/branch_score.py

@dataclass
class CBSMetric:
    """
    Normalized branch-score implementation of the TreeMetric protocol.

    ``distance`` = d = sqrt((1 - w) * CBS^2 / c_ref + w * RF / r_ref): the
    Kuhner-Felsenstein (1994) clade branch score combined with the Robinson-Foulds
    topological distance, weighted by ``w`` in [0, 1] (0 = pure branch score,
    1 = pure topology).  ``c_ref`` and ``r_ref`` normalize the two terms; the
    defaults 1.0 give the raw branch score at ``w = 0``.  Use :meth:`for_anchors`
    to derive the normalizers from an anchor set.
    ``quality`` = B_frac in [0, 1]; the length-overlap fraction.  Near 0 = branch-
    score saturation; near 1 = similar trees.

    kuhner1994branch, robinson1981rf.
    """

    w: float = METRIC_W
    c_ref: float = 1.0
    r_ref: float = 1.0

    def __call__(
        self, t1: HifukuTree, t2: HifukuTree
    ) -> tuple[float, float]:
        return _cbs_dist_quality(t1, t2, self.w, self.c_ref, self.r_ref)

    @classmethod
    def for_anchors(cls, anchors, taxon_table=None, w: float = METRIC_W) -> "CBSMetric":
        """
        Build a metric whose normalizers are the mean CBS^2 and mean RF over the
        anchor set.

        ``w`` then means the same across datasets, and every chart built on these
        anchors shares one metric, so the charts stitch into one common frame.
        The anchors must have ``global_leaf_indices`` set.
        """
        c_ref, r_ref = compute_normalizers(anchors, taxon_table)
        return cls(w=w, c_ref=c_ref, r_ref=r_ref)

for_anchors classmethod

for_anchors(anchors, taxon_table=None, w: float = METRIC_W) -> 'CBSMetric'

Build a metric whose normalizers are the mean CBS^2 and mean RF over the anchor set.

w then means the same across datasets, and every chart built on these anchors shares one metric, so the charts stitch into one common frame. The anchors must have global_leaf_indices set.

Source code in src/hifuku/metric/branch_score.py

@classmethod
def for_anchors(cls, anchors, taxon_table=None, w: float = METRIC_W) -> "CBSMetric":
    """
    Build a metric whose normalizers are the mean CBS^2 and mean RF over the
    anchor set.

    ``w`` then means the same across datasets, and every chart built on these
    anchors shares one metric, so the charts stitch into one common frame.
    The anchors must have ``global_leaf_indices`` set.
    """
    c_ref, r_ref = compute_normalizers(anchors, taxon_table)
    return cls(w=w, c_ref=c_ref, r_ref=r_ref)

compute_anchor_distances

compute_anchor_distances(r1: HifukuTree, r2: HifukuTree, r3: HifukuTree, metric: TreeMetric) -> tuple[np.ndarray, np.ndarray]

Compute pairwise distances and quality scores for three anchor trees.

Parameters:

Name	Type	Description	Default
r1	HifukuTree	Anchor trees, each restricted to the shared global taxon namespace with global_leaf_indices set.	required
r2	HifukuTree	Anchor trees, each restricted to the shared global taxon namespace with global_leaf_indices set.	required
r3	HifukuTree	Anchor trees, each restricted to the shared global taxon namespace with global_leaf_indices set.	required
metric	TreeMetric	The metric to use for the pairwise distances (see validate_anchors).	required

Returns:

Type	Description
(dist, qual) : two (3, 3) float64 arrays	dist[i,j] = the metric distance between anchor i and anchor j. qual[i,j] = B_frac quality score in [0, 1]; 1 = identical trees, 0 = metric saturation. Both arrays are symmetric with zero diagonal.

Source code in src/hifuku/anchor_triangle.py

def compute_anchor_distances(
    r1: HifukuTree,
    r2: HifukuTree,
    r3: HifukuTree,
    metric: TreeMetric,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Compute pairwise distances and quality scores for three anchor trees.

    Parameters
    ----------
    r1, r2, r3 : HifukuTree
        Anchor trees, each restricted to the shared global taxon namespace
        with ``global_leaf_indices`` set.
    metric : TreeMetric
        The metric to use for the pairwise distances (see ``validate_anchors``).

    Returns
    -------
    (dist, qual) : two (3, 3) float64 arrays
        ``dist[i,j]`` = the metric distance between anchor i and anchor j.
        ``qual[i,j]`` = B_frac quality score in [0, 1]; 1 = identical trees,
        0 = metric saturation.  Both arrays are symmetric with zero diagonal.
    """
    anchors = [r1, r2, r3]
    dist = np.zeros((3, 3), dtype=np.float64)
    qual = np.ones((3, 3), dtype=np.float64)  # diagonal = 1 (self-quality)

    for i in range(3):
        for j in range(i + 1, 3):
            s, q = metric(anchors[i], anchors[j])
            dist[i, j] = dist[j, i] = s
            qual[i, j] = qual[j, i] = q

    return dist, qual

embed_anchor_triangle

embed_anchor_triangle(d12: float, d13: float, d23: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]

Embed three anchor trees in a 2D plane from their pairwise tree distances.

This is the law-of-cosines (classical MDS for three points) embedding:

P1 = (0, 0)
P2 = (d12, 0)
x3 = (d12^2 + d13^2 - d23^2) / (2 * d12)
y3 = sqrt(d13^2 - x3^2)
P3 = (x3, y3)

The resulting positions satisfy |P1-P2| = d12, |P1-P3| = d13, |P2-P3| = d23, providing the local Euclidean chart of tree space.

Parameters:

Name	Type	Description	Default
d12	float	Pairwise CBS distances between anchors 1-2, 1-3, and 2-3.	required
d13	float	Pairwise CBS distances between anchors 1-2, 1-3, and 2-3.	required
d23	float	Pairwise CBS distances between anchors 1-2, 1-3, and 2-3.	required

Returns:

Type	Description
P1, P2, P3 : np.ndarray, shape (2,)	Anchor positions in the 2D chart.

Raises:

Type	Description
ValueError	If d12 == 0 (degenerate anchors 1 and 2), or if the three distances violate the triangle inequality (the square inside the sqrt is negative, indicating the distances are not Euclidean-embeddable). The latter is the Euclidean-embeddability (Schoenberg) condition for 3 points.

Notes

The 2D plane is a local Euclidean chart of tree space, which is CAT(0) / non-Euclidean. This chart is valid only where the sampled region is approximately flat. The out-of-plane residual z (Section 4.3/4.11) quantifies the a-posteriori deviation from flatness per sample.

kuhner1994branch : CBS split branch score; distances used here.

Source code in src/hifuku/anchor_triangle.py

def embed_anchor_triangle(
    d12: float,
    d13: float,
    d23: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Embed three anchor trees in a 2D plane from their pairwise tree distances.

    This is the law-of-cosines (classical MDS for three points) embedding:

        P1 = (0, 0)
        P2 = (d12, 0)
        x3 = (d12^2 + d13^2 - d23^2) / (2 * d12)
        y3 = sqrt(d13^2 - x3^2)
        P3 = (x3, y3)

    The resulting positions satisfy |P1-P2| = d12, |P1-P3| = d13,
    |P2-P3| = d23, providing the local Euclidean chart of tree space.

    Parameters
    ----------
    d12, d13, d23 : float
        Pairwise CBS distances between anchors 1-2, 1-3, and 2-3.

    Returns
    -------
    P1, P2, P3 : np.ndarray, shape (2,)
        Anchor positions in the 2D chart.

    Raises
    ------
    ValueError
        If d12 == 0 (degenerate anchors 1 and 2), or if the three distances
        violate the triangle inequality (the square inside the sqrt is negative,
        indicating the distances are not Euclidean-embeddable).  The latter is
        the Euclidean-embeddability (Schoenberg) condition for 3 points.

    Notes
    -----
    The 2D plane is a **local Euclidean chart** of tree space, which is
    CAT(0) / non-Euclidean.  This chart is valid only where the sampled region
    is approximately flat.  The out-of-plane residual z (Section 4.3/4.11)
    quantifies the a-posteriori deviation from flatness per sample.

    kuhner1994branch : CBS split branch score; distances used here.
    """
    if d12 == 0.0:
        raise ValueError(
            f"Anchors 1 and 2 have identical CBS distances (d12=0); "
            f"they may be the same tree."
        )

    P1 = np.array([0.0, 0.0])
    P2 = np.array([d12, 0.0])

    x3 = (d12 * d12 + d13 * d13 - d23 * d23) / (2.0 * d12)
    radicand = d13 * d13 - x3 * x3

    if radicand < 0.0:
        raise ValueError(
            f"Anchor distances violate the triangle inequality: "
            f"d12={d12:.6g}, d13={d13:.6g}, d23={d23:.6g} do not admit a "
            f"planar triangle (radicand = {radicand:.6g} < 0).  "
            f"The three anchor trees are not Euclidean-embeddable at these "
            f"distances; this is the Schoenberg chartability condition."
        )

    y3 = math.sqrt(radicand)
    P3 = np.array([x3, y3])

    return P1, P2, P3

triangle_quality

triangle_quality(P1: ndarray, P2: ndarray, P3: ndarray) -> tuple[float, float, float]

Compute the triangle quality Q = lambda_min / lambda_max.

Q is the ratio of the minimum to the maximum eigenvalue of the 2x2 population covariance of the three anchor positions:

centroid = (P1 + P2 + P3) / 3
Cov = (1/3) * sum_i (Pi - centroid)(Pi - centroid)^T
Q = lambda_min(Cov) / lambda_max(Cov)

Q in (0, 1]: 1 at the equilateral triangle, 0 at collinear anchors. Q is scale-invariant and pivot-independent. It is the a-priori proxy for geometric dilution of precision (GDOP): trilateration is most accurate when the query sees the three anchors ~120 degrees apart, which occurs at the centroid of an equilateral triangle.

Parameters:

Name	Type	Description	Default
P1	(ndarray, shape(2))	Anchor positions in the 2D Euclidean chart.	required
P2	(ndarray, shape(2))	Anchor positions in the 2D Euclidean chart.	required
P3	(ndarray, shape(2))	Anchor positions in the 2D Euclidean chart.	required

Returns:

Type	Description
(Q, lambda_min, lambda_max) : (float, float, float)	Q is the quality ratio; lambda_min and lambda_max are the covariance eigenvalues (for diagnostic reporting).

Raises:

Type	Description
ValueError	If lambda_max is zero (all three anchors are identical points).

Notes

parkinson1996gdop : GDOP, geometric dilution of precision in navigation; Q is the anchor-level GDOP proxy prior to any sample evaluation.

Source code in src/hifuku/anchor_triangle.py

def triangle_quality(
    P1: np.ndarray,
    P2: np.ndarray,
    P3: np.ndarray,
) -> tuple[float, float, float]:
    """
    Compute the triangle quality Q = lambda_min / lambda_max.

    Q is the ratio of the minimum to the maximum eigenvalue of the 2x2
    population covariance of the three anchor positions:

        centroid = (P1 + P2 + P3) / 3
        Cov = (1/3) * sum_i (Pi - centroid)(Pi - centroid)^T
        Q = lambda_min(Cov) / lambda_max(Cov)

    Q in (0, 1]: 1 at the equilateral triangle, 0 at collinear anchors.
    Q is scale-invariant and pivot-independent.  It is the a-priori proxy for
    geometric dilution of precision (GDOP): trilateration is most accurate when
    the query sees the three anchors ~120 degrees apart, which occurs at the
    centroid of an equilateral triangle.

    Parameters
    ----------
    P1, P2, P3 : np.ndarray, shape (2,)
        Anchor positions in the 2D Euclidean chart.

    Returns
    -------
    (Q, lambda_min, lambda_max) : (float, float, float)
        Q is the quality ratio; lambda_min and lambda_max are the covariance
        eigenvalues (for diagnostic reporting).

    Raises
    ------
    ValueError
        If lambda_max is zero (all three anchors are identical points).

    Notes
    -----
    parkinson1996gdop : GDOP, geometric dilution of precision in navigation;
    Q is the anchor-level GDOP proxy prior to any sample evaluation.
    """
    pts = np.array([P1, P2, P3], dtype=np.float64)  # (3, 2)
    centroid = pts.mean(axis=0)
    diffs = pts - centroid                            # (3, 2)
    cov = (diffs.T @ diffs) / 3.0                    # (2, 2)

    evals = np.linalg.eigvalsh(cov)   # ascending order
    lambda_min = float(evals[0])
    lambda_max = float(evals[1])

    if lambda_max == 0.0:
        raise ValueError(
            "All three anchor positions are identical; cannot compute Q."
        )

    Q = lambda_min / lambda_max
    return Q, lambda_min, lambda_max

barycentric

barycentric(s1: float, s2: float, s3: float, triangle: AnchorTriangle) -> tuple[np.ndarray, float]

Map a query tree's anchor distances to barycentric coordinates plus the out-of-plane residual z.

Given distances s1, s2, s3 from the query to the three anchors (obtained via the CBS TreeMetric), solve for the planar position X and convert to barycentric coordinates (lambda1, lambda2, lambda3).

Solve method: linear trilateration. Write:

f_i(X) = |X - Pi|^2 - si^2

Subtract f_1 from f_2 and f_3 to cancel |X|^2, giving a 2x2 linear system:

2 (P2 - P1) . X = (|P2|^2 - |P1|^2) - (s2^2 - s1^2)
2 (P3 - P1) . X = (|P3|^2 - |P1|^2) - (s3^2 - s1^2)

Solve for X (the in-plane position), then convert to barycentric via:

[lambda1, lambda2]^T = [P1-P3, P2-P3]^{-1} (X - P3)
lambda3 = 1 - lambda1 - lambda2

Out-of-plane residual (measure-only). After solving for X, the residual c = f_1(X) = f_2(X) = f_3(X) satisfies c = -z^2, where z is the distance of the true tree position from the anchor plane. When the three distances si are planar-consistent, z ~ 0. When the query lies off the plane (the 2D CV is lossy there), z > 0. No code in the loop may branch on z; it is stored as a diagnostic only.

Geometric dilution of precision (GDOP): the reliability of the in-plane position depends on the angles at which the query sees the three anchors. With unit line-of-sight vectors u_i = (X - Pi) / |X - Pi| and H = [u1^T; u2^T; u3^T], GDOP = sqrt(trace((H^T H)^{-1})). GDOP is smallest at the triangle centroid (equilateral, ~120 deg between anchors) and grows as the query moves outside the triangle. GDOP is the in-plane companion to z; both are measure-only diagnostics. parkinson1996gdop.

Parameters:

Name	Type	Description	Default
s1	float	CBS distances from the query tree to r1, r2, r3.	required
s2	float	CBS distances from the query tree to r1, r2, r3.	required
s3	float	CBS distances from the query tree to r1, r2, r3.	required
triangle	AnchorTriangle	Precomputed anchor geometry.	required

Returns:

Type	Description
(lambdas, z) : (np.ndarray shape (3,), float)	lambdas = [lambda1, lambda2, lambda3] with sum == 1. Trees inside the triangle have all lambda_i >= 0; trees outside have at least one lambda_i < 0. z >= 0 is the out-of-plane residual (measure-only, no code may branch on it in the loop).

Source code in src/hifuku/anchor_triangle.py

def barycentric(
    s1: float,
    s2: float,
    s3: float,
    triangle: AnchorTriangle,
) -> tuple[np.ndarray, float]:
    """
    Map a query tree's anchor distances to barycentric coordinates plus
    the out-of-plane residual z.

    Given distances s1, s2, s3 from the query to the three anchors (obtained
    via the CBS TreeMetric), solve for the planar position X and convert to
    barycentric coordinates (lambda1, lambda2, lambda3).

    **Solve method: linear trilateration.**  Write:

        f_i(X) = |X - Pi|^2 - si^2

    Subtract f_1 from f_2 and f_3 to cancel |X|^2, giving a 2x2 linear system:

        2 (P2 - P1) . X = (|P2|^2 - |P1|^2) - (s2^2 - s1^2)
        2 (P3 - P1) . X = (|P3|^2 - |P1|^2) - (s3^2 - s1^2)

    Solve for X (the in-plane position), then convert to barycentric via:

        [lambda1, lambda2]^T = [P1-P3, P2-P3]^{-1} (X - P3)
        lambda3 = 1 - lambda1 - lambda2

    **Out-of-plane residual (measure-only).**  After solving for X, the
    residual c = f_1(X) = f_2(X) = f_3(X) satisfies c = -z^2, where z is the
    distance of the true tree position from the anchor plane.  When the three
    distances si are planar-consistent, z ~ 0.  When the query lies off the
    plane (the 2D CV is lossy there), z > 0.  No code in the loop may branch
    on z; it is stored as a diagnostic only.

    Geometric dilution of precision (GDOP): the reliability of the in-plane
    position depends on the angles at which the query sees the three anchors.
    With unit line-of-sight vectors u_i = (X - Pi) / |X - Pi| and
    H = [u1^T; u2^T; u3^T], GDOP = sqrt(trace((H^T H)^{-1})).  GDOP is
    smallest at the triangle centroid (equilateral, ~120 deg between anchors)
    and grows as the query moves outside the triangle.  GDOP is the in-plane
    companion to z; both are measure-only diagnostics.
    parkinson1996gdop.

    Parameters
    ----------
    s1, s2, s3 : float
        CBS distances from the query tree to r1, r2, r3.
    triangle : AnchorTriangle
        Precomputed anchor geometry.

    Returns
    -------
    (lambdas, z) : (np.ndarray shape (3,), float)
        ``lambdas = [lambda1, lambda2, lambda3]`` with sum == 1.  Trees inside
        the triangle have all lambda_i >= 0; trees outside have at least one
        lambda_i < 0.
        ``z >= 0`` is the out-of-plane residual (measure-only, no code may
        branch on it in the loop).
    """
    P1, P2, P3 = triangle.P1, triangle.P2, triangle.P3

    # Set up the 2x2 linear system by subtracting the first sphere equation.
    rhs1 = (np.dot(P2, P2) - np.dot(P1, P1)) - (s2 * s2 - s1 * s1)
    rhs2 = (np.dot(P3, P3) - np.dot(P1, P1)) - (s3 * s3 - s1 * s1)

    A_mat = 2.0 * np.array([P2 - P1, P3 - P1])   # (2, 2)
    rhs = np.array([rhs1, rhs2])

    X = np.linalg.solve(A_mat, rhs)  # in-plane position (2,)

    # Barycentric coordinates: [P1-P3, P2-P3]^{-1} (X - P3).
    M = np.column_stack([P1 - P3, P2 - P3])  # (2, 2)
    lam12 = np.linalg.solve(M, X - P3)
    lambda1 = float(lam12[0])
    lambda2 = float(lam12[1])
    lambda3 = 1.0 - lambda1 - lambda2

    # Out-of-plane residual (measure-only).
    c = np.dot(X - P1, X - P1) - s1 * s1   # = -z^2
    z = math.sqrt(max(0.0, -c))

    return np.array([lambda1, lambda2, lambda3]), z

gdop

gdop(X: ndarray, triangle: AnchorTriangle) -> float

Geometric dilution of precision at in-plane query position X.

GDOP = sqrt(trace((H^T H)^{-1})) where H = [u1^T; u2^T; u3^T] and u_i = (X - Pi) / |X - Pi| are unit line-of-sight vectors to the three anchors. GDOP is a principled measure of in-plane trilateration reliability: it is smallest (2/sqrt(3) ~= 1.155) at the centroid of an equilateral triangle and grows without bound as the query moves far outside the triangle (lines of sight become near-parallel).

GDOP is the in-plane companion to the out-of-plane residual z: - z measures departure from the chart's plane (how non-Euclidean the local tree space is), - GDOP measures how weakly the in-plane position is constrained by the three anchor distances.

No code may branch on GDOP in the loop. It is a measure-only diagnostic. parkinson1996gdop : GDOP definition and navigation context.

Parameters:

Name	Type	Description	Default
X	(ndarray, shape(2))	In-plane position (e.g., from the Cartesian coordinates derived during the barycentric solve, or from lambda @ [P1, P2, P3]).	required
triangle	AnchorTriangle		required

Returns:

Type	Description
float	GDOP >= 0. Very large values indicate poor in-plane conditioning.

Source code in src/hifuku/anchor_triangle.py

def gdop(X: np.ndarray, triangle: AnchorTriangle) -> float:
    """
    Geometric dilution of precision at in-plane query position X.

    GDOP = sqrt(trace((H^T H)^{-1})) where H = [u1^T; u2^T; u3^T] and
    u_i = (X - Pi) / |X - Pi| are unit line-of-sight vectors to the three
    anchors.  GDOP is a principled measure of in-plane trilateration
    reliability: it is smallest (2/sqrt(3) ~= 1.155) at the centroid of an
    equilateral triangle and grows without bound as the query moves far outside
    the triangle (lines of sight become near-parallel).

    GDOP is the in-plane companion to the out-of-plane residual z:
      - z measures departure from the chart's plane (how non-Euclidean
        the local tree space is),
      - GDOP measures how weakly the in-plane position is constrained by
        the three anchor distances.

    No code may branch on GDOP in the loop.  It is a measure-only diagnostic.
    parkinson1996gdop : GDOP definition and navigation context.

    Parameters
    ----------
    X : np.ndarray, shape (2,)
        In-plane position (e.g., from the Cartesian coordinates derived
        during the barycentric solve, or from lambda @ [P1, P2, P3]).
    triangle : AnchorTriangle

    Returns
    -------
    float
        GDOP >= 0.  Very large values indicate poor in-plane conditioning.
    """
    H_rows = []
    for P in (triangle.P1, triangle.P2, triangle.P3):
        diff = X - P
        norm = float(np.linalg.norm(diff))
        if norm < 1e-12:
            return float("inf")
        H_rows.append(diff / norm)
    H = np.array(H_rows)  # (3, 2)
    HtH = H.T @ H         # (2, 2)
    det = float(HtH[0, 0] * HtH[1, 1] - HtH[0, 1] * HtH[1, 0])
    if abs(det) < 1e-30:
        return float("inf")
    # Trace of (H^T H)^{-1}
    inv_trace = (HtH[1, 1] + HtH[0, 0]) / det
    return math.sqrt(max(0.0, inv_trace))