Ordination

A survey illuminates one chart at a time. Synthesis puts the elite trees of every chart in one picture, which needs distances between trees that no single chart holds.

The tree metric is squared Euclidean by construction, so a tree can be written as a point whose ordinary distance to another point is the tree distance. split_coordinates is that construction, and it makes distances exact rather than an approximation. The whole cloud is therefore handled without subsampling: the coordinate width follows the number of distinct splits, which the taxon set bounds, not the number of trees.

faithfulness carries the scientific claim. A cloud of elite trees occupies more than two dimensions, so no plane keeps both its distances and its neighborhoods, and the measures state which was kept. neighbor_rank answers where a map is untrustworthy rather than only how much.

This module returns arrays and dataclasses. It imports no plotting library, and ordinate is a convenience whose t-SNE import is lazy. A caller may embed with any tool and hand the result to faithfulness, which scores every method on the same terms.

GOWER (1966), venna2001neighborhood, facco2017twonn. Each construction is verified in docs/theory/algebra/.

split_coordinates

split_coordinates(trees, metric) -> np.ndarray

Write each tree as a point whose ordinary distance is the tree distance.

The metric is

d^2 = (1 - w) * CBS^2 / c_ref + w * RF / r_ref

and both terms are squared Euclidean. The clade branch score is a squared Euclidean distance on the vector of split lengths, and the Robinson-Foulds count is the Hamming distance on the vector of split indicators, which is squared Euclidean because a binary difference squares to itself. A non-negative combination of squared Euclidean distances is squared Euclidean with the coordinates written side by side, each scaled by the square root of its weight:

x = concat( sqrt((1 - w) / c_ref) * lengths,
            sqrt(w / r_ref)       * indicators )

The columns run over the union of the splits of every tree, in ascending key order. A split absent from a tree contributes a length of zero and an indicator of zero, which is what makes the branch score of two trees with different split sets come out right.

Parameters:

Name	Type	Description	Default
trees	sequence of HifukuTree	Trees on one namespace, each with global_leaf_indices set.	required
metric	CBSMetric	The metric whose w, c_ref and r_ref scale the two blocks. Use :attr:hifuku.survey_io.Survey.metric to score a survey under the metric its charts were normalized against.	required

Returns:

Type	Description
ndarray	(n, 2S) float64, where S is the number of distinct splits over every tree given. The first S columns hold the scaled lengths and the rest hold the scaled indicators.

Notes

The array is dense in the splits, so its width grows with the diversity of the cloud rather than with the number of trees. A cloud whose trees share most of their splits stays narrow.

Source code in src/hifuku/ordination.py

def split_coordinates(trees, metric) -> np.ndarray:
    """Write each tree as a point whose ordinary distance is the tree distance.

    The metric is

        d^2 = (1 - w) * CBS^2 / c_ref + w * RF / r_ref

    and both terms are squared Euclidean.  The clade branch score is a squared
    Euclidean distance on the vector of split lengths, and the Robinson-Foulds
    count is the Hamming distance on the vector of split indicators, which is
    squared Euclidean because a binary difference squares to itself.  A
    non-negative combination of squared Euclidean distances is squared Euclidean
    with the coordinates written side by side, each scaled by the square root of
    its weight:

        x = concat( sqrt((1 - w) / c_ref) * lengths,
                    sqrt(w / r_ref)       * indicators )

    The columns run over the union of the splits of every tree, in ascending key
    order.  A split absent from a tree contributes a length of zero and an
    indicator of zero, which is what makes the branch score of two trees with
    different split sets come out right.

    Parameters
    ----------
    trees : sequence of HifukuTree
        Trees on one namespace, each with ``global_leaf_indices`` set.
    metric : CBSMetric
        The metric whose ``w``, ``c_ref`` and ``r_ref`` scale the two blocks.
        Use :attr:`hifuku.survey_io.Survey.metric` to score a survey under the
        metric its charts were normalized against.

    Returns
    -------
    numpy.ndarray
        ``(n, 2S)`` float64, where ``S`` is the number of distinct splits over
        every tree given.  The first ``S`` columns hold the scaled lengths and
        the rest hold the scaled indicators.

    Notes
    -----
    The array is dense in the splits, so its width grows with the diversity of
    the cloud rather than with the number of trees.  A cloud whose trees share
    most of their splits stays narrow.
    """
    maps = [split_lengths(t) for t in trees]
    union = sorted({key for m in maps for key in m})
    index = {key: col for col, key in enumerate(union)}
    n, s = len(maps), len(union)

    w = float(metric.w)
    length_scale = np.sqrt((1.0 - w) / float(metric.c_ref)) if w < 1.0 else 0.0
    indicator_scale = np.sqrt(w / float(metric.r_ref)) if w > 0.0 else 0.0

    X = np.zeros((n, 2 * s), dtype=np.float64)
    for row, m in enumerate(maps):
        for key, length in m.items():
            col = index[key]
            X[row, col] = length_scale * length
            X[row, s + col] = indicator_scale
    return X

