Chart Selection

A chart is three anchor trees embedded in a plane, and a survey runs genes on it. A gene charted on a poor triangle is measured through a badly conditioned coordinate frame, so the choice of charts matters.

Measurement is separated from policy. chart_table measures the whole candidate space once: the pairwise anchor graph, the quality of every triple, and how flat each gene sits on each triple. ChartTable.select is one policy over that table, and a caller who wants another may filter Q, read a fit row, or hand a chosen set straight to assign_genes.

Measurement needs no survey. Locating a gene's own tree on a candidate chart and reading its out-of-plane residual says how flat that chart is for that gene, and that is a distance calculation.

reconstruction_error scores a candidate set by re-embedding the pairs its triangles cover and measuring the result against every pair, including the ones no chart covers. Scoring only the covered pairs would reward a thin set for having fewer terms.

The selection thresholds default to the values AnchorTriangle.from_anchors enforces, so a chart this policy chooses is one the survey will accept.

Leeuw & Mair (2009), verified in docs/theory/algebra/deleeuw2009smacof-majorization.py.

chart_table

chart_table(trees, labels, metric, genes=None, gene_trees=None) -> ChartTable

Measure every chart a pool of anchors admits.

Parameters:

Name	Type	Description	Default
trees	sequence of HifukuTree	The anchor pool, on one namespace with global_leaf_indices set.	required
labels	sequence of str	One label per anchor.	required
metric	TreeMetric	The metric the charts will use.	required
genes	sequence of str	Gene labels to fit. Defaults to the anchor labels, which is the usual case: each alignment supplies both an anchor and a gene to survey.	None
gene_trees	sequence of HifukuTree	The tree of each gene. Defaults to the anchor trees.	None

Returns:

Type	Description
ChartTable

Source code in src/hifuku/charts.py

def chart_table(trees, labels, metric, genes=None, gene_trees=None) -> ChartTable:
    """Measure every chart a pool of anchors admits.

    Parameters
    ----------
    trees : sequence of HifukuTree
        The anchor pool, on one namespace with ``global_leaf_indices`` set.
    labels : sequence of str
        One label per anchor.
    metric : TreeMetric
        The metric the charts will use.
    genes : sequence of str, optional
        Gene labels to fit.  Defaults to the anchor labels, which is the usual
        case: each alignment supplies both an anchor and a gene to survey.
    gene_trees : sequence of HifukuTree, optional
        The tree of each gene.  Defaults to the anchor trees.

    Returns
    -------
    ChartTable
    """
    trees = list(trees)
    labels = list(labels)
    if len(trees) != len(labels):
        raise ValueError(
            f"{len(trees)} trees and {len(labels)} labels; they must match")
    if len(trees) < 3:
        raise ValueError(f"a chart needs three anchors, and the pool holds "
                         f"{len(trees)}")

    if genes is None:
        genes, gene_trees = list(labels), list(trees)
    else:
        genes = list(genes)
        if gene_trees is None:
            index = {lab: t for lab, t in zip(labels, trees)}
            gene_trees = [index[g] for g in genes]
        gene_trees = list(gene_trees)

    D, B = anchor_distances(trees, metric)
    triples = np.array(list(itertools.combinations(range(len(trees)), 3)),
                       dtype=int)
    Q = np.array([triple_quality(D, a, b, c) for a, b, c in triples])

    fit = np.full((len(genes), triples.shape[0]), np.nan, dtype=np.float64)
    for column, (a, b, c) in enumerate(triples):
        if not np.isfinite(Q[column]):
            continue
        triangle = _triangle_for(trees, D, B, a, b, c, metric)
        for row, tree in enumerate(gene_trees):
            fit[row, column] = _chart_misfit(tree, triangle, metric)

    return ChartTable(labels=labels, D=D, B=B, triples=triples, Q=Q,
                      genes=genes, fit=fit)

ChartTable dataclass

Every candidate chart a pool of anchors admits, measured.

Attributes:

