Add fusion moves implementation

Author: constantinpape

MIGRATION_GUIDE.md

@@ -408,6 +408,77 @@ Intentional differences vs. nifty:
   `GreedyAdditiveMulticut` and no fallthrough is given — the greedy solver
   already operates on each connected component internally.
 
+## Fusion Moves (Multicut)
+
+Nifty exposes the fusion-move multicut solver via the factory hierarchy with a
+chosen proposal generator and sub-solver factory.
+
+Nifty:
+
+```python
+import nifty.graph.opt.multicut as nmc
+
+objective = nmc.multicutObjective(graph, edge_costs)
+pgen = nmc.watershedProposals(sigma=1.0, numberOfSeeds=0.1)
+factory = nmc.fusionMoveBasedFactory(
+    proposalGenerator=pgen,
+    fusionMove=nmc.fusionMoveSettings(
+        mcFactory=nmc.greedyAdditiveFactory(),
+    ),
+    numberOfIterations=10,
+    stopIfNoImprovement=4,
+)
+labels = factory.create(objective).optimize()
+```
+
+bioimage-cpp:
+
+```python
+import bioimage_cpp as bic
+
+objective = bic.graph.MulticutObjective(graph, edge_costs)
+solver = bic.graph.FusionMoveMulticut(
+    proposal_generator=bic.graph.WatershedProposalGenerator(
+        sigma=1.0, n_seeds_fraction=0.1, seed=0,
+    ),
+    sub_solver=bic.graph.GreedyAdditiveMulticut(),
+    number_of_iterations=10,
+    stop_if_no_improvement=4,
+)
+labels = solver.optimize(objective)
+```
+
+Proposal generators:
+
+| nifty proposal generator | bioimage-cpp proposal generator |
+| --- | --- |
+| `watershedProposals(sigma=..., numberOfSeeds=...)` | `WatershedProposalGenerator(sigma=..., n_seeds_fraction=..., seed=...)` |
+| `greedyAdditiveProposals(sigma=..., weightStopCond=..., nodeNumStopCond=...)` | `GreedyAdditiveProposalGenerator(sigma=..., weight_stop=..., node_num_stop=..., seed=...)` |
+
+Sub-solvers: any built-in multicut solver (`GreedyAdditiveMulticut`,
+`GreedyFixationMulticut`, `KernighanLinMulticut`). If `sub_solver` is omitted,
+the default is `GreedyAdditiveMulticut` constructed with no-noise defaults.
+
+Intentional differences vs. nifty:
+
+- Single object construction: no separate factory / solver step.
+- Proposal generators are Python classes carrying their settings; the C++
+  proposal-generator object is built lazily when `optimize` is called.
+- The driver warm-starts from the trivial singleton labeling by running the
+  default greedy-additive sub-solver once before the proposal loop.
+- A best-of-three safety net (current, proposal, fused) keeps the running
+  energy monotonically non-increasing across iterations.
+- The current implementation is single-threaded and uses pairwise
+  (proposal, current) fuses. `number_of_threads` and
+  `number_of_parallel_proposals` are kept on the API for forward compatibility
+  but must be left at their defaults (`1` and `2` respectively).
+
+Notes:
+
+- Custom Python proposal generators are not yet supported; subclass
+  `ProposalGenerator` and provide your own `_build` returning a C++
+  proposal-generator object if you need to extend the set.
+
 ## Segmentation Overlaps
 
 Nifty:

include/bioimage_cpp/graph/detail/fusion_contract.hxx