distances

distances(source, metric=None) -> np.ndarray

The full matrix of tree distances, computed through the coordinates.

The result is exact. It is the same number :class:~hifuku.metric.branch_score.CBSMetric returns for a pair, computed for every pair at once.

Parameters:

Name	Type	Description	Default
source	Survey or sequence of HifukuTree	A survey uses its cloud, in the order :meth:Survey.cloud gives, and its own metric. A sequence of trees needs metric.	required
metric	CBSMetric	Required when source is a sequence of trees.	None

Returns:

Type	Description
ndarray	(n, n) float64, symmetric, with a zero diagonal.

Source code in src/hifuku/ordination.py

def distances(source, metric=None) -> np.ndarray:
    """The full matrix of tree distances, computed through the coordinates.

    The result is exact.  It is the same number
    :class:`~hifuku.metric.branch_score.CBSMetric` returns for a pair, computed
    for every pair at once.

    Parameters
    ----------
    source : Survey or sequence of HifukuTree
        A survey uses its cloud, in the order :meth:`Survey.cloud` gives, and
        its own metric.  A sequence of trees needs ``metric``.
    metric : CBSMetric, optional
        Required when ``source`` is a sequence of trees.

    Returns
    -------
    numpy.ndarray
        ``(n, n)`` float64, symmetric, with a zero diagonal.
    """
    from scipy.spatial.distance import pdist, squareform

    if hasattr(source, "cloud"):
        trees = [point.tree for point in source.cloud()]
        if metric is None:
            metric = source.metric
    else:
        trees = list(source)
        if metric is None:
            raise ValueError(
                "a sequence of trees needs a metric; pass one, or pass a Survey, "
                "which carries its own")
    return squareform(pdist(split_coordinates(trees, metric)))

mds

mds(D, dim: int = 2, *, report: bool = False)

Classical multidimensional scaling, also called principal coordinates.

Gower's section 3 gives the construction. Set a_ij = -d_ij^2 / 2, center the matrix by rows and columns, and take the latent roots and vectors. A vector scaled so that its sum of squares equals its root gives one column of the configuration, and the distances between the rows are the distances the input asked for.

A real configuration exists exactly when the centered matrix is positive semi-definite. The tree metric is squared Euclidean for any set of trees, so that condition always holds here and a negative root above rounding error means a defect in the implementation.

Parameters:

Name	Type	Description	Default
D	ndarray	(n, n) symmetric distances with a zero diagonal.	required
dim	int	Coordinates to keep. Fewer than the positive roots gives the best fit in that many dimensions, in the least-squares sense of Gower step (v).	2
report	bool	Return a :class:ScalingReport alongside the configuration.	False

Returns:

Type	Description
ndarray or tuple	(n, dim) float64, centered on the origin. With report set, a pair of the configuration and the report.

Source code in src/hifuku/ordination.py