Name	Type	Description
labels	list[str]	The N anchor labels, in the order the arrays index them.
D	ndarray	(N, N) pairwise anchor distance.
B	ndarray	(N, N) per-pair B_frac, the length-overlap quality. A value near zero means the metric can no longer tell that pair apart.
triples	ndarray	(T, 3) int, canonical with a < b < c.
Q	ndarray	(T,) triangle quality, or NaN where the triple admits no planar triangle.
genes	list[str]	The M gene labels.
fit	ndarray	(M, T) misfit of gene m on triple t, as a length in the units of the metric: the out-of-plane residual and the anisotropy of the triangle combined. Small is good. NaN where the triple is not a chart. See :func:_chart_misfit.

Notes

The triples are canonical rather than a dense (N, N, N) tensor. The dense form is symmetric under permutation of its three axes and undefined on the diagonal, so it is six times redundant. That is irrelevant at N = 10 and decides feasibility later: at N = M = 100 the fit array is 800 MB dense against 130 MB canonical. :meth:as_tensor builds the dense view when N is small enough to slice interactively.

Source code in src/hifuku/charts.py

@dataclass(frozen=True)
class ChartTable:
    """Every candidate chart a pool of anchors admits, measured.

    Attributes
    ----------
    labels : list[str]
        The ``N`` anchor labels, in the order the arrays index them.
    D : numpy.ndarray
        ``(N, N)`` pairwise anchor distance.
    B : numpy.ndarray
        ``(N, N)`` per-pair ``B_frac``, the length-overlap quality.  A value
        near zero means the metric can no longer tell that pair apart.
    triples : numpy.ndarray
        ``(T, 3)`` int, canonical with ``a < b < c``.
    Q : numpy.ndarray
        ``(T,)`` triangle quality, or NaN where the triple admits no planar
        triangle.
    genes : list[str]
        The ``M`` gene labels.
    fit : numpy.ndarray
        ``(M, T)`` misfit of gene ``m`` on triple ``t``, as a length in the
        units of the metric: the out-of-plane residual and the anisotropy of
        the triangle combined.  Small is good.  NaN where the triple is not a
        chart.  See :func:`_chart_misfit`.

    Notes
    -----
    The triples are canonical rather than a dense ``(N, N, N)`` tensor.  The
    dense form is symmetric under permutation of its three axes and undefined on
    the diagonal, so it is six times redundant.  That is irrelevant at ``N = 10``
    and decides feasibility later: at ``N = M = 100`` the ``fit`` array is
    800 MB dense against 130 MB canonical.  :meth:`as_tensor` builds the dense
    view when ``N`` is small enough to slice interactively.
    """

    labels: list
    D: np.ndarray
    B: np.ndarray
    triples: np.ndarray
    Q: np.ndarray
    genes: list
    fit: np.ndarray

    def as_tensor(self) -> np.ndarray:
        """The dense ``(N, N, N)`` view of ``Q``, for interactive slicing.

        Symmetric in its three axes and NaN wherever two indices repeat.  Build
        this only for a small pool: it holds ``N ** 3`` floats where the
        canonical form holds ``T``.
        """
        n = len(self.labels)
        dense = np.full((n, n, n), np.nan, dtype=np.float64)
        for row, (a, b, c) in enumerate(self.triples):
            for p in itertools.permutations((int(a), int(b), int(c))):
                dense[p] = self.Q[row]
        return dense

    def containing(self, label) -> np.ndarray:
        """Indices of the triples that use one anchor, by label or index."""
        index = self.labels.index(label) if isinstance(label, str) else int(label)
        return np.nonzero(np.any(self.triples == index, axis=1))[0]

    def reconstruction_error(self, triples) -> float:
        """How much of the anchor geometry a candidate set of charts keeps.

        A set of triangles covers some pairs of anchors and leaves others
        uncovered.  This re-embeds the anchors from the covered pairs alone,
        then scores that configuration against *every* pairwise distance,
        including the ones no chart covers.

        Scoring over every pair is the point.  Stress sums only over pairs with
        a positive weight, so dropping a pair deletes a non-negative term and a
        thinner candidate set scores lower for free.  That decrease measures the
        edge count, not the embedding.  The check is in
        ``docs/theory/algebra/deleeuw2009smacof-majorization.py``.

        This is a set function, not a table column.  What a triangle adds
        depends on which triangles are already chosen, so it cannot live in
        ``Q`` or ``fit``.

        The score is not monotone in the candidate set.  Adding a triangle
        usually lowers it and sometimes raises it, by a small amount: the
        re-embedding minimizes stress over the retained pairs, while the score
        measures stress over every pair, so a new pair changes the function
        being minimized and can move the result to a different local minimum.
        Measured on a pool of seven anchors, a rise appeared in four of six
        random orders, the largest being 4e-3 against scores near 0.1.
        :meth:`select` therefore takes the best candidate each round and stops
        on the size of the improvement, rather than assuming each step helps.

        Parameters
        ----------
        triples : array-like of int
            Row indices into :attr:`triples`.

        Returns
        -------
        float
            Normalized error over all pairs, zero when the retained set
            reproduces the whole geometry.
        """
        rows = np.atleast_1d(np.asarray(triples, dtype=int))
        n = len(self.labels)
        W = np.zeros((n, n), dtype=np.float64)
        for a, b, c in self.triples[rows]:
            for i, j in ((a, b), (a, c), (b, c)):
                W[i, j] = W[j, i] = 1.0

        X = smacof(self.D, W, dim=2)
        d = _pairwise(X)
        upper = np.triu_indices(n, k=1)
        total = float((self.D[upper] ** 2).sum())
        if total == 0.0:
            return 0.0
        return float(np.sqrt(((self.D[upper] - d[upper]) ** 2).sum() / total))

    def select(self, *, tol: float = 0.01, q_min: float = None,
               b_frac_min: float = None) -> np.ndarray:
        """Choose a covering set of charts, then add while it pays.

        The policy has two stages.  Cover every anchor, so that each one can
        host a chart, taking the triangle that lowers the reconstruction error
        most at each step.  Then keep adding triangles while each removes more
        than ``tol`` of the error.

        This is one policy over an inspectable table.  A caller who wants
        another may filter :attr:`Q`, read a :attr:`fit` row, or pass a chosen
        set straight to :func:`assign_genes`.

        Parameters
        ----------
        tol : float
            The least error a further chart must remove to earn its place.
        q_min : float
            Discard a triple whose triangle quality falls below this.  A triple
            with no quality at all is always discarded.  Defaults to the
            threshold the survey itself enforces.
        b_frac_min : float
            Discard a triple with a pair whose length overlap falls below this.
            Defaults to the threshold the survey itself enforces.

        Returns
        -------
        numpy.ndarray
            Row indices into :attr:`triples`.

        Notes
        -----
        The two defaults are the thresholds ``AnchorTriangle.from_anchors``
        applies, so a chart this policy chooses is one the survey will accept.
        A looser setting here emits a plan the survey then refuses.
        """
        from hifuku.params import SelectionParams

        defaults = SelectionParams()
        q_min = defaults.q_min if q_min is None else q_min
        b_frac_min = defaults.b_frac_min if b_frac_min is None else b_frac_min

        keeps_quality = np.isfinite(self.Q) & (self.Q >= q_min)
        rows = np.arange(self.triples.shape[0])
        keeps_overlap = np.array([
            min(self.B[a, b], self.B[a, c], self.B[b, c]) >= b_frac_min
            for a, b, c in self.triples])
        usable = rows[keeps_quality & keeps_overlap]
        if usable.size == 0:
            raise ValueError(
                "no triple of these anchors makes a usable chart: none has "
                f"triangle quality at least {q_min} with every pair overlapping "
                f"at least {b_frac_min}")

        n = len(self.labels)
        chosen: list[int] = []
        uncovered = set(range(n))
        while uncovered:
            candidates = [row for row in usable
                          if uncovered & set(self.triples[row].tolist())]
            if not candidates:
                break
            best = min(candidates,
                       key=lambda row: (
                           -len(uncovered & set(self.triples[row].tolist())),
                           self.reconstruction_error(chosen + [row])))
            chosen.append(best)
            uncovered -= set(self.triples[best].tolist())

        error = self.reconstruction_error(chosen)
        while True:
            rest = [row for row in usable if row not in chosen]
            if not rest:
                break
            scores = [(self.reconstruction_error(chosen + [row]), row)
                      for row in rest]
            best_error, best_row = min(scores)
            if error - best_error <= tol:
                break
            chosen.append(best_row)
            error = best_error

        return np.array(sorted(chosen), dtype=int)