@@ -0,0 +1,107 @@
+#pragma once
+
+#include "bioimage_cpp/detail/relabel.hxx"
+#include "bioimage_cpp/detail/union_find.hxx"
+#include "bioimage_cpp/graph/undirected_graph.hxx"
+
+#include <cstddef>
+#include <cstdint>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace bioimage_cpp::graph::detail {
+
+struct AgreementContraction {
+    // Contracted graph: one node per agreement component, dense ids in
+    // [0, number_of_components).
+    UndirectedGraph contracted_graph;
+    // For every original edge id: the contracted edge id it maps onto, or
+    // -1 if both endpoints collapsed into the same component.
+    std::vector<std::int64_t> contracted_edge_of_original;
+    // For every original node id: its dense agreement-component id.
+    std::vector<std::uint64_t> root_of_node;
+};
+
+// Build the agreement-projection contracted graph: merge `(u, v)` iff every
+// proposal labels `u` and `v` identically, then dense-relabel and assemble the
+// contracted graph with one edge per distinct surviving (root_u, root_v) pair.
+//
+// `proposals` is a row-major buffer of shape (n_proposals, number_of_nodes).
+// Passing 0 proposals collapses every edge (all proposals trivially agree)
+// and yields a single-node contracted graph; we reject that as a usage error.
+inline AgreementContraction contract_by_agreement(
+    const UndirectedGraph &graph,
+    const std::uint64_t *proposals,
+    const std::size_t n_proposals,
+    const std::size_t n_nodes_per_proposal
+) {
+    if (n_proposals == 0) {
+        throw std::invalid_argument("at least one proposal is required");
+    }
+    const auto number_of_nodes = graph.number_of_nodes();
+    const auto number_of_edges = graph.number_of_edges();
+    if (n_nodes_per_proposal != static_cast<std::size_t>(number_of_nodes)) {
+        throw std::invalid_argument(
+            "proposal width must equal number_of_nodes, got " +
+            std::to_string(n_nodes_per_proposal) + " for number_of_nodes=" +
+            std::to_string(number_of_nodes)
+        );
+    }
+
+    bioimage_cpp::detail::UnionFind sets(static_cast<std::size_t>(number_of_nodes));
+    for (std::uint64_t edge = 0; edge < number_of_edges; ++edge) {
+        const auto uv = graph.uv(edge);
+        const auto u = static_cast<std::size_t>(uv.first);
+        const auto v = static_cast<std::size_t>(uv.second);
+        bool agree = true;
+        for (std::size_t p = 0; p < n_proposals; ++p) {
+            const auto *row = proposals + p * n_nodes_per_proposal;
+            if (row[u] != row[v]) {
+                agree = false;
+                break;
+            }
+        }
+        if (agree) {
+            sets.merge(uv.first, uv.second);
+        }
+    }
+
+    std::vector<std::uint64_t> raw_root(static_cast<std::size_t>(number_of_nodes));
+    for (std::uint64_t node = 0; node < number_of_nodes; ++node) {
+        raw_root[static_cast<std::size_t>(node)] = sets.find(node);
+    }
+    auto root_of_node = bioimage_cpp::detail::dense_relabel(raw_root);
+
+    std::uint64_t number_of_components = 0;
+    for (const auto root : root_of_node) {
+        if (root + 1 > number_of_components) {
+            number_of_components = root + 1;
+        }
+    }
+
+    UndirectedGraph contracted_graph(number_of_components);
+    std::vector<std::int64_t> contracted_edge_of_original(
+        static_cast<std::size_t>(number_of_edges), -1
+    );
+
+    for (std::uint64_t edge = 0; edge < number_of_edges; ++edge) {
+        const auto uv = graph.uv(edge);
+        const auto ru = root_of_node[static_cast<std::size_t>(uv.first)];
+        const auto rv = root_of_node[static_cast<std::size_t>(uv.second)];
+        if (ru == rv) {
+            continue;
+        }
+        const auto inserted = contracted_graph.insert_edge(ru, rv);
+        contracted_edge_of_original[static_cast<std::size_t>(edge)] =
+            static_cast<std::int64_t>(inserted);
+    }
+
+    return AgreementContraction{
+        std::move(contracted_graph),
+        std::move(contracted_edge_of_original),
+        std::move(root_of_node),
+    };
+}
+
+} // namespace bioimage_cpp::graph::detail

include/bioimage_cpp/graph/multicut/fusion_move.hxx