def mds(D, dim: int = 2, *, report: bool = False):
    """Classical multidimensional scaling, also called principal coordinates.

    Gower's section 3 gives the construction.  Set ``a_ij = -d_ij^2 / 2``, center
    the matrix by rows and columns, and take the latent roots and vectors.  A
    vector scaled so that its sum of squares equals its root gives one column of
    the configuration, and the distances between the rows are the distances the
    input asked for.

    A real configuration exists exactly when the centered matrix is positive
    semi-definite.  The tree metric is squared Euclidean for any set of trees, so
    that condition always holds here and a negative root above rounding error
    means a defect in the implementation.

    Parameters
    ----------
    D : numpy.ndarray
        ``(n, n)`` symmetric distances with a zero diagonal.
    dim : int
        Coordinates to keep.  Fewer than the positive roots gives the best fit
        in that many dimensions, in the least-squares sense of Gower step (v).
    report : bool
        Return a :class:`ScalingReport` alongside the configuration.

    Returns
    -------
    numpy.ndarray or tuple
        ``(n, dim)`` float64, centered on the origin.  With ``report`` set, a
        pair of the configuration and the report.
    """
    D = np.asarray(D, dtype=np.float64)
    if D.ndim != 2 or D.shape[0] != D.shape[1]:
        raise ValueError(f"D must be a square matrix, got shape {D.shape}")
    n = D.shape[0]
    if dim < 1:
        raise ValueError(f"dim must be 1 or more, got {dim}")
    if not np.isfinite(D).all():
        # A non-finite entry propagates through the centering into the roots.
        # Sorting then puts it first, because a NaN compares below everything
        # and the order is reversed, so it would become the leading coordinate.
        raise ValueError("D must be finite; it holds a NaN or an infinity")

    # Gower equation 7: center a_ij = -d_ij^2 / 2 by rows and columns.  Equation
    # 8 shows this leaves every squared distance unchanged.
    A = -0.5 * (D ** 2)
    row = A.mean(axis=1, keepdims=True)
    col = A.mean(axis=0, keepdims=True)
    alpha = A - row - col + A.mean()
    alpha = 0.5 * (alpha + alpha.T)          # kill the asymmetry rounding leaves

    values, vectors = np.linalg.eigh(alpha)
    order = np.argsort(values)[::-1]
    values, vectors = values[order], vectors[:, order]

    keep = min(dim, n)
    positive = np.clip(values[:keep], 0.0, None)
    # Gower step (iii): scale each vector so its sum of squares is its root.
    X = vectors[:, :keep] * np.sqrt(positive)

    if not report:
        return X
    return X, ScalingReport(
        eigenvalues=values,
        positive_mass=float(values[values > 0].sum()),
        negative_mass=float(-values[values < 0].sum()),
        residual=float(max(values.sum() - positive.sum(), 0.0)),
    )

ScalingReport dataclass

What classical scaling found, beyond the configuration itself.

Attributes:

Name	Type	Description
eigenvalues	ndarray	The latent roots of the centered matrix, largest first.
positive_mass	float	Sum of the positive roots. The variance a real configuration holds.
negative_mass	float	Absolute sum of the negative roots. Zero for a squared Euclidean input. A value above rounding error means the distances are not of negative type, which for this metric means a defect in the implementation and not a property of the data.
residual	float	Gower step (v): the trace of the centered matrix less the roots kept. The sum of squares of the perpendiculars onto the reduced space.

Source code in src/hifuku/ordination.py

@dataclass(frozen=True)
class ScalingReport:
    """What classical scaling found, beyond the configuration itself.

    Attributes
    ----------
    eigenvalues : numpy.ndarray
        The latent roots of the centered matrix, largest first.
    positive_mass : float
        Sum of the positive roots.  The variance a real configuration holds.
    negative_mass : float
        Absolute sum of the negative roots.  Zero for a squared Euclidean input.
        A value above rounding error means the distances are not of negative
        type, which for this metric means a defect in the implementation and not
        a property of the data.
    residual : float
        Gower step (v): the trace of the centered matrix less the roots kept.
        The sum of squares of the perpendiculars onto the reduced space.
    """

    eigenvalues: np.ndarray
    positive_mass: float
    negative_mass: float
    residual: float

