diff --git a/AGENTS.md b/AGENTS.md index ed1d2c7..09b7ee7 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -1,459 +1,142 @@ # AGENTS.md — bioimage-cpp -This repository implements `bioimage-cpp`: a small C++/Python package for bioimage-analysis algorithms that are currently hard to install from existing C++-heavy libraries. The primary goal is to provide reliable PyPI wheels with a minimal dependency footprint. - -The package should be designed as a focused implementation of selected bioimage-analysis functionality, not as a compatibility layer for older libraries and not as a full reimplementation of large libraries such as `nifty` or `vigra`. +A small C++/Python package for bioimage-analysis algorithms that are otherwise hard to install. Goal: reliable PyPI wheels with a minimal dependency footprint. Not a compatibility layer for `nifty`/`vigra`/`affogato`, and not a full reimplementation — a focused implementation of selected functionality. ## Core principles -- Keep the package easy to build and easy to ship on PyPI. -- Prefer simple, dependency-light C++20 over large template libraries or external binary dependencies. -- Avoid dependencies that complicate wheels unless they are absolutely necessary. -- Do not introduce a C++ image container or multidimensional-array framework as a public dependency. -- Treat NumPy arrays as the Python-side data model and lightweight C++ array views as the C++-side data model. -- Keep I/O outside the C++ core. File formats such as TIFF, HDF5, zarr, N5, OME-NGFF, etc. should be handled by existing Python libraries, not by this package. -- Prefer explicit, narrow APIs over broad compatibility with old libraries. -- Keep the C++ algorithmic core independent of Python, NumPy, and nanobind. Nanobind should only be used in the binding layer. - -## Dependencies - -### Preferred dependencies - -- C++20 standard library. -- `nanobind` for Python bindings. -- `numpy` for array inputs and outputs. -- `scikit-build-core` as the Python build backend. -- `cmake` for configuring the C++ build. -- `pytest` for tests. -- `cibuildwheel` for wheel builds in CI. - -### Avoid in the core package - -Do not add these dependencies to the core package unless there is a strong reason and the wheel-building consequences have been considered carefully: - -- `xtensor` -- `Boost.Python` -- `pybind11` -- `vigra` -- `nifty` -- HDF5 -- z5 / N5 libraries -- libtiff, libpng, libjpeg, OpenEXR -- FFTW -- solver libraries with nontrivial binary dependencies -- conda-only dependencies -- OpenMP, unless the packaging implications are handled explicitly - -If optional functionality really requires a heavy dependency, isolate it behind an optional extension or a separate package. The default PyPI wheel should remain small and robust. - -## Suggested repository layout - -Use a `src` layout for the Python package and keep C++ headers, C++ implementation files, and Python bindings clearly separated. - -```text -bioimage-cpp/ -├── AGENTS.md -├── CMakeLists.txt -├── pyproject.toml -├── README.md -├── include/ -│ └── bioimage_cpp/ -│ ├── array_view.hxx -│ ├── dtype_dispatch.hxx -│ ├── errors.hxx -│ └── ... -├── src/ -│ ├── bioimage_cpp/ -│ │ ├── __init__.py -│ │ ├── _version.py -│ │ └── ... -│ ├── cpp/ -│ │ └── ... -│ └── bindings/ -│ ├── module.cxx -│ ├── array_conversions.cxx -│ └── ... -├── tests/ -│ ├── test_*.py -│ └── reference/ -│ └── ... -└── .github/ - └── workflows/ - ├── tests.yml - └── wheels.yml -``` - -The exact file names may change, but the separation should remain: - -- `include/bioimage_cpp/`: public or semi-public C++ headers. -- `src/cpp/`: C++ implementations. -- `src/bindings/`: Python binding code only. -- `src/bioimage_cpp/`: Python package code. -- `tests/`: Python-level tests and, where useful, C++ behavior tests exposed through Python. - -## Build system - -Use `scikit-build-core` with CMake. Avoid `setup.py`-based builds. +- Easy to build, easy to ship on PyPI. +- C++20 standard library only. Avoid template-heavy libraries and external binary dependencies. +- NumPy at the Python boundary; lightweight C++ array views internally. +- Keep I/O outside the C++ core (TIFF/HDF5/zarr/N5/OME-NGFF belong in Python). +- Keep the C++ algorithmic core independent of Python, NumPy, and nanobind — nanobind only in the binding layer. +- Explicit, narrow APIs. No compatibility namespaces for older libraries. -A minimal `pyproject.toml` should look conceptually like this: +## Build & test -```toml -[build-system] -requires = ["scikit-build-core", "nanobind"] -build-backend = "scikit_build_core.build" - -[project] -name = "bioimage-cpp" -requires-python = ">=3.10" -dependencies = ["numpy"] - -[tool.scikit-build] -cmake.version = ">=3.21" -wheel.packages = ["src/bioimage_cpp"] +```bash +pip install -e . --no-build-isolation # editable build (rebuilds C++) +python -m pytest tests/ -q # full Python test suite ``` -The package name on PyPI should be `bioimage-cpp`. The import name should be `bioimage_cpp`. +Required in the environment: `scikit-build-core`, `cmake>=3.21`, `ninja`, `nanobind`, `numpy`, `pytest`. The build produces one extension module: `bioimage_cpp._core`. Development scripts under `development/` compare against external references (nifty, affogato); they are not part of the test suite and must not be required by `pytest`. -If support for Python 3.10 and 3.11 is dropped later, the project may switch to Python >=3.12 and nanobind stable-ABI wheels. Do not enable nanobind `STABLE_ABI` while the package still supports Python <3.12. +## Reuse existing helpers — do not duplicate -## CMake guidelines +Before introducing union-finds, priority queues, edge hashing, stride math, threading helpers, or label relabeling, check `include/bioimage_cpp/detail/`: -- Build one main Python extension module, for example `bioimage_cpp._core`. -- Use C++20 consistently. -- Use nanobind, not pybind11. -- Keep CMake simple and portable. -- Avoid platform-specific compiler flags unless they are guarded carefully. -- Do not assume conda is available. -- Do not assume system libraries such as HDF5, TIFF, or Boost are available. -- Do not enable native architecture flags such as `-march=native` in wheel builds. -- Use hidden symbol visibility where appropriate. -- Keep the extension self-contained. +- `detail/union_find.hxx` — `UnionFind` (path compression + union by rank). `find`, `merge`, `merge_to`, `unite_roots`. +- `detail/indexed_heap.hxx` — addressable max-heap with mutable priorities. `DenseIndexedHeap` (vector-backed locator, integer keys in `[0, N)`) and `SparseIndexedHeap` (hashmap-backed, arbitrary keys). +- `detail/edge_hash.hxx` — `Edge`, `edge_key`, `EdgeHash` for hashing unordered node pairs. +- `detail/grid.hxx` — `c_order_strides`, `valid_offset_target`, `is_valid_grid_edge` for row-major grid offsets. +- `detail/relabel.hxx` — `dense_relabel` to map arbitrary labels onto `[0, k)` preserving first-occurrence order. +- `detail/threading.hxx` — `normalize_thread_count`, `parallel_for_chunks(n_threads, n_items, chunk)`. -A typical CMake structure should be: +If a needed helper does not exist but is generally useful, add it to `detail/` as a header-only utility with a focused API rather than inlining it into one algorithm. Keep the contract small and well-documented; over time `detail/` is what lets multiple modules stay small and consistent. -```cmake -cmake_minimum_required(VERSION 3.21) -project(bioimage_cpp LANGUAGES CXX) +## Dependencies -set(CMAKE_CXX_STANDARD 20) -set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_CXX_EXTENSIONS OFF) +**Allowed**: C++20 stdlib, `nanobind`, `numpy`, `scikit-build-core`, `cmake`, `pytest`, `cibuildwheel` (CI only). -find_package(Python 3.10 COMPONENTS Interpreter Development.Module REQUIRED) -find_package(nanobind CONFIG REQUIRED) +**Avoid in the core package** (require strong justification + wheel impact analysis): `xtensor`, `Boost.Python`, `pybind11`, `vigra`, `nifty`, HDF5, z5/N5, libtiff/libpng/libjpeg, OpenEXR, FFTW, heavy solver libraries, conda-only packages, and unguarded OpenMP. Isolate optional heavy features behind separate packages. -nanobind_add_module(_core - NB_STATIC - src/bindings/module.cxx - src/bindings/array_conversions.cxx - src/cpp/some_algorithm.cxx -) +## Repository layout -target_include_directories(_core PRIVATE include) -target_compile_features(_core PRIVATE cxx_std_20) ``` - -If the minimum Python version is raised to 3.12 or later, stable-ABI builds can be considered: - -```cmake -nanobind_add_module(_core - STABLE_ABI - NB_STATIC - src/bindings/module.cxx - src/bindings/array_conversions.cxx - src/cpp/some_algorithm.cxx -) +include/bioimage_cpp/ public/semi-public C++ headers (header-only algorithms) +include/bioimage_cpp/detail/ shared header-only utilities (see above) +src/cpp/ non-templated C++ implementation files +src/bindings/ nanobind binding code only (module.cxx + per-module .cxx) +src/bioimage_cpp/ Python package +tests/ pytest suite, runs against the installed package +development/ reference comparisons — not required by tests ``` -Adjust this as the project grows, but avoid turning the build into a large dependency-discovery system. - -## Wheel building - -Use `cibuildwheel` in GitHub Actions for binary wheels. - -The initial target platforms should be: - -- Linux x86_64, using manylinux wheels. -- macOS x86_64. -- macOS arm64. -- Windows x86_64. - -Later targets such as Linux aarch64 can be added once the core wheels are stable. - -Important wheel-building rules: - -- Wheels must not depend on conda. -- Wheels must not require users to have a local compiler at install time. -- Avoid external shared libraries unless they are bundled correctly and legally. -- Test the installed wheel, not only the source tree. -- Test importability with `python -c "import bioimage_cpp"`. -- Test core functions on small arrays after wheel installation. -- Do not use `-march=native`, unguarded OpenMP, or platform-specific compiler assumptions in release wheels. - -## Python API design - -The public Python API should be explicit and simple. Prefer functions such as: - -```python -import bioimage_cpp as bic - -result = bic.some_function(input_array, labels, ...) -``` - -Avoid exposing complex C++ classes unless they are genuinely needed. If a C++ object is exposed, its lifetime and ownership semantics must be obvious from Python. - -Do not implement a `vigra` or `nifty` compatibility namespace. The API should be designed for `bioimage-cpp` directly. - -## Multidimensional-array handling - -### General rule - -Use NumPy arrays at the Python boundary and lightweight C++ views internally. - -Do not use `xtensor`, `vigra::MultiArray`, or another multidimensional-array framework as a core dependency. +Build system is `scikit-build-core` + CMake. Build one extension module (`bioimage_cpp._core`). Do not use `setup.py`-style builds. Do not enable nanobind `STABLE_ABI` while Python <3.12 is supported. -### C++ array view - -Use a simple internal view type that stores: - -- data pointer -- shape -- strides -- number of dimensions -- dtype handled by template dispatch +## CMake guidelines -For example: +- C++20 throughout (`cxx_std_20`). +- Do not assume conda or system HDF5/TIFF/Boost. +- Do not enable `-march=native` in wheel builds. +- Use hidden symbol visibility where appropriate. +- Avoid platform-specific flags unless carefully guarded. -```cpp -#include -#include - -template -struct ArrayView { - T* data = nullptr; - std::vector shape; - std::vector strides; - - std::ptrdiff_t ndim() const { - return static_cast(shape.size()); - } -}; -``` +## Wheel building -For read-only inputs, use `const T*` or a separate `ConstArrayView`. +`cibuildwheel` in GitHub Actions. Targets: Linux x86_64 (manylinux), macOS x86_64 / arm64, Windows x86_64. Linux aarch64 once core wheels are stable. -The C++ algorithm implementation must not depend directly on Python, NumPy, or nanobind. +- No conda dependency. No compiler required at install time. +- No external shared libraries unless bundled correctly and legally. +- Test the installed wheel — importability and core functions on small arrays. +- No `-march=native`, no unguarded OpenMP, no platform-specific assumptions. -### Binding-layer responsibilities +## Python API -The Python binding layer is responsible for: +Prefer functions: `bic.some_function(input_array, labels, ...)`. Avoid exposing complex C++ classes unless genuinely needed; if exposed, lifetime semantics must be obvious. No `vigra`/`nifty` compatibility namespace. -- accepting `nanobind::ndarray` inputs -- checking dimensionality -- checking or converting dtype -- checking writability for output or in-place arguments -- checking shape compatibility between inputs -- checking C-contiguity if the C++ implementation requires it -- converting nanobind array metadata into `ArrayView` -- allocating output arrays -- dispatching to the correct C++ template instantiation +## Array handling -Use nanobind only at this boundary. Do not pass `nanobind::ndarray` or other nanobind types into the algorithmic core. +- NumPy at the Python boundary, `ArrayView` / `ConstArrayView` internally. The view carries `data`, `shape`, `strides`, `ndim()` only. +- Algorithm code must not depend on Python, NumPy, or nanobind. +- C-contiguous in C++ kernels; Python wrappers convert with `np.ascontiguousarray` when copying is acceptable. Support for arbitrary strides, if any, must be intentional and tested. +- NumPy axis order throughout. For shape `(z, y, x)`, coordinates are reported as `(z, y, x)` unless explicitly documented otherwise. +- Define supported dtypes explicitly per function. Typical: label arrays `uint32/uint64/int32/int64`; image arrays `uint8/uint16/uint32/float32/float64`; stats outputs `float64`; index outputs `int64/uint64`. Use explicit dtype dispatch in the binding layer; avoid large template instantiation matrices. -A typical binding file should use a local namespace alias: +### Memory ownership -```cpp -#include -#include +- Do not store raw pointers to array data past the call. Do not return views into temporaries. +- Prefer returning newly allocated arrays. For in-place outputs, require writable arrays and document mutation. +- If nanobind wraps memory owned by C++, attach an explicit capsule/owner so lifetimes are correct. -namespace nb = nanobind; -``` +### GIL -Example pattern: +Release the GIL only after all Python-object validation and metadata extraction: ```cpp -nb::ndarray, nb::c_contig> -some_function( - nb::ndarray, nb::c_contig> labels -) { - // Validate Python-facing arguments here. - // Convert labels to ConstArrayView. - // Release the GIL. - // Run the C++ implementation. - // Return a newly allocated output array. -} +{ nb::gil_scoped_release release; run_algorithm(...); } ``` -Prefer readable, explicit binding code over overly generic binding templates. - -### Contiguity policy - -Prefer a two-level policy: - -1. Fast C++ kernels assume C-contiguous arrays. -2. Python wrappers convert inputs with `np.ascontiguousarray` when copying is acceptable. - -For large arrays, avoid hidden expensive copies where possible. Public Python functions should make the copy policy clear, for example with an argument such as: - -```python -def some_function(array, *, copy=True): - ... -``` - -or by documenting that inputs are converted to contiguous arrays. - -If an algorithm supports arbitrary strides, this should be intentional and tested. Do not accidentally support some strided cases while failing on others. - -### Shape and axis conventions +Do not touch Python objects while the GIL is released. -Use NumPy conventions throughout: - -- arrays are indexed in row-major order -- dimensions are ordered as they appear in the NumPy array -- spatial axes should not be silently reordered -- coordinate outputs should be in NumPy axis order - -For an array with shape `(z, y, x)`, coordinates should be reported as `(z, y, x)` unless a function explicitly documents otherwise. - -### Dtype policy - -Do not silently accept every dtype. Each function should define supported dtypes explicitly. - -Common policies: - -- label arrays: `uint32`, `uint64`, `int32`, `int64` -- image arrays: `uint8`, `uint16`, `uint32`, `float32`, `float64` -- output statistics: usually `float64` for numeric stability, unless there is a reason to preserve dtype -- index arrays: prefer `int64` or `uint64` at the Python boundary - -Use explicit dtype dispatch in the binding layer. Avoid large, uncontrolled template instantiation matrices. - -### Memory ownership and lifetime - -- Never store raw pointers to array data beyond the duration of the function call unless the owning Python object is kept alive explicitly. -- Do not return C++ views into temporary arrays. -- Do not return references to internal buffers unless ownership is clear and tested. -- Prefer returning newly allocated arrays. -- For in-place operations, require writable arrays and document mutation clearly. -- If a nanobind array wraps memory owned by C++, attach an explicit owner/capsule so the lifetime is correct. - -### GIL handling +### Error handling -For long-running C++ kernels, release the Python GIL in the binding layer: +Throw `std::invalid_argument` / `std::runtime_error` from C++; nanobind translates them to Python exceptions. Messages must name the argument and report expected vs. actual shape/dtype. Example: `labels must have ndim >= 2, got ndim=1`. -```cpp -{ - nb::gil_scoped_release release; - run_algorithm(...); -} -``` +## Binding layer -Only release the GIL after all Python-object validation and array metadata extraction have completed. Do not touch Python objects while the GIL is released. +The binding layer (and only the binding layer) validates `nanobind::ndarray` inputs — dimensionality, dtype (or convert), writability for outputs, shape compatibility, C-contiguity when required — converts metadata into `ArrayView`, allocates outputs, dispatches to the right template instantiation, and releases the GIL. Do not pass `nanobind::ndarray` or other nanobind types into algorithm code. Use `namespace nb = nanobind;` at the top of binding files. Prefer readable, explicit binding code over overly generic templates. -### Error handling +## Testing -C++ code should throw standard exceptions such as `std::invalid_argument` or `std::runtime_error`. The binding layer should allow nanobind to translate these into Python exceptions. +Tests run against the installed Python package. For each public function, cover: small 2D / small 3D arrays, nontrivial labels, background-label handling, dtype variants, non-contiguous inputs (if accepted), invalid shapes, invalid dtypes, empty/degenerate inputs, and deterministic behavior. Reference comparisons against `nifty`/`affogato` are useful during development but must not be required by the default test suite, wheel tests, or source builds. -Error messages should include: +## Coding style -- the offending argument name -- expected dimensionality or dtype -- actual dimensionality or dtype -- expected shape relationship, if applicable +**C++**: clear, boring C++20. Small functions over template-heavy abstractions. `std::vector`/`std::array`/`std::span` where appropriate. Explicit integer types when overflow matters; careful with signed/unsigned conversions for shapes, strides, labels. Avoid clever metaprogramming. -Example: +**Python**: thin wrappers. Validate user-facing arguments in Python when this improves error messages. Use NumPy for lightweight preprocessing only. Avoid dependencies for small conveniences. -```text -labels must have ndim >= 2, got ndim=1 -``` +## Performance -## Testing guidelines +Correct first, fast second. Benchmark before adding complexity. Prefer algorithmic improvements over build-level tweaks. Release the GIL for expensive kernels. Add threading only after the single-threaded implementation is stable; it must be portable and user-controllable. -Tests should run against the installed Python package. +## Documentation -For each public function, test: +`README.md` covers: what `bioimage-cpp` is and isn't, install/build, minimal examples, design philosophy. Public functions need concise docstrings documenting input shapes, supported dtypes, output shapes/dtypes, copy behavior, background-label behavior (if relevant), and axis/coordinate conventions. -- small 2D arrays -- small 3D arrays where applicable -- nontrivial labels or regions -- background label behavior -- dtype variants -- non-contiguous inputs, if accepted -- incorrect shapes -- incorrect dtypes -- empty or degenerate inputs -- deterministic behavior +Update `MIGRATION_GUIDE.md` whenever public functionality changes, so users migrating from `nifty`/`affogato` see the corresponding `bioimage-cpp` API, behavioral differences, and intentional improvements. -Reference tests against older libraries may be useful during development, but they are not part of the public API contract. Do not require `vigra` or `nifty` for the default test suite, wheel tests, or source builds. +## Checklist when modifying code -## Coding style +1. Reuse existing `detail/` helpers; add new shared ones rather than duplicating. +2. Keep the build dependency-light. No external C++ dependencies without strong justification. +3. Algorithm code stays separate from binding code. nanobind only in the binding layer. +4. Validate arrays at the binding boundary. +5. Use internal array views — not external multidim-array libraries. +6. Tests for new behavior; tests run against the installed package. +7. Preserve PyPI wheel portability — no conda, no system libraries, no compiler at install time. +8. Clear error messages over permissive behavior. +9. No I/O format support in the C++ core. +10. No compatibility namespaces for older libraries. -### C++ - -- Use clear, boring C++20. -- Prefer small functions over large template-heavy abstractions. -- Keep algorithm code independent of nanobind. -- Use `std::vector`, `std::array`, `std::span` where available and appropriate. -- Avoid clever metaprogramming unless it clearly reduces maintenance burden. -- Use explicit integer types where overflow matters. -- Be careful with signed/unsigned conversions for shapes, strides, and labels. - -### Python - -- Keep Python wrappers thin. -- Validate user-facing arguments in Python when this improves error messages. -- Use NumPy for lightweight preprocessing only. -- Avoid adding dependencies for small convenience features. - -## Performance guidelines - -- Start with correct, portable implementations. -- Benchmark before introducing complexity. -- Prefer algorithmic improvements over build-level optimizations. -- Avoid `-march=native` in distributed wheels. -- Consider releasing the GIL for expensive kernels. -- Consider optional threading only after the single-threaded implementation is stable. -- If threading is added, make it portable and controllable by the user. - -## Documentation expectations - -The README should explain: - -- what `bioimage-cpp` is -- what it intentionally is not -- how to install it from PyPI -- how to build it from source -- minimal examples for the public API -- the dependency-light design philosophy - -Public functions should have concise docstrings documenting: - -- input shapes -- supported dtypes -- output shapes and dtypes -- copy behavior -- background-label behavior, if relevant -- axis and coordinate conventions - -When adding or changing public functionality, update `MIGRATION_GUIDE.md` so -users migrating from `nifty`, `affogato`, or related packages know the -corresponding `bioimage-cpp` API, important behavioral differences, and any -intentional improvements. - -## What an agent should do when modifying this repository - -When adding or changing code: - -1. Keep the build dependency-light. -2. Do not add external C++ dependencies without strong justification. -3. Keep algorithm code separate from binding code. -4. Validate arrays at the Python/binding boundary. -5. Use nanobind only in the binding layer. -6. Use internal C++ array views rather than external multidimensional-array libraries. -7. Add tests for new behavior. -8. Preserve PyPI wheel portability. -9. Prefer clear error messages over permissive behavior. -10. Do not add I/O format support to the C++ core. -11. Do not add compatibility namespaces for older libraries. - -When unsure, choose the simpler design that is easier to build, test, and ship as wheels. +When in doubt, choose the simpler design that is easier to build, test, and ship as wheels. diff --git a/CMakeLists.txt b/CMakeLists.txt index cdfa2c4..1ee5a7d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,5 +23,8 @@ nanobind_add_module(_core target_include_directories(_core PRIVATE include) target_compile_features(_core PRIVATE cxx_std_20) +target_compile_options(_core PRIVATE + $<$>:-O3> +) install(TARGETS _core LIBRARY DESTINATION bioimage_cpp) diff --git a/development/graph/multicut/_compatibility.py b/development/graph/multicut/_compatibility.py new file mode 100644 index 0000000..39812aa --- /dev/null +++ b/development/graph/multicut/_compatibility.py @@ -0,0 +1,134 @@ +from __future__ import annotations + +import argparse +from statistics import median +from time import perf_counter +from typing import Callable + +import numpy as np + + +def parser(description: str) -> argparse.ArgumentParser: + arg_parser = argparse.ArgumentParser(description=description) + arg_parser.add_argument( + "--path", + default=None, + help="Path to the external multicut problem. Defaults to the package cache.", + ) + arg_parser.add_argument( + "--repeats", + type=int, + default=3, + help="Number of timed repeats per implementation.", + ) + arg_parser.add_argument( + "--threads", + type=int, + default=1, + help="Number of threads for solvers that support it.", + ) + arg_parser.add_argument( + "--timeout", + type=float, + default=30.0, + help="Download timeout in seconds if the external problem is not cached.", + ) + arg_parser.add_argument( + "--energy-bound", + type=float, + default=-76900.0, + help="Maximum accepted energy for both implementations.", + ) + return arg_parser + + +def load_problem(path: str | None, *, timeout: float): + import bioimage_cpp as bic + import nifty + + uv_ids, costs = bic.graph.load_external_multicut_problem_data( + path, + timeout=timeout, + ) + bic_graph = bic.graph.UndirectedGraph.from_edges(int(uv_ids.max()) + 1, uv_ids) + nifty_graph = nifty.graph.undirectedGraph(int(uv_ids.max()) + 1) + nifty_graph.insertEdges(uv_ids) + return bic_graph, nifty_graph, costs + + +def time_call(function: Callable[[], np.ndarray], repeats: int): + timings = [] + result = None + for _ in range(repeats): + start = perf_counter() + result = function() + timings.append(perf_counter() - start) + assert result is not None + return timings, result + + +def optimize_bic_solver(make_bic_solver, objective): + objective.reset_labels() + return make_bic_solver().optimize(objective) + + +def bic_energy(graph, costs: np.ndarray, labels: np.ndarray) -> float: + import bioimage_cpp as bic + + return bic.graph.MulticutObjective(graph, costs).energy(labels) + + +def nifty_energy(graph, costs: np.ndarray, labels: np.ndarray) -> float: + import nifty.graph.opt.multicut as nmc + + return float(nmc.multicutObjective(graph, costs).evalNodeLabels(labels)) + + +def run_comparison( + name: str, + make_bic_solver, + make_nifty_solver, + args: argparse.Namespace, +) -> dict[str, float]: + import bioimage_cpp as bic + import nifty.graph.opt.multicut as nmc + + bic_graph, nifty_graph, costs = load_problem(args.path, timeout=args.timeout) + bic_objective = bic.graph.MulticutObjective(bic_graph, costs) + nifty_objective = nmc.multicutObjective(nifty_graph, costs) + + bic_timings, bic_labels = time_call( + lambda: optimize_bic_solver(make_bic_solver, bic_objective), + args.repeats, + ) + nifty_timings, nifty_labels = time_call( + lambda: make_nifty_solver(nifty_objective).create(nifty_objective).optimize(), + args.repeats, + ) + + bic_score = bic_energy(bic_graph, costs, bic_labels) + nifty_score = nifty_energy(nifty_graph, costs, nifty_labels) + if bic_score > args.energy_bound: + raise AssertionError( + f"bioimage-cpp {name} energy {bic_score:.6f} exceeds bound {args.energy_bound:.6f}" + ) + if nifty_score > args.energy_bound: + raise AssertionError( + f"nifty {name} energy {nifty_score:.6f} exceeds bound {args.energy_bound:.6f}" + ) + + result = { + "bioimage_cpp_energy": bic_score, + "nifty_energy": nifty_score, + "energy_difference": bic_score - nifty_score, + "bioimage_cpp_median_runtime": median(bic_timings), + "nifty_median_runtime": median(nifty_timings), + } + print(f"solver: {name}") + print(f"nodes: {bic_graph.number_of_nodes}, edges: {bic_graph.number_of_edges}") + print(f"bioimage-cpp energy: {bic_score:.6f}") + print(f"nifty energy: {nifty_score:.6f}") + print(f"energy difference: {bic_score - nifty_score:.6f}") + print(f"bioimage-cpp median runtime [s]: {median(bic_timings):.6f}") + print(f"nifty median runtime [s]: {median(nifty_timings):.6f}") + return result diff --git a/development/graph/multicut/check_chained.py b/development/graph/multicut/check_chained.py new file mode 100644 index 0000000..6832a45 --- /dev/null +++ b/development/graph/multicut/check_chained.py @@ -0,0 +1,29 @@ +from __future__ import annotations + +import bioimage_cpp as bic + +from _compatibility import parser, run_comparison + + +def main() -> None: + args = parser("Compare bioimage-cpp and nifty chained multicut solvers.").parse_args() + run_comparison( + "chained", + lambda: bic.graph.ChainedMulticutSolvers( + [ + bic.graph.GreedyAdditiveMulticut(), + bic.graph.KernighanLinMulticut(number_of_outer_iterations=5), + ] + ), + lambda objective: objective.chainedSolversFactory( + [ + objective.greedyAdditiveFactory(), + objective.kernighanLinFactory(numberOfOuterIterations=5), + ] + ), + args, + ) + + +if __name__ == "__main__": + main() diff --git a/development/graph/multicut/check_decomposer.py b/development/graph/multicut/check_decomposer.py new file mode 100644 index 0000000..01c24b6 --- /dev/null +++ b/development/graph/multicut/check_decomposer.py @@ -0,0 +1,23 @@ +from __future__ import annotations + +import bioimage_cpp as bic + +from _compatibility import parser, run_comparison + + +def main() -> None: + args = parser("Compare bioimage-cpp and nifty decomposer multicut.").parse_args() + run_comparison( + "decomposer", + lambda: bic.graph.MulticutDecomposer(bic.graph.GreedyAdditiveMulticut()), + lambda objective: objective.multicutDecomposerFactory( + submodelFactory=objective.greedyAdditiveFactory(), + fallthroughFactory=objective.greedyAdditiveFactory(), + numberOfThreads=args.threads, + ), + args, + ) + + +if __name__ == "__main__": + main() diff --git a/development/graph/multicut/check_greedy_additive.py b/development/graph/multicut/check_greedy_additive.py new file mode 100644 index 0000000..e9f7710 --- /dev/null +++ b/development/graph/multicut/check_greedy_additive.py @@ -0,0 +1,19 @@ +from __future__ import annotations + +import bioimage_cpp as bic + +from _compatibility import parser, run_comparison + + +def main() -> None: + args = parser("Compare bioimage-cpp and nifty greedy-additive multicut.").parse_args() + run_comparison( + "greedy_additive", + lambda: bic.graph.GreedyAdditiveMulticut(), + lambda objective: objective.greedyAdditiveFactory(), + args, + ) + + +if __name__ == "__main__": + main() diff --git a/development/graph/multicut/check_greedy_fixation.py b/development/graph/multicut/check_greedy_fixation.py new file mode 100644 index 0000000..29b033f --- /dev/null +++ b/development/graph/multicut/check_greedy_fixation.py @@ -0,0 +1,19 @@ +from __future__ import annotations + +import bioimage_cpp as bic + +from _compatibility import parser, run_comparison + + +def main() -> None: + args = parser("Compare bioimage-cpp and nifty greedy-fixation multicut.").parse_args() + run_comparison( + "greedy_fixation", + lambda: bic.graph.GreedyFixationMulticut(), + lambda objective: objective.greedyFixationFactory(), + args, + ) + + +if __name__ == "__main__": + main() diff --git a/development/graph/multicut/check_kernighan_lin.py b/development/graph/multicut/check_kernighan_lin.py new file mode 100644 index 0000000..843bcb2 --- /dev/null +++ b/development/graph/multicut/check_kernighan_lin.py @@ -0,0 +1,22 @@ +from __future__ import annotations + +import bioimage_cpp as bic + +from _compatibility import parser, run_comparison + + +def main() -> None: + args = parser("Compare bioimage-cpp and nifty Kernighan-Lin multicut.").parse_args() + run_comparison( + "kernighan_lin", + lambda: bic.graph.KernighanLinMulticut(number_of_outer_iterations=5), + lambda objective: objective.kernighanLinFactory( + warmStartGreedy=True, + numberOfOuterIterations=5, + ), + args, + ) + + +if __name__ == "__main__": + main() diff --git a/include/bioimage_cpp/detail/edge_hash.hxx b/include/bioimage_cpp/detail/edge_hash.hxx new file mode 100644 index 0000000..25ec4c2 --- /dev/null +++ b/include/bioimage_cpp/detail/edge_hash.hxx @@ -0,0 +1,27 @@ +#pragma once + +#include +#include +#include + +namespace bioimage_cpp::detail { + +using NodeId = std::uint64_t; +using Edge = std::pair; + +inline Edge edge_key(NodeId u, NodeId v) { + if (v < u) { + std::swap(u, v); + } + return {u, v}; +} + +struct EdgeHash { + std::size_t operator()(const Edge &edge) const { + const auto first = static_cast(edge.first); + const auto second = static_cast(edge.second); + return first ^ (second + 0x9e3779b97f4a7c15ULL + (first << 6U) + (first >> 2U)); + } +}; + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/grid.hxx b/include/bioimage_cpp/detail/grid.hxx new file mode 100644 index 0000000..4c0b4f8 --- /dev/null +++ b/include/bioimage_cpp/detail/grid.hxx @@ -0,0 +1,59 @@ +#pragma once + +#include +#include +#include + +namespace bioimage_cpp::detail { + +// C-order strides for a row-major array of the given shape, in units of array +// elements (not bytes). The innermost (last) axis has stride 1. +inline std::vector c_order_strides(const std::vector &shape) { + std::vector strides(shape.size(), 1); + for (std::ptrdiff_t axis = static_cast(shape.size()) - 2; axis >= 0; --axis) { + strides[static_cast(axis)] = + strides[static_cast(axis + 1)] * + shape[static_cast(axis + 1)]; + } + return strides; +} + +// Translate a flat node index by a per-axis offset on a row-major grid. +// +// Returns true when the offset keeps the result inside the grid, in which case +// `target_out` is set to the neighbor's flat index. Returns false otherwise and +// leaves `target_out` unchanged. `strides` must match `shape` and is typically +// `c_order_strides(shape)`. +inline bool valid_offset_target( + const std::uint64_t node, + const std::vector &offset, + const std::vector &shape, + const std::vector &strides, + std::uint64_t &target_out +) { + std::int64_t target_signed = static_cast(node); + for (std::size_t axis = 0; axis < shape.size(); ++axis) { + const auto coord = + static_cast(node / static_cast(strides[axis])) % + shape[axis]; + const auto neighbor = coord + offset[axis]; + if (neighbor < 0 || neighbor >= shape[axis]) { + return false; + } + target_signed += static_cast(offset[axis] * strides[axis]); + } + target_out = static_cast(target_signed); + return true; +} + +inline bool is_valid_grid_edge( + const std::uint64_t node, + const std::vector &offset, + const std::vector &shape, + const std::vector &strides +) { + std::uint64_t unused = 0; + return valid_offset_target(node, offset, shape, strides, unused); +} + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/indexed_heap.hxx b/include/bioimage_cpp/detail/indexed_heap.hxx new file mode 100644 index 0000000..82de272 --- /dev/null +++ b/include/bioimage_cpp/detail/indexed_heap.hxx @@ -0,0 +1,278 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::detail { + +// Sentinel used by both locators below to mark "this key is not in the heap". +inline constexpr std::size_t indexed_heap_not_in_heap = + std::numeric_limits::max(); + +// Locator concept used by IndexedHeap: +// std::size_t at(key) const — returns position in heap, or `not_in_heap`. +// void set(key, std::size_t pos) — records a key's position. `not_in_heap` clears it. + +// Locator for keys drawn from [0, capacity). Backed by a flat vector so the +// position lookup is a single load. +class DenseLocator { +public: + static constexpr std::size_t not_in_heap = indexed_heap_not_in_heap; + + explicit DenseLocator(const std::size_t capacity = 0) + : positions_(capacity, not_in_heap) {} + + [[nodiscard]] std::size_t at(const std::size_t key) const { + return positions_[key]; + } + + void set(const std::size_t key, const std::size_t pos) { + positions_[key] = pos; + } + + // Resize and reset every position to `not_in_heap`. Use between solver + // invocations on different graphs. + void reset_capacity(const std::size_t capacity) { + positions_.assign(capacity, not_in_heap); + } + +private: + std::vector positions_; +}; + +// Locator for keys with sparse / arbitrary identity. Position lookup is a +// hash-map probe. +template > +class SparseLocator { +public: + static constexpr std::size_t not_in_heap = indexed_heap_not_in_heap; + + [[nodiscard]] std::size_t at(const KeyT &key) const { + const auto it = positions_.find(key); + return it == positions_.end() ? not_in_heap : it->second; + } + + void set(const KeyT &key, const std::size_t pos) { + if (pos == not_in_heap) { + positions_.erase(key); + } else { + positions_.insert_or_assign(key, pos); + } + } + + void reserve(const std::size_t expected_keys) { + positions_.reserve(expected_keys); + } + + void clear() { + positions_.clear(); + } + +private: + std::unordered_map positions_; +}; + +// Addressable max-heap with mutable priorities. +// +// Each key maps to exactly one entry in the heap, tracked via the Locator. +// `push`, `change`, `push_or_change`, `erase`, and `pop` are O(log n) and keep +// the locator in sync with the heap permutation. `top`, `contains`, and +// `priority_of` are O(1). Compared to a `std::priority_queue` with lazy +// invalidation, the heap never carries stale entries — its size equals the +// number of currently active keys. +// +// The comparator `Compare` follows the standard "less-than → max-heap" +// convention used by `std::priority_queue`. Specialize with `std::greater<>` +// for a min-heap. +template < + class KeyT, + class PriorityT, + class Locator, + class Compare = std::less +> +class IndexedHeap { +public: + using key_type = KeyT; + using priority_type = PriorityT; + + struct Entry { + KeyT key; + PriorityT priority; + }; + + static constexpr std::size_t not_in_heap = Locator::not_in_heap; + + IndexedHeap() = default; + explicit IndexedHeap(Locator locator) : locator_(std::move(locator)) {} + + [[nodiscard]] std::size_t size() const { return heap_.size(); } + [[nodiscard]] bool empty() const { return heap_.empty(); } + + [[nodiscard]] bool contains(const KeyT &key) const { + return locator_.at(key) != not_in_heap; + } + + [[nodiscard]] const PriorityT &priority_of(const KeyT &key) const { + return heap_[locator_.at(key)].priority; + } + + [[nodiscard]] const Entry &top() const { + return heap_.front(); + } + + // Precondition: `key` is not currently in the heap. + void push(KeyT key, PriorityT priority) { + const auto pos = heap_.size(); + heap_.push_back({std::move(key), std::move(priority)}); + locator_.set(heap_.back().key, pos); + sift_up(pos); + } + + // Precondition: `key` is currently in the heap. + void change(const KeyT &key, PriorityT priority) { + const auto pos = locator_.at(key); + apply_change(pos, std::move(priority)); + } + + void push_or_change(KeyT key, PriorityT priority) { + const auto pos = locator_.at(key); + if (pos == not_in_heap) { + push(std::move(key), std::move(priority)); + } else { + apply_change(pos, std::move(priority)); + } + } + + // No-op when the key is not in the heap. This is the "remove if exists" + // semantic every caller needs. + void erase(const KeyT &key) { + const auto pos = locator_.at(key); + if (pos != not_in_heap) { + erase_at(pos); + } + } + + Entry pop() { + Entry top = heap_.front(); + erase_at(0); + return top; + } + + void clear() { + for (const auto &entry : heap_) { + locator_.set(entry.key, not_in_heap); + } + heap_.clear(); + } + + Locator &locator() { return locator_; } + const Locator &locator() const { return locator_; } + +private: + std::vector heap_; + Locator locator_; + Compare compare_; + + void swap_positions(const std::size_t a, const std::size_t b) { + std::swap(heap_[a], heap_[b]); + locator_.set(heap_[a].key, a); + locator_.set(heap_[b].key, b); + } + + void sift_up(std::size_t pos) { + while (pos > 0) { + const auto parent = (pos - 1) / 2; + if (!compare_(heap_[parent].priority, heap_[pos].priority)) { + break; + } + swap_positions(parent, pos); + pos = parent; + } + } + + void sift_down(std::size_t pos) { + const auto n = heap_.size(); + while (true) { + const auto left = 2 * pos + 1; + if (left >= n) { + break; + } + auto target = left; + const auto right = left + 1; + if (right < n && compare_(heap_[left].priority, heap_[right].priority)) { + target = right; + } + if (!compare_(heap_[pos].priority, heap_[target].priority)) { + break; + } + swap_positions(pos, target); + pos = target; + } + } + + void apply_change(const std::size_t pos, PriorityT new_priority) { + const bool increased = compare_(heap_[pos].priority, new_priority); + heap_[pos].priority = std::move(new_priority); + if (increased) { + sift_up(pos); + } else { + sift_down(pos); + } + } + + void erase_at(const std::size_t pos) { + const auto last = heap_.size() - 1; + locator_.set(heap_[pos].key, not_in_heap); + if (pos != last) { + heap_[pos] = std::move(heap_[last]); + locator_.set(heap_[pos].key, pos); + heap_.pop_back(); + // The replacement could need sifting in either direction. + sift_up(pos); + sift_down(pos); + } else { + heap_.pop_back(); + } + } +}; + +template > +class DenseIndexedHeap final + : public IndexedHeap { + using base = IndexedHeap; + +public: + DenseIndexedHeap() = default; + explicit DenseIndexedHeap(const std::size_t capacity) + : base(DenseLocator(capacity)) {} + + void reset_capacity(const std::size_t capacity) { + this->clear(); + this->locator().reset_capacity(capacity); + } +}; + +template < + class KeyT, + class PriorityT, + class Hash = std::hash, + class Compare = std::less +> +class SparseIndexedHeap final + : public IndexedHeap, Compare> { + using base = IndexedHeap, Compare>; + +public: + SparseIndexedHeap() = default; + + void reserve(const std::size_t expected_keys) { + this->locator().reserve(expected_keys); + } +}; + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/relabel.hxx b/include/bioimage_cpp/detail/relabel.hxx new file mode 100644 index 0000000..7edd3d1 --- /dev/null +++ b/include/bioimage_cpp/detail/relabel.hxx @@ -0,0 +1,25 @@ +#pragma once + +#include +#include +#include +#include + +namespace bioimage_cpp::detail { + +// Map an arbitrary labeling to the dense range [0, k) preserving the first +// occurrence order of each distinct input label. +inline std::vector dense_relabel(const std::vector &labels) { + std::unordered_map relabeling; + std::vector result(labels.size()); + for (std::size_t index = 0; index < labels.size(); ++index) { + auto found = relabeling.find(labels[index]); + if (found == relabeling.end()) { + found = relabeling.emplace(labels[index], static_cast(relabeling.size())).first; + } + result[index] = found->second; + } + return result; +} + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/threading.hxx b/include/bioimage_cpp/detail/threading.hxx new file mode 100644 index 0000000..79df7db --- /dev/null +++ b/include/bioimage_cpp/detail/threading.hxx @@ -0,0 +1,60 @@ +#pragma once + +#include +#include +#include +#include +#include + +namespace bioimage_cpp::detail { + +inline std::size_t normalize_thread_count( + const std::size_t requested, + const std::size_t number_of_work_items +) { + if (number_of_work_items == 0) { + return 1; + } + std::size_t n_threads = requested; + if (n_threads == 0) { + n_threads = std::thread::hardware_concurrency(); + if (n_threads == 0) { + n_threads = 1; + } + } + return std::max(1, std::min(n_threads, number_of_work_items)); +} + +// Split [0, number_of_work_items) into n_threads contiguous chunks and invoke +// `chunk(thread_id, begin, end)` on each chunk. Thread 0 runs on the calling +// thread; threads 1..n_threads-1 run in parallel std::threads. The caller is +// responsible for picking n_threads via normalize_thread_count and for the +// thread safety of `chunk`. +template +void parallel_for_chunks( + const std::size_t n_threads, + const std::size_t number_of_work_items, + Chunk &&chunk +) { + const auto bounds = [&](const std::size_t thread_id) { + const auto begin = thread_id * number_of_work_items / n_threads; + const auto end = (thread_id + 1) * number_of_work_items / n_threads; + return std::pair{begin, end}; + }; + + std::vector threads; + threads.reserve(n_threads > 0 ? n_threads - 1 : 0); + for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) { + const auto [begin, end] = bounds(thread_id); + threads.emplace_back([thread_id, begin, end, &chunk]() { + chunk(thread_id, begin, end); + }); + } + const auto [begin, end] = bounds(0); + chunk(0, begin, end); + for (auto &thread : threads) { + thread.join(); + } +} + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/union_find.hxx b/include/bioimage_cpp/detail/union_find.hxx new file mode 100644 index 0000000..f05fa4e --- /dev/null +++ b/include/bioimage_cpp/detail/union_find.hxx @@ -0,0 +1,78 @@ +#pragma once + +#include +#include +#include +#include +#include + +namespace bioimage_cpp::detail { + +// Disjoint-set / union-find with path compression and union-by-rank. +// +// Three merge variants are provided so both rank-driven and caller-driven +// merge directions are expressible: +// * merge(u, v): find both, then union by rank. +// * merge_to(stable, removed): find both, then force `stable` to be the +// parent regardless of rank. The caller picks +// which root survives. +// * unite_roots(a, b): preconditions: a and b are roots and a != b. Union +// by rank without an additional find(). +class UnionFind { +public: + explicit UnionFind(const std::size_t size) : parents_(size), ranks_(size, 0) { + std::iota(parents_.begin(), parents_.end(), std::uint64_t{0}); + } + + std::uint64_t find(const std::uint64_t node) { + auto current = static_cast(node); + while (parents_[current] != current) { + parents_[current] = parents_[parents_[current]]; + current = static_cast(parents_[current]); + } + return static_cast(current); + } + + std::uint64_t merge(std::uint64_t first, std::uint64_t second) { + first = find(first); + second = find(second); + if (first == second) { + return first; + } + return unite_roots(first, second); + } + + std::uint64_t merge_to(std::uint64_t stable, std::uint64_t removed) { + stable = find(stable); + removed = find(removed); + if (stable == removed) { + return stable; + } + parents_[static_cast(removed)] = stable; + if (ranks_[static_cast(stable)] <= ranks_[static_cast(removed)]) { + ranks_[static_cast(stable)] = ranks_[static_cast(removed)] + 1; + } + return stable; + } + + std::uint64_t unite_roots(std::uint64_t first, std::uint64_t second) { + if (ranks_[static_cast(first)] < ranks_[static_cast(second)]) { + std::swap(first, second); + } + parents_[static_cast(second)] = first; + if (ranks_[static_cast(first)] == ranks_[static_cast(second)]) { + ++ranks_[static_cast(first)]; + } + return first; + } + + [[nodiscard]] std::size_t size() const { + return parents_.size(); + } + +private: + std::vector parents_; + std::vector ranks_; +}; + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/graph/connected_components.hxx b/include/bioimage_cpp/graph/connected_components.hxx new file mode 100644 index 0000000..7f93fc5 --- /dev/null +++ b/include/bioimage_cpp/graph/connected_components.hxx @@ -0,0 +1,45 @@ +#pragma once + +#include "bioimage_cpp/detail/union_find.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" + +#include +#include +#include +#include + +namespace bioimage_cpp::graph { + +inline std::vector dense_labels_from_union_find( + detail::UnionFind &sets, + const std::uint64_t number_of_nodes +) { + std::unordered_map relabeling; + std::vector labels(static_cast(number_of_nodes)); + for (std::uint64_t node = 0; node < number_of_nodes; ++node) { + const auto root = sets.find(node); + auto found = relabeling.find(root); + if (found == relabeling.end()) { + found = relabeling.emplace(root, static_cast(relabeling.size())).first; + } + labels[static_cast(node)] = found->second; + } + return labels; +} + +inline std::vector connected_components( + const UndirectedGraph &graph, + const std::uint8_t *edge_mask = nullptr +) { + detail::UnionFind sets(static_cast(graph.number_of_nodes())); + for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) { + if (edge_mask != nullptr && edge_mask[static_cast(edge)] == 0) { + continue; + } + const auto uv = graph.uv(edge); + sets.merge(uv.first, uv.second); + } + return dense_labels_from_union_find(sets, graph.number_of_nodes()); +} + +} // namespace bioimage_cpp::graph diff --git a/include/bioimage_cpp/graph/feature_accumulation.hxx b/include/bioimage_cpp/graph/feature_accumulation.hxx index 1b98a0f..edf4297 100644 --- a/include/bioimage_cpp/graph/feature_accumulation.hxx +++ b/include/bioimage_cpp/graph/feature_accumulation.hxx @@ -1,6 +1,8 @@ #pragma once #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/detail/threading.hxx" #include "bioimage_cpp/graph/region_adjacency_graph.hxx" #include @@ -11,7 +13,6 @@ #include #include #include -#include #include #include @@ -123,15 +124,6 @@ inline std::int64_t edge_for_labels( return rag.find_edge(u, v); } -inline std::vector c_order_strides(const std::vector &shape) { - std::vector strides(shape.size(), 1); - for (std::ptrdiff_t axis = static_cast(shape.size()) - 2; axis >= 0; --axis) { - strides[static_cast(axis)] = - strides[static_cast(axis + 1)] * shape[static_cast(axis + 1)]; - } - return strides; -} - template void scan_edge_map_2d_chunk( const RegionAdjacencyGraph &rag, @@ -212,26 +204,6 @@ void scan_edge_map_3d_chunk( } } -inline bool valid_offset_target( - const std::uint64_t node, - const std::vector &offset, - const std::vector &shape, - const std::vector &strides, - std::uint64_t &target -) { - std::int64_t target_signed = static_cast(node); - for (std::size_t axis = 0; axis < shape.size(); ++axis) { - const auto coord = static_cast(node / static_cast(strides[axis])) % shape[axis]; - const auto neighbor = coord + offset[axis]; - if (neighbor < 0 || neighbor >= shape[axis]) { - return false; - } - target_signed += static_cast(offset[axis] * strides[axis]); - } - target = static_cast(target_signed); - return true; -} - template void scan_affinity_chunk( const RegionAdjacencyGraph &rag, @@ -243,13 +215,13 @@ void scan_affinity_chunk( const std::size_t node_end, std::vector &stats ) { - const auto spatial_strides = c_order_strides(shape); + const auto spatial_strides = bioimage_cpp::detail::c_order_strides(shape); const auto number_of_nodes = static_cast(number_of_pixels(shape)); for (std::size_t channel = 0; channel < offsets.size(); ++channel) { const auto channel_offset = static_cast(channel) * number_of_nodes; for (std::uint64_t node = node_begin; node < node_end; ++node) { std::uint64_t target = 0; - if (!valid_offset_target(node, offsets[channel], shape, spatial_strides, target)) { + if (!bioimage_cpp::detail::valid_offset_target(node, offsets[channel], shape, spatial_strides, target)) { continue; } const auto edge = edge_for_labels(rag, label_at(labels, node), label_at(labels, target)); @@ -355,76 +327,48 @@ void accumulate_edge_map_features( const auto work_items = static_cast(labels.shape[0]); const auto n_threads = detail::normalize_thread_count(number_of_threads, work_items); - std::vector threads; - threads.reserve(n_threads > 0 ? n_threads - 1 : 0); + const auto number_of_edges = static_cast(rag.number_of_edges()); + + const auto run_scan = [&](auto &per_thread_stats) { + bioimage_cpp::detail::parallel_for_chunks( + n_threads, + work_items, + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + if (labels.ndim() == 2) { + detail_features::scan_edge_map_2d_chunk( + rag, labels.data, edge_map.data, + static_cast(labels.shape[0]), + static_cast(labels.shape[1]), + begin, end, per_thread_stats[thread_id] + ); + } else { + detail_features::scan_edge_map_3d_chunk( + rag, labels.data, edge_map.data, + static_cast(labels.shape[0]), + static_cast(labels.shape[1]), + static_cast(labels.shape[2]), + begin, end, per_thread_stats[thread_id] + ); + } + } + ); + }; if (compute_complex_features) { std::vector> per_thread_stats( n_threads, - std::vector(static_cast(rag.number_of_edges())) + std::vector(number_of_edges) ); - const auto run_chunk = [&](const std::size_t thread_id) { - const auto begin = thread_id * work_items / n_threads; - const auto end = (thread_id + 1) * work_items / n_threads; - if (labels.ndim() == 2) { - detail_features::scan_edge_map_2d_chunk( - rag, labels.data, edge_map.data, - static_cast(labels.shape[0]), - static_cast(labels.shape[1]), - begin, end, per_thread_stats[thread_id] - ); - } else { - detail_features::scan_edge_map_3d_chunk( - rag, labels.data, edge_map.data, - static_cast(labels.shape[0]), - static_cast(labels.shape[1]), - static_cast(labels.shape[2]), - begin, end, per_thread_stats[thread_id] - ); - } - }; - for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) { - threads.emplace_back(run_chunk, thread_id); - } - run_chunk(0); - for (auto &thread : threads) { - thread.join(); - } - auto stats = detail_features::merge_stats(per_thread_stats, static_cast(rag.number_of_edges())); + run_scan(per_thread_stats); + auto stats = detail_features::merge_stats(per_thread_stats, number_of_edges); detail_features::write_complex_features(stats, out); } else { std::vector> per_thread_stats( n_threads, - std::vector(static_cast(rag.number_of_edges())) + std::vector(number_of_edges) ); - const auto run_chunk = [&](const std::size_t thread_id) { - const auto begin = thread_id * work_items / n_threads; - const auto end = (thread_id + 1) * work_items / n_threads; - if (labels.ndim() == 2) { - detail_features::scan_edge_map_2d_chunk( - rag, labels.data, edge_map.data, - static_cast(labels.shape[0]), - static_cast(labels.shape[1]), - begin, end, per_thread_stats[thread_id] - ); - } else { - detail_features::scan_edge_map_3d_chunk( - rag, labels.data, edge_map.data, - static_cast(labels.shape[0]), - static_cast(labels.shape[1]), - static_cast(labels.shape[2]), - begin, end, per_thread_stats[thread_id] - ); - } - }; - for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) { - threads.emplace_back(run_chunk, thread_id); - } - run_chunk(0); - for (auto &thread : threads) { - thread.join(); - } - auto stats = detail_features::merge_stats(per_thread_stats, static_cast(rag.number_of_edges())); + run_scan(per_thread_stats); + auto stats = detail_features::merge_stats(per_thread_stats, number_of_edges); detail_features::write_simple_features(stats, out); } } @@ -466,52 +410,36 @@ void accumulate_affinity_features( const auto number_of_nodes = detail_features::number_of_pixels(labels.shape); const auto n_threads = detail::normalize_thread_count(number_of_threads, number_of_nodes); - std::vector threads; - threads.reserve(n_threads > 0 ? n_threads - 1 : 0); + const auto number_of_edges = static_cast(rag.number_of_edges()); + + const auto run_scan = [&](auto &per_thread_stats) { + bioimage_cpp::detail::parallel_for_chunks( + n_threads, + number_of_nodes, + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + detail_features::scan_affinity_chunk( + rag, labels.data, affinities.data, offsets, labels.shape, + begin, end, per_thread_stats[thread_id] + ); + } + ); + }; if (compute_complex_features) { std::vector> per_thread_stats( n_threads, - std::vector(static_cast(rag.number_of_edges())) + std::vector(number_of_edges) ); - const auto run_chunk = [&](const std::size_t thread_id) { - const auto begin = thread_id * number_of_nodes / n_threads; - const auto end = (thread_id + 1) * number_of_nodes / n_threads; - detail_features::scan_affinity_chunk( - rag, labels.data, affinities.data, offsets, labels.shape, - begin, end, per_thread_stats[thread_id] - ); - }; - for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) { - threads.emplace_back(run_chunk, thread_id); - } - run_chunk(0); - for (auto &thread : threads) { - thread.join(); - } - auto stats = detail_features::merge_stats(per_thread_stats, static_cast(rag.number_of_edges())); + run_scan(per_thread_stats); + auto stats = detail_features::merge_stats(per_thread_stats, number_of_edges); detail_features::write_complex_features(stats, out); } else { std::vector> per_thread_stats( n_threads, - std::vector(static_cast(rag.number_of_edges())) + std::vector(number_of_edges) ); - const auto run_chunk = [&](const std::size_t thread_id) { - const auto begin = thread_id * number_of_nodes / n_threads; - const auto end = (thread_id + 1) * number_of_nodes / n_threads; - detail_features::scan_affinity_chunk( - rag, labels.data, affinities.data, offsets, labels.shape, - begin, end, per_thread_stats[thread_id] - ); - }; - for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) { - threads.emplace_back(run_chunk, thread_id); - } - run_chunk(0); - for (auto &thread : threads) { - thread.join(); - } - auto stats = detail_features::merge_stats(per_thread_stats, static_cast(rag.number_of_edges())); + run_scan(per_thread_stats); + auto stats = detail_features::merge_stats(per_thread_stats, number_of_edges); detail_features::write_simple_features(stats, out); } } diff --git a/include/bioimage_cpp/graph/multicut.hxx b/include/bioimage_cpp/graph/multicut.hxx new file mode 100644 index 0000000..1e75948 --- /dev/null +++ b/include/bioimage_cpp/graph/multicut.hxx @@ -0,0 +1,6 @@ +#pragma once + +#include "bioimage_cpp/graph/multicut/greedy_additive.hxx" +#include "bioimage_cpp/graph/multicut/greedy_fixation.hxx" +#include "bioimage_cpp/graph/multicut/kernighan_lin.hxx" +#include "bioimage_cpp/graph/multicut/objective.hxx" diff --git a/include/bioimage_cpp/graph/multicut/detail.hxx b/include/bioimage_cpp/graph/multicut/detail.hxx new file mode 100644 index 0000000..76ab5c9 --- /dev/null +++ b/include/bioimage_cpp/graph/multicut/detail.hxx @@ -0,0 +1,264 @@ +#pragma once + +#include "bioimage_cpp/detail/indexed_heap.hxx" +#include "bioimage_cpp/detail/union_find.hxx" +#include "bioimage_cpp/graph/connected_components.hxx" +#include "bioimage_cpp/graph/multicut/objective.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" + +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph::multicut::detail { + +inline constexpr std::size_t no_edge = std::numeric_limits::max(); + +// Per-super-node adjacency entry. +struct NeighborEntry { + std::size_t neighbor; + std::size_t edge_id; +}; + +// Per-edge data stored in a flat vector indexed by stable edge_id. Edge ids +// never grow past the original number of input edges — when two edges fold +// into one during a merge, one id survives and the other becomes orphaned +// (no adjacency entry or heap entry references it again). +struct DynamicEdge { + std::size_t u = 0; + std::size_t v = 0; + double weight = 0.0; + unsigned char is_constraint = 0; +}; + +// Heap keyed by dense stable edge_id. The vector-backed DenseLocator updates +// in O(1) per sift step, which was the missing ingredient for greedy_additive +// to beat the std::priority_queue baseline. +using EdgeHeap = bioimage_cpp::detail::DenseIndexedHeap; + +struct DynamicGraph { + explicit DynamicGraph(const UndirectedGraph &graph) + : adjacency(static_cast(graph.number_of_nodes())), + alive(static_cast(graph.number_of_nodes()), true), + alive_count(static_cast(graph.number_of_nodes())), + scratch_edge_id(static_cast(graph.number_of_nodes()), no_edge) { + for (std::uint64_t node = 0; node < graph.number_of_nodes(); ++node) { + const auto degree = graph.node_adjacency(node).size(); + adjacency[static_cast(node)].reserve(degree); + } + edges.resize(static_cast(graph.number_of_edges())); + } + + // O(degree(u)). Returns no_edge when (u, v) is not an edge. + [[nodiscard]] std::size_t find_edge(const std::size_t u, const std::size_t v) const { + if (u >= adjacency.size() || v >= adjacency.size()) { + return no_edge; + } + for (const auto &entry : adjacency[u]) { + if (entry.neighbor == v) { + return entry.edge_id; + } + } + return no_edge; + } + + [[nodiscard]] bool edge_exists(const std::size_t u, const std::size_t v) const { + return find_edge(u, v) != no_edge; + } + + [[nodiscard]] bool has_constraint(const std::size_t u, const std::size_t v) const { + const auto id = find_edge(u, v); + return id != no_edge && edges[id].is_constraint != 0; + } + + std::vector> adjacency; + std::vector edges; + std::vector alive; + std::size_t alive_count; + // Sized to #super-nodes, all entries == no_edge between merges. Used by + // merge_dynamic_nodes to find the existing stable-side edge for each of + // removed's neighbors in O(1) without hashing. + std::vector scratch_edge_id; +}; + +inline double priority_for(const double weight, const bool absolute_priority) { + return absolute_priority ? std::abs(weight) : weight; +} + +inline void initialize_dynamic_graph( + const UndirectedGraph &graph, + const std::vector &costs, + DynamicGraph &dynamic_graph, + EdgeHeap &heap, + const bool absolute_priority, + const bool add_noise = false, + const int seed = 42, + const double sigma = 1.0 +) { + const auto n_edges = static_cast(graph.number_of_edges()); + heap.reset_capacity(n_edges); + + std::mt19937 generator(seed); + std::normal_distribution noise(0.0, sigma); + for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) { + const auto uv = graph.uv(edge); + const auto u = static_cast(uv.first); + const auto v = static_cast(uv.second); + double weight = costs[static_cast(edge)]; + if (add_noise) { + weight += noise(generator); + } + const auto edge_id = static_cast(edge); + auto &e = dynamic_graph.edges[edge_id]; + e.u = u; + e.v = v; + e.weight = weight; + e.is_constraint = 0; + dynamic_graph.adjacency[u].push_back({v, edge_id}); + dynamic_graph.adjacency[v].push_back({u, edge_id}); + heap.push(edge_id, priority_for(weight, absolute_priority)); + } +} + +inline std::vector labels_from_sets( + bioimage_cpp::detail::UnionFind &sets, + const UndirectedGraph &graph +) { + return dense_labels_from_union_find(sets, graph.number_of_nodes()); +} + +inline std::size_t stop_node_count(const UndirectedGraph &graph, const double node_num_stop) { + return node_num_stop >= 1.0 + ? static_cast(node_num_stop) + : static_cast(double(graph.number_of_nodes()) * node_num_stop + 0.5); +} + +namespace internal { + +// Erase the first entry with neighbor == target from `list` (swap-with-back). +// Returns whether anything was erased. +inline bool erase_by_neighbor(std::vector &list, const std::size_t target) { + for (std::size_t i = 0; i < list.size(); ++i) { + if (list[i].neighbor == target) { + list[i] = list.back(); + list.pop_back(); + return true; + } + } + return false; +} + +// Rename the first entry whose neighbor == from_node to point to to_node. +inline void rename_neighbor( + std::vector &list, + const std::size_t from_node, + const std::size_t to_node +) { + for (auto &entry : list) { + if (entry.neighbor == from_node) { + entry.neighbor = to_node; + return; + } + } +} + +} // namespace internal + +// Contract the edge between u and v in the dynamic graph. The smaller-degree +// super-node is folded into the larger-degree one to keep the per-super-node +// adjacency growth amortized. The heap is kept in sync without staleness: each +// edge id appears at most once. +inline std::size_t merge_dynamic_nodes( + DynamicGraph &dynamic_graph, + bioimage_cpp::detail::UnionFind &sets, + EdgeHeap &heap, + std::size_t u, + std::size_t v, + const bool absolute_priority +) { + u = static_cast(sets.find(u)); + v = static_cast(sets.find(v)); + if (u == v) { + return u; + } + + auto stable = u; + auto removed = v; + if (dynamic_graph.adjacency[stable].size() < dynamic_graph.adjacency[removed].size()) { + std::swap(stable, removed); + } + sets.merge_to(stable, removed); + + // Stamp stable's neighbors so each removed-neighbor lookup is O(1). + for (const auto &entry : dynamic_graph.adjacency[stable]) { + dynamic_graph.scratch_edge_id[entry.neighbor] = entry.edge_id; + } + + // The contracted edge (stable, removed) must be in stable's adjacency. + const auto contracted_edge_id = dynamic_graph.scratch_edge_id[removed]; + heap.erase(contracted_edge_id); + dynamic_graph.scratch_edge_id[removed] = no_edge; + internal::erase_by_neighbor(dynamic_graph.adjacency[stable], removed); + + // Snapshot removed's neighbors before mutating its adjacency. + const auto removed_neighbors = dynamic_graph.adjacency[removed]; + + for (const auto &entry : removed_neighbors) { + const auto neighbor = entry.neighbor; + const auto removed_edge_id = entry.edge_id; + if (neighbor == stable) { + continue; + } + + const auto existing_id = dynamic_graph.scratch_edge_id[neighbor]; + if (existing_id == no_edge) { + // No existing stable-side edge to `neighbor`. Re-key the + // removed-side edge by replacing `removed` with `stable` on both + // sides. The edge_id, weight and constraint flag are preserved, + // so the heap entry remains valid without any priority change. + dynamic_graph.adjacency[stable].push_back({neighbor, removed_edge_id}); + dynamic_graph.scratch_edge_id[neighbor] = removed_edge_id; + internal::rename_neighbor(dynamic_graph.adjacency[neighbor], removed, stable); + auto &e = dynamic_graph.edges[removed_edge_id]; + if (e.u == removed) { + e.u = stable; + } else { + e.v = stable; + } + } else { + // Stable already has an edge to `neighbor`; fold removed's weight + // (and constraint flag) into it. Then drop the removed-side edge. + auto &keep = dynamic_graph.edges[existing_id]; + const auto &fold = dynamic_graph.edges[removed_edge_id]; + keep.weight += fold.weight; + const bool propagated_constraint = + keep.is_constraint != 0 || fold.is_constraint != 0; + keep.is_constraint = propagated_constraint ? 1 : 0; + + heap.erase(removed_edge_id); + internal::erase_by_neighbor(dynamic_graph.adjacency[neighbor], removed); + + if (propagated_constraint) { + heap.erase(existing_id); + } else { + heap.change(existing_id, priority_for(keep.weight, absolute_priority)); + } + } + } + + // Clear scratch via the updated stable adjacency (includes appended entries). + for (const auto &entry : dynamic_graph.adjacency[stable]) { + dynamic_graph.scratch_edge_id[entry.neighbor] = no_edge; + } + + dynamic_graph.adjacency[removed].clear(); + dynamic_graph.alive[removed] = false; + --dynamic_graph.alive_count; + return stable; +} + +} // namespace bioimage_cpp::graph::multicut::detail diff --git a/include/bioimage_cpp/graph/multicut/greedy_additive.hxx b/include/bioimage_cpp/graph/multicut/greedy_additive.hxx new file mode 100644 index 0000000..82f1b4f --- /dev/null +++ b/include/bioimage_cpp/graph/multicut/greedy_additive.hxx @@ -0,0 +1,80 @@ +#pragma once + +#include "bioimage_cpp/graph/multicut/detail.hxx" +#include "bioimage_cpp/graph/multicut/objective.hxx" + +#include +#include +#include + +namespace bioimage_cpp::graph::multicut { + +inline std::vector greedy_additive( + const UndirectedGraph &graph, + const std::vector &costs, + const double weight_stop, + const double node_num_stop, + const bool add_noise, + const int seed, + const double sigma +) { + validate_costs(graph, costs); + detail::DynamicGraph dynamic_graph(graph); + bioimage_cpp::detail::UnionFind sets(static_cast(graph.number_of_nodes())); + detail::EdgeHeap heap; + detail::initialize_dynamic_graph(graph, costs, dynamic_graph, heap, false, add_noise, seed, sigma); + + while (!heap.empty() && dynamic_graph.alive_count > 1) { + const auto top = heap.top(); + if (top.priority <= weight_stop) { + break; + } + if (node_num_stop > 0.0 + && dynamic_graph.alive_count <= detail::stop_node_count(graph, node_num_stop)) { + break; + } + const auto &edge = dynamic_graph.edges[top.key]; + detail::merge_dynamic_nodes(dynamic_graph, sets, heap, edge.u, edge.v, false); + } + return detail::labels_from_sets(sets, graph); +} + +class GreedyAdditiveSolver final : public SolverBase { +public: + GreedyAdditiveSolver( + const double weight_stop = 0.0, + const double node_num_stop = -1.0, + const bool add_noise = false, + const int seed = 42, + const double sigma = 1.0 + ) + : weight_stop_(weight_stop), + node_num_stop_(node_num_stop), + add_noise_(add_noise), + seed_(seed), + sigma_(sigma) { + } + + std::vector optimize(Objective &objective) const override { + auto labels = greedy_additive( + objective.graph(), + objective.costs(), + weight_stop_, + node_num_stop_, + add_noise_, + seed_, + sigma_ + ); + objective.set_labels(labels); + return labels; + } + +private: + double weight_stop_; + double node_num_stop_; + bool add_noise_; + int seed_; + double sigma_; +}; + +} // namespace bioimage_cpp::graph::multicut diff --git a/include/bioimage_cpp/graph/multicut/greedy_fixation.hxx b/include/bioimage_cpp/graph/multicut/greedy_fixation.hxx new file mode 100644 index 0000000..8795276 --- /dev/null +++ b/include/bioimage_cpp/graph/multicut/greedy_fixation.hxx @@ -0,0 +1,66 @@ +#pragma once + +#include "bioimage_cpp/graph/multicut/detail.hxx" +#include "bioimage_cpp/graph/multicut/objective.hxx" + +#include +#include +#include + +namespace bioimage_cpp::graph::multicut { + +inline std::vector greedy_fixation( + const UndirectedGraph &graph, + const std::vector &costs, + const double weight_stop, + const double node_num_stop +) { + validate_costs(graph, costs); + detail::DynamicGraph dynamic_graph(graph); + bioimage_cpp::detail::UnionFind sets(static_cast(graph.number_of_nodes())); + detail::EdgeHeap heap; + detail::initialize_dynamic_graph(graph, costs, dynamic_graph, heap, true); + + while (!heap.empty() && dynamic_graph.alive_count > 1) { + const auto top = heap.top(); + // Priority is |weight|, so this also handles weight == 0 (stop). + if (top.priority <= weight_stop) { + break; + } + if (node_num_stop > 0.0 + && dynamic_graph.alive_count <= detail::stop_node_count(graph, node_num_stop)) { + break; + } + const auto edge_id = top.key; + const auto &edge = dynamic_graph.edges[edge_id]; + if (edge.weight > 0.0) { + detail::merge_dynamic_nodes(dynamic_graph, sets, heap, edge.u, edge.v, true); + } else { + // weight < 0: forbid merging through this edge. merge_dynamic_nodes + // propagates the flag onto any merged successor edges. + heap.pop(); + dynamic_graph.edges[edge_id].is_constraint = 1; + } + } + return detail::labels_from_sets(sets, graph); +} + +class GreedyFixationSolver final : public SolverBase { +public: + GreedyFixationSolver(const double weight_stop = 0.0, const double node_num_stop = -1.0) + : weight_stop_(weight_stop), + node_num_stop_(node_num_stop) { + } + + std::vector optimize(Objective &objective) const override { + auto labels = greedy_fixation(objective.graph(), objective.costs(), weight_stop_, node_num_stop_); + objective.set_labels(labels); + return labels; + } + +private: + double weight_stop_; + double node_num_stop_; +}; + +} // namespace bioimage_cpp::graph::multicut diff --git a/include/bioimage_cpp/graph/multicut/kernighan_lin.hxx b/include/bioimage_cpp/graph/multicut/kernighan_lin.hxx new file mode 100644 index 0000000..69ae68e --- /dev/null +++ b/include/bioimage_cpp/graph/multicut/kernighan_lin.hxx @@ -0,0 +1,399 @@ +#pragma once + +#include "bioimage_cpp/detail/edge_hash.hxx" +#include "bioimage_cpp/detail/indexed_heap.hxx" +#include "bioimage_cpp/detail/union_find.hxx" +#include "bioimage_cpp/graph/multicut/objective.hxx" + +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph::multicut { + +namespace detail_kl { + +struct ClusterPair { + std::uint64_t a; + std::uint64_t b; + double cut_weight; +}; + +inline std::vector compute_cluster_pairs( + const UndirectedGraph &graph, + const std::vector &costs, + const std::vector &labels +) { + using bioimage_cpp::detail::Edge; + using bioimage_cpp::detail::EdgeHash; + using bioimage_cpp::detail::edge_key; + + std::unordered_map table; + table.reserve(static_cast(graph.number_of_edges())); + for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) { + const auto uv = graph.uv(edge); + const auto la = labels[static_cast(uv.first)]; + const auto lb = labels[static_cast(uv.second)]; + if (la == lb) { + continue; + } + table[edge_key(la, lb)] += costs[static_cast(edge)]; + } + + std::vector result; + result.reserve(table.size()); + for (const auto &entry : table) { + result.push_back({entry.first.first, entry.first.second, entry.second}); + } + std::sort(result.begin(), result.end(), [](const ClusterPair &lhs, const ClusterPair &rhs) { + return std::tie(lhs.a, lhs.b) < std::tie(rhs.a, rhs.b); + }); + return result; +} + +inline std::vector> build_cluster_to_nodes( + const std::vector &labels, + const std::uint64_t number_of_clusters +) { + std::vector> result(static_cast(number_of_clusters)); + for (std::uint64_t node = 0; node < labels.size(); ++node) { + result[static_cast(labels[static_cast(node)])].push_back(node); + } + return result; +} + +// Scratch flag used only during the initial gain scan, so the chain can tell +// whether each neighbor belongs to the current bipartition. After init the +// addressable heap is the source of truth (heap.contains(u) replaces the +// stale-tracking arrays the old hand-rolled heap required). +struct ChainBuffers { + std::vector in_pair; + + explicit ChainBuffers(const std::size_t n_nodes) : in_pair(n_nodes, 0) {} +}; + +struct ChainScratch { + std::vector queue_nodes; + bioimage_cpp::detail::DenseIndexedHeap heap; + + explicit ChainScratch(const std::size_t n_nodes) : heap(n_nodes) {} +}; + +// Run a Kernighan-Lin move-chain on the bipartition (cluster_a, cluster_b). +// +// Mutates `labels` and `cluster_to_nodes` if the chain commits any moves. +// `cluster_to_nodes[c]` is treated as an append-only list of "nodes ever in +// cluster c during this outer iteration"; the filter on `labels[v] == c` +// removes stale entries on the fly. Returns the committed cumulative gain +// (>= 0). +inline double run_chain( + const UndirectedGraph &graph, + const std::vector &costs, + std::vector &labels, + std::vector> &cluster_to_nodes, + ChainBuffers &bufs, + ChainScratch &scratch, + const std::uint64_t cluster_a, + const std::uint64_t cluster_b, + const double epsilon +) { + if (cluster_a == cluster_b) { + return 0.0; + } + + auto &queue_nodes = scratch.queue_nodes; + auto &heap = scratch.heap; + queue_nodes.clear(); + heap.clear(); + + const auto &stale_a = cluster_to_nodes[static_cast(cluster_a)]; + const auto &stale_b = cluster_to_nodes[static_cast(cluster_b)]; + queue_nodes.reserve(stale_a.size() + stale_b.size()); + + std::size_t live_a = 0; + std::size_t live_b = 0; + for (const auto v : stale_a) { + if (labels[static_cast(v)] == cluster_a) { + queue_nodes.push_back(v); + bufs.in_pair[static_cast(v)] = 1; + ++live_a; + } + } + for (const auto v : stale_b) { + if (labels[static_cast(v)] == cluster_b) { + queue_nodes.push_back(v); + bufs.in_pair[static_cast(v)] = 1; + ++live_b; + } + } + // No two-side swap is possible if either side is empty or both sides have + // a single node (the only non-trivial outcome is a join, handled by + // apply_joins). + if (live_a == 0 || live_b == 0 || (live_a == 1 && live_b == 1)) { + for (const auto v : queue_nodes) { + bufs.in_pair[static_cast(v)] = 0; + } + return 0.0; + } + + for (const auto v : queue_nodes) { + double w_to_a = 0.0; + double w_to_b = 0.0; + for (const auto adj : graph.node_adjacency(v)) { + if (!bufs.in_pair[static_cast(adj.node)]) { + continue; + } + const auto c = costs[static_cast(adj.edge)]; + if (labels[static_cast(adj.node)] == cluster_a) { + w_to_a += c; + } else { + w_to_b += c; + } + } + const auto v_label = labels[static_cast(v)]; + const double gain_v = (v_label == cluster_a) ? (w_to_b - w_to_a) : (w_to_a - w_to_b); + heap.push(static_cast(v), gain_v); + } + + struct Move { + std::uint64_t node; + std::uint64_t new_label; + }; + std::vector chain; + chain.reserve(queue_nodes.size()); + + double cumulative = 0.0; + double best_cumulative = 0.0; + std::size_t best_prefix = 0; + // The chain keeps running through negative moves because a later prefix + // can still recover. In practice deep negative runs almost never improve + // best_cumulative, so cap the lookahead at a small constant after each new + // best. + constexpr std::size_t max_steps_without_improvement = 32; + std::size_t steps_since_best = 0; + + while (!heap.empty()) { + const auto top = heap.pop(); + const auto v = static_cast(top.key); + const auto gain_v = top.priority; + const auto old_label = labels[static_cast(v)]; + const auto new_label = (old_label == cluster_a) ? cluster_b : cluster_a; + + cumulative += gain_v; + chain.push_back({v, new_label}); + + if (cumulative > best_cumulative + epsilon) { + best_cumulative = cumulative; + best_prefix = chain.size(); + steps_since_best = 0; + } else { + ++steps_since_best; + if (steps_since_best > max_steps_without_improvement) { + break; + } + } + + for (const auto adj : graph.node_adjacency(v)) { + const auto u_key = static_cast(adj.node); + if (!heap.contains(u_key)) { + continue; + } + // u is in the heap, so it's an unmoved member of the current + // bipartition. labels[u] is therefore u's actual side (chain + // moves are tentative — they only mutate `labels` on commit). + const auto c = costs[static_cast(adj.edge)]; + const auto u_label = labels[static_cast(adj.node)]; + const double delta = (u_label == old_label) ? 2.0 * c : -2.0 * c; + heap.change(u_key, heap.priority_of(u_key) + delta); + } + } + + for (const auto v : queue_nodes) { + bufs.in_pair[static_cast(v)] = 0; + } + + if (best_cumulative > epsilon) { + for (std::size_t i = 0; i < best_prefix; ++i) { + const auto v = chain[i].node; + const auto new_label = chain[i].new_label; + labels[static_cast(v)] = new_label; + cluster_to_nodes[static_cast(new_label)].push_back(v); + } + return best_cumulative; + } + return 0.0; +} + +// Single-node polish pass. For each node, find the adjacent cluster (if any) +// for which the unilateral move-out gain is positive, then move. The chain +// already considers these moves implicitly, but its greedy ordering can lock +// a node into a suboptimal commit prefix; this cheap O(|V| + |E|) sweep +// recovers those decisions without altering converged optima. +inline bool single_node_polish( + const UndirectedGraph &graph, + const std::vector &costs, + std::vector &labels, + const double epsilon +) { + bool improved = false; + // Hoist scratch across nodes; clear via the touched-keys list per node. + std::unordered_map sums; + std::vector touched; + for (std::uint64_t v = 0; v < graph.number_of_nodes(); ++v) { + sums.clear(); + touched.clear(); + for (const auto adj : graph.node_adjacency(v)) { + const auto label = labels[static_cast(adj.node)]; + auto it = sums.find(label); + if (it == sums.end()) { + sums.emplace(label, costs[static_cast(adj.edge)]); + touched.push_back(label); + } else { + it->second += costs[static_cast(adj.edge)]; + } + } + const auto cur = labels[static_cast(v)]; + const auto cur_it = sums.find(cur); + const auto cur_sum = (cur_it == sums.end()) ? 0.0 : cur_it->second; + + double best_gain = 0.0; + std::uint64_t best = cur; + for (const auto candidate : touched) { + if (candidate == cur) { + continue; + } + const auto gain = sums[candidate] - cur_sum; + if (gain > best_gain + epsilon) { + best_gain = gain; + best = candidate; + } + } + if (best != cur) { + labels[static_cast(v)] = best; + improved = true; + } + } + return improved; +} + +// Apply beneficial cluster joins in a single pass using union-find over cluster +// ids. Merging two clusters across cut_weight > 0 edges removes those edges +// from the multicut energy. +inline bool apply_joins( + std::vector &labels, + const std::vector &pairs, + const std::uint64_t number_of_clusters, + const double epsilon +) { + if (number_of_clusters == 0) { + return false; + } + bioimage_cpp::detail::UnionFind sets(static_cast(number_of_clusters)); + bool any_join = false; + for (const auto &pair : pairs) { + if (pair.cut_weight <= epsilon) { + continue; + } + const auto root_a = sets.find(pair.a); + const auto root_b = sets.find(pair.b); + if (root_a == root_b) { + continue; + } + sets.merge(root_a, root_b); + any_join = true; + } + if (any_join) { + for (auto &label : labels) { + label = sets.find(label); + } + } + return any_join; +} + +} // namespace detail_kl + +inline std::vector kernighan_lin( + const UndirectedGraph &graph, + const std::vector &costs, + std::vector labels, + const std::uint64_t number_of_outer_iterations, + const double epsilon +) { + validate_costs(graph, costs); + validate_labels(graph, labels); + labels = dense_relabel(labels); + + const auto n_nodes = static_cast(graph.number_of_nodes()); + detail_kl::ChainBuffers bufs(n_nodes); + detail_kl::ChainScratch scratch(n_nodes); + + for (std::uint64_t iteration = 0; iteration < number_of_outer_iterations; ++iteration) { + bool improved = false; + + const auto pairs_for_chain = detail_kl::compute_cluster_pairs(graph, costs, labels); + const auto number_of_clusters = labels.empty() + ? std::uint64_t{0} + : (*std::max_element(labels.begin(), labels.end()) + 1); + auto cluster_to_nodes = detail_kl::build_cluster_to_nodes(labels, number_of_clusters); + + for (const auto &pair : pairs_for_chain) { + const auto delta = detail_kl::run_chain( + graph, costs, labels, cluster_to_nodes, bufs, scratch, pair.a, pair.b, epsilon + ); + if (delta > epsilon) { + improved = true; + } + } + + const auto pairs_for_join = detail_kl::compute_cluster_pairs(graph, costs, labels); + const auto current_number_of_clusters = labels.empty() + ? std::uint64_t{0} + : (*std::max_element(labels.begin(), labels.end()) + 1); + if (detail_kl::apply_joins(labels, pairs_for_join, current_number_of_clusters, epsilon)) { + improved = true; + } + + if (detail_kl::single_node_polish(graph, costs, labels, epsilon)) { + improved = true; + } + + labels = dense_relabel(labels); + if (!improved) { + break; + } + } + return labels; +} + +class KernighanLinSolver final : public SolverBase { +public: + KernighanLinSolver( + const std::uint64_t number_of_outer_iterations = 100, + const double epsilon = 1.0e-6 + ) + : number_of_outer_iterations_(number_of_outer_iterations), + epsilon_(epsilon) { + } + + std::vector optimize(Objective &objective) const override { + auto labels = kernighan_lin( + objective.graph(), + objective.costs(), + objective.labels(), + number_of_outer_iterations_, + epsilon_ + ); + objective.set_labels(labels); + return labels; + } + +private: + std::uint64_t number_of_outer_iterations_; + double epsilon_; +}; + +} // namespace bioimage_cpp::graph::multicut diff --git a/include/bioimage_cpp/graph/multicut/objective.hxx b/include/bioimage_cpp/graph/multicut/objective.hxx new file mode 100644 index 0000000..6a67492 --- /dev/null +++ b/include/bioimage_cpp/graph/multicut/objective.hxx @@ -0,0 +1,114 @@ +#pragma once + +#include "bioimage_cpp/detail/relabel.hxx" +#include "bioimage_cpp/graph/undirected_graph.hxx" + +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp::graph::multicut { + +inline void validate_costs(const UndirectedGraph &graph, const std::vector &costs) { + if (costs.size() != static_cast(graph.number_of_edges())) { + throw std::invalid_argument("edge costs length must match graph number_of_edges"); + } +} + +inline void validate_labels(const UndirectedGraph &graph, const std::vector &labels) { + if (labels.size() != static_cast(graph.number_of_nodes())) { + throw std::invalid_argument("labels length must match graph number_of_nodes"); + } +} + +inline std::vector singleton_labels(const UndirectedGraph &graph) { + std::vector labels(static_cast(graph.number_of_nodes())); + std::iota(labels.begin(), labels.end(), std::uint64_t{0}); + return labels; +} + +using bioimage_cpp::detail::dense_relabel; + +inline double energy( + const UndirectedGraph &graph, + const std::vector &costs, + const std::vector &labels +) { + validate_costs(graph, costs); + validate_labels(graph, labels); + double result = 0.0; + for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) { + const auto uv = graph.uv(edge); + if (labels[static_cast(uv.first)] != labels[static_cast(uv.second)]) { + result += costs[static_cast(edge)]; + } + } + return result; +} + +class Objective { +public: + Objective(const UndirectedGraph &graph, std::vector costs) + : graph_(graph), + costs_(std::move(costs)), + labels_(singleton_labels(graph)) { + validate_costs(graph_, costs_); + } + + Objective( + const UndirectedGraph &graph, + std::vector costs, + std::vector labels + ) + : graph_(graph), + costs_(std::move(costs)), + labels_(std::move(labels)) { + validate_costs(graph_, costs_); + validate_labels(graph_, labels_); + } + + [[nodiscard]] const UndirectedGraph &graph() const { + return graph_; + } + + [[nodiscard]] const std::vector &costs() const { + return costs_; + } + + [[nodiscard]] const std::vector &labels() const { + return labels_; + } + + void set_labels(std::vector labels) { + validate_labels(graph_, labels); + labels_ = dense_relabel(labels); + } + + void reset_labels() { + labels_ = singleton_labels(graph_); + } + + [[nodiscard]] double eval() const { + return energy(graph_, costs_, labels_); + } + + [[nodiscard]] double eval(const std::vector &labels) const { + return energy(graph_, costs_, labels); + } + +private: + const UndirectedGraph &graph_; + std::vector costs_; + std::vector labels_; +}; + +class SolverBase { +public: + virtual ~SolverBase() = default; + virtual std::vector optimize(Objective &objective) const = 0; +}; + +} // namespace bioimage_cpp::graph::multicut diff --git a/include/bioimage_cpp/graph/region_adjacency_graph.hxx b/include/bioimage_cpp/graph/region_adjacency_graph.hxx index 7248376..58e9f1c 100644 --- a/include/bioimage_cpp/graph/region_adjacency_graph.hxx +++ b/include/bioimage_cpp/graph/region_adjacency_graph.hxx @@ -1,6 +1,8 @@ #pragma once #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/edge_hash.hxx" +#include "bioimage_cpp/detail/threading.hxx" #include "bioimage_cpp/graph/undirected_graph.hxx" #include @@ -10,7 +12,6 @@ #include #include #include -#include #include #include #include @@ -40,22 +41,10 @@ private: namespace detail { -using Edge = UndirectedGraph::Edge; - -struct EdgeHash { - std::size_t operator()(const Edge &edge) const { - const auto first = static_cast(edge.first); - const auto second = static_cast(edge.second); - return first ^ (second + 0x9e3779b97f4a7c15ULL + (first << 6U) + (first >> 2U)); - } -}; - -inline Edge edge_key(std::uint64_t u, std::uint64_t v) { - if (v < u) { - std::swap(u, v); - } - return {u, v}; -} +using bioimage_cpp::detail::Edge; +using bioimage_cpp::detail::EdgeHash; +using bioimage_cpp::detail::edge_key; +using bioimage_cpp::detail::normalize_thread_count; template std::uint64_t checked_label_to_node(const T value) { @@ -82,23 +71,6 @@ std::uint64_t max_label(const ConstArrayView &labels) { return max_value; } -inline std::size_t normalize_thread_count( - const std::size_t requested, - const std::size_t number_of_work_items -) { - if (number_of_work_items == 0) { - return 1; - } - std::size_t n_threads = requested; - if (n_threads == 0) { - n_threads = std::thread::hardware_concurrency(); - if (n_threads == 0) { - n_threads = 1; - } - } - return std::max(1, std::min(n_threads, number_of_work_items)); -} - template void add_edge_if_different( const T *data, @@ -213,41 +185,33 @@ RegionAdjacencyGraph region_adjacency_graph( const auto work_items = static_cast(labels.shape[0]); const auto n_threads = detail::normalize_thread_count(number_of_threads, work_items); std::vector> per_thread_edges(n_threads); - std::vector threads; - threads.reserve(n_threads > 0 ? n_threads - 1 : 0); - const auto run_chunk = [&](const std::size_t thread_id) { - const auto begin = thread_id * work_items / n_threads; - const auto end = (thread_id + 1) * work_items / n_threads; - if (labels.ndim() == 2) { - detail::scan_2d_chunk( - labels.data, - static_cast(labels.shape[0]), - static_cast(labels.shape[1]), - begin, - end, - per_thread_edges[thread_id] - ); - } else { - detail::scan_3d_chunk( - labels.data, - static_cast(labels.shape[0]), - static_cast(labels.shape[1]), - static_cast(labels.shape[2]), - begin, - end, - per_thread_edges[thread_id] - ); + bioimage_cpp::detail::parallel_for_chunks( + n_threads, + work_items, + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + if (labels.ndim() == 2) { + detail::scan_2d_chunk( + labels.data, + static_cast(labels.shape[0]), + static_cast(labels.shape[1]), + begin, + end, + per_thread_edges[thread_id] + ); + } else { + detail::scan_3d_chunk( + labels.data, + static_cast(labels.shape[0]), + static_cast(labels.shape[1]), + static_cast(labels.shape[2]), + begin, + end, + per_thread_edges[thread_id] + ); + } } - }; - - for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) { - threads.emplace_back(run_chunk, thread_id); - } - run_chunk(0); - for (auto &thread : threads) { - thread.join(); - } + ); auto graph = RegionAdjacencyGraph(max_node + 1, std::move(shape)); for (const auto edge : detail::merge_edge_sets(per_thread_edges)) { diff --git a/include/bioimage_cpp/graph/undirected_graph.hxx b/include/bioimage_cpp/graph/undirected_graph.hxx index a949965..dc34605 100644 --- a/include/bioimage_cpp/graph/undirected_graph.hxx +++ b/include/bioimage_cpp/graph/undirected_graph.hxx @@ -1,5 +1,7 @@ #pragma once +#include "bioimage_cpp/detail/edge_hash.hxx" + #include #include #include @@ -21,7 +23,7 @@ class UndirectedGraph { public: using NodeId = std::uint64_t; using EdgeId = std::uint64_t; - using Edge = std::pair; + using Edge = detail::Edge; using AdjacencyList = std::vector; explicit UndirectedGraph( @@ -108,7 +110,7 @@ public: throw std::invalid_argument("self edges are not supported"); } - const auto key = edge_key(u, v); + const auto key = detail::edge_key(u, v); const auto found = edge_lookup_.find(key); if (found != edge_lookup_.end()) { return found->second; @@ -123,7 +125,7 @@ public: return -1; } - const auto key = edge_key(u, v); + const auto key = detail::edge_key(u, v); const auto found = edge_lookup_.find(key); if (found == edge_lookup_.end()) { return -1; @@ -193,25 +195,10 @@ protected: } private: - struct EdgeHash { - std::size_t operator()(const Edge &edge) const { - const auto first = static_cast(edge.first); - const auto second = static_cast(edge.second); - return first ^ (second + 0x9e3779b97f4a7c15ULL + (first << 6U) + (first >> 2U)); - } - }; - - static Edge edge_key(NodeId u, NodeId v) { - if (v < u) { - std::swap(u, v); - } - return {u, v}; - } - NodeId number_of_nodes_; std::vector edges_; std::vector adjacency_; - std::unordered_map edge_lookup_; + std::unordered_map edge_lookup_; }; } // namespace bioimage_cpp::graph diff --git a/include/bioimage_cpp/segmentation/mutex_watershed.hxx b/include/bioimage_cpp/segmentation/mutex_watershed.hxx index 47cf445..263779e 100644 --- a/include/bioimage_cpp/segmentation/mutex_watershed.hxx +++ b/include/bioimage_cpp/segmentation/mutex_watershed.hxx @@ -1,6 +1,8 @@ #pragma once #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/detail/union_find.hxx" #include #include @@ -13,35 +15,6 @@ namespace bioimage_cpp { -class DisjointSets { -public: - explicit DisjointSets(const std::size_t size) : parents_(size), ranks_(size, 0) { - std::iota(parents_.begin(), parents_.end(), std::uint64_t{0}); - } - - std::uint64_t find(const std::uint64_t node) { - if (parents_[node] != node) { - parents_[node] = find(parents_[node]); - } - return parents_[node]; - } - - std::uint64_t unite_roots(std::uint64_t first, std::uint64_t second) { - if (ranks_[first] < ranks_[second]) { - std::swap(first, second); - } - parents_[second] = first; - if (ranks_[first] == ranks_[second]) { - ++ranks_[first]; - } - return first; - } - -private: - std::vector parents_; - std::vector ranks_; -}; - using MutexStorage = std::vector>; template @@ -93,15 +66,6 @@ inline void merge_mutexes( mutexes_from.clear(); } -inline std::vector c_order_strides(const std::vector &shape) { - std::vector strides(shape.size(), 1); - for (std::ptrdiff_t axis = static_cast(shape.size()) - 2; axis >= 0; --axis) { - strides[static_cast(axis)] = - strides[static_cast(axis + 1)] * shape[static_cast(axis + 1)]; - } - return strides; -} - template void mutex_watershed_grid( const ConstArrayView &affinities, @@ -157,7 +121,7 @@ void mutex_watershed_grid( std::ptrdiff_t{1}, [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; } )); - const auto spatial_strides = c_order_strides(spatial_shape); + const auto spatial_strides = detail::c_order_strides(spatial_shape); std::vector offset_strides(number_of_channels, 0); for (std::size_t channel = 0; channel < number_of_channels; ++channel) { @@ -181,7 +145,7 @@ void mutex_watershed_grid( return first.weight > second.weight; }); - DisjointSets sets(static_cast(number_of_nodes)); + detail::UnionFind sets(static_cast(number_of_nodes)); MutexStorage mutexes(static_cast(number_of_nodes)); for (const auto &edge : edge_order) { diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx index b26436f..b4979db 100644 --- a/src/bindings/graph.cxx +++ b/src/bindings/graph.cxx @@ -1,7 +1,9 @@ #include "graph.hxx" #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/graph/connected_components.hxx" #include "bioimage_cpp/graph/feature_accumulation.hxx" +#include "bioimage_cpp/graph/multicut.hxx" #include "bioimage_cpp/graph/region_adjacency_graph.hxx" #include "bioimage_cpp/graph/undirected_graph.hxx" @@ -25,6 +27,7 @@ namespace { using Graph = graph::UndirectedGraph; using Rag = graph::RegionAdjacencyGraph; using UInt64Array = nb::ndarray; +using ConstUInt8Array = nb::ndarray; using ConstUInt64Array = nb::ndarray; using Int64Array = nb::ndarray; using DoubleArray = nb::ndarray; @@ -71,6 +74,46 @@ DoubleArray make_double_array(const std::vector &shape) { return DoubleArray(data, shape.size(), shape.data(), owner); } +UInt64Array vector_to_uint64_array(const std::vector &values) { + auto result = make_uint64_array({values.size()}); + std::copy(values.begin(), values.end(), result.data()); + return result; +} + +std::vector double_array_to_vector( + ConstDoubleArray array, + const char *argument_name, + const std::uint64_t expected_size +) { + if (array.ndim() != 1) { + throw std::invalid_argument(std::string(argument_name) + " must be a 1D float64 array"); + } + if (array.shape(0) != static_cast(expected_size)) { + throw std::invalid_argument( + std::string(argument_name) + " length must match expected size" + ); + } + const auto *data = array.data(); + return std::vector(data, data + array.shape(0)); +} + +std::vector uint64_array_to_vector( + ConstUInt64Array array, + const char *argument_name, + const std::uint64_t expected_size +) { + if (array.ndim() != 1) { + throw std::invalid_argument(std::string(argument_name) + " must be a 1D uint64 array"); + } + if (array.shape(0) != static_cast(expected_size)) { + throw std::invalid_argument( + std::string(argument_name) + " length must match expected size" + ); + } + const auto *data = array.data(); + return std::vector(data, data + array.shape(0)); +} + UInt64Array graph_nodes(const Graph &graph) { auto result = make_uint64_array({static_cast(graph.number_of_nodes())}); auto *data = result.data(); @@ -217,6 +260,100 @@ UInt64Array graph_edges_from_node_list(const Graph &graph, ConstUInt64Array node return extracted.first; } +UInt64Array graph_connected_components(const Graph &graph) { + std::vector labels; + { + nb::gil_scoped_release release; + labels = graph::connected_components(graph); + } + return vector_to_uint64_array(labels); +} + +UInt64Array graph_connected_components_masked(const Graph &graph, ConstUInt8Array edge_mask) { + if (edge_mask.ndim() != 1) { + throw std::invalid_argument("edge_mask must be a 1D uint8 array"); + } + if (edge_mask.shape(0) != static_cast(graph.number_of_edges())) { + throw std::invalid_argument("edge_mask length must match graph number_of_edges"); + } + std::vector labels; + { + nb::gil_scoped_release release; + labels = graph::connected_components(graph, edge_mask.data()); + } + return vector_to_uint64_array(labels); +} + +double multicut_energy(const Graph &graph, ConstDoubleArray costs, ConstUInt64Array labels) { + const auto cost_vector = double_array_to_vector(costs, "edge_costs", graph.number_of_edges()); + const auto label_vector = uint64_array_to_vector(labels, "labels", graph.number_of_nodes()); + nb::gil_scoped_release release; + return graph::multicut::energy(graph, cost_vector, label_vector); +} + +UInt64Array multicut_greedy_additive( + const Graph &graph, + ConstDoubleArray costs, + const double weight_stop, + const double node_num_stop, + const bool add_noise, + const int seed, + const double sigma +) { + const auto cost_vector = double_array_to_vector(costs, "edge_costs", graph.number_of_edges()); + std::vector labels; + { + nb::gil_scoped_release release; + labels = graph::multicut::greedy_additive( + graph, + cost_vector, + weight_stop, + node_num_stop, + add_noise, + seed, + sigma + ); + } + return vector_to_uint64_array(labels); +} + +UInt64Array multicut_greedy_fixation( + const Graph &graph, + ConstDoubleArray costs, + const double weight_stop, + const double node_num_stop +) { + const auto cost_vector = double_array_to_vector(costs, "edge_costs", graph.number_of_edges()); + std::vector labels; + { + nb::gil_scoped_release release; + labels = graph::multicut::greedy_fixation(graph, cost_vector, weight_stop, node_num_stop); + } + return vector_to_uint64_array(labels); +} + +UInt64Array multicut_kernighan_lin( + const Graph &graph, + ConstDoubleArray costs, + ConstUInt64Array initial_labels, + const std::uint64_t number_of_outer_iterations, + const double epsilon +) { + const auto cost_vector = double_array_to_vector(costs, "edge_costs", graph.number_of_edges()); + auto label_vector = uint64_array_to_vector(initial_labels, "initial_labels", graph.number_of_nodes()); + { + nb::gil_scoped_release release; + label_vector = graph::multicut::kernighan_lin( + graph, + cost_vector, + std::move(label_vector), + number_of_outer_iterations, + epsilon + ); + } + return vector_to_uint64_array(label_vector); +} + template Rag region_adjacency_graph_t( LabelArray labels, @@ -416,6 +553,49 @@ void bind_graph(nb::module_ &m) { nb::class_(m, "RegionAdjacencyGraph") .def_prop_ro("shape", &Rag::shape); + m.def("_connected_components", &graph_connected_components, nb::arg("graph")); + m.def( + "_connected_components_masked", + &graph_connected_components_masked, + nb::arg("graph"), + nb::arg("edge_mask") + ); + m.def( + "_multicut_energy", + &multicut_energy, + nb::arg("graph"), + nb::arg("edge_costs"), + nb::arg("labels") + ); + m.def( + "_multicut_greedy_additive", + &multicut_greedy_additive, + nb::arg("graph"), + nb::arg("edge_costs"), + nb::arg("weight_stop"), + nb::arg("node_num_stop"), + nb::arg("add_noise"), + nb::arg("seed"), + nb::arg("sigma") + ); + m.def( + "_multicut_greedy_fixation", + &multicut_greedy_fixation, + nb::arg("graph"), + nb::arg("edge_costs"), + nb::arg("weight_stop"), + nb::arg("node_num_stop") + ); + m.def( + "_multicut_kernighan_lin", + &multicut_kernighan_lin, + nb::arg("graph"), + nb::arg("edge_costs"), + nb::arg("initial_labels"), + nb::arg("number_of_outer_iterations"), + nb::arg("epsilon") + ); + m.def( "_region_adjacency_graph_uint32", ®ion_adjacency_graph_t, diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index a3fd91d..98d01b9 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -2,9 +2,18 @@ from __future__ import annotations +from abc import ABC, abstractmethod + import numpy as np from .. import _core +from ._external import ( + DEFAULT_EXTERNAL_MULTICUT_PROBLEM_PATH, + EXTERNAL_MULTICUT_PROBLEM_URL, + external_multicut_problem_path, + load_external_multicut_problem, + load_external_multicut_problem_data, +) _REGION_ADJACENCY_GRAPH_BY_DTYPE = { np.dtype("uint32"): _core._region_adjacency_graph_uint32, @@ -118,6 +127,47 @@ def _as_serialization_array(serialization) -> np.ndarray: return np.ascontiguousarray(array) +def _copy_graph(graph: UndirectedGraph | RegionAdjacencyGraph) -> UndirectedGraph: + copied = UndirectedGraph(int(graph.number_of_nodes), int(graph.number_of_edges)) + if graph.number_of_edges: + copied.insert_edges(graph.uv_ids()) + return copied + + +def _as_edge_costs(edge_costs, graph: UndirectedGraph | RegionAdjacencyGraph) -> np.ndarray: + array = np.asarray(edge_costs, dtype=np.float64) + if array.ndim != 1: + raise ValueError("edge_costs must be a 1D array") + if array.shape[0] != graph.number_of_edges: + raise ValueError("edge_costs length must match graph number_of_edges") + return np.ascontiguousarray(array) + + +def _as_node_labels(labels, graph: UndirectedGraph | RegionAdjacencyGraph) -> np.ndarray: + array = np.asarray(labels, dtype=np.uint64) + if array.ndim != 1: + raise ValueError("labels must be a 1D array") + if array.shape[0] != graph.number_of_nodes: + raise ValueError("labels length must match graph number_of_nodes") + return np.ascontiguousarray(array) + + +def _dense_labels(labels) -> np.ndarray: + labels = np.asarray(labels, dtype=np.uint64) + _, dense = np.unique(labels, return_inverse=True) + return np.ascontiguousarray(dense.astype(np.uint64, copy=False)) + + +def _subproblem_from_edges(number_of_nodes: int, nodes, uvs, edge_costs): + local_ids = np.full(int(number_of_nodes), -1, dtype=np.int64) + local_ids[nodes] = np.arange(nodes.size, dtype=np.int64) + local_uvs = local_ids[np.asarray(uvs, dtype=np.uint64)] + sub_graph = UndirectedGraph(int(nodes.size), int(len(edge_costs))) + if local_uvs.size: + sub_graph.insert_edges(np.ascontiguousarray(local_uvs.astype(np.uint64, copy=False))) + return sub_graph, np.ascontiguousarray(np.asarray(edge_costs, dtype=np.float64)) + + def undirected_graph(number_of_nodes: int) -> UndirectedGraph: """Create an :class:`UndirectedGraph`.""" return UndirectedGraph(number_of_nodes) @@ -126,6 +176,240 @@ def undirected_graph(number_of_nodes: int) -> UndirectedGraph: RegionAdjacencyGraph = _core.RegionAdjacencyGraph +def connected_components( + graph: UndirectedGraph | RegionAdjacencyGraph, + edge_mask: np.ndarray | None = None, +) -> np.ndarray: + """Compute dense connected-component labels for graph nodes. + + If ``edge_mask`` is given, only edges with a true mask value contribute to + the connected components. + """ + if edge_mask is None: + return _core._connected_components(graph) + + mask = np.asarray(edge_mask) + if mask.dtype != np.dtype("bool"): + raise TypeError(f"edge_mask must have dtype bool, got dtype={mask.dtype}") + if mask.ndim != 1: + raise ValueError("edge_mask must be a 1D array") + if mask.shape[0] != graph.number_of_edges: + raise ValueError("edge_mask length must match graph number_of_edges") + return _core._connected_components_masked( + graph, np.ascontiguousarray(mask.astype(np.uint8, copy=False)) + ) + + +class MulticutObjective: + """Multicut objective for an undirected graph and edge costs.""" + + def __init__( + self, + graph: UndirectedGraph | RegionAdjacencyGraph, + edge_costs, + initial_labels=None, + ): + self._graph = _copy_graph(graph) + self._edge_costs = _as_edge_costs(edge_costs, self._graph) + if initial_labels is None: + self._labels = np.arange(self._graph.number_of_nodes, dtype=np.uint64) + else: + self._labels = _as_node_labels(initial_labels, self._graph) + + @property + def graph(self) -> UndirectedGraph: + return self._graph + + @property + def edge_costs(self) -> np.ndarray: + return self._edge_costs + + @property + def labels(self) -> np.ndarray: + return self._labels + + @labels.setter + def labels(self, labels) -> None: + self._labels = _as_node_labels(labels, self._graph) + + def set_labels(self, labels) -> None: + self.labels = labels + + def reset_labels(self) -> None: + self._labels = np.arange(self._graph.number_of_nodes, dtype=np.uint64) + + def energy(self, labels=None) -> float: + label_array = self._labels if labels is None else _as_node_labels(labels, self._graph) + return float(_core._multicut_energy(self._graph, self._edge_costs, label_array)) + + +class MulticutSolver(ABC): + """Base class for multicut solvers.""" + + @abstractmethod + def optimize(self, objective: MulticutObjective) -> np.ndarray: + """Optimize ``objective`` and return the node labeling.""" + + +class GreedyAdditiveMulticut(MulticutSolver): + def __init__( + self, + *, + weight_stop: float = 0.0, + node_num_stop: float = -1.0, + add_noise: bool = False, + seed: int = 42, + sigma: float = 1.0, + ): + self.weight_stop = float(weight_stop) + self.node_num_stop = float(node_num_stop) + self.add_noise = bool(add_noise) + self.seed = int(seed) + self.sigma = float(sigma) + + def optimize(self, objective: MulticutObjective) -> np.ndarray: + labels = _core._multicut_greedy_additive( + objective.graph, + objective.edge_costs, + self.weight_stop, + self.node_num_stop, + self.add_noise, + self.seed, + self.sigma, + ) + objective.labels = labels + return objective.labels + + +class GreedyFixationMulticut(MulticutSolver): + def __init__(self, *, weight_stop: float = 0.0, node_num_stop: float = -1.0): + self.weight_stop = float(weight_stop) + self.node_num_stop = float(node_num_stop) + + def optimize(self, objective: MulticutObjective) -> np.ndarray: + labels = _core._multicut_greedy_fixation( + objective.graph, + objective.edge_costs, + self.weight_stop, + self.node_num_stop, + ) + objective.labels = labels + return objective.labels + + +class KernighanLinMulticut(MulticutSolver): + def __init__( + self, + *, + number_of_outer_iterations: int = 100, + number_of_inner_iterations: int | None = None, + epsilon: float = 1.0e-6, + ): + self.number_of_outer_iterations = int(number_of_outer_iterations) + if self.number_of_outer_iterations < 0: + raise ValueError("number_of_outer_iterations must be non-negative") + self.number_of_inner_iterations = ( + None if number_of_inner_iterations is None else int(number_of_inner_iterations) + ) + self.epsilon = float(epsilon) + + def optimize(self, objective: MulticutObjective) -> np.ndarray: + initial_labels = objective.labels + if np.array_equal( + initial_labels, + np.arange(objective.graph.number_of_nodes, dtype=np.uint64), + ): + initial_labels = _core._multicut_greedy_additive( + objective.graph, + objective.edge_costs, + 0.0, + -1.0, + False, + 42, + 1.0, + ) + labels = _core._multicut_kernighan_lin( + objective.graph, + objective.edge_costs, + initial_labels, + self.number_of_outer_iterations, + self.epsilon, + ) + objective.labels = labels + return objective.labels + + +class ChainedMulticutSolvers(MulticutSolver): + def __init__(self, solvers): + self.solvers = list(solvers) + if len(self.solvers) == 0: + raise ValueError("solvers must contain at least one solver") + if not all(isinstance(solver, MulticutSolver) for solver in self.solvers): + raise TypeError("all solvers must inherit from MulticutSolver") + + def optimize(self, objective: MulticutObjective) -> np.ndarray: + labels = objective.labels + for solver in self.solvers: + labels = solver.optimize(objective) + return labels + + +class MulticutDecomposer(MulticutSolver): + def __init__( + self, + sub_solver: MulticutSolver, + *, + fallthrough_solver: MulticutSolver | None = None, + number_of_threads: int = 0, + ): + if not isinstance(sub_solver, MulticutSolver): + raise TypeError("sub_solver must inherit from MulticutSolver") + if fallthrough_solver is not None and not isinstance(fallthrough_solver, MulticutSolver): + raise TypeError("fallthrough_solver must inherit from MulticutSolver") + self.sub_solver = sub_solver + self.fallthrough_solver = fallthrough_solver + self.number_of_threads = _normalize_number_of_threads(number_of_threads) + + def optimize(self, objective: MulticutObjective) -> np.ndarray: + if self.fallthrough_solver is None and isinstance(self.sub_solver, GreedyAdditiveMulticut): + return self.sub_solver.optimize(objective) + + component_labels = connected_components( + objective.graph, + edge_mask=objective.edge_costs > 0.0, + ) + number_of_components = int(component_labels.max()) + 1 if component_labels.size else 0 + if number_of_components <= 1: + solver = self.fallthrough_solver or self.sub_solver + return solver.optimize(objective) + + global_labels = np.empty(objective.graph.number_of_nodes, dtype=np.uint64) + label_offset = 0 + all_uvs = objective.graph.uv_ids() + for component in range(number_of_components): + nodes = np.flatnonzero(component_labels == component).astype(np.uint64) + if nodes.size == 1: + global_labels[int(nodes[0])] = label_offset + label_offset += 1 + continue + + edge_ids = objective.graph.edges_from_node_list(nodes) + sub_graph, sub_costs = _subproblem_from_edges( + objective.graph.number_of_nodes, + nodes, + all_uvs[edge_ids], + objective.edge_costs[edge_ids], + ) + sub_objective = MulticutObjective(sub_graph, sub_costs) + sub_labels = self.sub_solver.optimize(sub_objective) + sub_labels = _dense_labels(sub_labels) + global_labels[nodes] = sub_labels + label_offset + label_offset += int(sub_labels.max()) + 1 + + objective.labels = _dense_labels(global_labels) + return objective.labels + + def region_adjacency_graph( labels: np.ndarray, *, @@ -327,14 +611,27 @@ def _normalize_number_of_threads(number_of_threads: int) -> int: __all__ = [ + "ChainedMulticutSolvers", "COMPLEX_EDGE_FEATURE_NAMES", + "DEFAULT_EXTERNAL_MULTICUT_PROBLEM_PATH", + "EXTERNAL_MULTICUT_PROBLEM_URL", + "GreedyAdditiveMulticut", + "GreedyFixationMulticut", + "KernighanLinMulticut", + "MulticutDecomposer", + "MulticutObjective", + "MulticutSolver", "RegionAdjacencyGraph", "SIMPLE_EDGE_FEATURE_NAMES", "UndirectedGraph", "affinity_features", "affinity_features_complex", + "connected_components", "edge_map_features", "edge_map_features_complex", + "external_multicut_problem_path", + "load_external_multicut_problem", + "load_external_multicut_problem_data", "region_adjacency_graph", "undirected_graph", ] diff --git a/src/bioimage_cpp/graph/_external.py b/src/bioimage_cpp/graph/_external.py new file mode 100644 index 0000000..5401281 --- /dev/null +++ b/src/bioimage_cpp/graph/_external.py @@ -0,0 +1,90 @@ +"""External graph problems used for development and regression checks.""" + +from __future__ import annotations + +import os +import urllib.request +from pathlib import Path +from urllib.error import URLError + +import numpy as np + + +EXTERNAL_MULTICUT_PROBLEM_URL = "https://oc.embl.de/index.php/s/yVKwyQ8VoPXYkft/download" +DEFAULT_EXTERNAL_MULTICUT_PROBLEM_PATH = Path("/tmp/bioimage_cpp_external_multicut_problem.txt") + + +def external_multicut_problem_path( + path: str | os.PathLike | None = None, + *, + download: bool = True, + timeout: float = 30.0, +) -> Path: + """Return the local path for the external multicut problem. + + The path can be supplied explicitly, via + ``BIOIMAGE_CPP_EXTERNAL_MULTICUT_PATH``, or via the cache path + ``BIOIMAGE_CPP_EXTERNAL_MULTICUT_CACHE``. If no existing file is found and + ``download`` is true, the problem is downloaded into the cache path. + """ + explicit_path = path or os.environ.get("BIOIMAGE_CPP_EXTERNAL_MULTICUT_PATH") + if explicit_path is not None: + resolved = Path(explicit_path) + if not resolved.exists(): + raise FileNotFoundError(f"external multicut problem does not exist: {resolved}") + return resolved + + cache_path = Path( + os.environ.get( + "BIOIMAGE_CPP_EXTERNAL_MULTICUT_CACHE", + str(DEFAULT_EXTERNAL_MULTICUT_PROBLEM_PATH), + ) + ) + if cache_path.exists(): + return cache_path + if not download: + raise FileNotFoundError(f"external multicut problem does not exist: {cache_path}") + + cache_path.parent.mkdir(parents=True, exist_ok=True) + try: + with urllib.request.urlopen(EXTERNAL_MULTICUT_PROBLEM_URL, timeout=timeout) as response: + cache_path.write_bytes(response.read()) + except URLError as error: + raise RuntimeError(f"could not download external multicut problem: {error}") from error + return cache_path + + +def load_external_multicut_problem_data( + path: str | os.PathLike | None = None, + *, + download: bool = True, + timeout: float = 30.0, +) -> tuple[np.ndarray, np.ndarray]: + """Load external multicut problem edge ids and costs.""" + problem_path = external_multicut_problem_path( + path, + download=download, + timeout=timeout, + ) + problem = np.genfromtxt(problem_path) + uv_ids = np.ascontiguousarray(problem[:, :2].astype(np.uint64, copy=False)) + costs = np.ascontiguousarray(problem[:, -1].astype(np.float64, copy=False)) + return uv_ids, costs + + +def load_external_multicut_problem( + path: str | os.PathLike | None = None, + *, + download: bool = True, + timeout: float = 30.0, +): + """Load the external multicut problem as a bioimage-cpp graph and costs.""" + from . import UndirectedGraph + + uv_ids, costs = load_external_multicut_problem_data( + path, + download=download, + timeout=timeout, + ) + graph = UndirectedGraph.from_edges(int(uv_ids.max()) + 1, uv_ids) + return graph, costs diff --git a/tests/graph/multicut/__init__.py b/tests/graph/multicut/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/tests/graph/multicut/__init__.py @@ -0,0 +1 @@ + diff --git a/tests/graph/multicut/_helpers.py b/tests/graph/multicut/_helpers.py new file mode 100644 index 0000000..250f049 --- /dev/null +++ b/tests/graph/multicut/_helpers.py @@ -0,0 +1,16 @@ +import numpy as np + + +def same_partition(labels, expected): + labels = np.asarray(labels) + expected = np.asarray(expected) + assert labels.shape == expected.shape + np.testing.assert_array_equal( + labels[:, None] == labels[None, :], + expected[:, None] == expected[None, :], + ) + + +def edge_cut_labels(graph, labels): + uv_ids = graph.uv_ids() + return labels[uv_ids[:, 0]] != labels[uv_ids[:, 1]] diff --git a/tests/graph/multicut/conftest.py b/tests/graph/multicut/conftest.py new file mode 100644 index 0000000..3d0ffee --- /dev/null +++ b/tests/graph/multicut/conftest.py @@ -0,0 +1,55 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +@pytest.fixture +def external_toy_problem(): + graph = bic.graph.UndirectedGraph.from_edges( + 6, + [ + [0, 1], + [0, 3], + [1, 2], + [1, 4], + [2, 5], + [3, 4], + [4, 5], + ], + ) + costs = np.array([5, -20, 5, 5, -20, 5, 5], dtype=np.float64) + expected_cut_edges = np.array([False, True, False, False, True, False, False]) + return graph, costs, expected_cut_edges + + +@pytest.fixture +def chain_problem(): + graph = bic.graph.UndirectedGraph.from_edges(4, [[0, 1], [1, 2], [2, 3]]) + costs = np.array([1.0, 2.0, -5.0], dtype=np.float64) + return graph, costs + + +@pytest.fixture +def frustrated_triangle(): + graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2], [0, 2]]) + costs = np.array([2.0, 2.0, -5.0], dtype=np.float64) + return graph, costs + + +@pytest.fixture +def grid_problem(): + edges = [] + costs = [] + shape = (5, 5) + for y in range(shape[0]): + for x in range(shape[1]): + node = y * shape[1] + x + if x + 1 < shape[1]: + edges.append([node, node + 1]) + costs.append(1.5 if x != 2 else -4.0) + if y + 1 < shape[0]: + edges.append([node, node + shape[1]]) + costs.append(1.0 if y != 2 else -3.0) + graph = bic.graph.UndirectedGraph.from_edges(shape[0] * shape[1], edges) + return graph, np.array(costs, dtype=np.float64) diff --git a/tests/graph/multicut/test_chain_decomposer.py b/tests/graph/multicut/test_chain_decomposer.py new file mode 100644 index 0000000..a7b428d --- /dev/null +++ b/tests/graph/multicut/test_chain_decomposer.py @@ -0,0 +1,83 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + +from ._helpers import same_partition + + +def test_chained_multicut_solvers(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs) + solver = bic.graph.ChainedMulticutSolvers( + [ + bic.graph.GreedyAdditiveMulticut(), + bic.graph.KernighanLinMulticut(number_of_outer_iterations=5), + ] + ) + + labels = solver.optimize(objective) + + same_partition(labels, [0, 0, 0, 1]) + assert objective.energy(labels) == pytest.approx(-5.0) + + +def test_chained_solver_rejects_empty_chain(): + with pytest.raises(ValueError, match="at least one"): + bic.graph.ChainedMulticutSolvers([]) + + +def test_multicut_decomposer_solves_positive_components(): + graph = bic.graph.UndirectedGraph.from_edges(4, [[0, 1], [1, 2], [2, 3]]) + objective = bic.graph.MulticutObjective(graph, [1.0, -5.0, 1.0]) + solver = bic.graph.MulticutDecomposer(bic.graph.GreedyAdditiveMulticut()) + + labels = solver.optimize(objective) + + same_partition(labels, [0, 0, 1, 1]) + assert objective.energy() == pytest.approx(-5.0) + + +def test_multicut_decomposer_uses_fallthrough_solver_for_single_component(): + class SingletonSolver(bic.graph.MulticutSolver): + def optimize(self, objective): + objective.labels = np.arange(objective.graph.number_of_nodes, dtype=np.uint64) + return objective.labels + + graph = bic.graph.UndirectedGraph.from_edges(2, [[0, 1]]) + objective = bic.graph.MulticutObjective(graph, [1.0]) + solver = bic.graph.MulticutDecomposer( + bic.graph.GreedyAdditiveMulticut(), + fallthrough_solver=SingletonSolver(), + ) + + labels = solver.optimize(objective) + + same_partition(labels, [0, 1]) + + +def test_decomposer_on_external_toy_problem(external_toy_problem): + graph, costs, _ = external_toy_problem + objective = bic.graph.MulticutObjective(graph, costs) + solver = bic.graph.MulticutDecomposer( + bic.graph.ChainedMulticutSolvers( + [ + bic.graph.GreedyAdditiveMulticut(), + bic.graph.KernighanLinMulticut(number_of_outer_iterations=10), + ] + ) + ) + + labels = solver.optimize(objective) + + assert objective.energy(labels) <= -35.0 + + +def test_decomposer_energy_bound_on_grid_problem(grid_problem): + graph, costs = grid_problem + objective = bic.graph.MulticutObjective(graph, costs) + solver = bic.graph.MulticutDecomposer(bic.graph.GreedyAdditiveMulticut()) + + labels = solver.optimize(objective) + + assert objective.energy(labels) <= -20.0 diff --git a/tests/graph/multicut/test_external_problem.py b/tests/graph/multicut/test_external_problem.py new file mode 100644 index 0000000..78a9f7f --- /dev/null +++ b/tests/graph/multicut/test_external_problem.py @@ -0,0 +1,29 @@ +import pytest + +import bioimage_cpp as bic + + +def test_solvers_on_full_external_problem(tmp_path): + try: + graph, costs = bic.graph.load_external_multicut_problem(timeout=30) + except (FileNotFoundError, RuntimeError) as error: + pytest.skip(str(error)) + + solvers = [ + bic.graph.GreedyAdditiveMulticut(), + bic.graph.GreedyFixationMulticut(), + bic.graph.KernighanLinMulticut(number_of_outer_iterations=5), + bic.graph.ChainedMulticutSolvers( + [ + bic.graph.GreedyAdditiveMulticut(), + bic.graph.KernighanLinMulticut(number_of_outer_iterations=5), + ] + ), + bic.graph.MulticutDecomposer(bic.graph.GreedyAdditiveMulticut()), + ] + + for solver in solvers: + objective = bic.graph.MulticutObjective(graph, costs) + labels = solver.optimize(objective) + assert labels.shape == (graph.number_of_nodes,) + assert objective.energy(labels) <= -76900.0 diff --git a/tests/graph/multicut/test_greedy_additive.py b/tests/graph/multicut/test_greedy_additive.py new file mode 100644 index 0000000..610e3d0 --- /dev/null +++ b/tests/graph/multicut/test_greedy_additive.py @@ -0,0 +1,44 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + +from ._helpers import same_partition + + +def test_greedy_additive_merges_positive_components(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyAdditiveMulticut().optimize(objective) + + same_partition(labels, [0, 0, 0, 1]) + np.testing.assert_array_equal(objective.labels, labels) + assert objective.energy() == pytest.approx(-5.0) + + +def test_greedy_additive_respects_weight_stop(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyAdditiveMulticut(weight_stop=1.5).optimize(objective) + + same_partition(labels, [0, 1, 1, 2]) + + +def test_greedy_additive_on_external_toy_problem(external_toy_problem): + graph, costs, _ = external_toy_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyAdditiveMulticut().optimize(objective) + + assert objective.energy(labels) <= -35.0 + + +def test_greedy_additive_energy_bound_on_grid_problem(grid_problem): + graph, costs = grid_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyAdditiveMulticut().optimize(objective) + + assert objective.energy(labels) <= -20.0 diff --git a/tests/graph/multicut/test_greedy_fixation.py b/tests/graph/multicut/test_greedy_fixation.py new file mode 100644 index 0000000..c1f65a5 --- /dev/null +++ b/tests/graph/multicut/test_greedy_fixation.py @@ -0,0 +1,38 @@ +import bioimage_cpp as bic + + +def test_greedy_fixation_respects_negative_constraints(frustrated_triangle): + graph, costs = frustrated_triangle + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyFixationMulticut().optimize(objective) + + assert labels[0] != labels[2] + assert objective.energy(labels) <= -3.0 + + +def test_greedy_fixation_node_num_stop(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyFixationMulticut(node_num_stop=3).optimize(objective) + + assert len(set(labels.tolist())) == 3 + + +def test_greedy_fixation_on_external_toy_problem(external_toy_problem): + graph, costs, _ = external_toy_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyFixationMulticut().optimize(objective) + + assert objective.energy(labels) <= -35.0 + + +def test_greedy_fixation_energy_bound_on_grid_problem(grid_problem): + graph, costs = grid_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.GreedyFixationMulticut().optimize(objective) + + assert objective.energy(labels) <= -20.0 diff --git a/tests/graph/multicut/test_kernighan_lin.py b/tests/graph/multicut/test_kernighan_lin.py new file mode 100644 index 0000000..b1cd562 --- /dev/null +++ b/tests/graph/multicut/test_kernighan_lin.py @@ -0,0 +1,54 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def test_kernighan_lin_improves_or_preserves_energy(frustrated_triangle): + graph, costs = frustrated_triangle + objective = bic.graph.MulticutObjective(graph, costs) + before = objective.energy() + + labels = bic.graph.KernighanLinMulticut(number_of_outer_iterations=10).optimize(objective) + + assert objective.energy(labels) <= before + assert labels.dtype == np.uint64 + + +def test_kernighan_lin_finds_frustrated_triangle_optimum(frustrated_triangle): + graph, costs = frustrated_triangle + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.KernighanLinMulticut(number_of_outer_iterations=10).optimize(objective) + + assert objective.energy(labels) == pytest.approx(-3.0) + + +def test_kernighan_lin_accepts_initial_labels(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs, initial_labels=[0, 0, 0, 1]) + before = objective.energy() + + labels = bic.graph.KernighanLinMulticut(number_of_outer_iterations=5).optimize(objective) + + assert objective.energy(labels) <= before + + +def test_kernighan_lin_on_external_toy_problem(external_toy_problem): + graph, costs, _ = external_toy_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.KernighanLinMulticut(number_of_outer_iterations=20).optimize(objective) + + # The move-chain implementation reaches the optimum on this instance. + assert objective.energy(labels) == pytest.approx(-35.0) + + +def test_kernighan_lin_energy_bound_on_grid_problem(grid_problem): + graph, costs = grid_problem + objective = bic.graph.MulticutObjective(graph, costs) + + labels = bic.graph.KernighanLinMulticut(number_of_outer_iterations=10).optimize(objective) + + # Regression guard pinned to the move-chain implementation's converged energy. + assert objective.energy(labels) <= -34.0 diff --git a/tests/graph/multicut/test_objective.py b/tests/graph/multicut/test_objective.py new file mode 100644 index 0000000..fde684d --- /dev/null +++ b/tests/graph/multicut/test_objective.py @@ -0,0 +1,40 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def test_multicut_objective_energy_and_validation(frustrated_triangle): + graph, costs = frustrated_triangle + objective = bic.graph.MulticutObjective(graph, costs) + + assert objective.energy([0, 0, 1]) == pytest.approx(-3.0) + objective.labels = [0, 0, 1] + assert objective.energy() == pytest.approx(-3.0) + + with pytest.raises(ValueError, match="edge_costs"): + bic.graph.MulticutObjective(graph, [1.0, 2.0]) + with pytest.raises(ValueError, match="labels"): + objective.labels = [0, 1] + + +def test_multicut_objective_copies_graph_topology(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs) + + graph.insert_edge(0, 3) + + assert objective.graph.number_of_edges == 3 + np.testing.assert_array_equal( + objective.graph.uv_ids(), + np.array([[0, 1], [1, 2], [2, 3]], dtype=np.uint64), + ) + + +def test_multicut_objective_reset_labels(chain_problem): + graph, costs = chain_problem + objective = bic.graph.MulticutObjective(graph, costs, initial_labels=[0, 0, 1, 1]) + + objective.reset_labels() + + np.testing.assert_array_equal(objective.labels, np.arange(4, dtype=np.uint64)) diff --git a/tests/graph/test_connected_components.py b/tests/graph/test_connected_components.py new file mode 100644 index 0000000..cbccc07 --- /dev/null +++ b/tests/graph/test_connected_components.py @@ -0,0 +1,54 @@ +import numpy as np +import pytest + +import bioimage_cpp as bic + + +def test_connected_components_without_edges(): + graph = bic.graph.UndirectedGraph(3) + + np.testing.assert_array_equal( + bic.graph.connected_components(graph), + np.array([0, 1, 2], dtype=np.uint64), + ) + + +def test_connected_components_for_undirected_graph(): + graph = bic.graph.UndirectedGraph.from_edges(5, [[0, 1], [2, 3], [1, 4]]) + + np.testing.assert_array_equal( + bic.graph.connected_components(graph), + np.array([0, 0, 1, 1, 0], dtype=np.uint64), + ) + + +def test_connected_components_with_edge_mask(): + graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2]]) + + np.testing.assert_array_equal( + bic.graph.connected_components(graph, edge_mask=np.array([True, False])), + np.array([0, 0, 1], dtype=np.uint64), + ) + + +def test_connected_components_accepts_region_adjacency_graph(): + labels = np.array([[1, 1, 2], [3, 3, 2]], dtype=np.uint64) + rag = bic.graph.region_adjacency_graph(labels) + + components = bic.graph.connected_components(rag) + + assert components.dtype == np.uint64 + assert components.shape == (4,) + assert components[0] != components[1] + np.testing.assert_array_equal(components[[1, 2, 3]], np.repeat(components[1], 3)) + + +def test_connected_components_rejects_invalid_mask(): + graph = bic.graph.UndirectedGraph.from_edges(3, [[0, 1], [1, 2]]) + + with pytest.raises(TypeError, match="edge_mask must have dtype bool"): + bic.graph.connected_components(graph, edge_mask=np.array([1, 0], dtype=np.uint8)) + with pytest.raises(ValueError, match="1D"): + bic.graph.connected_components(graph, edge_mask=np.ones((1, 2), dtype=bool)) + with pytest.raises(ValueError, match="number_of_edges"): + bic.graph.connected_components(graph, edge_mask=np.array([True]))