as_tensor

as_tensor() -> np.ndarray

The dense (N, N, N) view of Q, for interactive slicing.

Symmetric in its three axes and NaN wherever two indices repeat. Build this only for a small pool: it holds N ** 3 floats where the canonical form holds T.

Source code in src/hifuku/charts.py

def as_tensor(self) -> np.ndarray:
    """The dense ``(N, N, N)`` view of ``Q``, for interactive slicing.

    Symmetric in its three axes and NaN wherever two indices repeat.  Build
    this only for a small pool: it holds ``N ** 3`` floats where the
    canonical form holds ``T``.
    """
    n = len(self.labels)
    dense = np.full((n, n, n), np.nan, dtype=np.float64)
    for row, (a, b, c) in enumerate(self.triples):
        for p in itertools.permutations((int(a), int(b), int(c))):
            dense[p] = self.Q[row]
    return dense

containing

containing(label) -> np.ndarray

Indices of the triples that use one anchor, by label or index.

Source code in src/hifuku/charts.py

def containing(self, label) -> np.ndarray:
    """Indices of the triples that use one anchor, by label or index."""
    index = self.labels.index(label) if isinstance(label, str) else int(label)
    return np.nonzero(np.any(self.triples == index, axis=1))[0]

reconstruction_error

reconstruction_error(triples) -> float