faithfulness

faithfulness(D, X, k: int = 15) -> Faithfulness

Score an embedding against the distances it came from.

Parameters:

Name	Type	Description	Default
D	ndarray	(n, n) distances in the original space.	required
X	ndarray	(n, dim) embedding to score. Any tool may have produced it.	required
k	int	Neighborhood size. Venna and Kaski state the scaling for k < n / 2, which this function enforces.	15

Returns:

Type	Description
Faithfulness

Source code in src/hifuku/ordination.py

def faithfulness(D, X, k: int = 15) -> Faithfulness:
    """Score an embedding against the distances it came from.

    Parameters
    ----------
    D : numpy.ndarray
        ``(n, n)`` distances in the original space.
    X : numpy.ndarray
        ``(n, dim)`` embedding to score.  Any tool may have produced it.
    k : int
        Neighborhood size.  Venna and Kaski state the scaling for ``k < n / 2``,
        which this function enforces.

    Returns
    -------
    Faithfulness
    """
    from scipy.spatial.distance import pdist, squareform

    D = np.asarray(D, dtype=np.float64)
    X = np.asarray(X, dtype=np.float64)
    n = D.shape[0]
    if X.shape[0] != n:
        raise ValueError(
            f"the embedding has {X.shape[0]} rows and the distances {n}")
    if not 1 <= k < n / 2:
        raise ValueError(
            f"k must be at least 1 and below n / 2 = {n / 2}, got {k}; the "
            "scaling of the measures is stated for that range only")

    rank_true = _ranks(D)
    rank_map = _ranks(squareform(pdist(X)))

    # The k nearest in each space, self excluded: ranks 1 .. k.
    near_true = (rank_true >= 1) & (rank_true <= k)
    near_map = (rank_map >= 1) & (rank_map <= k)

    entered = near_map & ~near_true          # U_k: on the map but not in truth
    left = near_true & ~near_map             # V_k: in truth but not on the map

    scale = 2.0 / (n * k * (2 * n - 3 * k - 1))
    trust = 1.0 - scale * float((rank_true[entered] - k).sum())
    cont = 1.0 - scale * float((rank_map[left] - k).sum())

    # Where the map is untrustworthy, not only how much.
    neighbor_rank = np.where(near_map, rank_true, 0).sum(axis=1) / k
    false_rate = float(np.count_nonzero(entered)) / (n * k)

    return Faithfulness(
        k=int(k),
        trustworthiness=float(trust),
        continuity=float(cont),
        neighbor_rank=neighbor_rank.astype(np.float64),
        ideal_rank=(k + 1) / 2.0,
        false_neighbor_rate=false_rate,
    )

Faithfulness dataclass

How much of a cloud's neighborhood structure a map keeps.

A cloud of elite trees occupies more than two dimensions, so no plane keeps both its distances and its neighborhoods. These measures state which was kept and by how much.

trustworthiness and continuity are the two measures of Venna and Kaski, equations 1 and 2. The penalty is weighted by rank rather than counted, so a point drawn onto the map from true rank 1000 costs far more than one from rank 16. A map that moves nothing scores 1 and the worst possible map scores 0; the scaling is derived in docs/theory/algebra/venna2001neighborhood-trustworthiness.py.

neighbor_rank is not from the paper. It answers where a map is untrustworthy rather than how much: for each point, the mean true rank of the neighbors the map puts around it.

Attributes:

Name	Type	Description
k	int	Neighborhood size the measures used.
trustworthiness	float	Penalizes points the map drew in from far away. A reader has no way to tell an introduced proximity from a real one, which is why Venna and Kaski call this the more damaging error.
continuity	float	Penalizes true neighbors the map pushed away. The paper calls this preservation of the original neighborhoods.
neighbor_rank	ndarray	(n,). Per point, the mean true rank of its map neighbors.
ideal_rank	float	(k + 1) / 2, the value a faithful map returns.
false_neighbor_rate	float	Share of map neighbors whose true rank is above k.

