diff --git a/cpp/Makefile b/cpp/Makefile new file mode 100644 index 0000000..67df584 --- /dev/null +++ b/cpp/Makefile @@ -0,0 +1,12 @@ + +CXX = g++ +CXXFLAGS = -std=c++11 -Wall -O2 -fsanitize=address -fsanitize=undefined +LDLIBS = -l:libboost_unit_test_framework.a + +test_mwmatching: test_mwmatching.cpp mwmatching.h + $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(LDFLAGS) -o $@ $< $(LDLIBS) + +.PHONY: clean +clean: + $(RM) test_mwmatching + diff --git a/cpp/mwmatching.h b/cpp/mwmatching.h new file mode 100644 index 0000000..a011a05 --- /dev/null +++ b/cpp/mwmatching.h @@ -0,0 +1,2161 @@ +/** + * Algorithm for finding a maximum weight matching in general graphs. + * + * Depends on Boost.Hash. + * Tested with Boost v1.74, available from https://www.boost.org/ + */ + +#ifndef MWMATCHING_H_ +#define MWMATCHING_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +namespace mwmatching { + + +/** Type representing the unique ID of a vertex. */ +using VertexId = unsigned int; + +/** Type representing a pair of vertices. */ +using VertexPair = std::pair; + + +/** Type representing a weighted edge. + * + * @tparam WeightType Numeric type used to represent edge weights. + */ +template +struct Edge +{ + static_assert(std::numeric_limits::is_specialized, + "Edge weight must be a numeric type"); + + /** Incident vertices. */ + VertexPair vt; + + /** Edge weight. */ + WeightType weight; + + Edge(VertexPair vt, WeightType w) + : vt(vt), weight(w) + { } + + Edge(VertexId x, VertexId y, WeightType w) + : vt(x, y), weight(w) + { } +}; + + +namespace impl { + + +/** Value used to mark an invalid or undefined vertex. */ +constexpr VertexId NO_VERTEX = std::numeric_limits::max(); + + +/** Top-level blossoms may be labeled "S" or "T" or unlabeled. */ +enum BlossomLabel { LABEL_NONE = 0, LABEL_S = 1, LABEL_T = 2 }; + + +/** Return a pair of vertices in flipped order. */ +inline VertexPair flip_vertex_pair(const VertexPair& vt) +{ + return std::make_pair(vt.second, vt.first); +} + + +/** + * Check that the input is a valid graph. + * + * The graph may not have self-edges or multiple edges between + * the same pair of vertices. Edge weights must be within a valid range. + * + * This function takes time O(m). + * + * @param edges Vector of weighted edges defining the graph. + * @throw std::invalid_argument If the input graph is not valid. + */ +template +void check_input_graph(const std::vector>& edges) +{ + const VertexId max_num_vertex = std::numeric_limits::max(); + + for (const Edge& edge : edges) { + + // Check that vertex IDs are valid. + if ((edge.vt.first >= max_num_vertex) + || (edge.vt.second >= max_num_vertex)) { + throw std::invalid_argument("Vertex ID out of range"); + } + + // Check that edge weight is a finite number. + if (! std::numeric_limits::is_integer) { + if (! std::isfinite(edge.weight)) { + throw std::invalid_argument( + "Edge weights must be finite numbers"); + } + } + + // Check that edge weight will not cause overflow. + if (edge.weight > std::numeric_limits::max() / 4) { + throw std::invalid_argument( + "Edge weight exceeds maximum supported value"); + } + } + + // Check that the graph has no self-edges. + for (const Edge& edge : edges) { + if (edge.vt.first == edge.vt.second) { + throw std::invalid_argument("Self-edges are not supported"); + } + } + + // Check that the graph does not contain duplicate edges. + std::unordered_set> edges_seen; + for (const Edge& edge : edges) { + VertexPair vt = edge.vt; + if (vt.first > vt.second) { + vt = flip_vertex_pair(vt); + } + auto status = edges_seen.insert(vt); + if (!status.second) { + throw std::invalid_argument("Duplicate edges are not supported"); + } + } +} + + +/** + * Copy edges and remove edges with negative weight. + * + * This function takes time O(m). + * + * @param edges Vector of weighted edges defining the graph. + * @return Vector of edges with non-negative weight. + */ +template +std::vector> remove_negative_weight_edges( + const std::vector>& edges) +{ + std::vector> real_edges; + real_edges.reserve(edges.size()); + + for (const Edge& edge : edges) { + if (edge.weight >= 0) { + real_edges.push_back(edge); + } + } + + return real_edges; +} + + +// Forward declaration. +template struct NonTrivialBlossom; + + +/** + * Represents a blossom. + * + * A blossom is an odd-length alternating cycle over sub-blossoms. + * An alternating path consists of alternating matched and unmatched edges. + * An alternating cycle is an alternating path that starts and ends in + * the same sub-blossom. + * + * A single vertex by itself is also a blossom: a "trivial blossom". + * + * An instance of this class represents either a trivial blossom, + * or a non-trivial blossom which consists of multiple sub-blossoms. + * + * Blossoms are recursive structures: A non-trivial blossoms contains + * sub-blossoms, which may themselves contain sub-blossoms etc. + * + * Each blossom contains exactly one vertex that is not matched to another + * vertex in the same blossom. This is the "base vertex" of the blossom. + */ +template +struct Blossom { + + /** Parent of this blossom, or "nullptr" if this blossom is top-level. */ + NonTrivialBlossom* parent; + + /** Index of the base vertex of this blossom. */ + VertexId base_vertex; + + /** Label S or T if this is a top-level blossom in an alternating tree. */ + BlossomLabel label; + + /** Optional edge that attaches this blossom to the alternating tree. */ + VertexPair tree_edge; + + /** + * In case of a top-level S-blossom, "best_edge" is the least-slack edge + * that links to a different S-blossom, or "nullptr" if no such edge + * has been found. + */ +// TODO : consider storing a copy of the edge instead of pointer + const Edge* best_edge; + + /** Initialize a trivial (single-vertex) blossom. */ + explicit Blossom(VertexId base_vertex) + : parent(nullptr), + base_vertex(base_vertex), + label(LABEL_NONE), + best_edge(nullptr) + { } + + // Dummy virtual method to allow the use of dynamic_cast<>. + virtual ~Blossom() = default; +}; + + +/** + * Represents a non-trivial blossom. + * + * A non-trivial blossom is a blossom that contains multiple sub-blossoms + * (at least 3 sub-blossoms, since all blossoms have odd length). + * + * Non-trivial blossoms maintain a list of their sub-blossoms and the edges + * between their subblossoms. + * + * Unlike trivial blossoms, each non-trivial blossom is associated with + * a variable in the dual LPP problem. + * + * Non-trivial blossoms are created and destroyed by the matching algorithm. + * This implies that not every odd-length alternating cycle is a blossom; + * it only becomes a blossom through an explicit action of the algorithm. + * An existing blossom may change when the matching is augmented along + * a path that runs through the blossom. + */ +template +struct NonTrivialBlossom : public Blossom +{ + struct SubBlossom { + /** Pointer to sub-blossom. */ + Blossom* blossom; + + /** + * Ordered pair of vertices (x, y) where "x" is a vertex in "blossom" + * and "y" is a vertex in the next sub-blossom. + */ + VertexPair edge; + }; + + /** + * List of sub-blossoms. + * Ordered by their appearance in the alternating cycle. + * First blossom in the list is the start/end of the alternating cycle + * and contains the base vertex. + */ + std::list subblossoms; + + /** Dual LPP variable for this blossom. */ + WeightType dual_var; + + /** + * In case of a top-level S-blossom, "best_edge_set" is a list of + * least-slack edges between this blossom and other S-blossoms. + */ +// TODO : consider storing a copy of the edge instead of pointer + std::list*> best_edge_set; + + /** Initialize a non-trivial blossom. */ + NonTrivialBlossom( + const std::vector*>& subblossoms, + const std::deque& edges) + : Blossom(subblossoms.front()->base_vertex), + dual_var(0) + { + assert(subblossoms.size() == edges.size()); + assert(subblossoms.size() % 2 == 1); + assert(subblossoms.size() >= 3); + + auto blossom_it = subblossoms.begin(); + auto blossom_end = subblossoms.end(); + auto edge_it = edges.begin(); + while (blossom_it != blossom_end) { + this->subblossoms.emplace_back(); + this->subblossoms.back().blossom = *blossom_it; + this->subblossoms.back().edge = *edge_it; + ++blossom_it; + ++edge_it; + } + } + + /** Find the position of the specified subblossom. */ + std::pair::iterator> + find_subblossom(Blossom* subblossom) + { + VertexId pos = 0; + auto it = subblossoms.begin(); + while (it->blossom != subblossom) { + ++it; + ++pos; + assert(it != subblossoms.end()); + } + return make_pair(pos, it); + } +}; + + +/** + * Represents a list of edges forming an alternating path or + * an alternating cycle. + * + * The path is defined over top-level blossoms; it skips parts of the path + * that are internal to blossoms. Vertex pairs are oriented to match the + * direction of the path. + */ +struct AlternatingPath +{ +// TODO : consider avoiding dynamic allocations for the alternating path + std::deque edges; +}; + + +/** + * This class holds all data used by the matching algorithm. + * + * It contains the input graph, a partial solution, and several + * auxiliary data structures. + */ +template +struct MatchingContext +{ + // Alias for edge type. + using EdgeT = Edge; + using BlossomT = Blossom; + using NonTrivialBlossomT = NonTrivialBlossom; + + /** Specification of a delta step. */ + struct DeltaStep + { + /** Type of delta step: 1, 2, 3 or 4. */ + int kind; + + /** Delta value. */ + WeightType value; + + /** Edge on which the minimum delta occurs (for delta type 2 or 3). */ + VertexPair edge; + + /** Blossom on which the minimum delta occurs (for delta type 4). */ + NonTrivialBlossomT* blossom; + }; + + /** Scale integer edge weights to enable integer calculations. */ + static constexpr WeightType weight_factor = + std::numeric_limits::is_integer ? 2 : 1; + + /** Graph edges. */ + const std::vector edges; + + /** Number of vertices. */ + const VertexId num_vertex; + + /** For each vertex, a vector of pointers to its incident edges. */ +// TODO : consider copying the Edge instead of pointers + const std::vector> adjacent_edges; + + /** + * Current state of the matching. + * + * vertex_mate[x] == y if vertex "x" is matched to vertex "y". + * vertex_mate[x] == NO_VERTEX if vertex "x" is unmatched. + */ +// TODO : consider merging all per-vertex info into a single "struct Vertex" + std::vector vertex_mate; + + /** + * Each vertex is associated with a trivial blossom. + * + * "trivial_blossom[x]" is the trivial blossom that contains only + * vertex "x". + */ + std::vector trivial_blossom; + + /** + * Non-trivial blossoms may be created and destroyed during + * the course of the algorithm. + */ +// NOTE - this MUST be a list, because we delete items from it while keeping pointers to other items + std::list nontrivial_blossom; + + /** + * Every vertex is contained in exactly one top-level blossom + * (possibly the trivial blossom that contains just that vertex). + * + * "vertex_top_blossom[x]" is the top-level blossom that contains + * vertex "x". + */ + std::vector vertex_top_blossom; + + /** + * Every vertex has a variable in the dual LPP. + * + * "vertex_dual[x]" is the dual variable of vertex "x". + */ + std::vector vertex_dual; + + /** + * In case of T-vertex or unlabeled vertex "x", + * "vertex_best_edge[x]" is the least-slack edge to any S-vertex, + * or "nullptr" if no such edge has been found. + */ +// TODO - consider storing a copy of the edge instead of a pointer + std::vector vertex_best_edge; + + /** Queue of S-vertices to be scanned. */ + std::deque queue; + + /** Markers placed while tracing an alternating path. */ + std::vector vertex_marker; + + /** Vertices marked while tracing an alternating path. */ + std::vector marked_vertex; + + /** + * Initialize the matching algorithm. + * + * This function takes time O(n + m). + */ + explicit MatchingContext(const std::vector& edges_in) + : edges(remove_negative_weight_edges(edges_in)), + num_vertex(count_num_vertex(edges)), + adjacent_edges(build_adjacent_edges(edges, num_vertex)) + { + // Initially all vertices are unmatched. + vertex_mate.resize(num_vertex, NO_VERTEX); + + // Create a trivial blossom for each vertex. + trivial_blossom.reserve(num_vertex); + for (VertexId x = 0; x < num_vertex; ++x) { + trivial_blossom.emplace_back(x); + } + + // Initially all vertices are trivial top-level blossoms. + vertex_top_blossom.reserve(num_vertex); + for (VertexId x = 0; x < num_vertex; ++x) { + vertex_top_blossom.push_back(&trivial_blossom[x]); + } + + // Vertex duals are initialized to half the maximum edge weight. + WeightType max_weight = 0; + for (const EdgeT& edge : edges) { + max_weight = std::max(max_weight, edge.weight); + } + WeightType init_vertex_dual = max_weight * (weight_factor / 2); + vertex_dual.resize(num_vertex, init_vertex_dual); + + // Initialize "vertex_best_edge". + vertex_best_edge.resize(num_vertex, nullptr); + + // Allocate temporary arrays for path tracing. + vertex_marker.resize(num_vertex); + marked_vertex.reserve(num_vertex); + } + + // Prevent copying. + MatchingContext(const MatchingContext&) = delete; + MatchingContext& operator=(const MatchingContext&) = delete; + + /* + * Static methods used by the constructor: + */ + + /** Count the number of vertices in the graph. */ + static VertexId count_num_vertex(const std::vector& edges) + { + const VertexId max_num_vertex = std::numeric_limits::max(); + + VertexId num_vertex = 0; + for (const Edge& edge : edges) { + VertexId m = std::max(edge.vt.first, edge.vt.second); + assert(m < max_num_vertex); + num_vertex = std::max(num_vertex, m + 1); + } + + return num_vertex; + } + + /** + * Build adjacency lists for a graph. + * + * Note that the adjacency lists contain pointers to edge instances, + * and therefore become invalid if the input edge vector is modified + * in any way that invalidates its iterators. + * + * This function takes time O(m). + * + * @param num_vertex Number of vertices in the graph. + * @param edges List of edges in the graph. + * @return Vector of incident edges for each vertex. + */ + static std::vector> build_adjacent_edges( + const std::vector& edges, + VertexId num_vertex) + { + // Count the number of incident edges per vertex. + std::vector vertex_degree(num_vertex); + for (const EdgeT& edge : edges) { + assert(edge.vt.first < num_vertex); + assert(edge.vt.second < num_vertex); + vertex_degree[edge.vt.first]++; + vertex_degree[edge.vt.second]++; + } + + // Build adjacency lists. + std::vector> adjacent_edges(num_vertex); + for (VertexId i = 0; i < num_vertex; ++i) { + adjacent_edges[i].reserve(vertex_degree[i]); + } + for (const EdgeT& edge : edges) { + adjacent_edges[edge.vt.first].push_back(&edge); + adjacent_edges[edge.vt.second].push_back(&edge); + } + + return adjacent_edges; + } + + /* + * General support routines: + */ + + /** Calculate edge slack. */ + WeightType edge_slack(const EdgeT& edge) const + { + VertexId x = edge.vt.first; + VertexId y = edge.vt.second; + assert(vertex_top_blossom[x] != vertex_top_blossom[y]); + return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight; + } + + /** Call a function for every vertex inside the specified blossom. */ + template + void for_vertices_in_blossom(BlossomT* blossom, Func func) + { + NonTrivialBlossomT* ntb = dynamic_cast(blossom); + if (ntb) { + // Visit all vertices in the non-trivial blossom. + // Use an explicit stack to avoid deep call chains. + // TODO : re-use a global stack instance + std::vector stack; + stack.push_back(ntb); + + while (! stack.empty()) { + NonTrivialBlossomT* b = stack.back(); + stack.pop_back(); + + for (const auto& sub : b->subblossoms) { + ntb = dynamic_cast(sub.blossom); + if (ntb) { + stack.push_back(ntb); + } else { + func(sub.blossom->base_vertex); + } + } + } + + } else { + // A trivial blossom contains just one vertex. + func(blossom->base_vertex); + } + } + + /* + * Least-slack edge tracking: + * + * To calculate delta steps, the matching algorithm needs to find + * - the least-slack edge between any S-vertex and an unlabeled vertex; + * - the least-slack edge between any pair of top-level S-blossoms. + * + * For each unlabeled vertex and each T-vertex, we keep track of the + * least-slack edge to any S-vertex. Tracking for unlabeled vertices + * serves to provide the least-slack edge for the delta step. + * Tracking for T-vertices is done because such vertices can turn into + * unlabeled vertices if they are part of a T-blossom that gets expanded. + * + * For each top-level S-blossom, we keep track of the least-slack edge + * to any S-vertex not in the same blossom. + * + * Furthermore, for each top-level S-blossom, we keep a list of + * least-slack edges to other top-level S-blossoms. For any pair of + * top-level S-blossoms, the least-slack edge between them is contained + * in the edge list of at least one of the blossoms. An edge list may + * contain multiple edges to the same S-blossom. Such redundant edges are + * pruned during blossom merging to limit the number of tracked edges. + * + * Note: For a given vertex or blossom, the identity of the least-slack + * edge to any S-blossom remains unchanged during a delta step. + * Although the delta step changes edge slacks, it changes the slack + * of every edge to an S-vertex by the same amount. Therefore the edge + * that had least slack before the delta step, will still have least slack + * after the delta step. + */ + + /** + * Reset least-slack edge tracking. + * + * This function takes time O(n + m). + */ + void lset_reset() + { + for (VertexId x = 0; x < num_vertex; ++x) { + vertex_best_edge[x] = nullptr; + } + + for (BlossomT& blossom : trivial_blossom) { + blossom.best_edge = nullptr; + } + + for (NonTrivialBlossomT& blossom : nontrivial_blossom) { + blossom.best_edge = nullptr; + blossom.best_edge_set.clear(); + } + } + + /** + * Add edge "e" from an S-vertex to unlabeled vertex or T-vertex "y". + * + * This function takes time O(1) per call. + * This function is called O(m) times per stage. + */ + void lset_add_vertex_edge(VertexId y, const EdgeT* edge, WeightType slack) + { + const EdgeT* cur_best_edge = vertex_best_edge[y]; + if ((cur_best_edge == nullptr) + || (slack < edge_slack(*cur_best_edge))) { + vertex_best_edge[y] = edge; + } + } + + /** + * Return the index and slack of the least-slack edge between + * any S-vertex and unlabeled vertex. + * + * This function takes time O(n) per call. + * This function takes total time O(n**2) per stage. + * + * @return Pair (edge, slack) if there is a least-slack edge, + * or (nullptr, 0) if there is no suitable edge. + */ + std::pair lset_get_best_vertex_edge() + { + const EdgeT* best_edge = nullptr; + WeightType best_slack = 0; + + for (VertexId x = 0; x < num_vertex; ++x) { + if (vertex_top_blossom[x]->label == LABEL_NONE) { + const EdgeT* edge = vertex_best_edge[x]; + if (edge != nullptr) { + WeightType slack = edge_slack(*edge); + if ((best_edge == nullptr) || (slack < best_slack)) { + best_edge = edge; + best_slack = slack; + } + } + } + } + + return std::make_pair(best_edge, best_slack); + } + + /** Start tracking edges for a new S-blossom. */ + void lset_new_blossom(BlossomT* blossom) + { + assert(blossom->best_edge == nullptr); + NonTrivialBlossomT* ntb = dynamic_cast(blossom); + if (ntb) { + assert(ntb->best_edge_set.empty()); + } + } + + /** + * Add edge "e" between the specified S-blossom and another S-blossom. + * + * This function takes time O(1) per call. + * This function is called O(m) times per stage. + */ + void lset_add_blossom_edge(BlossomT* blossom, + const EdgeT* edge, + WeightType slack) + { + const EdgeT* cur_best_edge = blossom->best_edge; + if ((cur_best_edge == nullptr) + || (slack < edge_slack(*cur_best_edge))) { + blossom->best_edge = edge; + } + + NonTrivialBlossomT* ntb = dynamic_cast(blossom); + if (ntb) { + ntb->best_edge_set.push_back(edge); + } + } + + /** + * Update least-slack edge tracking after merging sub-blossoms + * into a new S-blossom. + * + * This function takes total time O(n**2) per stage. + */ + void lset_merge_blossoms(NonTrivialBlossomT* blossom) + { + assert(blossom->best_edge == nullptr); + assert(blossom->best_edge_set.empty()); + + // Collect edges from the sub-blossoms that used to be S-blossoms. + for (auto& subblossom_item : blossom->subblossoms) { + BlossomT* sub = subblossom_item.blossom; + if (sub->label == LABEL_S) { + NonTrivialBlossomT* ntb = + dynamic_cast(sub); + if (ntb) { + // Take least-slack edge set from this subblossom. + blossom->best_edge_set.splice( + blossom->best_edge_set.end(), + ntb->best_edge_set); + } else { + // Trivial blossoms don't maintain a least-slack edge set. + // Just consider all incident edges. + for (const EdgeT* edge : + adjacent_edges[sub->base_vertex]) { + // Only take edges between different S-blossoms. + VertexId x = edge->vt.first; + VertexId y = edge->vt.second; + BlossomT* bx = vertex_top_blossom[x]; + BlossomT* by = vertex_top_blossom[y]; + if ((bx != by) + && (bx->label == LABEL_S) + && (by->label == LABEL_S)) { + blossom->best_edge_set.push_back(edge); + } + } + } + } + } + + // Build a temporary array holding the least-slack edge index to + // each top-level S-blossom. This array is indexed by the base vertex + // of the blossoms. + std::vector> + best_edge_to_blossom(num_vertex); + + for (const EdgeT* edge : blossom->best_edge_set) { + BlossomT* bx = vertex_top_blossom[edge->vt.first]; + BlossomT* by = vertex_top_blossom[edge->vt.second]; + assert(bx == blossom || by == blossom); + + // Ignore internal edges. + if (bx != by) { + bx = (bx == blossom) ? by : bx; + + // Only consider edges to S-blossoms. + if (bx->label == LABEL_S) { + + // Keep only the least-slack edge to blossom "bx". + WeightType slack = edge_slack(*edge); + VertexId bx_base = bx->base_vertex; + auto& best_edge_item = best_edge_to_blossom[bx_base]; + if ((best_edge_item.first == nullptr) + || (slack < best_edge_item.second)) { + best_edge_item.first = edge; + best_edge_item.second = slack; + } + } + } + }; + + // Rebuild a compact list of least-slack edges. + // Also find the overall least-slack edge to any other S-blossom. + WeightType best_slack = 0; + blossom->best_edge_set.clear(); + for (auto& best_edge_item : best_edge_to_blossom) { + const EdgeT* edge = best_edge_item.first; + if (edge != nullptr) { + blossom->best_edge_set.push_back(edge); + WeightType slack = best_edge_item.second; + if ((blossom->best_edge == nullptr) || (slack < best_slack)) { + blossom->best_edge = edge; + best_slack = slack; + } + } + } + } + + /** + * Return the index and slack of the least-slack edge between + * any pair of top-level S-blossoms. + * + * This function takes time O(n) per call. + * This function takes total time O(n**2) per stage. + * + * @return Pair (edge, slack) if there is a least-slack edge, + * or (nullptr, 0) if there is no suitable edge. + */ + std::pair lset_get_best_blossom_edge() + { + const EdgeT* best_edge = nullptr; + WeightType best_slack = 0; + +// TODO - clean this up when we have a linked list of top-level blossoms + auto consider_blossom = + [this,&best_edge,&best_slack](BlossomT* blossom) + { + if ((blossom->parent == nullptr) && (blossom->label == LABEL_S)) { + const EdgeT* edge = blossom->best_edge; + if (edge != nullptr) { + WeightType slack = edge_slack(*edge); + if ((best_edge == nullptr) || (slack < best_slack)) { + best_edge = edge; + best_slack = slack; + } + } + } + }; + + for (BlossomT& blossom : trivial_blossom) { + consider_blossom(&blossom); + } + for (BlossomT& blossom : nontrivial_blossom) { + consider_blossom(&blossom); + } + + return std::make_pair(best_edge, best_slack); + } + + /* + * Creating and expanding blossoms: + */ + + /** + * Trace back through the alternating trees from vertices "x" and "y". + * + * If both vertices are part of the same alternating tree, this function + * discovers a new blossom. In this case it returns an alternating path + * through the blossom that starts and ends in the same sub-blossom. + * + * If the vertices are part of different alternating trees, this function + * discovers an augmenting path. In this case it returns an alternating + * path that starts and ends in an unmatched vertex. + * + * This function takes time O(k) to discover a blossom, where "k" is the + * number of sub-blossoms, or time O(n) to discover an augmenting path. + */ + AlternatingPath trace_alternating_paths(VertexId x, VertexId y) + { + assert(vertex_top_blossom[x] != vertex_top_blossom[y]); + + // Initialize a path containing only the edge (x, y). + AlternatingPath path; + path.edges.emplace_back(x, y); + + // "first_common" is the first common ancestor of "x" and "y" + // in the alternating tree, or "nullptr" if no common ancestor + // has been found. + BlossomT* first_common = nullptr; + + // Alternate between tracing the path from "x" and the path from "y". + // This ensures that the search time is bounded by the size of any + // newly found blossom. + while ((x != NO_VERTEX) || (y != NO_VERTEX)) { + + // Trace path from vertex "x". + if (x != NO_VERTEX) { + + // Stop if we found a common ancestor. + BlossomT* bx = vertex_top_blossom[x]; + if (vertex_marker[bx->base_vertex]) { + first_common = bx; + break; + } + + // Mark blossom as potential common ancestor. + vertex_marker[bx->base_vertex] = true; + marked_vertex.push_back(bx->base_vertex); + + // Trace back to the parent in the alternating tree. + x = bx->tree_edge.first; + if (x != NO_VERTEX) { + path.edges.push_front(bx->tree_edge); + } + } + + // Trace path from vertex "y". + if (y != NO_VERTEX) { + + // Stop if we found a common ancestor. + BlossomT* by = vertex_top_blossom[y]; + if (vertex_marker[by->base_vertex]) { + first_common = by; + break; + } + + // Mark blossom as potential common ancestor. + vertex_marker[by->base_vertex] = true; + marked_vertex.push_back(by->base_vertex); + + // Trace back to the parent in the alternating tree. + y = by->tree_edge.first; + if (y != NO_VERTEX) { + path.edges.emplace_back(by->tree_edge.second, y); + } + } + } + + // Remove all markers we placed. + for (VertexId k : marked_vertex) { + vertex_marker[k] = false; + } + marked_vertex.clear(); + + // If we found a common ancestor, trim the paths so they end there. + if (first_common) { + assert(first_common->label == LABEL_S); + while (vertex_top_blossom[path.edges.front().first] + != first_common) { + path.edges.pop_front(); + } + while (vertex_top_blossom[path.edges.back().second] + != first_common) { + path.edges.pop_back(); + } + } + + // Any alternating path between S-blossoms must have odd length. + assert(path.edges.size() % 2 == 1); + + return path; + } + + /** + * Create a new blossom from an alternating cycle. + * + * Assign label S to the new blossom. + * Relabel all T-sub-blossoms as S and add their vertices to the queue. + * + * This function takes total time O(n**2) per stage. + */ + void make_blossom(const AlternatingPath& path) + { + assert(path.edges.size() % 2 == 1); + assert(path.edges.size() >= 3); + + // Construct the list of sub-blossoms. + std::vector subblossoms; + subblossoms.reserve(path.edges.size()); + for (VertexPair edge : path.edges) { + subblossoms.push_back(vertex_top_blossom[edge.first]); + } + + // Check that the path is cyclic. + VertexId pos = 0; + for (VertexPair edge : path.edges) { + pos = (pos + 1) % subblossoms.size(); + assert(vertex_top_blossom[edge.second] == subblossoms[pos]); + } + + // Create the new blossom object. + nontrivial_blossom.emplace_back(subblossoms, path.edges); + NonTrivialBlossomT* blossom = &nontrivial_blossom.back(); + + // Link the subblossoms to the their new parent. + for (BlossomT* sub : subblossoms) { + sub->parent = blossom; + } + + // Mark vertices as belonging to the new blossom. + for_vertices_in_blossom(blossom, [this,blossom](VertexId x) { + vertex_top_blossom[x] = blossom; + }); + + // Assign label S to the new blossom and link to the alternating tree. + assert(subblossoms.front()->label == LABEL_S); + blossom->label = LABEL_S; + blossom->tree_edge = subblossoms.front()->tree_edge; + + // Consider vertices inside former T-sub-blossoms which now + // became S-vertices; add them to the queue. + for (BlossomT* sub : subblossoms) { + if (sub->label == LABEL_T) { + for_vertices_in_blossom(sub, + [this](VertexId x) { + queue.push_back(x); + }); + } + } + + // Merge least-slack edges for the S-sub-blossoms. + lset_merge_blossoms(blossom); + } + + /** Erase the specified non-trivial blossom. */ + void erase_nontrivial_blossom(NonTrivialBlossomT* blossom) + { + auto blossom_it = std::find_if( + nontrivial_blossom.begin(), + nontrivial_blossom.end(), + [blossom](const NonTrivialBlossomT& b) { + return (&b == blossom); + }); + assert(blossom_it != nontrivial_blossom.end()); + nontrivial_blossom.erase(blossom_it); + } + + /** + * Expand the specified T-blossom. + * + * This function takes time O(n). + */ + void expand_t_blossom(NonTrivialBlossomT* blossom) + { + assert(blossom->parent == nullptr); + assert(blossom->label == LABEL_T); + + // Convert sub-blossoms into top-level blossoms. + for (const auto& sub : blossom->subblossoms) { + BlossomT* sub_blossom = sub.blossom; + assert(sub_blossom->parent == blossom); + assert(sub_blossom->label == LABEL_NONE); + sub_blossom->parent = nullptr; + for_vertices_in_blossom(sub_blossom, + [this,sub_blossom](VertexId x) { + vertex_top_blossom[x] = sub_blossom; + }); + } + + // The expanded blossom was part of an alternating tree. + // We must now reconstruct the part of the alternating tree + // that ran through this blossom by linking some of the sub-blossoms + // into the tree. + + // Find the sub-blossom that was attached to the parent node + // in the alternating tree. + BlossomT* entry = vertex_top_blossom[blossom->tree_edge.second]; + + // Assign label T to that blossom and link to the alternating tree. + entry->label = LABEL_T; + entry->tree_edge = blossom->tree_edge; + + // Find the position of this sub-blossom within the expanding blossom. + auto subblossom_loc = blossom->find_subblossom(entry); + VertexId entry_pos = subblossom_loc.first; + auto entry_it = subblossom_loc.second; + + // Walk around the blossom from "entry" to the base + // in an even number of steps. + auto sub_it = entry_it; + if (entry_pos % 2 == 0) { + + // Walk backward to the base. + auto sub_begin = blossom->subblossoms.begin(); + while (sub_it != sub_begin) { + // Assign label S to the next node on the path. + --sub_it; + assign_label_s(sub_it->edge.first); + + // Assign label T to the next node on the path. + assert(sub_it != sub_begin); + --sub_it; + sub_it->blossom->label = LABEL_T; + sub_it->blossom->tree_edge = flip_vertex_pair(sub_it->edge); + } + + } else { + + // Walk forward to the base. + auto sub_end = blossom->subblossoms.end(); + while (sub_it != sub_end) { + // Assign label S to the next node on the path. + assign_label_s(sub_it->edge.second); + ++sub_it; + + // Assign label T to the next node on the path. + // We may wrap past the end of the subblossom list. + assert(sub_it != sub_end); + VertexPair& tree_edge = sub_it->edge; + ++sub_it; + + BlossomT *sub_blossom = (sub_it == sub_end) ? + blossom->subblossoms.front().blossom : + sub_it->blossom; + sub_blossom->label = LABEL_T; + sub_blossom->tree_edge = tree_edge; + } + } + + // Delete the expanded blossom. + erase_nontrivial_blossom(blossom); + } + + /** + * Expand the specified unlabeled blossom. + * + * This function takes time O(n). + */ + void expand_unlabeled_blossom(NonTrivialBlossomT* blossom) + { + assert(blossom->parent == nullptr); + assert(blossom->label == LABEL_NONE); + + // Convert sub-blossoms into top-level blossoms. + for (const auto& sub : blossom->subblossoms) { + BlossomT* sub_blossom = sub.blossom; + assert(sub_blossom->parent == blossom); + assert(sub_blossom->label == LABEL_NONE); + sub_blossom->parent = nullptr; + for_vertices_in_blossom(sub_blossom, + [this,sub_blossom](VertexId x) { + vertex_top_blossom[x] = sub_blossom; + }); + } + + // Delete the expanded blossom. + erase_nontrivial_blossom(blossom); + } + + /* + * Augmenting: + */ + + /** + * Augment along an alternating path through the specified blossom, + * from sub-blossom "entry" to the base vertex of the blossom. + * + * Modify the blossom to reflect that sub-blossom "entry" contains + * the base vertex after augmenting. + */ + void augment_blossom_rec( + NonTrivialBlossomT* blossom, + BlossomT* entry, + std::stack>& rec_stack) + { + // Find the position of "entry" within "blossom". + auto subblossom_loc = blossom->find_subblossom(entry); + VertexId entry_pos = subblossom_loc.first; + auto entry_it = subblossom_loc.second; + + // Walk around the blossom from "entry" to the base + // in an even number of steps. + auto sub_begin = blossom->subblossoms.begin(); + auto sub_end = blossom->subblossoms.end(); + auto sub_it = entry_it; + while ((sub_it != sub_begin) && (sub_it != sub_end)) { + VertexId x, y; + BlossomT* bx; + BlossomT* by; + + if (entry_pos % 2 == 0) { + // Walk backward to the base. + --sub_it; + by = sub_it->blossom; + assert(sub_it != sub_begin); + --sub_it; + bx = sub_it->blossom; + std::tie(x, y) = sub_it->edge; + } else { + // Walk forward to the base. + // We may wrap past the end of the subblossom list. + ++sub_it; + assert(sub_it != sub_end); + std::tie(x, y) = sub_it->edge; + bx = sub_it->blossom; + ++sub_it; + by = (sub_it == sub_end) ? + blossom->subblossoms.front().blossom : + sub_it->blossom; + } + + // Pull this edge into the matching. + vertex_mate[x] = y; + vertex_mate[y] = x; + + // Augment through any non-trivial subblossoms touching this edge. + NonTrivialBlossomT* bx_nt = dynamic_cast(bx); + if (bx_nt != nullptr) { + rec_stack.emplace(bx_nt, &trivial_blossom[x]); + } + + NonTrivialBlossomT* by_nt = dynamic_cast(by); + if (by_nt != nullptr) { + rec_stack.emplace(by_nt, &trivial_blossom[y]); + } + } + + // Re-orient the blossom. + if (entry_it != sub_begin) { + // Rotate the subblossom list. + blossom->subblossoms.splice(sub_begin, + blossom->subblossoms, + entry_it, + sub_end); + + // Update the base vertex. + // We can pull the new base vertex from the entry sub-blossom + // since its augmentation has already finished. + blossom->base_vertex = entry->base_vertex; + } + } + + /** + * Augment along an alternating path through the specified blossom, + * from sub-blossom "entry" to the base vertex of the blossom. + * + * Recursively handle sub-blossoms as needed. + * + * This function takes time O(n). + */ + void augment_blossom(NonTrivialBlossomT* blossom, BlossomT* entry) + { + // Use an explicit stack to avoid deep recursion. + std::stack> rec_stack; + rec_stack.emplace(blossom, entry); + + while (!rec_stack.empty()) { + NonTrivialBlossomT* outer_blossom; + BlossomT* inner_entry; + std::tie(outer_blossom, inner_entry) = rec_stack.top(); + + NonTrivialBlossomT* inner_blossom = inner_entry->parent; + assert(inner_blossom != nullptr); + + if (inner_blossom != outer_blossom) { + // After augmenting "inner_blossom", + // continue by augmenting its parent. + rec_stack.top() = std::make_pair(outer_blossom, + inner_blossom); + } else { + // After augmenting "inner_blossom", + // this entire "outer_blossom" will be finished. + rec_stack.pop(); + } + + // Augment "inner_blossom". + augment_blossom_rec(inner_blossom, inner_entry, rec_stack); + } + } + + /** + * Augment the matching through the specified augmenting path. + * + * This function takes time O(n). + */ + void augment_matching(const AlternatingPath& path) + { + // Check that the path starts and ends in an unmatched blossom. + assert(path.edges.size() % 2 == 1); + assert(vertex_mate[ + vertex_top_blossom[path.edges.front().first]->base_vertex] + == NO_VERTEX); + assert(vertex_mate[ + vertex_top_blossom[path.edges.back().second]->base_vertex] + == NO_VERTEX); + + // Process the unmatched edges on the augmenting path. + auto edge_it = path.edges.begin(); + auto edge_end = path.edges.end(); + while (edge_it != edge_end) { + VertexId x = edge_it->first; + VertexId y = edge_it->second; + + // Augment any non-trivial blossoms that touch this edge. + BlossomT* bx = vertex_top_blossom[x]; + NonTrivialBlossomT* bx_nt = + dynamic_cast(bx); + if (bx_nt != nullptr) { + augment_blossom(bx_nt, &trivial_blossom[x]); + } + + BlossomT* by = vertex_top_blossom[y]; + NonTrivialBlossomT* by_nt = + dynamic_cast(by); + if (by_nt != nullptr) { + augment_blossom(by_nt, &trivial_blossom[y]); + } + + // Pull this edge into the matching. + vertex_mate[x] = y; + vertex_mate[y] = x; + + // Move forward through the augmenting path to + // the next edge to be matched. + ++edge_it; + if (edge_it == edge_end) { + break; + } + ++edge_it; + } + } + + /* + * Labels and alternating tree: + */ + + /** + * Assign label S to the unlabeled blossom that contains vertex "x". + * + * If vertex "x" is matched, it becomes attached to the alternating tree + * via its matched edge. If vertex "x" is unmatched, it becomes the root + * of an alternating tree. + * + * All vertices in the newly labeled blossom are added to the scan queue. + * + * @pre "x" is an unlabeled vertex, either unmatched or matched to + * a T-vertex via a tight edge. + */ + void assign_label_s(VertexId x) + { + // Assign label S to the blossom that contains vertex "x". + BlossomT* bx = vertex_top_blossom[x]; + assert(bx->label == LABEL_NONE); // TODO - this assertion sometimes fails + bx->label = LABEL_S; + + VertexId y = vertex_mate[x]; + if (y == NO_VERTEX) { + // Vertex "x" is unmatched. + // It must be either a top-level vertex or the base vertex of + // a top-level blossom. + assert(bx->base_vertex == x); + + // Mark the blossom as root of an alternating tree. + bx->tree_edge = std::make_pair(NO_VERTEX, x); + + } else { + // Vertex "x" is matched to T-vertex "y". + assert(vertex_top_blossom[y]->label == LABEL_T); + + // Attach the blossom to the alternating tree via vertex "y". + bx->tree_edge = std::make_pair(y, x); + } + + // Start least-slack edge tracking for the S-blossom. + lset_new_blossom(bx); + + // Add all vertices inside the newly labeled S-blossom to the queue. + for_vertices_in_blossom(bx, + [this](VertexId x) { + queue.push_back(x); + }); + } + + /** + * Assign label T to the unlabeled blossom that contains vertex "y". + * + * Attach it to the alternating tree via edge (x, y). + * Then immediately assign label S to the mate of vertex "y". + * + * Note that this function may expand blossoms that contain vertex "y". + * + * @pre "x" is an S-vertex. + * @pre "y" is an unlabeled, matched vertex. + * @pre There is a tight edge between vertices "x" and "y". + */ + void assign_label_t(VertexId x, VertexId y) + { + assert(vertex_top_blossom[x]->label == LABEL_S); + + BlossomT* by = vertex_top_blossom[y]; + assert(by->label == LABEL_NONE); + + // If "y" is part of a zero-dual blossom, expand it. + // This would otherwise likely happen through a zero-delta4 step, + // so we can just do it now and avoid a substage. + NonTrivialBlossomT* ntb = dynamic_cast(by); + while (ntb != nullptr && ntb->dual_var == 0) { + expand_unlabeled_blossom(ntb); + by = vertex_top_blossom[y]; + assert(by->label == LABEL_NONE); + ntb = dynamic_cast(by); + } + + // Assign label T to the top-level blossom that contains vertex "y". + by->label = LABEL_T; + by->tree_edge = std::make_pair(x, y); + + // Assign label S to the blossom that is mated to the T-blossom. + VertexId z = vertex_mate[by->base_vertex]; + assert(z != NO_VERTEX); + assign_label_s(z); + } + + /** + * Add the edge between S-vertices "x" and "y". + * + * If the edge connects blossoms that are part of the same alternating + * tree, this function creates a new S-blossom and returns false. + * + * If the edge connects two different alternating trees, an augmenting + * path has been discovered. In this case the matching is augmented. + * + * @return True if the matching was augmented; otherwise false. + */ + bool add_s_to_s_edge(VertexId x, VertexId y) + { + assert(vertex_top_blossom[x]->label == LABEL_S); + assert(vertex_top_blossom[y]->label == LABEL_S); + + // Trace back through the alternating trees from "x" and "y". + AlternatingPath path = trace_alternating_paths(x, y); + + // If the path is a cycle, create a new blossom. + // Otherwise the path is an augmenting path. + // Note that an alternating starts and ends in the same blossom, + // but not necessarily in the same vertex within that blossom. + VertexId p = path.edges.front().first; + VertexId q = path.edges.back().second; + if (vertex_top_blossom[p] == vertex_top_blossom[q]) { + make_blossom(path); + return false; + } else { + augment_matching(path); + return true; + } + } + + /** + * Scan queued S-vertices to expand the alternating trees. + * + * The scan proceeds until either an augmenting path is found, + * or the queue of S-vertices becomes empty. + * If an augmenting path is found, it is used to augment the matching. + * + * New blossoms may be created during the scan. + * + * @return True if the matching was augmented; otherwise false. + */ + bool substage_scan() + { + // Process the queue of S-vertices to be scanned. + // This loop runs through O(n) iterations per stage. + while (! queue.empty()) { + + // Take one vertex from the queue. + VertexId x = queue.front(); + queue.pop_front(); + + assert(vertex_top_blossom[x]->label == LABEL_S); + + // Scan the edges that are incident on "x". + // This loop runs through O(m) iterations per stage. + for (const EdgeT* edge : adjacent_edges[x]) { + VertexId y = (edge->vt.first != x) ? edge->vt.first + : edge->vt.second; + + // Note: The top-level blossom of vertex "x" may change + // during this loop, so we need to refresh it in each pass. + BlossomT* bx = vertex_top_blossom[x]; + + // Ignore edges that are internal to a blossom. + if (bx == vertex_top_blossom[y]) { + continue; + } + + BlossomLabel ylabel = vertex_top_blossom[y]->label; + + // Check whether this edge is tight (has zero slack). + // Only tight edges may be part of an alternating tree. + WeightType slack = edge_slack(*edge); + if (slack <= 0) { + if (ylabel == LABEL_NONE) { + // Found a tight edge to an unlabeled blossom. + // Assign label T to the blossom that contains "y". + assign_label_t(x, y); + } else if (ylabel == LABEL_S) { + // Found a tight edge between two S-blossoms. + // Find either a new blossom or an augmenting path. + bool augmented = add_s_to_s_edge(x, y); + if (augmented) { + return true; + } + } + } else if (ylabel == LABEL_S) { + // Found a non-tight edge between two S-blossoms. + // Pass it to the least-slack edge tracker. + lset_add_blossom_edge(bx, edge, slack); + } + + if (ylabel != LABEL_S) { + // Found an to a T-vertex or unlabeled vertex "y". + // Pass it to the least-slack edge tracker. + // Tight edges must also be tracked in this way. + lset_add_vertex_edge(y, edge, slack); + } + } + } + + // No further S vertices to scan, and no augmenting path found. + return false; + } + + /* + * Delta steps: + */ + + /** + * Calculate a delta step in the dual LPP problem. + * + * This function returns the minimum of the 4 types of delta values, + * and the type of delta which obtain the minimum, and the edge or + * blossom that produces the minimum delta, if applicable. + * + * This function takes time O(n). + * + * @pre There is at least one S-vertex. + * + * @return Tuple (delta_type, delta_value, delta_edge, delta_blossom) + */ + DeltaStep substage_calc_dual_delta() + { + DeltaStep delta; + delta.blossom = nullptr; + + // Compute delta1: minimum dual variable of any S-vertex. + delta.kind = 1; + delta.value = std::numeric_limits::max(); + for (VertexId x = 0; x < num_vertex; ++x) { + if (vertex_top_blossom[x]->label == LABEL_S) { + delta.value = std::min(delta.value, vertex_dual[x]); + } + } + + // Compute delta2: minimum slack of any edge between an S-vertex and + // an unlabeled vertex. + const EdgeT* edge; + WeightType slack; + std::tie(edge, slack) = lset_get_best_vertex_edge(); + if ((edge != nullptr) && (slack <= delta.value)) { + delta.kind = 2; + delta.value = slack; + delta.edge = edge->vt; + } + + // Compute delta3: half minimum slack of any edge between two + // top-level S-blossoms. + std::tie(edge, slack) = lset_get_best_blossom_edge(); + if ((edge != nullptr) && (slack / 2 <= delta.value)) { + delta.kind = 3; + delta.value = slack / 2; + delta.edge = edge->vt; + } + + // Compute delta4: half minimum dual of a top-level T-blossom. + for (NonTrivialBlossomT& blossom : nontrivial_blossom) { + if ((! blossom.parent) && (blossom.label == LABEL_T)) { + if (blossom.dual_var / 2 <= delta.value) { + delta.kind = 4; + delta.value = blossom.dual_var / 2; + delta.blossom = &blossom; + } + } + } + + return delta; + } + + /** Apply a delta step to the dual LPP variables. */ + void substage_apply_delta_step(WeightType delta) + { + // Apply delta to dual variables of all vertices. + for (VertexId x = 0; x < num_vertex; ++x) { + BlossomLabel xlabel = vertex_top_blossom[x]->label; + if (xlabel == LABEL_S) { + // S-vertex: subtract delta from dual variable. + vertex_dual[x] -= delta; + } else if (xlabel == LABEL_T) { + // T-vertex: add delta to dual variable. + vertex_dual[x] += delta; + } + } + + // Apply delta to dual variables of top-level non-trivial blossoms. + for (NonTrivialBlossomT& blossom : nontrivial_blossom) { + if (blossom.parent == nullptr) { + if (blossom.label == LABEL_S) { + // S-blossom: add 2*delta to dual variable. + blossom.dual_var += 2 * delta; + } else if (blossom.label == LABEL_T) { + // T-blossom: subtract 2*delta from dual variable. + blossom.dual_var -= 2 * delta; + } + } + } + } + + /* + * Main algorithm: + */ + + /** + * Reset data which are only valid during a stage. + * + * Mark all blossoms as unlabeled, clear the queue, + * reset least-slack edge tracking. + */ + void reset_stage() + { + // Remove blossom labels. + for (BlossomT& blossom : trivial_blossom) { + blossom.label = LABEL_NONE; + } + for (BlossomT& blossom : nontrivial_blossom) { + blossom.label = LABEL_NONE; + } + + // Clear the S-vertex queue. + queue.clear(); + + // Reset least-slack edge tracking. + lset_reset(); + } + + /** + * Run one stage of the matching algorithm. + * + * The stage searches a maximum-weight augmenting path. + * If this path is found, it is used to augment the matching, + * thereby increasing the number of matched edges by 1. + * If no such path is found, the matching must already be optimal. + * + * This function takes time O(n**2). + * + * @return True if the matching was successfully augmented; + * false if no further improvement is possible. + */ + bool run_stage() + { + // Assign label S to all unmatched vertices and put them in the queue. + for (VertexId x = 0; x < num_vertex; ++x) { + if (vertex_mate[x] == NO_VERTEX) { + assign_label_s(x); + } + } + + // Stop if all vertices are matched. + // No further improvement is possible in that case. + // This avoids messy calculations of delta steps without any S-vertex. + if (queue.empty()) { + return false; + } + + // Each pass through the following loop is a "substage". + // The substage tries to find an augmenting path. + // If an augmenting path is found, we augment the matching and end + // the stage. Otherwise we update the dual LPP problem and enter the + // next substage, or stop if no further improvement is possible. + // + // This loop runs through at most O(n) iterations per stage. + bool augmented = false; + while (true) { + + // Grow alternating trees. + // End the stage if an augmenting path is found. + augmented = substage_scan(); + if (augmented) { + break; + } + + // Calculate delta step in the dual LPP problem. + DeltaStep delta = substage_calc_dual_delta(); + + // Apply the delta step to the dual variables. + substage_apply_delta_step(delta.value); + + if (delta.kind == 2) { + // Use the edge from S-vertex to unlabeled vertex that got + // unlocked through the delta update. + VertexId x = delta.edge.first; + VertexId y = delta.edge.second; + if (vertex_top_blossom[x]->label != LABEL_S) { + std::swap(x, y); + } + assign_label_t(x, y); + + } else if (delta.kind == 3) { + // Use the S-to-S edge that got unlocked by the delta update. + // This may reveal an augmenting path. + VertexId x = delta.edge.first; + VertexId y = delta.edge.second; + augmented = add_s_to_s_edge(x, y); + if (augmented) { + break; + } + + } else if (delta.kind == 4) { + // Expand the T-blossom that reached dual value 0 through + // the delta update. + assert(delta.blossom); + expand_t_blossom(delta.blossom); + + } else { + // No further improvement possible. End the stage. + assert(delta.kind == 1); + break; + } + } + + // Remove all labels, clear queue. + reset_stage(); + + // Return True if the matching was augmented. + return augmented; + } + + /** Run the matching algorithm. */ + void run() + { + // Improve the solution until no further improvement is possible. + // + // Each successful pass through this loop increases the number + // of matched edges by 1. + // + // This loop runs through at most (n/2 + 1) iterations. + // Each iteration takes time O(n**2). + while (run_stage()) ; + } + + /* + * Verify optimal solution: + */ + + /** + * Helper function for verifying the solution. + * + * Descend down the blossom tree to find edges that are contained + * in blossoms. + * + * Adjust the slack of all contained edges to account for the dual + * variables of its containing blossoms. + * + * On the way down, keep track of the sum of dual variables of + * the containing blossoms. + * + * On the way up, keep track of the total number of matched edges + * in the subblossoms. Then check that all blossoms with non-zero + * dual variable are "full". + * + * @return True if successful; + * false if a blossom with non-zero dual is not full. + */ + bool verify_blossom_edges(NonTrivialBlossomT* blossom, + std::unordered_map>& edge_duals) + { +// TODO : fix line length +// TODO : simplify code +// TODO : optimize +// TODO : try to think of a simpler way + + // For each vertex "x", + // "vertex_depth[x]" is the depth of the smallest blossom on + // the current descent path that contains "x". + std::vector vertex_depth(num_vertex); + + // At each depth, keep track of the sum of blossom duals + // along the current descent path. + std::vector path_sum_dual = {0}; + + // At each depth, keep track of the number of matched edges + // along the current ascent path. + std::vector path_num_matched = {0}; + + // Use an explicit stack to avoid deep recursion. + std::stack::iterator>> stack; + stack.emplace(blossom, blossom->subblossoms.begin()); + + while (! stack.empty()) { + VertexId depth = stack.size(); + auto& stack_elem = stack.top(); + + blossom = stack_elem.first; + auto subblossom_it = stack_elem.second; + + if (subblossom_it == blossom->subblossoms.begin()) { + // We just entered this sub-blossom. + // Update the depth of all vertices in this sub-blossom. + for_vertices_in_blossom(blossom, + [&vertex_depth,depth](VertexId x) { + vertex_depth[x] = depth; + }); + + // Calculate the sum of blossom duals at the new depth. + path_sum_dual.push_back(path_sum_dual.back() + + blossom->dual_var); + + // Initialize the number of matched edges at the new depth. + path_num_matched.push_back(0); + + if (blossom->subblossoms.size() < 3) { + return false; + } + } + + if (subblossom_it != blossom->subblossoms.end()) { + + // Update the sub-blossom pointer at the current depth. + ++stack_elem.second; + + // Examine the current sub-blossom. + BlossomT* sub = subblossom_it->blossom; + NonTrivialBlossomT* ntb = dynamic_cast(sub); + if (ntb) { + // Prepare to descend into this sub-blossom. + stack.emplace(std::make_pair(ntb, ntb->subblossoms.begin())); + } else { + // Handle single vertex. + // For each incident edge, find the smallest blossom + // that contains it. + VertexId x = sub->base_vertex; + for (const EdgeT* edge : adjacent_edges[x]) { + // Only consider edges pointing out from "x". + if (edge->vt.first == x) { + VertexId y = edge->vt.second; + VertexId edge_depth = vertex_depth[y]; + if (edge_depth > 0) { + // This edge is contained in a blossom. + // Update its slack. + edge_duals[edge->vt] += path_sum_dual[edge_depth]; + + // Update the number of matched edges in the blossom. + if (vertex_mate[x] == y) { + path_num_matched[edge_depth] += 1; + } + } + } + } + } + + } else { + // We are leaving the current sub-blossom. + + // Count the number of vertices inside this blossom. + VertexId blossom_num_vertex = 0; + for_vertices_in_blossom(blossom, + [&blossom_num_vertex](VertexId x) { + ++blossom_num_vertex; + }); + + // Check that the blossom is "full". + // A blossom is full if all except one of its vertices + // are matched to another vertex within the blossom. + VertexId blossom_num_matched = path_num_matched[depth]; + if (blossom_num_vertex != 2 * blossom_num_matched + 1) { + return false; + } + + // Update the number of matched edges in the parent blossom. + path_num_matched[depth - 1] += path_num_matched[depth]; + + // Push vertices in this sub-blossom back to the parent. + for_vertices_in_blossom(blossom, + [&vertex_depth,depth](VertexId x) { + vertex_depth[x] = depth - 1; + }); + + // Trim the descending path. + path_sum_dual.pop_back(); + path_num_matched.pop_back(); + + // Remove the current blossom from the stack. + stack.pop(); + } + } + + return true; + } + + /** + * Verify that the optimum solution has been found. + * + * This function takes time O(n**2). + * + * @return True if the solution is optimal; otherwise false. + */ + bool verify_optimum() + { + // Count matched vertices and check symmetry of "vertex_mate". + VertexId num_matched_vertex = 0; + for (VertexId x = 0; x < num_vertex; ++x) { + VertexId y = vertex_mate[x]; + if (y != NO_VERTEX) { + ++num_matched_vertex; + if (vertex_mate[y] != x) { + return false; + } + } + } + + // Check that all matched edges exist in the graph. + VertexId num_matched_edge = 0; + for (const EdgeT& edge : edges) { + if (vertex_mate[edge.vt.first] == edge.vt.second) { + ++num_matched_edge; + } + } + + if (num_matched_vertex != 2 * num_matched_edge) { + return false; + } + + // Check that all dual variables are non-negative. + for (VertexId x = 0; x < num_vertex; ++x) { + if (vertex_dual[x] < 0) { + return false; + } + } + + for (NonTrivialBlossomT& blossom : nontrivial_blossom) { + if (blossom.dual_var < 0) { + return false; + } + } + + // Check that all unmatched vertices have zero dual. + for (VertexId x = 0; x < num_vertex; ++x) { + if ((vertex_mate[x] == NO_VERTEX) && (vertex_dual[x] != 0)) { + return false; + } + } + + // Calculate the slack of each edge. + // A correction will be needed for edges inside blossoms. + std::unordered_map> edge_duals; + for (const EdgeT& edge : edges) { +// TODO : consider potential numeric overflow + edge_duals[edge.vt] = vertex_dual[edge.vt.first] + + vertex_dual[edge.vt.second]; + } + + // Descend down each top-level blossom. + // Adjust edge slacks to account for blossom dual. + // Also check that all blossoms are full. + // This takes total time O(n**2). + for (NonTrivialBlossomT& blossom : nontrivial_blossom) { + if (blossom.parent == nullptr) { + verify_blossom_edges(&blossom, edge_duals); + } + } + + // Check that all edges have non-negative slack. + // Check that all matched edges have zero slack. +// TODO : fix potential numeric overflow - WeightType may be unsigned + for (const EdgeT& edge : edges) { + WeightType slack = edge_duals[edge.vt] + - weight_factor * edge.weight; + if (slack < 0) { + return false; + } + if (vertex_mate[edge.vt.first] == edge.vt.second) { + if (slack != 0) { + return false; + } + } + } + + // Optimum solution confirmed. + return true; + } +}; + + +} // namespace impl + + +/** + * Compute a maximum-weighted matching in the general undirected weighted + * graph given by "edges". + * + * The graph is specified as a list of edges, each edge specified as a tuple + * of its two vertices and the edge weight. + * There may be at most one edge between any pair of vertices. + * No vertex may have an edge to itself. + * The graph may be non-connected (i.e. contain multiple components). + * + * Vertices are indexed by consecutive, non-negative integers, such that + * the first vertex has index 0 and the last vertex has index (n-1). + * Edge weights may be integers or floating point numbers. + * + * Isolated vertices (not incident to any edge) are allowed, but not + * recommended since such vertices consume time and memory but have + * no effect on the maximum-weight matching. + * Edges with negative weight are ignored. + * + * This function takes time O(n**3), where "n" is the number of vertices. + * This function uses O(n + m) memory, where "m" is the number of edges. + * + * @tparam WeightType Type used to represent edge weights. + * For example "long" or "double". + * @param edges List of weighted edges defining the graph. + * + * @return Vector of pairs of matched vertex indices. + * This is a subset of the edges in the graph. + * + * @throw std::invalid_argument If the input graph is not valid. + */ +template +std::vector maximum_weight_matching( + const std::vector>& edges) +{ + // Check that the input meets all constraints. + impl::check_input_graph(edges); + + // Run matching algorithm. + impl::MatchingContext matching(edges); + matching.run(); + + // Verify that the solution is optimal (works only for integer weights). + if (std::numeric_limits::is_integer) { + assert(matching.verify_optimum()); + } + + // Extract the matched edges. + std::vector solution; + for (const Edge& edge : matching.edges) { + if (matching.vertex_mate[edge.vt.first] == edge.vt.second) { + solution.push_back(edge.vt); + } + } + + return solution; +} + + +/** + * Adjust edge weights to compute a maximum cardinality matching. + * + * This function adjusts edge weights in the graph such that the + * maximum-weight matching of the adjusted graph is a maximum-cardinality + * matching, equal to a matching in the original graph that has maximum weight + * out of all matchings with maximum cardinality. + * + * The graph is specified as a vector of edges. + * Negative edge weights are allowed. + * + * This function increases all edge weights by an equal amount such that + * the adjusted weights satisfy the following conditions: + * - All edge weights are positive; + * - The minimum edge weight is at least "n" times the difference between + * maximum and minimum edge weight. + * + * This function increases edge weights by an amount that is proportional + * to the product of the unadjusted weight range and the number of vertices + * in the graph. This may fail if it would cause adjusted edge weights + * to exceed the supported numeric range. In case of floating point weights, + * the weight adjustment may also cause increased rounding errors. + * + * This function takes time O(m), where "m" is the number of edges. + * + * @tparam WeightType Type used to represent edge weights. + * @param edges Vector of weighted edges defining the graph. + * + * @return Vector of edges with adjusted weights. + * If no adjustments are necessary, this will be a copy of the + * input vector. + * + * @throw std::invalid_argument If the graph is invalid or edge weights + * exceed the supported range. + */ +template +std::vector> adjust_weights_for_maximum_cardinality_matching( + const std::vector>& edges_in) +{ + const WeightType min_safe_weight = + std::numeric_limits::min() / 4; + const WeightType max_safe_weight = + std::numeric_limits::max() / 4; + + // Copy edges. + std::vector> edges(edges_in); + + // Don't worry about empty graphs. + if (edges.empty()) { + return edges; + } + + // Count number of vertices. + // Determine minimum and maximum edge weight. + VertexId num_vertex = 0; + WeightType min_weight = edges.first().weight; + WeightType max_weight = min_weight; + + const VertexId max_num_vertex = std::numeric_limits::max(); + for (const Edge& edge : edges) { + VertexId m = std::max(edge.vt.first, edge.vt.second); + if (m >= max_num_vertex) { + throw std::invalid_argument("Vertex ID out of range"); + } + num_vertex = std::max(num_vertex, m + 1); + + if (! std::numeric_limits::is_integer) { + if (! std::isfinite(edge.weight)) { + throw std::invalid_argument( + "Edge weights must be finite numbers"); + } + } + + min_weight = std::min(min_weight, edge.weight); + max_weight = std::max(max_weight, edge.weight); + } + + // Calculate weight range and required weight adjustment. + if ((min_weight < min_safe_weight) || (max_weight > max_safe_weight)) { + throw std::invalid_argument( + "Edge weight exceeds maximum supported value"); + } + + WeightType weight_range = max_weight - min_weight; + + if (weight_range > max_safe_weight / num_vertex) { + throw std::invalid_argument( + "Adjusted edge weight exceeds maximum supported value"); + } + + // Do nothing if the weights already ensure maximum-cardinality. + if ((min_weight > 0) && (min_weight >= num_vertex * weight_range)) { + return edges; + } + + WeightType delta; + if (weight_range > 0) { + // Increase weights to make minimum edge weight large enough + // to improve any non-maximum-cardinality matching. + delta = num_vertex * weight_range - min_weight; + } else { + // All weights are the same. Increase weights to make them positive. + delta = 1 - min_weight; + } + + assert(delta >= 0); + + if (delta > max_safe_weight - max_weight) { + throw std::invalid_argument( + "Adjusted edge weight exceeds maximum supported value"); + } + + // Increase all edge weights by "delta". + for (Edge& edge: edges) { + edge.weight += delta; + } + + return edges; +} + + +} // namespace mwmatching + +#endif // MWMATCHING_H_ diff --git a/cpp/test_mwmatching.cpp b/cpp/test_mwmatching.cpp new file mode 100644 index 0000000..19d14cb --- /dev/null +++ b/cpp/test_mwmatching.cpp @@ -0,0 +1,618 @@ +/* + * Unit tests for maximum weight matching. + * + * Depends on the Boost.Test unit test framework. + * Tested with Boost v1.74, available from https://www.boost.org/ + */ + +#include +#include +#include +#include + +#define BOOST_TEST_MODULE mwmatching +#include + +#include "mwmatching.h" + + +using EdgeVectorLong = std::vector>; +using EdgeVectorDouble = std::vector>; +using Matching = std::vector; + + +BOOST_AUTO_TEST_SUITE(test_maximum_weight_matching) + +BOOST_AUTO_TEST_CASE(test10_empty) +{ + EdgeVectorLong edges = {}; + Matching expect = {}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test11_single_edge) +{ + EdgeVectorLong edges = {{0, 1, 1}}; + Matching expect = {{0, 1}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test12_simple) +{ + EdgeVectorLong edges = {{1, 2, 10}, {2, 3, 11}}; + Matching expect = {{2, 3}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test13_simple) +{ + EdgeVectorLong edges = {{1, 2, 5}, {2, 3, 11}, {3, 4, 5}}; + Matching expect = {{2, 3}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test15_float) +{ + EdgeVectorDouble edges = {{1, 2, 3.1415}, {2, 3, 2.7183}, {1, 3, 3.0}, {1, 4, 1.4142}}; + Matching expect = {{2, 3}, {1, 4}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test16_negative) +{ + EdgeVectorLong edges = {{1, 2, 2}, {1, 3, -2}, {2, 3, 1}, {2, 4, -1}, {3, 4, -6}}; + Matching expect = {{1, 2}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test20_sblossom) +{ + EdgeVectorLong edges = {{1, 2, 8}, {1, 3, 9}, {2, 3, 10}, {3, 4, 7}}; + Matching expect = {{1, 2}, {3, 4}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test20a_sblossom) +{ + EdgeVectorLong edges = {{1, 2, 8}, {1, 3, 9}, {2, 3, 10}, {3, 4, 7}, {1, 6, 5}, {4, 5, 6}}; + Matching expect = {{2, 3}, {1, 6}, {4, 5}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test21_tblossom) +{ + EdgeVectorLong edges = {{1, 2, 9}, {1, 3, 8}, {2, 3, 10}, {1, 4, 5}, {4, 5, 4}, {1, 6, 3}}; + Matching expect = {{2, 3}, {4, 5}, {1, 6}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test21a_tblossom) +{ + EdgeVectorLong edges = {{1, 2, 9}, {1, 3, 8}, {2, 3, 10}, {1, 4, 5}, {4, 5, 3}, {1, 6, 4}}; + Matching expect = {{2, 3}, {4, 5}, {1, 6}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test21b_tblossom) +{ + EdgeVectorLong edges = {{1, 2, 9}, {1, 3, 8}, {2, 3, 10}, {1, 4, 5}, {4, 5, 3}, {3, 6, 4}}; + Matching expect = {{1, 2}, {4, 5}, {3, 6}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test22_s_nest) +{ + EdgeVectorLong edges = {{1, 2, 9}, {1, 3, 9}, {2, 3, 10}, {2, 4, 8}, {3, 5, 8}, {4, 5, 10}, {5, 6, 6}}; + Matching expect = {{1, 3}, {2, 4}, {5, 6}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test23_s_relabel_nest) +{ + EdgeVectorLong edges = { + {1, 2, 10}, {1, 7, 10}, {2, 3, 12}, {3, 4, 20}, {3, 5, 20}, {4, 5, 25}, {5, 6, 10}, {6, 7, 10}, {7, 8, 8}}; + Matching expect = {{1, 2}, {3, 4}, {5, 6}, {7, 8}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test24_s_nest_expand) +{ + EdgeVectorLong edges = { + {1, 2, 8}, {1, 3, 8}, {2, 3, 10}, {2, 4, 12}, {3, 5, 12}, {4, 5, 14}, {4, 6, 12}, {5, 7, 12}, {6, 7, 14}, + {7, 8, 12}}; + Matching expect = {{1, 2}, {3, 5}, {4, 6}, {7, 8}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test25_s_t_expand) +{ + EdgeVectorLong edges = { + {1, 2, 23}, {1, 5, 22}, {1, 6, 15}, {2, 3, 25}, {3, 4, 22}, {4, 5, 25}, {4, 8, 14}, {5, 7, 13}}; + Matching expect = {{1, 6}, {2, 3}, {4, 8}, {5, 7}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test26_s_nest_t_expand) +{ + EdgeVectorLong edges = { + {1, 2, 19}, {1, 3, 20}, {1, 8, 8}, {2, 3, 25}, {2, 4, 18}, {3, 5, 18}, {4, 5, 13}, {4, 7, 7}, {5, 6, 7}}; + Matching expect = {{1, 8}, {2, 3}, {4, 7}, {5, 6}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test30_tnasty_expand) +{ + EdgeVectorLong edges = { + {1, 2, 45}, {1, 5, 45}, {2, 3, 50}, {3, 4, 45}, {4, 5, 50}, {1, 6, 30}, {3, 9, 35}, {4, 8, 35}, + {5, 7, 26}, {9, 10, 5}}; + Matching expect = {{2, 3}, {1, 6}, {4, 8}, {5, 7}, {9, 10}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test31_tnasty_expand) +{ + EdgeVectorLong edges = { + {1, 2, 45}, {1, 5, 45}, {2, 3, 50}, {3, 4, 45}, {4, 5, 50}, {1, 6, 30}, {3, 9, 35}, {4, 8, 26}, + {5, 7, 40}, {9, 10, 5}}; + Matching expect = {{2, 3}, {1, 6}, {4, 8}, {5, 7}, {9, 10}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test32_t_expand_leastslack) +{ + EdgeVectorLong edges = { + {1, 2, 45}, {1, 5, 45}, {2, 3, 50}, {3, 4, 45}, {4, 5, 50}, {1, 6, 30}, {3, 9, 35}, {4, 8, 28}, + {5, 7, 26}, {9, 10, 5}}; + Matching expect = {{2, 3}, {1, 6}, {4, 8}, {5, 7}, {9, 10}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test33_nest_tnasty_expand) +{ + EdgeVectorLong edges = { + {1, 2, 45}, {1, 7, 45}, {2, 3, 50}, {3, 4, 45}, {4, 5, 95}, {4, 6, 94}, {5, 6, 94}, {6, 7, 50}, + {1, 8, 30}, {3, 11, 35}, {5, 9, 36}, {7, 10, 26}, {11, 12, 5}}; + Matching expect = {{2, 3}, {4, 6}, {1, 8}, {5, 9}, {7, 10}, {11, 12}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test34_nest_relabel_expand) +{ + EdgeVectorLong edges = { + {1, 2, 40}, {1, 3, 40}, {2, 3, 60}, {2, 4, 55}, {3, 5, 55}, {4, 5, 50}, {1, 8, 15}, {5, 7, 30}, + {7, 6, 10}, {8, 10, 10}, {4, 9, 30}}; + Matching expect = {{1, 2}, {3, 5}, {7, 6}, {8, 10}, {4, 9}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test41_nonmax_card) +{ + EdgeVectorLong edges = {{0, 1, 2}, {0, 4, 3}, {1, 2, 7}, {1, 5, 2}, {2, 3, 9}, {2, 5, 4}, {3, 4, 8}, {3, 5, 4}}; + Matching expect = {{1, 2}, {3, 4}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test42_s_nest_partial_expand) +{ + /* + * [0]--8--[1]--6--[3]--5--[5] + * \ | | + * \ 9 8 + * 8 | | + * \--[2]--7--[4] + */ + EdgeVectorLong edges = {{0, 1, 8}, {0, 2, 8}, {1, 2, 9}, {1, 3, 6}, {2, 4, 7}, {3, 4, 8}, {3, 5, 5}}; + Matching expect = {{0, 1}, {2, 4}, {3, 5}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test43_s_nest_noexpand) +{ + /* + * [1]--9--[2] + * | / + * 7 ___7 + * | / + * [0] [5]--2--[6] + * | \___ + * 7 7 + * | \ + * [3]--9--[4] + */ + EdgeVectorLong edges = {{0, 1, 7}, {0, 2, 7}, {1, 2, 9}, {0, 3, 7}, {0, 4, 7}, {3, 4, 9}, {5, 6, 2}}; + Matching expect = {{1, 2}, {3, 4}, {5, 6}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test44_blossom_redundant_edge) +{ + /* + * [1]----9---[2] + * / | \ + * 7 8 \ + * / | 1 + * [0]--6--[4]--9--[3] | + * \ | + * \----1----[5] + */ + EdgeVectorLong edges = {{0, 1, 7}, {0, 4, 6}, {1, 2, 9}, {2, 3, 8}, {3, 4, 9}, {3, 5, 1}, {4, 5, 1}}; + Matching expect = {{1, 2}, {3, 4}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test45_flip_order) +{ + EdgeVectorLong edges = { + {0, 1, 2}, {0, 4, 3}, {1, 2, 7}, {1, 5, 1}, {2, 3, 9}, {2, 5, 4}, {3, 4, 8}, {3, 5, 4}}; + Matching expect = {{1, 2}, {3, 4}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); + + EdgeVectorLong redges(edges.rbegin(), edges.rend()); + Matching rexpect(expect.rbegin(), expect.rend()); + BOOST_TEST(mwmatching::maximum_weight_matching(redges) == rexpect); + + for (mwmatching::Edge& edge : edges) { + std::swap(edge.vt.first, edge.vt.second); + } + for (mwmatching::VertexPair& vt : expect) { + std::swap(vt.first, vt.second); + } + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); + + edges = {{0, 2, 4}, {0, 3, 4}, {0, 4, 1}, {1, 2, 8}, {1, 5, 3}, {2, 3, 9}, {3, 4, 7}, {4, 5, 2}}; + expect = {{1, 2}, {3, 4}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test46_expand_unlabeled_blossom) +{ + /* + * 3--[1] + * / | + * [2]--1--[4] 7 + * \ | + * 5--[3]--5--[5] + */ + EdgeVectorLong edges = {{1, 3, 7}, {1, 4, 3}, {2, 4, 1}, {3, 4, 5}, {3, 5, 5}}; + Matching expect1 = {{1, 3}, {2, 4}}; + Matching expect2 = {{1, 4}, {3, 5}}; + Matching result = mwmatching::maximum_weight_matching(edges); + BOOST_TEST((result == expect1 || result == expect2)); +} + +BOOST_AUTO_TEST_CASE(test47_expand_unlabeled_outer) +{ + /* + * [3]--10--[1]--15--[2]--12--[5] + * _/ \_ | | + * 11 16_ 8 15 + * / \ | | + * [4] [6]---7--[7] + */ + EdgeVectorLong edges = { + {1, 2, 15}, {1, 3, 10}, {1, 4, 11}, {1, 6, 17}, {2, 5, 12}, {2, 6, 8}, {5, 7, 15}, {6, 7, 7}}; + Matching expect = {{1, 4}, {2, 6}, {5, 7}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test48_expand_unlabeled_nested) +{ + /* + * [5]--11--[1]--11 14--[3] + * | \ / | + * 12 [4] 14 + * | / \ | + * [6]--11--[2]--12 11--[7] + */ + EdgeVectorLong edges = { + {1, 2, 12}, {1, 4, 11}, {1, 5, 11}, {2, 4, 12}, {2, 6, 11}, {3, 4, 14}, {3, 7, 14}, {4, 7, 11}}; + Matching expect = {{1, 5}, {2, 4}, {3, 7}}; + BOOST_TEST(mwmatching::maximum_weight_matching(edges) == expect); +} + +BOOST_AUTO_TEST_CASE(test_fail_bad_input) +{ + EdgeVectorDouble inf_weight = {{1, 2, std::numeric_limits::infinity()}}; + BOOST_CHECK_THROW(mwmatching::maximum_weight_matching(inf_weight), std::invalid_argument); + + EdgeVectorDouble huge_weight = {{1, 2, std::numeric_limits::max() / 2}}; + BOOST_CHECK_THROW(mwmatching::maximum_weight_matching(huge_weight), std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE(test_fail_bad_graph) +{ + EdgeVectorLong self_edge = {{0, 1, 2}, {1, 1, 1}}; + BOOST_CHECK_THROW(mwmatching::maximum_weight_matching(self_edge), std::invalid_argument); + + EdgeVectorLong duplicate_edge = {{0, 1, 2}, {1, 2, 1}, {2, 1, 1}}; + BOOST_CHECK_THROW(mwmatching::maximum_weight_matching(duplicate_edge), std::invalid_argument); +} + + +// TODO : test with unsigned weight type + + +BOOST_AUTO_TEST_SUITE_END() + + +/* + +class TestCornerCases(unittest.TestCase): + """Test cases that would catch certain errors in the algorithm. + + These graphs were generated semi-automatically to fail when + specific bugs are introduced in the code. + """ + + def test1(self): + pairs = mwm([(0,4,26), (1,3,31), (1,4,49)]) + self.assertEqual(pairs, [(0,4), (1,3)]) + + def test2(self): + pairs = mwm([(0,2,42), (0,4,36), (2,3,26)]) + self.assertEqual(pairs, [(0,4), (2,3)]) + + def test3(self): + pairs = mwm([(0,4,43), (1,4,28), (2,4,38)]) + self.assertEqual(pairs, [(0,4)]) + + def test4(self): + pairs = mwm([(0,1,50), (0,3,46), (0,4,45)]) + self.assertEqual(pairs, [(0,1)]) + + def test5(self): + pairs = mwm([(0,1,35), (0,3,36), (0,4,46)]) + self.assertEqual(pairs, [(0,4)]) + + def test6(self): + pairs = mwm([(0,1,50), (0,4,51), (0,5,34), (1,2,43), (1,4,57), (2,5,47), (3,4,17)]) + self.assertEqual(pairs, [(0,1), (2,5), (3,4)]) + + def test7(self): + pairs = mwm([(0,1,34), (0,3,19), (1,2,45), (1,3,30), (1,4,37), (2,4,36)]) + self.assertEqual(pairs, [(0,1), (2,4)]) + + def test8(self): + pairs = mwm([(0,1,48), (0,3,42), (0,4,57), (1,3,51), (1,5,36), (2,3,23), (4,5,46)]) + self.assertEqual(pairs, [(0,1), (2,3), (4,5)]) + + def test9(self): + pairs = mwm([(0,1,21), (0,2,25), (0,5,42), (1,4,40), (2,3,10), (2,5,40), (3,5,31), (4,5,58)]) + self.assertEqual(pairs, [(0,2), (1,4), (3,5)]) + + def test10(self): + pairs = mwm([(0,2,7), (0,5,20), (1,2,50), (1,4,46), (2,3,35), (2,4,8), (2,5,25), (3,5,47)]) + self.assertEqual(pairs, [(0,5), (1,4), (2,3)]) + + def test11(self): + pairs = mwm([(0,1,42), (0,2,60), (1,3,34), (1,4,58), (1,5,52), (2,5,60), (3,5,34), (4,5,57)]) + self.assertEqual(pairs, [(0,2), (1,4), (3,5)]) + + def test12(self): + pairs = mwm([(0,1,23), (0,2,26), (0,3,22), (0,4,41), (2,4,36)]) + self.assertEqual(pairs, [(0,1), (2,4)]) + + def test13(self): + pairs = mwm([(0,3,58), (0,4,49), (1,5,34), (2,3,22), (2,5,42), (4,5,36)]) + self.assertEqual(pairs, [(0,4), (1,5), (2,3)]) + + def test14(self): + pairs = mwm([(0,1,29), (0,3,35), (0,4,42), (1,2,12), (2,4,29), (3,4,44)]) + self.assertEqual(pairs, [(0,1), (3,4)]) + + def test15(self): + pairs = mwm([(0,4,53), (0,5,42), (1,4,45), (2,4,59), (2,6,39), (4,5,69), (4,6,52)]) + self.assertEqual(pairs, [(0,5), (1,4), (2,6)]) + + def test16(self): + pairs = mwm([(0,2,13), (1,2,11), (2,3,39), (2,4,17), (3,4,35)]) + self.assertEqual(pairs, [(0,2), (3,4)]) + + def test17(self): + pairs = mwm([(0,1,48), (0,2,44), (0,4,48), (1,4,36), (3,4,31)]) + self.assertEqual(pairs, [(0,2), (1,4)]) + + +class TestAdjustWeightForMaxCardinality(unittest.TestCase): + """Test adjust_weights_for_maximum_cardinality_matching() function.""" + + def test_empty(self): + self.assertEqual(adj([]), []) + + def test_chain(self): + self.assertEqual( + adj([(0,1,2), (1,2,8), (2,3,3), (3,4,9), (4,5,1), (5,6,7), (6,7,4)]), + [(0,1,65), (1,2,71), (2,3,66), (3,4,72), (4,5,64), (5,6,70), (6,7,67)]) + + def test_chain_preadjusted(self): + self.assertEqual( + adj([(0,1,65), (1,2,71), (2,3,66), (3,4,72), (4,5,64), (5,6,70), (6,7,67)]), + [(0,1,65), (1,2,71), (2,3,66), (3,4,72), (4,5,64), (5,6,70), (6,7,67)]) + + def test_flat(self): + self.assertEqual( + adj([(0,1,0), (0,4,0), (1,2,0), (1,5,0), (2,3,0), (2,5,0), (3,4,0), (3,5,0)]), + [(0,1,1), (0,4,1), (1,2,1), (1,5,1), (2,3,1), (2,5,1), (3,4,1), (3,5,1)]) + + def test14_maxcard(self): + self.assertEqual( + adj([(1,2,5), (2,3,11), (3,4,5)]), + [(1,2,30), (2,3,36), (3,4,30)]) + + def test16_negative(self): + self.assertEqual( + adj([(1,2,2), (1,3,-2), (2,3,1), (2,4,-1), (3,4,-6)]), + [(1,2,48), (1,3,44), (2,3,47), (2,4,45), (3,4,40)]) + + +class TestMaximumCardinalityMatching(unittest.TestCase): + """Test maximum cardinality matching.""" + + def test14_maxcard(self): + """maximum cardinality""" + self.assertEqual( + mwm(adj([(1,2,5), (2,3,11), (3,4,5)])), + [(1,2), (3,4)]) + + def test16_negative(self): + """negative weights""" + self.assertEqual( + mwm(adj([(1,2,2), (1,3,-2), (2,3,1), (2,4,-1), (3,4,-6)])), + [(1,3), (2,4)]) + + def test43_maxcard(self): + """maximum cardinality""" + self.assertIn( + mwm(adj([(0,1,2), (0,4,3), (1,2,7), (1,5,2), (2,3,9), (2,5,4), (3,4,8), (3,5,4)])), + ([(0,1), (2,5), (3,4)], + [(0,4), (1,2), (3,5)])) + + +class TestGraphInfo(unittest.TestCase): + """Test _GraphInfo helper class.""" + + def test_empty(self): + graph = mwmatching._GraphInfo([]) + self.assertEqual(graph.num_vertex, 0) + self.assertEqual(graph.edges, []) + self.assertEqual(graph.adjacent_edges, []) + + +class TestVerificationFail(unittest.TestCase): + """Test failure handling in verification routine.""" + + def _make_context( + self, + edges, + vertex_mate, + vertex_dual_2x, + nontrivial_blossom): + ctx = Mock(spec=mwmatching._MatchingContext) + ctx.graph = mwmatching._GraphInfo(edges) + ctx.vertex_mate = vertex_mate + ctx.vertex_dual_2x = vertex_dual_2x + ctx.nontrivial_blossom = nontrivial_blossom + return ctx + + def test_success(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[-1, 2, 1], + vertex_dual_2x=[0, 20, 2], + nontrivial_blossom=[]) + mwmatching._verify_optimum(ctx) + + def test_asymmetric_matching(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[-1, 2, 0], + vertex_dual_2x=[0, 20, 2], + nontrivial_blossom=[]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_nonexistent_matched_edge(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[2, -1, 0], + vertex_dual_2x=[11, 11, 11], + nontrivial_blossom=[]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_negative_vertex_dual(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[-1, 2, 1], + vertex_dual_2x=[-2, 22, 0], + nontrivial_blossom=[]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_unmatched_nonzero_dual(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[-1, 2, 1], + vertex_dual_2x=[9, 11, 11], + nontrivial_blossom=[]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_negative_edge_slack(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[-1, 2, 1], + vertex_dual_2x=[0, 11, 11], + nontrivial_blossom=[]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_matched_edge_slack(self): + edges = [(0,1,10), (1,2,11)] + ctx = self._make_context( + edges, + vertex_mate=[-1, 2, 1], + vertex_dual_2x=[0, 20, 11], + nontrivial_blossom=[]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_negative_blossom_dual(self): + # + # [0]--7--[1]--9--[2]--6--[3] + # \ / + # \----8-----/ + # + edges = [(0,1,7), (0,2,8), (1,2,9), (2,3,6)] + blossom = mwmatching._NonTrivialBlossom( + subblossoms=[ + mwmatching._Blossom(0), + mwmatching._Blossom(1), + mwmatching._Blossom(2)], + edges=[0,2,1]) + for sub in blossom.subblossoms: + sub.parent = blossom + blossom.dual_var = -1 + ctx = self._make_context( + edges, + vertex_mate=[1, 0, 3, 2], + vertex_dual_2x=[4, 6, 8, 4], + nontrivial_blossom=[blossom]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + + def test_blossom_not_full(self): + # + # [3] [4] + # | | + # 8 8 + # | | + # [0]--7--[1]--5--[2] + # \ / + # \----2-----/ + # + edges = [(0,1,7), (0,2,2), (1,2,5), (0,3,8), (1,4,8)] + blossom = mwmatching._NonTrivialBlossom( + subblossoms=[ + mwmatching._Blossom(0), + mwmatching._Blossom(1), + mwmatching._Blossom(2)], + edges=[0,2,1]) + for sub in blossom.subblossoms: + sub.parent = blossom + blossom.dual_var = 2 + ctx = self._make_context( + edges, + vertex_mate=[3, 4, -1, 0, 1], + vertex_dual_2x=[4, 10, 0, 12, 6], + nontrivial_blossom=[blossom]) + with self.assertRaises(mwmatching.MatchingError): + mwmatching._verify_optimum(ctx) + +*/