How much of the anchor geometry a candidate set of charts keeps.

A set of triangles covers some pairs of anchors and leaves others uncovered. This re-embeds the anchors from the covered pairs alone, then scores that configuration against every pairwise distance, including the ones no chart covers.

Scoring over every pair is the point. Stress sums only over pairs with a positive weight, so dropping a pair deletes a non-negative term and a thinner candidate set scores lower for free. That decrease measures the edge count, not the embedding. The check is in docs/theory/algebra/deleeuw2009smacof-majorization.py.

This is a set function, not a table column. What a triangle adds depends on which triangles are already chosen, so it cannot live in Q or fit.

The score is not monotone in the candidate set. Adding a triangle usually lowers it and sometimes raises it, by a small amount: the re-embedding minimizes stress over the retained pairs, while the score measures stress over every pair, so a new pair changes the function being minimized and can move the result to a different local minimum. Measured on a pool of seven anchors, a rise appeared in four of six random orders, the largest being 4e-3 against scores near 0.1. :meth:select therefore takes the best candidate each round and stops on the size of the improvement, rather than assuming each step helps.

Parameters:

Name	Type	Description	Default
triples	array-like of int	Row indices into :attr:triples.	required

Returns:

Type	Description
float	Normalized error over all pairs, zero when the retained set reproduces the whole geometry.

Source code in src/hifuku/charts.py