Source code in src/hifuku/ordination.py

@dataclass(frozen=True)
class Faithfulness:
    """How much of a cloud's neighborhood structure a map keeps.

    A cloud of elite trees occupies more than two dimensions, so no plane keeps
    both its distances and its neighborhoods.  These measures state which was
    kept and by how much.

    ``trustworthiness`` and ``continuity`` are the two measures of Venna and
    Kaski, equations 1 and 2.  The penalty is weighted by rank rather than
    counted, so a point drawn onto the map from true rank 1000 costs far more
    than one from rank 16.  A map that moves nothing scores 1 and the worst
    possible map scores 0; the scaling is derived in
    ``docs/theory/algebra/venna2001neighborhood-trustworthiness.py``.

    ``neighbor_rank`` is not from the paper.  It answers where a map is
    untrustworthy rather than how much: for each point, the mean true rank of
    the neighbors the map puts around it.

    Attributes
    ----------
    k : int
        Neighborhood size the measures used.
    trustworthiness : float
        Penalizes points the map drew in from far away.  A reader has no way to
        tell an introduced proximity from a real one, which is why Venna and
        Kaski call this the more damaging error.
    continuity : float
        Penalizes true neighbors the map pushed away.  The paper calls this
        preservation of the original neighborhoods.
    neighbor_rank : numpy.ndarray
        ``(n,)``.  Per point, the mean true rank of its map neighbors.
    ideal_rank : float
        ``(k + 1) / 2``, the value a faithful map returns.
    false_neighbor_rate : float
        Share of map neighbors whose true rank is above ``k``.
    """

    k: int
    trustworthiness: float
    continuity: float
    neighbor_rank: np.ndarray
    ideal_rank: float
    false_neighbor_rate: float

intrinsic_dimension

intrinsic_dimension(D, discard: float = 0.1) -> float

Estimate the intrinsic dimension from first and second neighbors.

The TWO-NN estimator of Facco et al. For each point take the distances to its first and second nearest neighbors and form mu = r2 / r1. For a locally uniform density the distribution of mu depends on the intrinsic dimension and not on the density, which cancels, so

-log(1 - F(mu)) / log(mu) = d.

The points (log mu, -log(1 - F)) therefore lie on a line through the origin whose slope is the dimension. F comes from the sample: sort mu upward and set F(mu_i) = i / n.

Using only the two nearest neighbors is what makes the estimate usable on a cloud of elite trees, which is not uniformly distributed. Local uniformity is needed only at the scale of the second neighbor.

Parameters:

Name	Type	Description	Default
D	ndarray	(n, n) distances.	required
discard	float	Fraction of the largest mu to drop before fitting. The paper drops a tenth, because in a heavy-tailed cloud a few very large ratios dominate the slope. On well-behaved data the estimate barely moves.	0.1

Returns:

Type	Description
float	The estimated dimension.

Notes

Report this alongside the classical scaling spectrum rather than instead of it. This measures the local dimension of the manifold, and the spectrum measures how many linear axes the distances occupy. A disagreement between them says something about the cloud.

Source code in src/hifuku/ordination.py