@@ -0,0 +1,188 @@
+#pragma once
+
+#include "bioimage_cpp/graph/detail/fusion_contract.hxx"
+#include "bioimage_cpp/graph/multicut/greedy_additive.hxx"
+#include "bioimage_cpp/graph/multicut/objective.hxx"
+#include "bioimage_cpp/graph/proposal_generator.hxx"
+#include "bioimage_cpp/graph/undirected_graph.hxx"
+
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <memory>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+namespace bioimage_cpp::graph::multicut {
+
+class FusionMoveSolver final : public SolverBase {
+public:
+    // The proposal generator is borrowed; the caller owns its lifetime.
+    // The sub-solver pointer is optional: when null, the driver uses a default
+    // greedy-additive sub-solver internally.
+    //
+    // `number_of_threads` and `number_of_parallel_proposals` are reserved for
+    // future use; v1 only supports the single-threaded pairwise
+    // (proposal, current) fuse and rejects other values.
+    FusionMoveSolver(
+        ProposalGeneratorBase &proposal_generator,
+        const SolverBase *sub_solver = nullptr,
+        const std::size_t number_of_iterations = 10,
+        const std::size_t stop_if_no_improvement = 4,
+        const std::size_t number_of_threads = 1,
+        const std::size_t number_of_parallel_proposals = 2
+    )
+        : proposal_generator_(proposal_generator),
+          sub_solver_(sub_solver),
+          number_of_iterations_(number_of_iterations),
+          stop_if_no_improvement_(stop_if_no_improvement) {
+        if (number_of_threads != 1) {
+            throw std::invalid_argument(
+                "FusionMoveSolver currently supports number_of_threads=1 only"
+            );
+        }
+        if (number_of_parallel_proposals != 2) {
+            throw std::invalid_argument(
+                "FusionMoveSolver currently supports number_of_parallel_proposals=2 only"
+            );
+        }
+    }
+
+    std::vector<std::uint64_t> optimize(Objective &objective) const override {
+        const auto &graph = objective.graph();
+        const auto &costs = objective.costs();
+        const auto number_of_nodes = graph.number_of_nodes();
+
+        std::vector<std::uint64_t> current = objective.labels();
+        if (number_of_nodes == 0 || graph.number_of_edges() == 0) {
+            objective.set_labels(current);
+            return objective.labels();
+        }
+
+        const auto default_sub_solver = GreedyAdditiveSolver();
+        const auto &sub_solver = sub_solver_ != nullptr
+            ? *sub_solver_
+            : static_cast<const SolverBase &>(default_sub_solver);
+
+        // Warm start from greedy-additive if the caller passed the trivial
+        // singleton labeling.
+        if (is_singleton_labeling(current)) {
+            Objective warm_objective(graph, costs);
+            current = default_sub_solver.optimize(warm_objective);
+        }
+
+        double current_energy = energy(graph, costs, current);
+
+        std::vector<std::uint64_t> proposal(static_cast<std::size_t>(number_of_nodes));
+        std::size_t iterations_without_improvement = 0;
+
+        for (std::size_t iteration = 0; iteration < number_of_iterations_; ++iteration) {
+            proposal_generator_.generate(current, proposal);
+
+            const auto fused = fuse_pair(graph, costs, current, proposal, sub_solver);
+            const auto fused_energy = energy(graph, costs, fused);
+            const auto proposal_energy = energy(graph, costs, proposal);
+
+            // Best-of safety net across the three candidates so an iteration
+            // can never raise the running energy.
+            double best_energy = current_energy;
+            const std::vector<std::uint64_t> *best = &current;
+            if (fused_energy < best_energy) {
+                best_energy = fused_energy;
+                best = &fused;
+            }
+            if (proposal_energy < best_energy) {
+                best_energy = proposal_energy;
+                best = &proposal;
+            }
+
+            if (best_energy < current_energy) {
+                current = *best;
+                current_energy = best_energy;
+                iterations_without_improvement = 0;
+            } else {
+                ++iterations_without_improvement;
+                if (iterations_without_improvement >= stop_if_no_improvement_) {
+                    break;
+                }
+            }
+        }
+
+        objective.set_labels(current);
+        return objective.labels();
+    }
+
+private:
+    static bool is_singleton_labeling(const std::vector<std::uint64_t> &labels) {
+        for (std::size_t index = 0; index < labels.size(); ++index) {
+            if (labels[index] != static_cast<std::uint64_t>(index)) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    static std::vector<std::uint64_t> fuse_pair(
+        const UndirectedGraph &graph,
+        const std::vector<double> &costs,
+        const std::vector<std::uint64_t> &current,
+        const std::vector<std::uint64_t> &proposal,
+        const SolverBase &sub_solver
+    ) {
+        const auto number_of_nodes = static_cast<std::size_t>(graph.number_of_nodes());
+        std::vector<std::uint64_t> stacked(2 * number_of_nodes);
+        std::copy(current.begin(), current.end(), stacked.begin());
+        std::copy(proposal.begin(), proposal.end(), stacked.begin() + number_of_nodes);
+
+        auto contraction = ::bioimage_cpp::graph::detail::contract_by_agreement(
+            graph, stacked.data(), 2, number_of_nodes
+        );
+
+        const auto &contracted_graph = contraction.contracted_graph;
+        const auto number_of_contracted_edges = contracted_graph.number_of_edges();
+
+        std::vector<double> contracted_costs(
+            static_cast<std::size_t>(number_of_contracted_edges), 0.0
+        );
+        for (std::uint64_t edge = 0; edge < graph.number_of_edges(); ++edge) {
+            const auto target = contraction.contracted_edge_of_original[
+                static_cast<std::size_t>(edge)
+            ];
+            if (target < 0) {
+                continue;
+            }
+            contracted_costs[static_cast<std::size_t>(target)] +=
+                costs[static_cast<std::size_t>(edge)];
+        }
+
+        if (number_of_contracted_edges == 0) {
+            std::vector<std::uint64_t> result(number_of_nodes);
+            for (std::uint64_t node = 0; node < graph.number_of_nodes(); ++node) {
+                result[static_cast<std::size_t>(node)] = contraction.root_of_node[
+                    static_cast<std::size_t>(node)
+                ];
+            }
+            return result;
+        }
+
+        Objective sub_objective(contracted_graph, std::move(contracted_costs));
+        const auto sub_labels = sub_solver.optimize(sub_objective);
+
+        std::vector<std::uint64_t> result(number_of_nodes);
+        for (std::uint64_t node = 0; node < graph.number_of_nodes(); ++node) {
+            const auto root = contraction.root_of_node[static_cast<std::size_t>(node)];
+            result[static_cast<std::size_t>(node)] = sub_labels[
+                static_cast<std::size_t>(root)
+            ];
+        }
+        return result;
+    }
+
+    ProposalGeneratorBase &proposal_generator_;
+    const SolverBase *sub_solver_;
+    std::size_t number_of_iterations_;
+    std::size_t stop_if_no_improvement_;
+};
+
+} // namespace bioimage_cpp::graph::multicut