def reconstruction_error(self, triples) -> float:
    """How much of the anchor geometry a candidate set of charts keeps.

    A set of triangles covers some pairs of anchors and leaves others
    uncovered.  This re-embeds the anchors from the covered pairs alone,
    then scores that configuration against *every* pairwise distance,
    including the ones no chart covers.

    Scoring over every pair is the point.  Stress sums only over pairs with
    a positive weight, so dropping a pair deletes a non-negative term and a
    thinner candidate set scores lower for free.  That decrease measures the
    edge count, not the embedding.  The check is in
    ``docs/theory/algebra/deleeuw2009smacof-majorization.py``.

    This is a set function, not a table column.  What a triangle adds
    depends on which triangles are already chosen, so it cannot live in
    ``Q`` or ``fit``.

    The score is not monotone in the candidate set.  Adding a triangle
    usually lowers it and sometimes raises it, by a small amount: the
    re-embedding minimizes stress over the retained pairs, while the score
    measures stress over every pair, so a new pair changes the function
    being minimized and can move the result to a different local minimum.
    Measured on a pool of seven anchors, a rise appeared in four of six
    random orders, the largest being 4e-3 against scores near 0.1.
    :meth:`select` therefore takes the best candidate each round and stops
    on the size of the improvement, rather than assuming each step helps.

    Parameters
    ----------
    triples : array-like of int
        Row indices into :attr:`triples`.

    Returns
    -------
    float
        Normalized error over all pairs, zero when the retained set
        reproduces the whole geometry.
    """
    rows = np.atleast_1d(np.asarray(triples, dtype=int))
    n = len(self.labels)
    W = np.zeros((n, n), dtype=np.float64)
    for a, b, c in self.triples[rows]:
        for i, j in ((a, b), (a, c), (b, c)):
            W[i, j] = W[j, i] = 1.0

    X = smacof(self.D, W, dim=2)
    d = _pairwise(X)
    upper = np.triu_indices(n, k=1)
    total = float((self.D[upper] ** 2).sum())
    if total == 0.0:
        return 0.0
    return float(np.sqrt(((self.D[upper] - d[upper]) ** 2).sum() / total))

select

select(*, tol: float = 0.01, q_min: float = None, b_frac_min: float = None) -> np.ndarray

Choose a covering set of charts, then add while it pays.

The policy has two stages. Cover every anchor, so that each one can host a chart, taking the triangle that lowers the reconstruction error most at each step. Then keep adding triangles while each removes more than tol of the error.

This is one policy over an inspectable table. A caller who wants another may filter :attr:Q, read a :attr:fit row, or pass a chosen set straight to :func:assign_genes.

Parameters:

Name	Type	Description	Default
tol	float	The least error a further chart must remove to earn its place.	0.01
q_min	float	Discard a triple whose triangle quality falls below this. A triple with no quality at all is always discarded. Defaults to the threshold the survey itself enforces.	None
b_frac_min	float	Discard a triple with a pair whose length overlap falls below this. Defaults to the threshold the survey itself enforces.	None

Returns:

Type	Description
ndarray	Row indices into :attr:triples.

Notes

The two defaults are the thresholds AnchorTriangle.from_anchors applies, so a chart this policy chooses is one the survey will accept. A looser setting here emits a plan the survey then refuses.

Source code in src/hifuku/charts.py

