joris
/
maximum-weight-matching
1
0
Fork 0
Code Pull Requests Activity

Add struct Graph to represent input graph

This commit is contained in:
Joris van Rantwijk 2023-05-27 23:29:05 +02:00
parent 9b064de4d8
commit 4b7fab3f43
1 changed files with 171 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

287
cpp/mwmatching.h
View File

@ -148,6 +148,42 @@ void check_input_graph(const std::vector<Edge<WeightType>>& edges)
} }
/* **************************************************
* ** struct Graph **
* ************************************************** */
/**
* Representation of the input graph.
*/
template <typename WeightType>
struct Graph
{
using EdgeT = Edge<WeightType>;
/** Graph edges. */
const std::vector<EdgeT> edges;
/** Number of vertices. */
const VertexId num_vertex;
/** For each vertex, a vector of pointers to its incident edges. */
const std::vector<std::vector<const EdgeT*>> adjacent_edges;
/**
* Initialize the graph representation and prepare adjacent edge lists.
*
* This function takes time O(n + m).
*/
explicit Graph(const std::vector<EdgeT>& edges_in)
: edges(remove_negative_weight_edges(edges_in)),
num_vertex(count_num_vertex(edges)),
adjacent_edges(build_adjacent_edges(edges, num_vertex))
{ }
// Prevent copying.
Graph(const Graph&) = delete;
Graph& operator=(const Graph&) = delete;
/** /**
* Copy edges and remove edges with negative weight. * Copy edges and remove edges with negative weight.
* *
@ -156,14 +192,13 @@ void check_input_graph(const std::vector<Edge<WeightType>>& edges)
* @param edges Vector of weighted edges defining the graph. * @param edges Vector of weighted edges defining the graph.
* @return Vector of edges with non-negative weight. * @return Vector of edges with non-negative weight.
*/ */
template <typename WeightType> static std::vector<EdgeT> remove_negative_weight_edges(
std::vector<Edge<WeightType>> remove_negative_weight_edges( const std::vector<EdgeT>& edges_in)
const std::vector<Edge<WeightType>>& edges)
{ {
std::vector<Edge<WeightType>> real_edges; std::vector<EdgeT> real_edges;
real_edges.reserve(edges.size()); real_edges.reserve(edges_in.size());
for (const Edge<WeightType>& edge : edges) { for (const Edge<WeightType>& edge : edges_in) {
if (edge.weight >= 0) { if (edge.weight >= 0) {
real_edges.push_back(edge); real_edges.push_back(edge);
} }
@ -172,6 +207,59 @@ std::vector<Edge<WeightType>> remove_negative_weight_edges(
return real_edges; return real_edges;
} }
/** Count the number of vertices in the graph. */
static VertexId count_num_vertex(const std::vector<EdgeT>& edges)
{
VertexId num_vertex = 0;
for (const Edge<WeightType>& edge : edges) {
VertexId m = std::max(edge.vt.first, edge.vt.second);
assert(m < std::numeric_limits<VertexId>::max());
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<std::vector<const EdgeT*>> build_adjacent_edges(
const std::vector<EdgeT>& edges,
VertexId num_vertex)
{
// Count the number of incident edges per vertex.
std::vector<VertexId> 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<std::vector<const EdgeT*>> 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;
}
};
/* ************************************************** /* **************************************************
* ** struct Blossom ** * ** struct Blossom **
@ -410,8 +498,10 @@ struct AlternatingPath
* auxiliary data structures. * auxiliary data structures.
*/ */
template <typename WeightType> template <typename WeightType>
struct MatchingContext class MatchingContext
{ {
public:
using EdgeT = Edge<WeightType>; using EdgeT = Edge<WeightType>;
using BlossomT = Blossom<WeightType>; using BlossomT = Blossom<WeightType>;
using NonTrivialBlossomT = NonTrivialBlossom<WeightType>; using NonTrivialBlossomT = NonTrivialBlossom<WeightType>;
@ -436,14 +526,8 @@ struct MatchingContext
static constexpr WeightType weight_factor = static constexpr WeightType weight_factor =
std::numeric_limits<WeightType>::is_integer ? 2 : 1; std::numeric_limits<WeightType>::is_integer ? 2 : 1;
/** Graph edges. */ /** Input graph. */
const std::vector<EdgeT> edges; const Graph<WeightType> graph;
/** Number of vertices. */
const VertexId num_vertex;
/** For each vertex, a vector of pointers to its incident edges. */
const std::vector<std::vector<const EdgeT*>> adjacent_edges;
/** /**
* Current state of the matching. * Current state of the matching.
@ -506,99 +590,43 @@ struct MatchingContext
* This function takes time O(n + m). * This function takes time O(n + m).
*/ */
explicit MatchingContext(const std::vector<EdgeT>& edges_in) explicit MatchingContext(const std::vector<EdgeT>& edges_in)
: edges(remove_negative_weight_edges(edges_in)), : graph(edges_in)
num_vertex(count_num_vertex(edges)),
adjacent_edges(build_adjacent_edges(edges, num_vertex))
{ {
// Initially all vertices are unmatched. // Initially all vertices are unmatched.
vertex_mate.resize(num_vertex, NO_VERTEX); vertex_mate.resize(graph.num_vertex, NO_VERTEX);
// Create a trivial blossom for each vertex. // Create a trivial blossom for each vertex.
trivial_blossom.reserve(num_vertex); trivial_blossom.reserve(graph.num_vertex);
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
trivial_blossom.emplace_back(x); trivial_blossom.emplace_back(x);
} }
// Initially all vertices are trivial top-level blossoms. // Initially all vertices are trivial top-level blossoms.
vertex_top_blossom.reserve(num_vertex); vertex_top_blossom.reserve(graph.num_vertex);
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
vertex_top_blossom.push_back(&trivial_blossom[x]); vertex_top_blossom.push_back(&trivial_blossom[x]);
} }
// Vertex duals are initialized to half the maximum edge weight. // Vertex duals are initialized to half the maximum edge weight.
WeightType max_weight = 0; WeightType max_weight = 0;
for (const EdgeT& edge : edges) { for (const EdgeT& edge : graph.edges) {
max_weight = std::max(max_weight, edge.weight); max_weight = std::max(max_weight, edge.weight);
} }
WeightType init_vertex_dual = max_weight * (weight_factor / 2); WeightType init_vertex_dual = max_weight * (weight_factor / 2);
vertex_dual.resize(num_vertex, init_vertex_dual); vertex_dual.resize(graph.num_vertex, init_vertex_dual);
// Initialize "vertex_best_edge". // Initialize "vertex_best_edge".
vertex_best_edge.resize(num_vertex, nullptr); vertex_best_edge.resize(graph.num_vertex, nullptr);
// Allocate temporary arrays for path tracing. // Allocate temporary arrays for path tracing.
vertex_marker.resize(num_vertex); vertex_marker.resize(graph.num_vertex);
marked_vertex.reserve(num_vertex); marked_vertex.reserve(graph.num_vertex);
} }
// Prevent copying. // Prevent copying.
MatchingContext(const MatchingContext&) = delete; MatchingContext(const MatchingContext&) = delete;
MatchingContext& operator=(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<EdgeT>& edges)
{
VertexId num_vertex = 0;
for (const Edge<WeightType>& edge : edges) {
VertexId m = std::max(edge.vt.first, edge.vt.second);
assert(m < std::numeric_limits<VertexId>::max());
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<std::vector<const EdgeT*>> build_adjacent_edges(
const std::vector<EdgeT>& edges,
VertexId num_vertex)
{
// Count the number of incident edges per vertex.
std::vector<VertexId> 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<std::vector<const EdgeT*>> 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: ********** */ /* ********** General support routines: ********** */
/** Calculate edge slack. */ /** Calculate edge slack. */
@ -648,7 +676,7 @@ struct MatchingContext
*/ */
void lset_reset() void lset_reset()
{ {
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
vertex_best_edge[x] = nullptr; vertex_best_edge[x] = nullptr;
} }
@ -692,7 +720,7 @@ struct MatchingContext
const EdgeT* best_edge = nullptr; const EdgeT* best_edge = nullptr;
WeightType best_slack = 0; WeightType best_slack = 0;
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
if (vertex_top_blossom[x]->label == LABEL_NONE) { if (vertex_top_blossom[x]->label == LABEL_NONE) {
const EdgeT* edge = vertex_best_edge[x]; const EdgeT* edge = vertex_best_edge[x];
if (edge != nullptr) { if (edge != nullptr) {
@ -762,7 +790,7 @@ struct MatchingContext
} else { } else {
// Trivial blossoms don't maintain a least-slack edge set. // Trivial blossoms don't maintain a least-slack edge set.
// Just consider all incident edges. // Just consider all incident edges.
for (const EdgeT* edge : adjacent_edges[sub->base_vertex]) { for (const EdgeT* edge : graph.adjacent_edges[sub->base_vertex]) {
// Only take edges between different S-blossoms. // Only take edges between different S-blossoms.
VertexId x = edge->vt.first; VertexId x = edge->vt.first;
VertexId y = edge->vt.second; VertexId y = edge->vt.second;
@ -782,7 +810,7 @@ struct MatchingContext
// each top-level S-blossom. This array is indexed by the base vertex // each top-level S-blossom. This array is indexed by the base vertex
// of the blossoms. // of the blossoms.
std::vector<std::pair<const EdgeT*, WeightType>> std::vector<std::pair<const EdgeT*, WeightType>>
best_edge_to_blossom(num_vertex); best_edge_to_blossom(graph.num_vertex);
for (const EdgeT* edge : blossom->best_edge_set) { for (const EdgeT* edge : blossom->best_edge_set) {
BlossomT* bx = vertex_top_blossom[edge->vt.first]; BlossomT* bx = vertex_top_blossom[edge->vt.first];
@ -1467,7 +1495,7 @@ struct MatchingContext
// Scan the edges that are incident on "x". // Scan the edges that are incident on "x".
// This loop runs through O(m) iterations per stage. // This loop runs through O(m) iterations per stage.
for (const EdgeT* edge : adjacent_edges[x]) { for (const EdgeT* edge : graph.adjacent_edges[x]) {
VertexId y = (edge->vt.first != x) ? edge->vt.first VertexId y = (edge->vt.first != x) ? edge->vt.first
: edge->vt.second; : edge->vt.second;
@ -1540,7 +1568,7 @@ struct MatchingContext
// Compute delta1: minimum dual variable of any S-vertex. // Compute delta1: minimum dual variable of any S-vertex.
delta.kind = 1; delta.kind = 1;
delta.value = std::numeric_limits<WeightType>::max(); delta.value = std::numeric_limits<WeightType>::max();
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
if (vertex_top_blossom[x]->label == LABEL_S) { if (vertex_top_blossom[x]->label == LABEL_S) {
delta.value = std::min(delta.value, vertex_dual[x]); delta.value = std::min(delta.value, vertex_dual[x]);
} }
@ -1584,7 +1612,7 @@ struct MatchingContext
void substage_apply_delta_step(WeightType delta) void substage_apply_delta_step(WeightType delta)
{ {
// Apply delta to dual variables of all vertices. // Apply delta to dual variables of all vertices.
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
BlossomLabel xlabel = vertex_top_blossom[x]->label; BlossomLabel xlabel = vertex_top_blossom[x]->label;
if (xlabel == LABEL_S) { if (xlabel == LABEL_S) {
// S-vertex: subtract delta from dual variable. // S-vertex: subtract delta from dual variable.
@ -1650,7 +1678,7 @@ struct MatchingContext
bool run_stage() bool run_stage()
{ {
// Assign label S to all unmatched vertices and put them in the queue. // Assign label S to all unmatched vertices and put them in the queue.
for (VertexId x = 0; x < num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
if (vertex_mate[x] == NO_VERTEX) { if (vertex_mate[x] == NO_VERTEX) {
assign_label_s(x); assign_label_s(x);
} }
@ -1747,24 +1775,36 @@ struct MatchingContext
/** Helper class to verify that an optimal solution has been found. */ /** Helper class to verify that an optimal solution has been found. */
template <typename WeightType> template <typename WeightType>
struct MatchingVerifier class MatchingVerifier
{ {
public:
MatchingVerifier(const MatchingContext<WeightType>& ctx)
: ctx(ctx),
graph(ctx.graph),
edge_duals(ctx.graph.edges.size())
{ }
/**
* 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()
{
return (verify_vertex_mate()
&& verify_vertex_duals()
&& verify_blossom_duals()
&& verify_blossoms_and_calc_edge_duals()
&& verify_edge_slack());
}
private:
using EdgeT = Edge<WeightType>; using EdgeT = Edge<WeightType>;
using BlossomT = Blossom<WeightType>; using BlossomT = Blossom<WeightType>;
using NonTrivialBlossomT = NonTrivialBlossom<WeightType>; using NonTrivialBlossomT = NonTrivialBlossom<WeightType>;
// Reference to the MatchingContext instance.
const MatchingContext<WeightType>& ctx;
// For each edge, the sum of duals of its incident vertices
// and duals of all blossoms that contain the edge.
std::vector<WeightType> edge_duals;
MatchingVerifier(const MatchingContext<WeightType>& ctx)
: ctx(ctx),
edge_duals(ctx.edges.size())
{ }
static bool checked_add(WeightType& result, WeightType a, WeightType b) static bool checked_add(WeightType& result, WeightType a, WeightType b)
{ {
if (a > std::numeric_limits<WeightType>::max() - b) { if (a > std::numeric_limits<WeightType>::max() - b) {
@ -1778,7 +1818,7 @@ struct MatchingVerifier
/** Convert edge pointer to its index in the vector "edges". */ /** Convert edge pointer to its index in the vector "edges". */
std::size_t edge_index(const EdgeT* edge) std::size_t edge_index(const EdgeT* edge)
{ {
return edge - ctx.edges.data(); return edge - graph.edges.data();
} }
/** Check that the array "vertex_mate" is consistent. */ /** Check that the array "vertex_mate" is consistent. */
@ -1786,7 +1826,7 @@ struct MatchingVerifier
{ {
// Count matched vertices and check symmetry of "vertex_mate". // Count matched vertices and check symmetry of "vertex_mate".
VertexId num_matched_vertex = 0; VertexId num_matched_vertex = 0;
for (VertexId x = 0; x < ctx.num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
VertexId y = ctx.vertex_mate[x]; VertexId y = ctx.vertex_mate[x];
if (y != NO_VERTEX) { if (y != NO_VERTEX) {
++num_matched_vertex; ++num_matched_vertex;
@ -1798,7 +1838,7 @@ struct MatchingVerifier
// Count matched edges. // Count matched edges.
VertexId num_matched_edge = 0; VertexId num_matched_edge = 0;
for (const EdgeT& edge : ctx.edges) { for (const EdgeT& edge : graph.edges) {
if (ctx.vertex_mate[edge.vt.first] == edge.vt.second) { if (ctx.vertex_mate[edge.vt.first] == edge.vt.second) {
++num_matched_edge; ++num_matched_edge;
} }
@ -1814,7 +1854,7 @@ struct MatchingVerifier
*/ */
bool verify_vertex_duals() bool verify_vertex_duals()
{ {
for (VertexId x = 0; x < ctx.num_vertex; ++x) { for (VertexId x = 0; x < graph.num_vertex; ++x) {
if (ctx.vertex_dual[x] < 0) { if (ctx.vertex_dual[x] < 0) {
return false; return false;
} }
@ -1859,7 +1899,7 @@ struct MatchingVerifier
// For each vertex "x", // For each vertex "x",
// "vertex_depth[x]" is the depth of the smallest blossom on // "vertex_depth[x]" is the depth of the smallest blossom on
// the current descent path that contains "x". // the current descent path that contains "x".
std::vector<VertexId> vertex_depth(ctx.num_vertex); std::vector<VertexId> vertex_depth(graph.num_vertex);
// At each depth, keep track of the sum of blossom duals // At each depth, keep track of the sum of blossom duals
// along the current descent path. // along the current descent path.
@ -1922,7 +1962,7 @@ struct MatchingVerifier
// For each incident edge, find the smallest blossom // For each incident edge, find the smallest blossom
// that contains it. // that contains it.
VertexId x = sub->base_vertex; VertexId x = sub->base_vertex;
for (const EdgeT* edge : ctx.adjacent_edges[x]) { for (const EdgeT* edge : graph.adjacent_edges[x]) {
// Only consider edges pointing out from "x". // Only consider edges pointing out from "x".
if (edge->vt.first == x) { if (edge->vt.first == x) {
VertexId y = edge->vt.second; VertexId y = edge->vt.second;
@ -1991,7 +2031,7 @@ struct MatchingVerifier
bool verify_blossoms_and_calc_edge_duals() bool verify_blossoms_and_calc_edge_duals()
{ {
// For each edge, calculate the sum of its vertex duals. // For each edge, calculate the sum of its vertex duals.
for (const EdgeT& edge : ctx.edges) { for (const EdgeT& edge : graph.edges) {
if (checked_add(edge_duals[edge_index(&edge)], if (checked_add(edge_duals[edge_index(&edge)],
ctx.vertex_dual[edge.vt.first], ctx.vertex_dual[edge.vt.first],
ctx.vertex_dual[edge.vt.second])) { ctx.vertex_dual[edge.vt.second])) {
@ -2022,7 +2062,7 @@ struct MatchingVerifier
*/ */
bool verify_edge_slack() bool verify_edge_slack()
{ {
for (const EdgeT& edge : ctx.edges) { for (const EdgeT& edge : graph.edges) {
WeightType duals = edge_duals[edge_index(&edge)]; WeightType duals = edge_duals[edge_index(&edge)];
WeightType weight = ctx.weight_factor * edge.weight; WeightType weight = ctx.weight_factor * edge.weight;
@ -2040,21 +2080,18 @@ struct MatchingVerifier
return true; return true;
} }
private:
/** Reference to the MatchingContext instance. */
const MatchingContext<WeightType>& ctx;
/** Reference to the input graph. */
const Graph<WeightType>& graph;
/** /**
* Verify that the optimum solution has been found. * For each edge, the sum of duals of its incident vertices
* * and duals of all blossoms that contain the edge.
* This function takes time O(n**2).
*
* @return True if the solution is optimal; otherwise false.
*/ */
bool verify() std::vector<WeightType> edge_duals;
{
return (verify_vertex_mate()
&& verify_vertex_duals()
&& verify_blossom_duals()
&& verify_blossoms_and_calc_edge_duals()
&& verify_edge_slack());
}
}; };
@ -2114,7 +2151,7 @@ std::vector<VertexPair> maximum_weight_matching(
// Extract the matched edges. // Extract the matched edges.
std::vector<VertexPair> solution; std::vector<VertexPair> solution;
for (const Edge<WeightType>& edge : matching.edges) { for (const Edge<WeightType>& edge : matching.graph.edges) {
if (matching.vertex_mate[edge.vt.first] == edge.vt.second) { if (matching.vertex_mate[edge.vt.first] == edge.vt.second) {
solution.push_back(edge.vt); solution.push_back(edge.vt);
} }