From 2af75e8177b223e219a997bdf669632c8e61c020 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 16 May 2026 20:40:52 +0200 Subject: [PATCH 1/3] Implement mutex watershed clustering --- include/bioimage_cpp/detail/mutex_storage.hxx | 54 +++++ .../bioimage_cpp/graph/mutex_watershed.hxx | 139 ++++++++++++ .../segmentation/mutex_watershed.hxx | 47 +--- src/bindings/graph.cxx | 40 ++++ src/bioimage_cpp/graph/__init__.py | 52 +++++ tests/graph/test_mutex_watershed.py | 205 ++++++++++++++++++ 6 files changed, 491 insertions(+), 46 deletions(-) create mode 100644 include/bioimage_cpp/detail/mutex_storage.hxx create mode 100644 include/bioimage_cpp/graph/mutex_watershed.hxx create mode 100644 tests/graph/test_mutex_watershed.py diff --git a/include/bioimage_cpp/detail/mutex_storage.hxx b/include/bioimage_cpp/detail/mutex_storage.hxx new file mode 100644 index 0000000..b22c9fd --- /dev/null +++ b/include/bioimage_cpp/detail/mutex_storage.hxx @@ -0,0 +1,54 @@ +#pragma once + +#include +#include +#include + +namespace bioimage_cpp { + +using MutexStorage = std::vector>; + +inline bool check_mutex( + const std::uint64_t first, + const std::uint64_t second, + const MutexStorage &mutexes +) { + const auto &first_mutexes = mutexes[first]; + const auto &second_mutexes = mutexes[second]; + if (first_mutexes.size() < second_mutexes.size()) { + return first_mutexes.find(second) != first_mutexes.end(); + } + return second_mutexes.find(first) != second_mutexes.end(); +} + +inline void insert_mutex( + const std::uint64_t first, + const std::uint64_t second, + MutexStorage &mutexes +) { + mutexes[first].insert(second); + mutexes[second].insert(first); +} + +inline void merge_mutexes( + const std::uint64_t root_from, + const std::uint64_t root_to, + MutexStorage &mutexes +) { + auto &mutexes_from = mutexes[root_from]; + auto &mutexes_to = mutexes[root_to]; + + for (const auto other_root : mutexes_from) { + auto &other_mutexes = mutexes[other_root]; + other_mutexes.erase(root_from); + if (other_root != root_to) { + other_mutexes.insert(root_to); + mutexes_to.insert(other_root); + } + } + mutexes_to.erase(root_from); + mutexes_to.erase(root_to); + mutexes_from.clear(); +} + +} // namespace bioimage_cpp diff --git a/include/bioimage_cpp/graph/mutex_watershed.hxx b/include/bioimage_cpp/graph/mutex_watershed.hxx new file mode 100644 index 0000000..c2f4847 --- /dev/null +++ b/include/bioimage_cpp/graph/mutex_watershed.hxx @@ -0,0 +1,139 @@ +#pragma once + +#include "bioimage_cpp/detail/mutex_storage.hxx" +#include "bioimage_cpp/detail/relabel.hxx" +#include "bioimage_cpp/detail/union_find.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" + +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph { + +// Mutex watershed clustering on an arbitrary undirected graph. +// +// `graph` defines the attractive edges with one cost per edge in `edge_costs`. +// `mutex_uvs` defines the long-range repulsive (mutex) edges as (u, v) pairs +// with one cost per edge in `mutex_costs`. Higher costs win — they are +// processed first when sorting jointly by descending weight. +// +// Returns dense node labels in [0, k) following first-occurrence order. +// +// This is a port of `compute_mws_clustering` from the affogato library, +// adapted to bioimage-cpp's UndirectedGraph and detail/ primitives. +inline std::vector mutex_watershed_clustering( + const UndirectedGraph &graph, + const std::vector &edge_costs, + const std::vector> &mutex_uvs, + const std::vector &mutex_costs +) { + const auto number_of_edges = static_cast(graph.number_of_edges()); + if (edge_costs.size() != number_of_edges) { + throw std::invalid_argument( + "edge_costs size must match graph.number_of_edges(), got edge_costs size=" + + std::to_string(edge_costs.size()) + + ", number_of_edges=" + std::to_string(number_of_edges) + ); + } + if (mutex_costs.size() != mutex_uvs.size()) { + throw std::invalid_argument( + "mutex_costs size must match mutex_uvs size, got mutex_costs size=" + + std::to_string(mutex_costs.size()) + + ", mutex_uvs size=" + std::to_string(mutex_uvs.size()) + ); + } + + const auto number_of_nodes = static_cast(graph.number_of_nodes()); + const auto number_of_mutex = mutex_uvs.size(); + + for (std::size_t index = 0; index < number_of_mutex; ++index) { + const auto u = mutex_uvs[index][0]; + const auto v = mutex_uvs[index][1]; + if (u >= number_of_nodes || v >= number_of_nodes) { + throw std::invalid_argument( + "mutex_uvs endpoints must be < number_of_nodes, got u=" + + std::to_string(u) + ", v=" + std::to_string(v) + + ", number_of_nodes=" + std::to_string(number_of_nodes) + ); + } + } + + struct WeightedEdge { + double weight; + std::uint64_t index; + bool is_mutex; + }; + + std::vector edge_order; + edge_order.reserve(number_of_edges + number_of_mutex); + for (std::size_t index = 0; index < number_of_edges; ++index) { + edge_order.push_back( + WeightedEdge{edge_costs[index], static_cast(index), false} + ); + } + for (std::size_t index = 0; index < number_of_mutex; ++index) { + edge_order.push_back( + WeightedEdge{mutex_costs[index], static_cast(index), true} + ); + } + + std::sort(edge_order.begin(), edge_order.end(), [](const auto &first, const auto &second) { + if (first.weight != second.weight) { + return first.weight > second.weight; + } + if (first.is_mutex != second.is_mutex) { + return !first.is_mutex; + } + return first.index < second.index; + }); + + bioimage_cpp::detail::UnionFind sets(number_of_nodes); + MutexStorage mutexes(number_of_nodes); + + for (const auto &edge : edge_order) { + std::uint64_t u; + std::uint64_t v; + if (edge.is_mutex) { + const auto &pair = mutex_uvs[static_cast(edge.index)]; + u = pair[0]; + v = pair[1]; + } else { + const auto uv = graph.uv(edge.index); + u = uv.first; + v = uv.second; + } + if (u == v) { + continue; + } + + const auto root_u = sets.find(u); + const auto root_v = sets.find(v); + if (root_u == root_v) { + continue; + } + + if (edge.is_mutex) { + insert_mutex(root_u, root_v, mutexes); + } else { + if (check_mutex(root_u, root_v, mutexes)) { + continue; + } + const auto new_root = sets.unite_roots(root_u, root_v); + const auto old_root = (new_root == root_u) ? root_v : root_u; + merge_mutexes(old_root, new_root, mutexes); + } + } + + std::vector roots(number_of_nodes); + for (std::size_t node = 0; node < number_of_nodes; ++node) { + roots[node] = sets.find(static_cast(node)); + } + return bioimage_cpp::detail::dense_relabel(roots); +} + +} // namespace bioimage_cpp::graph diff --git a/include/bioimage_cpp/segmentation/mutex_watershed.hxx b/include/bioimage_cpp/segmentation/mutex_watershed.hxx index 263779e..18f888a 100644 --- a/include/bioimage_cpp/segmentation/mutex_watershed.hxx +++ b/include/bioimage_cpp/segmentation/mutex_watershed.hxx @@ -2,6 +2,7 @@ #include "bioimage_cpp/array_view.hxx" #include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/detail/mutex_storage.hxx" #include "bioimage_cpp/detail/union_find.hxx" #include @@ -10,62 +11,16 @@ #include #include #include -#include #include namespace bioimage_cpp { -using MutexStorage = std::vector>; - template struct WeightedGridEdge { T weight; std::uint64_t id; }; -inline bool check_mutex( - const std::uint64_t first, - const std::uint64_t second, - const MutexStorage &mutexes -) { - const auto &first_mutexes = mutexes[first]; - const auto &second_mutexes = mutexes[second]; - if (first_mutexes.size() < second_mutexes.size()) { - return first_mutexes.find(second) != first_mutexes.end(); - } - return second_mutexes.find(first) != second_mutexes.end(); -} - -inline void insert_mutex( - const std::uint64_t first, - const std::uint64_t second, - MutexStorage &mutexes -) { - mutexes[first].insert(second); - mutexes[second].insert(first); -} - -inline void merge_mutexes( - const std::uint64_t root_from, - const std::uint64_t root_to, - MutexStorage &mutexes -) { - auto &mutexes_from = mutexes[root_from]; - auto &mutexes_to = mutexes[root_to]; - - for (const auto other_root : mutexes_from) { - auto &other_mutexes = mutexes[other_root]; - other_mutexes.erase(root_from); - if (other_root != root_to) { - other_mutexes.insert(root_to); - mutexes_to.insert(other_root); - } - } - mutexes_to.erase(root_from); - mutexes_to.erase(root_to); - mutexes_from.clear(); -} - template void mutex_watershed_grid( const ConstArrayView &affinities, diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx index 82d516b..77e74da 100644 --- a/src/bindings/graph.cxx +++ b/src/bindings/graph.cxx @@ -10,6 +10,7 @@ #include "bioimage_cpp/graph/lifted_from_affinities.hxx" #include "bioimage_cpp/graph/lifted_multicut.hxx" #include "bioimage_cpp/graph/multicut.hxx" +#include "bioimage_cpp/graph/mutex_watershed.hxx" #include "bioimage_cpp/graph/multicut/fusion_move.hxx" #include "bioimage_cpp/graph/multicut/greedy_additive.hxx" #include "bioimage_cpp/graph/multicut/greedy_fixation.hxx" @@ -790,6 +791,36 @@ UInt64Array lifted_multicut_kernighan_lin( return vector_to_uint64_array(label_vector); } +UInt64Array mutex_watershed_clustering( + const Graph &graph, + ConstDoubleArray edge_costs, + ConstUInt64Array mutex_uvs, + ConstDoubleArray mutex_costs +) { + const auto edge_cost_vector = + double_array_to_vector(edge_costs, "edge_costs", graph.number_of_edges()); + require_uv_array(mutex_uvs, "mutex_uvs"); + const auto n_mutex = mutex_uvs.shape(0); + const auto mutex_cost_vector = + double_array_to_vector(mutex_costs, "mutex_costs", static_cast(n_mutex)); + + std::vector> mutex_uv_vector(n_mutex); + const auto *uv_data = mutex_uvs.data(); + for (std::size_t index = 0; index < n_mutex; ++index) { + mutex_uv_vector[index][0] = uv_data[2 * index]; + mutex_uv_vector[index][1] = uv_data[2 * index + 1]; + } + + std::vector labels; + { + nb::gil_scoped_release release; + labels = graph::mutex_watershed_clustering( + graph, edge_cost_vector, mutex_uv_vector, mutex_cost_vector + ); + } + return vector_to_uint64_array(labels); +} + UInt64Array multicut_fusion_move( const Graph &graph, ConstDoubleArray costs, @@ -1429,6 +1460,15 @@ void bind_graph(nb::module_ &m) { nb::arg("epsilon") ); + m.def( + "_mutex_watershed_clustering", + &mutex_watershed_clustering, + nb::arg("graph"), + nb::arg("edge_costs"), + nb::arg("mutex_uvs"), + nb::arg("mutex_costs") + ); + // Lifted multicut sub-solver hierarchy. Same shape as the multicut sub- // solver bindings — opaque to Python, used by future fusion-move drivers. nb::class_(m, "_LiftedMulticutSolverBase"); diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index bdc310a..3558148 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -599,6 +599,57 @@ def edge_weighted_watershed( return run(graph, weight_array, seed_array) +def mutex_watershed_clustering( + graph: UndirectedGraph | RegionAdjacencyGraph, + edge_costs, + mutex_uvs, + mutex_costs, +) -> np.ndarray: + """Mutex watershed clustering on an undirected graph. + + Attractive edges come from ``graph`` (one cost per edge in + ``edge_costs``); repulsive long-range edges are supplied separately as + ``mutex_uvs`` with weights ``mutex_costs``. All edges are jointly sorted + by descending weight and processed in a single pass: an attractive edge + merges its two components unless a mutex constraint already separates + them; a mutex edge installs a constraint between the two current + components. + + The input format matches :class:`LiftedMulticutObjective` — the same + ``(graph, edge_costs, lifted_uvs, lifted_costs)`` arrays can be passed + here as ``(graph, edge_costs, mutex_uvs, mutex_costs)`` — though the + algorithms differ (mutex constraints are hard; lifted costs are soft). + + Parameters + ---------- + graph: + :class:`UndirectedGraph` or :class:`RegionAdjacencyGraph` defining + the attractive edges. + edge_costs: + 1D float64 array of length ``graph.number_of_edges``. Higher + values are more attractive. + mutex_uvs: + ``(n_mutex, 2)`` uint64 array of (u, v) pairs for the mutex edges. + mutex_costs: + 1D float64 array of length ``n_mutex``. Higher values are stronger + repulsions. + + Returns + ------- + np.ndarray + ``(graph.number_of_nodes,)`` uint64 array. Dense node labels in + ``[0, k)`` in first-occurrence order. + """ + edge_cost_array = _as_edge_costs(edge_costs, graph) + mutex_uv_array = _as_uv_array(mutex_uvs, "mutex_uvs") + mutex_cost_array = _as_1d_array( + mutex_costs, np.float64, "mutex_costs", int(mutex_uv_array.shape[0]) + ) + return _core._mutex_watershed_clustering( + graph, edge_cost_array, mutex_uv_array, mutex_cost_array + ) + + class MulticutObjective: """Multicut objective for an undirected graph and edge costs.""" @@ -1818,6 +1869,7 @@ def _normalize_number_of_threads(number_of_threads: int) -> int: "load_multicut_problem", "load_multicut_problem_data", "multicut_problem_path", + "mutex_watershed_clustering", "project_node_labels_to_pixels", "region_adjacency_graph", "undirected_graph", diff --git a/tests/graph/test_mutex_watershed.py b/tests/graph/test_mutex_watershed.py new file mode 100644 index 0000000..aee3d72 --- /dev/null +++ b/tests/graph/test_mutex_watershed.py @@ -0,0 +1,205 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def _graph_from_edges(n: int, uvs): + uvs = np.asarray(uvs, dtype=np.uint64) + return bic.graph.UndirectedGraph.from_edges(n, uvs) + + +def _canonicalize(labels): + # Map labels to first-occurrence dense ids so two segmentations agree iff + # they produce the same partition (regardless of integer values). + array = np.asarray(labels, dtype=np.uint64) + _, inverse = np.unique(array, return_inverse=True) + # `np.unique`'s `return_inverse` is sorted by value, not by first + # occurrence — remap one more time to get first-occurrence order. + seen = {} + out = np.empty_like(inverse) + for index, value in enumerate(array): + key = int(value) + if key not in seen: + seen[key] = len(seen) + out[index] = seen[key] + return out + + +def test_all_attractive_merges_to_one_component(): + graph = _graph_from_edges(5, [[0, 1], [1, 2], [2, 3], [3, 4]]) + edge_costs = np.array([1.0, 1.0, 1.0, 1.0], dtype=np.float64) + mutex_uvs = np.zeros((0, 2), dtype=np.uint64) + mutex_costs = np.zeros(0, dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + assert labels.dtype == np.uint64 + assert labels.shape == (5,) + np.testing.assert_array_equal(labels, np.zeros(5, dtype=np.uint64)) + + +def test_mutex_only_keeps_singletons(): + graph = bic.graph.UndirectedGraph(4) + edge_costs = np.zeros(0, dtype=np.float64) + mutex_uvs = np.array( + [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]], dtype=np.uint64 + ) + mutex_costs = np.ones(6, dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + np.testing.assert_array_equal(_canonicalize(labels), np.array([0, 1, 2, 3])) + + +def test_mutex_beats_attractive_when_higher_weight(): + # Triangle a-b, b-c attractive (weights 1.0, 1.0) + mutex a-c (weight 2.0). + # Mutex arrives first (highest weight), then a-b merges, then b-c is + # blocked because root(a)=root(b) carries a mutex with c. + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + edge_costs = np.array([1.0, 1.0], dtype=np.float64) + mutex_uvs = np.array([[0, 2]], dtype=np.uint64) + mutex_costs = np.array([2.0], dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + # {0,1} together, {2} alone. + np.testing.assert_array_equal(_canonicalize(labels), np.array([0, 0, 1])) + + +def test_mutex_constraint_propagates_through_merge(): + # Nodes 0,1,2 chained attractively (weights 10, 9). Mutex 0-3 with + # weight 5. Attractive 2-3 with weight 4. After 0,1,2 merge into one + # root, the mutex propagated via merge_mutexes must block 2-3. + graph = _graph_from_edges(4, [[0, 1], [1, 2], [2, 3]]) + edge_costs = np.array([10.0, 9.0, 4.0], dtype=np.float64) + mutex_uvs = np.array([[0, 3]], dtype=np.uint64) + mutex_costs = np.array([5.0], dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + np.testing.assert_array_equal(_canonicalize(labels), np.array([0, 0, 0, 1])) + + +def test_attractive_higher_than_mutex_still_merges(): + # Mirror of test_mutex_beats_attractive but with the order swapped: + # attractive triangle wins so all three nodes end up together. + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + edge_costs = np.array([5.0, 5.0], dtype=np.float64) + mutex_uvs = np.array([[0, 2]], dtype=np.uint64) + mutex_costs = np.array([0.1], dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + np.testing.assert_array_equal(_canonicalize(labels), np.array([0, 0, 0])) + + +def test_lifted_multicut_inputs_accepted_unchanged(): + # The (graph, edge_costs, lifted_uvs, lifted_costs) shape used to build + # a LiftedMulticutObjective must also be a valid call to + # mutex_watershed_clustering — same arrays, no reshape. We do not assert + # that the algorithms agree on the labels (they don't, in general), + # only that the input format is compatible. + graph = _graph_from_edges(4, [[0, 1], [1, 2], [2, 3]]) + edge_costs = np.array([1.0, -0.5, 1.0], dtype=np.float64) + lifted_uvs = np.array([[0, 2], [0, 3], [1, 3]], dtype=np.uint64) + lifted_costs = np.array([-1.0, -1.0, 0.5], dtype=np.float64) + + objective = bic.graph.LiftedMulticutObjective( + graph, + edge_costs, + lifted_uvs=lifted_uvs, + lifted_costs=lifted_costs, + ) + assert objective is not None + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, lifted_uvs, lifted_costs + ) + assert labels.shape == (4,) + assert labels.dtype == np.uint64 + + +def test_deterministic_across_runs(): + rng = np.random.default_rng(7) + n = 30 + uvs = [] + for u in range(n): + for v in range(u + 1, min(u + 4, n)): + uvs.append([u, v]) + uvs = np.array(uvs, dtype=np.uint64) + graph = _graph_from_edges(n, uvs) + edge_costs = rng.uniform(-1.0, 1.0, size=int(graph.number_of_edges)).astype(np.float64) + + mutex_uvs = rng.integers(0, n, size=(40, 2), dtype=np.uint64) + mutex_uvs = mutex_uvs[mutex_uvs[:, 0] != mutex_uvs[:, 1]] + mutex_costs = rng.uniform(0.0, 1.0, size=mutex_uvs.shape[0]).astype(np.float64) + + first = bic.graph.mutex_watershed_clustering(graph, edge_costs, mutex_uvs, mutex_costs) + second = bic.graph.mutex_watershed_clustering(graph, edge_costs, mutex_uvs, mutex_costs) + np.testing.assert_array_equal(first, second) + + +def test_dense_label_range(): + graph = _graph_from_edges(6, [[0, 1], [2, 3], [4, 5]]) + edge_costs = np.array([1.0, 1.0, 1.0], dtype=np.float64) + mutex_uvs = np.zeros((0, 2), dtype=np.uint64) + mutex_costs = np.zeros(0, dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + # Three components → labels must be {0, 1, 2}, dense. + assert set(int(value) for value in labels) == {0, 1, 2} + assert labels.max() == 2 + + +def test_invalid_edge_costs_length_raises(): + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + bad_costs = np.array([1.0], dtype=np.float64) # expected size 2 + mutex_uvs = np.zeros((0, 2), dtype=np.uint64) + mutex_costs = np.zeros(0, dtype=np.float64) + + with pytest.raises(ValueError): + bic.graph.mutex_watershed_clustering(graph, bad_costs, mutex_uvs, mutex_costs) + + +def test_invalid_mutex_uvs_shape_raises(): + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + edge_costs = np.array([1.0, 1.0], dtype=np.float64) + bad_mutex_uvs = np.array([0, 2], dtype=np.uint64) # 1D, not (n, 2) + mutex_costs = np.array([1.0], dtype=np.float64) + + with pytest.raises(ValueError): + bic.graph.mutex_watershed_clustering(graph, edge_costs, bad_mutex_uvs, mutex_costs) + + +def test_mismatched_mutex_costs_raises(): + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + edge_costs = np.array([1.0, 1.0], dtype=np.float64) + mutex_uvs = np.array([[0, 2]], dtype=np.uint64) + bad_mutex_costs = np.array([1.0, 2.0], dtype=np.float64) # expected size 1 + + with pytest.raises(ValueError): + bic.graph.mutex_watershed_clustering(graph, edge_costs, mutex_uvs, bad_mutex_costs) + + +def test_out_of_range_mutex_endpoint_raises(): + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + edge_costs = np.array([1.0, 1.0], dtype=np.float64) + bad_mutex_uvs = np.array([[0, 99]], dtype=np.uint64) + mutex_costs = np.array([1.0], dtype=np.float64) + + with pytest.raises((ValueError, RuntimeError)): + bic.graph.mutex_watershed_clustering(graph, edge_costs, bad_mutex_uvs, mutex_costs) From b4055f2e227983e4330098928ed387c3da4b7d8d Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 16 May 2026 21:00:20 +0200 Subject: [PATCH 2/3] Update mutex impl --- development/graph/check_mutex_clustering.py | 341 ++++++++++++++++++ .../bioimage_cpp/graph/mutex_watershed.hxx | 13 +- src/bindings/graph.cxx | 33 +- src/bioimage_cpp/graph/__init__.py | 60 ++- ...rshed.py => test_mutex_watershed_graph.py} | 60 +++ 5 files changed, 480 insertions(+), 27 deletions(-) create mode 100644 development/graph/check_mutex_clustering.py rename tests/graph/{test_mutex_watershed.py => test_mutex_watershed_graph.py} (77%) diff --git a/development/graph/check_mutex_clustering.py b/development/graph/check_mutex_clustering.py new file mode 100644 index 0000000..48db1ad --- /dev/null +++ b/development/graph/check_mutex_clustering.py @@ -0,0 +1,341 @@ +"""Compare bioimage-cpp `mutex_watershed_clustering` against affogato's +`compute_mws_clustering` reference, using the three registered lifted multicut +problems as inputs (`local_uvs`/`local_costs` as attractive edges, +`lifted_uvs`/`lifted_costs` as mutex edges). + +Two modes: + +* ``check`` (default): run a single problem (``--size 2d|3d|grid``), report + partition equivalence + runtimes. +* ``evaluate``: run all registered sizes and print a markdown comparison + table — mirrors the layout of + ``development/graph/multicut/evaluate_solvers.py`` and + ``development/graph/lifted_multicut/evaluate_solvers.py``. + +Not part of the pytest suite (per AGENTS.md). Run manually with affogato +installed. +""" + +from __future__ import annotations + +import argparse +from statistics import median +from time import perf_counter +from typing import Callable + +import numpy as np + + +PROBLEMS = ("2d", "3d", "grid") +DTYPES = ("float32", "float64") + + +def load_problem(size: str, *, timeout: float): + import bioimage_cpp as bic + + problem = bic.graph.load_lifted_multicut_problem(size, timeout=timeout) + bic_graph = bic.graph.UndirectedGraph.from_edges(problem.n_nodes, problem.local_uvs) + return bic_graph, problem + + +def run_bioimage_cpp(bic_graph, problem, *, dtype: np.dtype) -> np.ndarray: + import bioimage_cpp as bic + + return bic.graph.mutex_watershed_clustering( + bic_graph, + problem.local_costs.astype(dtype, copy=False), + problem.lifted_uvs, + problem.lifted_costs.astype(dtype, copy=False), + ) + + +def run_affogato_reference(problem) -> np.ndarray: + # affogato's `compute_mws_clustering` uses float32 weights. + from affogato.segmentation import compute_mws_clustering + + return compute_mws_clustering( + int(problem.n_nodes), + problem.local_uvs.astype(np.uint64, copy=False), + problem.lifted_uvs.astype(np.uint64, copy=False), + problem.local_costs.astype(np.float32, copy=False), + problem.lifted_costs.astype(np.float32, copy=False), + ) + + +def _load_validation_metrics(): + try: + from elf.validation import rand_index, variation_of_information + + return "elf.validation", rand_index, variation_of_information + except ImportError: + from elf.evaluation import rand_index, variation_of_information + + return "elf.evaluation", rand_index, variation_of_information + + +def _canonical_labels(labels: np.ndarray) -> np.ndarray: + # Map to dense ids in first-occurrence order, so two partitions compare + # equal iff they induce the same node grouping (independent of which + # integer happened to be assigned to which cluster). + array = np.asarray(labels) + _, first_index, inverse = np.unique( + array, return_index=True, return_inverse=True + ) + order = np.argsort(first_index) + remap = np.empty_like(order) + remap[order] = np.arange(order.size) + return remap[inverse].astype(np.uint64, copy=False) + + +def compare_partitions( + candidate: np.ndarray, + reference: np.ndarray, +) -> dict: + source, rand_index, variation_of_information = _load_validation_metrics() + vi_split, vi_merge = variation_of_information(candidate, reference) + adapted_rand_error, ri = rand_index(candidate, reference) + partition_equal = bool( + np.array_equal(_canonical_labels(candidate), _canonical_labels(reference)) + ) + return { + "validation_source": source, + "vi_split": float(vi_split), + "vi_merge": float(vi_merge), + "adapted_rand_error": float(adapted_rand_error), + "rand_index": float(ri), + "partition_equal": partition_equal, + "n_clusters_bic": int(np.unique(candidate).size), + "n_clusters_reference": int(np.unique(reference).size), + } + + +def time_function_interleaved( + bic_run: Callable[[], np.ndarray], + reference_run: Callable[[], np.ndarray], + repeats: int, +) -> tuple[list[float], np.ndarray, list[float], np.ndarray]: + # Warm up both implementations so JIT / first-call allocation costs do + # not contaminate the timed runs. + bic_result = bic_run() + ref_result = reference_run() + + bic_timings: list[float] = [] + ref_timings: list[float] = [] + for repeat in range(repeats): + if repeat % 2 == 0: + start = perf_counter() + bic_result = bic_run() + bic_timings.append(perf_counter() - start) + start = perf_counter() + ref_result = reference_run() + ref_timings.append(perf_counter() - start) + else: + start = perf_counter() + ref_result = reference_run() + ref_timings.append(perf_counter() - start) + start = perf_counter() + bic_result = bic_run() + bic_timings.append(perf_counter() - start) + return bic_timings, bic_result, ref_timings, ref_result + + +def run_size( + size: str, *, repeats: int, timeout: float, dtype: np.dtype +) -> dict: + bic_graph, problem = load_problem(size, timeout=timeout) + + bic_timings, bic_labels, ref_timings, ref_labels = time_function_interleaved( + lambda: run_bioimage_cpp(bic_graph, problem, dtype=dtype), + lambda: run_affogato_reference(problem), + repeats, + ) + + metrics = compare_partitions(bic_labels, ref_labels) + bic_median = median(bic_timings) + ref_median = median(ref_timings) + return { + "problem": size, + "dtype": np.dtype(dtype).name, + "nodes": int(problem.n_nodes), + "local_edges": int(problem.local_uvs.shape[0]), + "lifted_edges": int(problem.lifted_uvs.shape[0]), + "bic_runtime_s": bic_median, + "affogato_runtime_s": ref_median, + "runtime_ratio": ref_median / bic_median if bic_median > 0 else float("inf"), + **metrics, + } + + +def print_check_report(result: dict) -> None: + print(f"problem: size={result['problem']}, dtype={result['dtype']}, " + f"nodes={result['nodes']}, local edges={result['local_edges']}, " + f"lifted edges={result['lifted_edges']}") + print(f"validation metrics: {result['validation_source']}") + print( + "VI split/merge: " + f"{result['vi_split']:.6g} / {result['vi_merge']:.6g}" + ) + print( + "adapted rand error / rand index: " + f"{result['adapted_rand_error']:.6g} / {result['rand_index']:.12g}" + ) + print(f"partition equality (after canonical relabel): {result['partition_equal']}") + print(f"clusters (bic / affogato): " + f"{result['n_clusters_bic']} / {result['n_clusters_reference']}") + print(f"bioimage-cpp median runtime [s]: {result['bic_runtime_s']:.6f}") + print(f"affogato median runtime [s]: {result['affogato_runtime_s']:.6f}") + print(f"affogato / bioimage-cpp runtime ratio: {result['runtime_ratio']:.3f}x") + + +def format_float(value: float) -> str: + return f"{value:.6g}" + + +def print_markdown_table(rows: list[dict]) -> None: + headers = [ + "problem", + "dtype", + "nodes", + "local_edges", + "lifted_edges", + "n_clusters_bic", + "n_clusters_affogato", + "vi_split", + "vi_merge", + "adapted_rand_error", + "rand_index", + "partition_equal", + "bic_runtime_s", + "affogato_runtime_s", + "runtime_ratio_affogato_over_bic", + ] + print("| " + " | ".join(headers) + " |") + print("| " + " | ".join(["---"] * len(headers)) + " |") + for row in rows: + values = [ + row["problem"], + row["dtype"], + str(row["nodes"]), + str(row["local_edges"]), + str(row["lifted_edges"]), + str(row["n_clusters_bic"]), + str(row["n_clusters_reference"]), + format_float(row["vi_split"]), + format_float(row["vi_merge"]), + format_float(row["adapted_rand_error"]), + format_float(row["rand_index"]), + str(row["partition_equal"]), + format_float(row["bic_runtime_s"]), + format_float(row["affogato_runtime_s"]), + format_float(row["runtime_ratio"]), + ] + print("| " + " | ".join(values) + " |") + + +def build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + description=( + "Compare bioimage-cpp `mutex_watershed_clustering` against " + "affogato's `compute_mws_clustering` on the registered lifted " + "multicut problems." + ) + ) + subparsers = parser.add_subparsers(dest="mode", required=False) + + check_parser = subparsers.add_parser( + "check", + help=( + "Run a single problem and print a partition-equivalence + " + "runtime report (default mode)." + ), + ) + check_parser.add_argument( + "--size", + choices=PROBLEMS, + default="3d", + help="Lifted multicut problem instance to load (default: 3d).", + ) + check_parser.add_argument("--repeats", type=int, default=3) + check_parser.add_argument("--timeout", type=float, default=60.0) + check_parser.add_argument( + "--dtype", + choices=DTYPES, + default="float32", + help=( + "Weight dtype for the bioimage-cpp call (default: float32, " + "matching the precision affogato's reference uses internally)." + ), + ) + + evaluate_parser = subparsers.add_parser( + "evaluate", + help=( + "Run all registered problem sizes and print a markdown " + "comparison table." + ), + ) + evaluate_parser.add_argument( + "--problems", + nargs="+", + choices=PROBLEMS, + default=PROBLEMS, + help="Problems to evaluate. Defaults to all.", + ) + evaluate_parser.add_argument("--n-repeats", type=int, default=1) + evaluate_parser.add_argument("--timeout", type=float, default=60.0) + evaluate_parser.add_argument( + "--dtypes", + nargs="+", + choices=DTYPES, + default=("float32",), + help=( + "Weight dtype(s) for the bioimage-cpp call. Each (problem, " + "dtype) pair becomes a row. Defaults to float32 (matches " + "affogato's internal precision)." + ), + ) + + return parser + + +def main() -> None: + parser = build_parser() + args = parser.parse_args() + + # Default to `check` if no subcommand was passed, matching the existing + # `check_*` scripts in this directory. + if args.mode is None or args.mode == "check": + if args.mode is None: + args.size = "3d" + args.repeats = 3 + args.timeout = 60.0 + args.dtype = "float32" + if args.repeats < 1: + raise ValueError("--repeats must be at least 1") + result = run_size( + args.size, + repeats=args.repeats, + timeout=args.timeout, + dtype=np.dtype(args.dtype), + ) + print_check_report(result) + return + + if args.n_repeats < 1: + raise ValueError("--n-repeats must be at least 1") + rows = [] + for problem in args.problems: + for dtype in args.dtypes: + rows.append( + run_size( + problem, + repeats=args.n_repeats, + timeout=args.timeout, + dtype=np.dtype(dtype), + ) + ) + print_markdown_table(rows) + + +if __name__ == "__main__": + main() diff --git a/include/bioimage_cpp/graph/mutex_watershed.hxx b/include/bioimage_cpp/graph/mutex_watershed.hxx index c2f4847..bf929be 100644 --- a/include/bioimage_cpp/graph/mutex_watershed.hxx +++ b/include/bioimage_cpp/graph/mutex_watershed.hxx @@ -26,11 +26,16 @@ namespace bioimage_cpp::graph { // // This is a port of `compute_mws_clustering` from the affogato library, // adapted to bioimage-cpp's UndirectedGraph and detail/ primitives. -inline std::vector mutex_watershed_clustering( +// +// Templated on the weight type. Concrete instantiations for `float` and +// `double` are provided by the binding layer; other floating types are +// supported but must be instantiated explicitly. +template +std::vector mutex_watershed_clustering( const UndirectedGraph &graph, - const std::vector &edge_costs, + const std::vector &edge_costs, const std::vector> &mutex_uvs, - const std::vector &mutex_costs + const std::vector &mutex_costs ) { const auto number_of_edges = static_cast(graph.number_of_edges()); if (edge_costs.size() != number_of_edges) { @@ -64,7 +69,7 @@ inline std::vector mutex_watershed_clustering( } struct WeightedEdge { - double weight; + WeightT weight; std::uint64_t index; bool is_mutex; }; diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx index 77e74da..675c60b 100644 --- a/src/bindings/graph.cxx +++ b/src/bindings/graph.cxx @@ -791,18 +791,19 @@ UInt64Array lifted_multicut_kernighan_lin( return vector_to_uint64_array(label_vector); } -UInt64Array mutex_watershed_clustering( +template +UInt64Array mutex_watershed_clustering_t( const Graph &graph, - ConstDoubleArray edge_costs, + ConstArray1D edge_costs, ConstUInt64Array mutex_uvs, - ConstDoubleArray mutex_costs + ConstArray1D mutex_costs ) { const auto edge_cost_vector = - double_array_to_vector(edge_costs, "edge_costs", graph.number_of_edges()); + array_1d_to_vector(edge_costs, "edge_costs", graph.number_of_edges()); require_uv_array(mutex_uvs, "mutex_uvs"); const auto n_mutex = mutex_uvs.shape(0); const auto mutex_cost_vector = - double_array_to_vector(mutex_costs, "mutex_costs", static_cast(n_mutex)); + array_1d_to_vector(mutex_costs, "mutex_costs", static_cast(n_mutex)); std::vector> mutex_uv_vector(n_mutex); const auto *uv_data = mutex_uvs.data(); @@ -814,7 +815,7 @@ UInt64Array mutex_watershed_clustering( std::vector labels; { nb::gil_scoped_release release; - labels = graph::mutex_watershed_clustering( + labels = graph::mutex_watershed_clustering( graph, edge_cost_vector, mutex_uv_vector, mutex_cost_vector ); } @@ -1460,14 +1461,18 @@ void bind_graph(nb::module_ &m) { nb::arg("epsilon") ); - m.def( - "_mutex_watershed_clustering", - &mutex_watershed_clustering, - nb::arg("graph"), - nb::arg("edge_costs"), - nb::arg("mutex_uvs"), - nb::arg("mutex_costs") - ); + const auto register_mutex_watershed_clustering = [&m](const char *name) { + m.def( + name, + &mutex_watershed_clustering_t, + nb::arg("graph"), + nb::arg("edge_costs"), + nb::arg("mutex_uvs"), + nb::arg("mutex_costs") + ); + }; + register_mutex_watershed_clustering.operator()("_mutex_watershed_clustering_float32"); + register_mutex_watershed_clustering.operator()("_mutex_watershed_clustering_float64"); // Lifted multicut sub-solver hierarchy. Same shape as the multicut sub- // solver bindings — opaque to Python, used by future fusion-move drivers. diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index 3558148..be0c11a 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -599,6 +599,30 @@ def edge_weighted_watershed( return run(graph, weight_array, seed_array) +_MUTEX_WATERSHED_CLUSTERING_BY_DTYPE = { + np.dtype("float32"): _core._mutex_watershed_clustering_float32, + np.dtype("float64"): _core._mutex_watershed_clustering_float64, +} + + +def _resolve_weight_dtype(array, name: str) -> np.ndarray: + """Coerce a weights input to a supported floating dtype. + + ``float32`` and ``float64`` pass through unchanged. Other floating + dtypes are cast to ``float32`` (matches the convention used by + :func:`edge_weighted_watershed`). Non-floating dtypes raise. + """ + array = np.asarray(array) + if array.dtype in (np.dtype("float32"), np.dtype("float64")): + return array + if np.issubdtype(array.dtype, np.floating): + return array.astype(np.float32, copy=False) + raise TypeError( + f"{name} must have a floating dtype (float32 or float64), got " + f"dtype={array.dtype}" + ) + + def mutex_watershed_clustering( graph: UndirectedGraph | RegionAdjacencyGraph, edge_costs, @@ -626,13 +650,15 @@ def mutex_watershed_clustering( :class:`UndirectedGraph` or :class:`RegionAdjacencyGraph` defining the attractive edges. edge_costs: - 1D float64 array of length ``graph.number_of_edges``. Higher - values are more attractive. + 1D array of length ``graph.number_of_edges``. Supported dtypes are + ``float32`` and ``float64``; other floating dtypes are cast to + ``float32``. Higher values are more attractive. mutex_uvs: ``(n_mutex, 2)`` uint64 array of (u, v) pairs for the mutex edges. mutex_costs: - 1D float64 array of length ``n_mutex``. Higher values are stronger - repulsions. + 1D array of length ``n_mutex``. Same dtype rules as ``edge_costs``; + if the two dtypes differ both are promoted to ``float64``. Higher + values are stronger repulsions. Returns ------- @@ -640,14 +666,30 @@ def mutex_watershed_clustering( ``(graph.number_of_nodes,)`` uint64 array. Dense node labels in ``[0, k)`` in first-occurrence order. """ - edge_cost_array = _as_edge_costs(edge_costs, graph) + edge_cost_array = _resolve_weight_dtype(edge_costs, "edge_costs") + mutex_cost_array = _resolve_weight_dtype(mutex_costs, "mutex_costs") + # Use a single instantiation for both arrays. If the user supplied + # mismatched dtypes (one float32, one float64) we promote both to + # float64 — the wider type — rather than silently downcasting. + if edge_cost_array.dtype != mutex_cost_array.dtype: + edge_cost_array = edge_cost_array.astype(np.float64, copy=False) + mutex_cost_array = mutex_cost_array.astype(np.float64, copy=False) + + edge_cost_array = _as_1d_array( + edge_cost_array, + edge_cost_array.dtype, + "edge_costs", + int(graph.number_of_edges), + ) mutex_uv_array = _as_uv_array(mutex_uvs, "mutex_uvs") mutex_cost_array = _as_1d_array( - mutex_costs, np.float64, "mutex_costs", int(mutex_uv_array.shape[0]) - ) - return _core._mutex_watershed_clustering( - graph, edge_cost_array, mutex_uv_array, mutex_cost_array + mutex_cost_array, + mutex_cost_array.dtype, + "mutex_costs", + int(mutex_uv_array.shape[0]), ) + run = _MUTEX_WATERSHED_CLUSTERING_BY_DTYPE[edge_cost_array.dtype] + return run(graph, edge_cost_array, mutex_uv_array, mutex_cost_array) class MulticutObjective: diff --git a/tests/graph/test_mutex_watershed.py b/tests/graph/test_mutex_watershed_graph.py similarity index 77% rename from tests/graph/test_mutex_watershed.py rename to tests/graph/test_mutex_watershed_graph.py index aee3d72..44765af 100644 --- a/tests/graph/test_mutex_watershed.py +++ b/tests/graph/test_mutex_watershed_graph.py @@ -195,6 +195,66 @@ def test_mismatched_mutex_costs_raises(): bic.graph.mutex_watershed_clustering(graph, edge_costs, mutex_uvs, bad_mutex_costs) +def test_float32_inputs_supported(): + graph = _graph_from_edges(4, [[0, 1], [1, 2], [2, 3]]) + edge_costs = np.array([10.0, 9.0, 4.0], dtype=np.float32) + mutex_uvs = np.array([[0, 3]], dtype=np.uint64) + mutex_costs = np.array([5.0], dtype=np.float32) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + + assert labels.dtype == np.uint64 + np.testing.assert_array_equal(_canonicalize(labels), np.array([0, 0, 0, 1])) + + +def test_float32_and_float64_agree_on_simple_problem(): + # On inputs where no weights tie, float32 and float64 must produce the + # exact same partition. + graph = _graph_from_edges(5, [[0, 1], [1, 2], [2, 3], [3, 4]]) + edge_costs = np.array([1.0, -0.25, 0.75, -0.5], dtype=np.float64) + mutex_uvs = np.array([[0, 4], [1, 3]], dtype=np.uint64) + mutex_costs = np.array([2.0, 0.1], dtype=np.float64) + + labels_64 = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + labels_32 = bic.graph.mutex_watershed_clustering( + graph, + edge_costs.astype(np.float32), + mutex_uvs, + mutex_costs.astype(np.float32), + ) + np.testing.assert_array_equal(_canonicalize(labels_32), _canonicalize(labels_64)) + + +def test_mismatched_dtypes_are_promoted(): + # Passing float32 edge_costs alongside float64 mutex_costs is accepted — + # the wrapper promotes both to float64 before dispatching. + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + edge_costs = np.array([1.0, 1.0], dtype=np.float32) + mutex_uvs = np.array([[0, 2]], dtype=np.uint64) + mutex_costs = np.array([2.0], dtype=np.float64) + + labels = bic.graph.mutex_watershed_clustering( + graph, edge_costs, mutex_uvs, mutex_costs + ) + np.testing.assert_array_equal(_canonicalize(labels), np.array([0, 0, 1])) + + +def test_integer_edge_costs_rejected(): + graph = _graph_from_edges(3, [[0, 1], [1, 2]]) + bad_edge_costs = np.array([1, 1], dtype=np.int32) + mutex_uvs = np.zeros((0, 2), dtype=np.uint64) + mutex_costs = np.zeros(0, dtype=np.float64) + + with pytest.raises(TypeError): + bic.graph.mutex_watershed_clustering( + graph, bad_edge_costs, mutex_uvs, mutex_costs + ) + + def test_out_of_range_mutex_endpoint_raises(): graph = _graph_from_edges(3, [[0, 1], [1, 2]]) edge_costs = np.array([1.0, 1.0], dtype=np.float64) From 6249766bb062716d3911df67333d28235d2115c4 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Sat, 16 May 2026 21:02:48 +0200 Subject: [PATCH 3/3] Update migration guide --- MIGRATION_GUIDE.md | 74 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 5dbcd45..6ac442b 100644 --- a/MIGRATION_GUIDE.md +++ b/MIGRATION_GUIDE.md @@ -811,6 +811,13 @@ Notes: ## Mutex Watershed +`bioimage-cpp` ships two mutex-watershed entry points, mirroring the two +affogato APIs: one that consumes a dense affinity grid, and one that +consumes an arbitrary graph with a separate list of mutex (long-range +repulsive) edges. + +### Grid-based mutex watershed (affinity volumes) + Affogato: ```python @@ -853,6 +860,73 @@ Important migration notes: masked pixels are set to label `0`. - Output labels are `uint64`, consecutive, and 1-based for foreground pixels. +### Mutex watershed on a generic graph + +For mutex watershed on an arbitrary undirected graph (region adjacency graph +or otherwise) with a separate list of long-range repulsive edges, +`bioimage-cpp` provides `bic.graph.mutex_watershed_clustering`. This is a +port of affogato's `compute_mws_clustering` using the same input format as +`LiftedMulticutObjective`: a base graph carries the attractive edges, and +long-range (called *mutex* here) edges are supplied alongside as a `(M, 2)` +node-pair array. The same `(graph, edge_costs, lifted_uvs, lifted_costs)` +tuple used to build a lifted multicut problem can be passed to the mutex +watershed clustering without any reshaping. + +Affogato: + +```python +from affogato.segmentation import compute_mws_clustering + +labels = compute_mws_clustering( + number_of_nodes, + uvs.astype(np.uint64), + mutex_uvs.astype(np.uint64), + weights.astype(np.float32), + mutex_weights.astype(np.float32), +) +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +graph = bic.graph.UndirectedGraph.from_edges(number_of_nodes, uvs) +labels = bic.graph.mutex_watershed_clustering( + graph, + weights, + mutex_uvs, + mutex_weights, +) +``` + +Notes: + +- The attractive edges are the edges of the base graph; the count and the + ordering of `weights` must match `graph.number_of_edges`. Mutex edges + are supplied separately as `(M, 2)` `uint64` pairs with matching + `mutex_weights`. +- Both `weights` and `mutex_weights` accept `float32` and `float64`. The + wrapper dispatches to a templated C++ instantiation per dtype; other + floating dtypes are cast to `float32`. If the two arrays' dtypes do not + match, both are promoted to `float64` rather than silently downcast. +- Higher weights are processed first (in descending order) — the same + convention affogato uses. +- The implementation reuses the union-find and per-root mutex-set helpers + shared with the grid-based mutex watershed (`detail/mutex_storage.hxx`), + so behavior is consistent between the two entry points. +- Output labels are dense `uint64` ids in `0 .. number_of_clusters - 1`, + assigned in first-occurrence order (matches the convention of the graph + multicut solvers, *not* the 1-based foreground labels produced by the + grid-based variant). +- The function accepts both `UndirectedGraph` and `RegionAdjacencyGraph`. +- Tie-breaking is deterministic: when weights are equal, attractive edges + are processed before mutex edges, then by index. Affogato's reference + uses a non-stable `std::sort`, so on inputs with many ties the two + implementations may produce slightly different (but very similar) + partitions. See `development/graph/check_mutex_clustering.py` for a + comparison harness. + ## Dictionary-Based Relabeling If you used a small helper to apply a dictionary to an integer label array, use