def select(self, *, tol: float = 0.01, q_min: float = None,
           b_frac_min: float = None) -> np.ndarray:
    """Choose a covering set of charts, then add while it pays.

    The policy has two stages.  Cover every anchor, so that each one can
    host a chart, taking the triangle that lowers the reconstruction error
    most at each step.  Then keep adding triangles while each removes more
    than ``tol`` of the error.

    This is one policy over an inspectable table.  A caller who wants
    another may filter :attr:`Q`, read a :attr:`fit` row, or pass a chosen
    set straight to :func:`assign_genes`.

    Parameters
    ----------
    tol : float
        The least error a further chart must remove to earn its place.
    q_min : float
        Discard a triple whose triangle quality falls below this.  A triple
        with no quality at all is always discarded.  Defaults to the
        threshold the survey itself enforces.
    b_frac_min : float
        Discard a triple with a pair whose length overlap falls below this.
        Defaults to the threshold the survey itself enforces.

    Returns
    -------
    numpy.ndarray
        Row indices into :attr:`triples`.

    Notes
    -----
    The two defaults are the thresholds ``AnchorTriangle.from_anchors``
    applies, so a chart this policy chooses is one the survey will accept.
    A looser setting here emits a plan the survey then refuses.
    """
    from hifuku.params import SelectionParams

    defaults = SelectionParams()
    q_min = defaults.q_min if q_min is None else q_min
    b_frac_min = defaults.b_frac_min if b_frac_min is None else b_frac_min

    keeps_quality = np.isfinite(self.Q) & (self.Q >= q_min)
    rows = np.arange(self.triples.shape[0])
    keeps_overlap = np.array([
        min(self.B[a, b], self.B[a, c], self.B[b, c]) >= b_frac_min
        for a, b, c in self.triples])
    usable = rows[keeps_quality & keeps_overlap]
    if usable.size == 0:
        raise ValueError(
            "no triple of these anchors makes a usable chart: none has "
            f"triangle quality at least {q_min} with every pair overlapping "
            f"at least {b_frac_min}")

    n = len(self.labels)
    chosen: list[int] = []
    uncovered = set(range(n))
    while uncovered:
        candidates = [row for row in usable
                      if uncovered & set(self.triples[row].tolist())]
        if not candidates:
            break
        best = min(candidates,
                   key=lambda row: (
                       -len(uncovered & set(self.triples[row].tolist())),
                       self.reconstruction_error(chosen + [row])))
        chosen.append(best)
        uncovered -= set(self.triples[best].tolist())

    error = self.reconstruction_error(chosen)
    while True:
        rest = [row for row in usable if row not in chosen]
        if not rest:
            break
        scores = [(self.reconstruction_error(chosen + [row]), row)
                  for row in rest]
        best_error, best_row = min(scores)
        if error - best_error <= tol:
            break
        chosen.append(best_row)
        error = best_error

    return np.array(sorted(chosen), dtype=int)

assign_genes

assign_genes(table: ChartTable, triples) -> list

Put every gene on the chart it sits flattest on.

Parameters:

Name	Type	Description	Default
table	ChartTable		required
triples	array-like of int	The chosen charts, as row indices into table.triples.	required

Returns:

Type	Description
list[ChartSpec]	Ready for :class:~hifuku.survey_io.SurveyPlan. A chart that no gene chose is left out, so the plan holds no chart with nothing to survey.

Source code in src/hifuku/charts.py

def assign_genes(table: ChartTable, triples) -> list:
    """Put every gene on the chart it sits flattest on.

    Parameters
    ----------
    table : ChartTable
    triples : array-like of int
        The chosen charts, as row indices into ``table.triples``.

    Returns
    -------
    list[ChartSpec]
        Ready for :class:`~hifuku.survey_io.SurveyPlan`.  A chart that no gene
        chose is left out, so the plan holds no chart with nothing to survey.
    """
    from hifuku.survey_io import ChartSpec

    rows = np.atleast_1d(np.asarray(triples, dtype=int))
    if rows.size == 0:
        raise ValueError("no charts were given to assign genes to")

    placed: dict[int, list] = {int(row): [] for row in rows}
    for gene_index, gene in enumerate(table.genes):
        candidates = table.fit[gene_index, rows]
        if np.all(np.isnan(candidates)):
            raise ValueError(
                f"gene {gene!r} does not fit any of the chosen charts")
        placed[int(rows[int(np.nanargmin(candidates))])].append(gene)

    return [ChartSpec(anchors=tuple(table.labels[i] for i in table.triples[row]),
                      genes=names)
            for row, names in placed.items() if names]

smacof

smacof(D, W, dim: int = 2, *, seed: int = 0, iters: int = _SMACOF_ITERS, tol: float = _SMACOF_TOL, history: bool = False)

Weighted multidimensional scaling by majorization.

Minimizes sum_{i<j} w_ij (delta_ij - d_ij(X))^2. A weight of zero drops a pair, which is the missing-value structure the paper describes and the way chart selection expresses a candidate set: a pair of anchors that no chosen triangle covers is simply not fitted.

Weighted stress has no closed form, so each step builds a quadratic surrogate that sits above stress and touches it at the current configuration, then jumps to the surrogate's minimum. That minimum is the Guttman transform X = V^+ B(Y) Y. The sandwich inequality makes every step non-increasing in stress.