def intrinsic_dimension(D, discard: float = 0.1) -> float:
    """Estimate the intrinsic dimension from first and second neighbors.

    The TWO-NN estimator of Facco et al.  For each point take the distances to
    its first and second nearest neighbors and form ``mu = r2 / r1``.  For a
    locally uniform density the distribution of ``mu`` depends on the intrinsic
    dimension and not on the density, which cancels, so

        -log(1 - F(mu)) / log(mu) = d.

    The points ``(log mu, -log(1 - F))`` therefore lie on a line through the
    origin whose slope is the dimension.  ``F`` comes from the sample: sort
    ``mu`` upward and set ``F(mu_i) = i / n``.

    Using only the two nearest neighbors is what makes the estimate usable on a
    cloud of elite trees, which is not uniformly distributed.  Local uniformity
    is needed only at the scale of the second neighbor.

    Parameters
    ----------
    D : numpy.ndarray
        ``(n, n)`` distances.
    discard : float
        Fraction of the largest ``mu`` to drop before fitting.  The paper drops
        a tenth, because in a heavy-tailed cloud a few very large ratios
        dominate the slope.  On well-behaved data the estimate barely moves.

    Returns
    -------
    float
        The estimated dimension.

    Notes
    -----
    Report this alongside the classical scaling spectrum rather than instead of
    it.  This measures the local dimension of the manifold, and the spectrum
    measures how many linear axes the distances occupy.  A disagreement between
    them says something about the cloud.
    """
    D = np.asarray(D, dtype=np.float64)
    n = D.shape[0]
    if n < 3:
        raise ValueError(f"the estimator needs at least 3 points, got {n}")
    if not 0.0 <= discard < 1.0:
        raise ValueError(f"discard must be in [0, 1), got {discard}")

    # The two smallest distances from each point, with the self distance out.
    part = np.partition(D + np.diag(np.full(n, np.inf)), 1, axis=1)[:, :2]
    r1, r2 = part[:, 0], part[:, 1]
    usable = r1 > 0.0                       # coincident points give no ratio
    mu = np.sort(r2[usable] / r1[usable])
    if mu.size < 2:
        raise ValueError("too few distinct neighbor ratios to fit a dimension")

    # The cumulate runs over every ratio, which is what fixes its scale.  Points
    # with mu of exactly 1 are dropped from the fit below and not from here:
    # removing them first would change the count and shift F for all the rest.
    F = np.arange(1, mu.size + 1, dtype=np.float64) / mu.size
    keep = max(2, int(round(mu.size * (1.0 - discard))))
    keep = min(keep, mu.size - 1)           # F = 1 makes -log(1 - F) infinite
    x = np.log(mu[:keep])
    y = -np.log1p(-F[:keep])

    # A ratio of 1 sits at the origin of the fit and says nothing about a slope.
    on_the_line = x > 0.0
    if on_the_line.sum() < 2:
        raise ValueError("too few distinct neighbor ratios to fit a dimension")
    x, y = x[on_the_line], y[on_the_line]

    return float((x @ y) / (x @ x))

project_landmarks

project_landmarks(coords, X, landmarks, sigma=None) -> np.ndarray

Place points that are not in the cloud into an existing embedding.

An anchor tree is a coordinate reference, not a member of the cloud, so it has no row in the embedding. This puts one on the map by averaging the embedded positions of the cloud, weighted by closeness in the original space.

Parameters:

Name	Type	Description	Default
coords	ndarray	(n, dim) embedding of the cloud.	required
X	ndarray	(n, p) cloud in the original space, for example the output of :func:split_coordinates.	required
landmarks	ndarray	(m, p) points to place, in the same space as X.	required
sigma	float	Width of the weighting kernel. Defaults to the median distance from a landmark to its nearest cloud point, which follows the scale of the data. A small width places a landmark on its nearest cloud point, and a large one pulls every landmark to the centroid.	None

Returns:

Type	Description
ndarray	(m, dim) coordinates in the same frame as coords.

Source code in src/hifuku/ordination.py