Leeuw & Mair (2009), equations 4 to 14, verified in docs/theory/algebra/deleeuw2009smacof-majorization.py.

Parameters:

Name	Type	Description	Default
D	ndarray	(n, n) dissimilarities, symmetric and hollow.	required
W	ndarray	(n, n) weights, symmetric, non-negative and hollow.	required
dim	int	Dimensions of the configuration.	2
seed	int	Seed for the starting configuration, used when classical scaling cannot supply one.	0
iters	int	Maximum majorization steps.	_SMACOF_ITERS
tol	float	Stop when a step lowers stress by less than this.	_SMACOF_TOL
history	bool	Also return the stress after each step.	False

Returns:

Type	Description
ndarray or tuple	(n, dim) configuration, or that with the stress history.

Source code in src/hifuku/charts.py

def smacof(D, W, dim: int = 2, *, seed: int = 0, iters: int = _SMACOF_ITERS,
           tol: float = _SMACOF_TOL, history: bool = False):
    """Weighted multidimensional scaling by majorization.

    Minimizes ``sum_{i<j} w_ij (delta_ij - d_ij(X))^2``.  A weight of zero drops
    a pair, which is the missing-value structure the paper describes and the way
    chart selection expresses a candidate set: a pair of anchors that no chosen
    triangle covers is simply not fitted.

    Weighted stress has no closed form, so each step builds a quadratic
    surrogate that sits above stress and touches it at the current
    configuration, then jumps to the surrogate's minimum.  That minimum is the
    Guttman transform ``X = V^+ B(Y) Y``.  The sandwich inequality makes every
    step non-increasing in stress.

    Leeuw & Mair (2009), equations 4 to 14, verified in
    ``docs/theory/algebra/deleeuw2009smacof-majorization.py``.

    Parameters
    ----------
    D : numpy.ndarray
        ``(n, n)`` dissimilarities, symmetric and hollow.
    W : numpy.ndarray
        ``(n, n)`` weights, symmetric, non-negative and hollow.
    dim : int
        Dimensions of the configuration.
    seed : int
        Seed for the starting configuration, used when classical scaling cannot
        supply one.
    iters : int
        Maximum majorization steps.
    tol : float
        Stop when a step lowers stress by less than this.
    history : bool
        Also return the stress after each step.

    Returns
    -------
    numpy.ndarray or tuple
        ``(n, dim)`` configuration, or that with the stress history.
    """
    from hifuku.ordination import mds

    D = np.asarray(D, dtype=np.float64)
    W = np.asarray(W, dtype=np.float64).copy()
    n = D.shape[0]
    np.fill_diagonal(W, 0.0)

    upper = np.triu_indices(n, k=1)

    def stress(X):
        d = _pairwise(X)
        return float((W[upper] * (D[upper] - d[upper]) ** 2).sum())

    # V from equation 6, summed over i < j.  Its pseudo-inverse solves the
    # Guttman step.  V is singular by construction, because shifting every point
    # leaves stress alone, so the pseudo-inverse is the right tool.
    V = np.diag(W.sum(axis=1)) - W
    V_pinv = np.linalg.pinv(V)

    X = mds(D, dim=dim)
    if not np.all(np.isfinite(X)) or np.allclose(X, 0.0):
        X = np.random.default_rng(seed).normal(size=(n, dim))

    trace = [stress(X)]
    for _ in range(iters):
        d = _pairwise(X)
        with np.errstate(divide="ignore", invalid="ignore"):
            ratio = np.where(d > 0.0, D / np.where(d > 0.0, d, 1.0), 0.0)
        # B(Y) from equation 8, with the same i < j convention as V.
        BW = W * ratio
        B = np.diag(BW.sum(axis=1)) - BW
        X = V_pinv @ (B @ X)
        trace.append(stress(X))
        if trace[-2] - trace[-1] < tol:
            break

    return (X, trace) if history else X