def project_landmarks(coords, X, landmarks, sigma=None) -> np.ndarray:
    """Place points that are not in the cloud into an existing embedding.

    An anchor tree is a coordinate reference, not a member of the cloud, so it
    has no row in the embedding.  This puts one on the map by averaging the
    embedded positions of the cloud, weighted by closeness in the original
    space.

    Parameters
    ----------
    coords : numpy.ndarray
        ``(n, dim)`` embedding of the cloud.
    X : numpy.ndarray
        ``(n, p)`` cloud in the original space, for example the output of
        :func:`split_coordinates`.
    landmarks : numpy.ndarray
        ``(m, p)`` points to place, in the same space as ``X``.
    sigma : float, optional
        Width of the weighting kernel.  Defaults to the median distance from a
        landmark to its nearest cloud point, which follows the scale of the
        data.  A small width places a landmark on its nearest cloud point, and a
        large one pulls every landmark to the centroid.

    Returns
    -------
    numpy.ndarray
        ``(m, dim)`` coordinates in the same frame as ``coords``.
    """
    from scipy.spatial.distance import cdist

    coords = np.asarray(coords, dtype=np.float64)
    X = np.asarray(X, dtype=np.float64)
    landmarks = np.atleast_2d(np.asarray(landmarks, dtype=np.float64))
    if X.shape[0] != coords.shape[0]:
        raise ValueError(
            f"the cloud has {X.shape[0]} rows and the embedding "
            f"{coords.shape[0]}")

    D = cdist(landmarks, X)
    if sigma is None:
        sigma = float(np.median(D.min(axis=1)))
    if sigma <= 0.0:
        sigma = np.finfo(np.float64).tiny

    # Shift by the row minimum before the exponential, so a landmark far from
    # every cloud point still gets usable weights instead of underflowing.
    W = np.exp(-0.5 * ((D - D.min(axis=1, keepdims=True)) / sigma) ** 2)
    total = W.sum(axis=1, keepdims=True)
    W = np.where(total > 0.0, W / np.where(total > 0.0, total, 1.0), 0.0)
    return W @ coords

ordinate

ordinate(D, method: str = 'tsne', **kw) -> np.ndarray

Embed a distance matrix with an outside tool.

Barnes-Hut t-SNE is the one step this package does not own, so this is a convenience behind the analysis extra with a lazy import, in the same way :meth:Survey.frame reaches for pandas. A caller may skip it, embed with any tool, and hand the result to :func:faithfulness, which scores every method on the same terms.

Parameters:

Name	Type	Description	Default
D	ndarray	(n, n) distances.	required
method	str	"tsne" or "mds". "mds" runs :func:mds and needs nothing beyond the hard dependencies.	'tsne'
**kw		Passed to the underlying estimator. random_state fixes the seed.	{}

Returns:

Type	Description
ndarray	(n, dim) coordinates.

Source code in src/hifuku/ordination.py

def ordinate(D, method: str = "tsne", **kw) -> np.ndarray:
    """Embed a distance matrix with an outside tool.

    Barnes-Hut t-SNE is the one step this package does not own, so this is a
    convenience behind the ``analysis`` extra with a lazy import, in the same
    way :meth:`Survey.frame` reaches for pandas.  A caller may skip it, embed
    with any tool, and hand the result to :func:`faithfulness`, which scores
    every method on the same terms.

    Parameters
    ----------
    D : numpy.ndarray
        ``(n, n)`` distances.
    method : str
        ``"tsne"`` or ``"mds"``.  ``"mds"`` runs :func:`mds` and needs nothing
        beyond the hard dependencies.
    **kw
        Passed to the underlying estimator.  ``random_state`` fixes the seed.

    Returns
    -------
    numpy.ndarray
        ``(n, dim)`` coordinates.
    """
    D = np.asarray(D, dtype=np.float64)
    if method == "mds":
        return mds(D, dim=int(kw.pop("dim", 2)))
    if method != "tsne":
        raise ValueError(f"unknown method {method!r}; use 'tsne' or 'mds'")
    try:
        from sklearn.manifold import TSNE
    except ImportError as exc:                              # pragma: no cover
        raise ImportError(
            "ordinate(method='tsne') needs scikit-learn, which is in the "
            "'analysis' extra: pip install 'hifuku[analysis]'.  Or embed with "
            "any tool and score the result with faithfulness()."
        ) from exc
    kw.setdefault("metric", "precomputed")
    kw.setdefault("init", "random")
    return np.asarray(TSNE(**kw).fit_transform(D), dtype=np.float64)