From 8dce659bde812db966122a3eb0dd72d5f0c50e04 Mon Sep 17 00:00:00 2001 From: Fabian Murariu <2404621+fabianmurariu@users.noreply.github.com> Date: Tue, 12 Nov 2024 09:25:54 +0000 Subject: [PATCH 1/2] remove snmalloc (#1856) --- Cargo.lock | 28 ---------------------------- Cargo.toml | 1 - raphtory/Cargo.toml | 3 --- raphtory/src/lib.rs | 7 ------- 4 files changed, 39 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index a2fe9477e..df9ba1a1b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1117,15 +1117,6 @@ version = "0.7.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1462739cb27611015575c0c11df5df7601141071f07518d56fcc1be504cbec97" -[[package]] -name = "cmake" -version = "0.1.51" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fb1e43aa7fd152b1f968787f7dbcdeb306d1867ff373c69955211876c053f91a" -dependencies = [ - "cc", -] - [[package]] name = "colorchoice" version = "1.0.3" @@ -4875,7 +4866,6 @@ dependencies = [ "rustc-hash 2.0.0", "serde", "serde_json", - "snmalloc-rs", "sorted_vector_map", "streaming-stats", "tantivy", @@ -5758,24 +5748,6 @@ version = "1.1.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1b6b67fb9a61334225b5b790716f609cd58395f895b3fe8b328786812a40bc3b" -[[package]] -name = "snmalloc-rs" -version = "0.3.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2504c9edd7ca7a1cfe637296dc0d263ce1e9975c4ec43f3652616ebce9d1df1c" -dependencies = [ - "snmalloc-sys", -] - -[[package]] -name = "snmalloc-sys" -version = "0.3.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8d448599db5c3263b35d67ab26a2399e74ca0265211f5f5dd4cb9f4c3ccada6a" -dependencies = [ - "cmake", -] - [[package]] name = "socket2" version = "0.5.7" diff --git a/Cargo.toml b/Cargo.toml index b48edd4c0..0cb4fff4e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,7 +47,6 @@ async-graphql = { version = ">=7.0.5, <7.0.8", features = [ "dynamic-schema", ] } # 7.0.8+ is borked, see https://github.com/async-graphql/async-graphql/issues/1586 bincode = "1.3.3" -snmalloc-rs = { version = "0.3.6" } async-graphql-poem = "7.0.5" dynamic-graphql = "0.9.0" reqwest = { version = "0.12.8", default-features = false, features = [ diff --git a/raphtory/Cargo.toml b/raphtory/Cargo.toml index 727b0a6fa..b21b95e78 100644 --- a/raphtory/Cargo.toml +++ b/raphtory/Cargo.toml @@ -78,9 +78,6 @@ pometry-storage = { workspace = true, optional = true } prost = { workspace = true, optional = true } prost-types = { workspace = true, optional = true } -[target.'cfg(target_os = "macos")'.dependencies] -snmalloc-rs = { workspace = true } - [dev-dependencies] csv = { workspace = true } pretty_assertions = { workspace = true } diff --git a/raphtory/src/lib.rs b/raphtory/src/lib.rs index f5e1ff2aa..860c952ed 100644 --- a/raphtory/src/lib.rs +++ b/raphtory/src/lib.rs @@ -88,13 +88,6 @@ pub mod core; pub mod db; pub mod graphgen; -#[cfg(target_os = "macos")] -use snmalloc_rs; - -#[cfg(target_os = "macos")] -#[global_allocator] -static ALLOC: snmalloc_rs::SnMalloc = snmalloc_rs::SnMalloc; - #[cfg(feature = "storage")] pub mod disk_graph; From b6f020caca98fa243f81a6bce2f04163996e789b Mon Sep 17 00:00:00 2001 From: Ben Steer Date: Wed, 13 Nov 2024 14:26:36 +0000 Subject: [PATCH 2/2] max weight matching (#1602) * Added max_weight_matching * Made weight optional * fmt * mostly working Mapping output * add python wrappers for Matching and improve the Matching output * cleanup, add tests and docs * fix non-deterministic order --------- Co-authored-by: Lucas Jeub --- python/python/raphtory/__init__.pyi | 8 +- .../python/raphtory/algorithms/__init__.pyi | 112 ++ python/python/raphtory/graphql/__init__.pyi | 20 +- python/python/raphtory/vectors/__init__.pyi | 2 +- python/scripts/gen-stubs.py | 2 + python/tests/test_algorithms.py | 27 + raphtory-api/src/core/entities/mod.rs | 26 + .../bipartite/max_weight_matching.rs | 1588 +++++++++++++++++ raphtory/src/algorithms/bipartite/mod.rs | 1 + .../algorithms/community_detection/louvain.rs | 4 +- raphtory/src/algorithms/mod.rs | 1 + raphtory/src/db/graph/graph.rs | 7 +- .../python/algorithm/max_weight_matching.rs | 121 ++ raphtory/src/python/algorithm/mod.rs | 1 + raphtory/src/python/packages/algorithms.rs | 49 + raphtory/src/python/packages/base_modules.rs | 11 +- 16 files changed, 1957 insertions(+), 23 deletions(-) create mode 100644 raphtory/src/algorithms/bipartite/max_weight_matching.rs create mode 100644 raphtory/src/algorithms/bipartite/mod.rs create mode 100644 raphtory/src/python/algorithm/max_weight_matching.rs diff --git a/python/python/raphtory/__init__.pyi b/python/python/raphtory/__init__.pyi index c7badbc32..7a6ad44fc 100644 --- a/python/python/raphtory/__init__.pyi +++ b/python/python/raphtory/__init__.pyi @@ -1162,7 +1162,7 @@ class Graph(GraphView): num_shards (int, optional): The number of locks to use in the storage to allow for multithreaded updates. """ - def __new__(self, num_shards: Optional[int] = None): + def __new__(self, num_shards: Optional[int] = None) -> Graph: """Create and return a new object. See help(type) for accurate signature.""" def __reduce__(self): ... @@ -3426,7 +3426,7 @@ class Nodes(object): class PersistentGraph(GraphView): """A temporal graph that allows edges and nodes to be deleted.""" - def __new__(self): + def __new__(self) -> PersistentGraph: """Create and return a new object. See help(type) for accurate signature.""" def __reduce__(self): ... @@ -4009,7 +4009,7 @@ class Prop(object): def __ne__(self, value): """Return self!=value.""" - def __new__(self, name): + def __new__(self, name) -> Prop: """Create and return a new object. See help(type) for accurate signature.""" def any(self, values): @@ -4103,7 +4103,7 @@ class PyGraphEncoder(object): """Call self as a function.""" def __getstate__(self): ... - def __new__(self): + def __new__(self) -> PyGraphEncoder: """Create and return a new object. See help(type) for accurate signature.""" def __setstate__(self): ... diff --git a/python/python/raphtory/algorithms/__init__.pyi b/python/python/raphtory/algorithms/__init__.pyi index 272bb6c48..880a5c91e 100644 --- a/python/python/raphtory/algorithms/__init__.pyi +++ b/python/python/raphtory/algorithms/__init__.pyi @@ -15,6 +15,78 @@ from raphtory.typing import * from datetime import datetime from pandas import DataFrame +class Matching(object): + """A Matching (i.e., a set of edges that do not share any nodes)""" + + def __bool__(self): + """True if self else False""" + + def __contains__(self, key): + """Return bool(key in self).""" + + def __iter__(self): + """Implement iter(self).""" + + def __len__(self): + """Return len(self).""" + + def __repr__(self): + """Return repr(self).""" + + def dst(self, src: InputNode) -> Optional[Node]: + """ + Get the matched destination node for a source node + + Arguments: + src (InputNode): The source node + + Returns: + Optional[Node]: The matched destination node if it exists + + """ + + def edge_for_dst(self, dst: InputNode) -> Optional[Edge]: + """ + Get the matched edge for a destination node + + Arguments: + dst (InputNode): The source node + + Returns: + Optional[Edge]: The matched edge if it exists + """ + + def edge_for_src(self, src: InputNode) -> Optional[Edge]: + """ + Get the matched edge for a source node + + Arguments: + src (InputNode): The source node + + Returns: + Optional[Edge]: The matched edge if it exists + """ + + def edges(self) -> Edges: + """ + Get a view of the matched edges + + Returns: + Edges: The edges in the matching + """ + + def src(self, dst: InputNode) -> Optional[Node]: + """ + Get the matched source node for a destination node + + Arguments: + dst (InputNode): The destination node + + Returns: + Optional[Node]: The matched source node if it exists + + """ + def all_local_reciprocity(g: GraphView): """ Local reciprocity - measure of the symmetry of relationships associated with a node @@ -405,6 +477,46 @@ def max_out_degree(g: GraphView): int : value of the largest outdegree """ +def max_weight_matching( + graph: GraphView, + weight_prop: Optional[str] = None, + max_cardinality: bool = True, + verify_optimum_flag: bool = False, +) -> Matching: + """ + Compute a maximum-weighted matching in the general undirected weighted + graph given by "edges". If `max_cardinality` is true, only + maximum-cardinality matchings are considered as solutions. + + The algorithm is based on "Efficient Algorithms for Finding Maximum + Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986. + + Based on networkx implementation + + + With reference to the standalone protoype implementation from: + + + + + The function takes time O(n**3) + + Arguments: + graph (GraphView): The graph to compute the maximum weight matching for + weight_prop (str, optional): The property on the edge to use for the weight. If not + provided, + max_cardinality (bool): If set to true compute the maximum-cardinality matching + with maximum weight among all maximum-cardinality matchings. Defaults to True. + verify_optimum_flag (bool): If true prior to returning an additional routine + to verify the optimal solution was found will be run after computing + the maximum weight matching. If it's true and the found matching is not + an optimal solution this function will panic. This option should + normally be only set true during testing. Defaults to False. + + Returns: + Matching: The matching + """ + def min_degree(g: GraphView) -> int: """ Returns the smallest degree found in the graph diff --git a/python/python/raphtory/graphql/__init__.pyi b/python/python/raphtory/graphql/__init__.pyi index 95904281d..28827d487 100644 --- a/python/python/raphtory/graphql/__init__.pyi +++ b/python/python/raphtory/graphql/__init__.pyi @@ -29,7 +29,7 @@ class GraphServer(object): otlp_agent_port=None, otlp_tracing_service_name=None, config_path=None, - ): + ) -> GraphServer: """Create and return a new object. See help(type) for accurate signature.""" def run(self, port: int = 1736, timeout_ms: int = 180000): @@ -147,7 +147,7 @@ class GraphqlGraphs(object): class RaphtoryClient(object): """A client for handling GraphQL operations in the context of Raphtory.""" - def __new__(self, url): + def __new__(self, url) -> RaphtoryClient: """Create and return a new object. See help(type) for accurate signature.""" def copy_graph(self, path, new_path): @@ -268,7 +268,7 @@ class RaphtoryClient(object): """ class RemoteEdge(object): - def __new__(self, path, client, src, dst): + def __new__(self, path, client, src, dst) -> RemoteEdge: """Create and return a new object. See help(type) for accurate signature.""" def add_constant_properties( @@ -323,11 +323,13 @@ class RemoteEdge(object): """ class RemoteEdgeAddition(object): - def __new__(self, src, dst, layer=None, constant_properties=None, updates=None): + def __new__( + self, src, dst, layer=None, constant_properties=None, updates=None + ) -> RemoteEdgeAddition: """Create and return a new object. See help(type) for accurate signature.""" class RemoteGraph(object): - def __new__(self, path, client): + def __new__(self, path, client) -> RemoteGraph: """Create and return a new object. See help(type) for accurate signature.""" def add_constant_properties(self, properties: dict): @@ -456,7 +458,7 @@ class RemoteGraph(object): """ class RemoteNode(object): - def __new__(self, path, client, id): + def __new__(self, path, client, id) -> RemoteNode: """Create and return a new object. See help(type) for accurate signature.""" def add_constant_properties(self, properties: Dict[str, Prop]): @@ -501,11 +503,13 @@ class RemoteNode(object): """ class RemoteNodeAddition(object): - def __new__(self, name, node_type=None, constant_properties=None, updates=None): + def __new__( + self, name, node_type=None, constant_properties=None, updates=None + ) -> RemoteNodeAddition: """Create and return a new object. See help(type) for accurate signature.""" class RemoteUpdate(object): - def __new__(self, time, properties=None): + def __new__(self, time, properties=None) -> RemoteUpdate: """Create and return a new object. See help(type) for accurate signature.""" class RunningGraphServer(object): diff --git a/python/python/raphtory/vectors/__init__.pyi b/python/python/raphtory/vectors/__init__.pyi index 255046f0c..efc95d211 100644 --- a/python/python/raphtory/vectors/__init__.pyi +++ b/python/python/raphtory/vectors/__init__.pyi @@ -16,7 +16,7 @@ from datetime import datetime from pandas import DataFrame class Document(object): - def __new__(self, content, life=None): + def __new__(self, content, life=None) -> Document: """Create and return a new object. See help(type) for accurate signature.""" def __repr__(self): diff --git a/python/scripts/gen-stubs.py b/python/scripts/gen-stubs.py index 8b5e6d2ae..336535dcb 100755 --- a/python/scripts/gen-stubs.py +++ b/python/scripts/gen-stubs.py @@ -279,6 +279,8 @@ def gen_class(cls: type, name) -> str: # Get __init__ signature from class info signature = cls_signature(cls) if signature is not None: + if obj_name == "__new__": + signature = signature.replace(return_annotation=name) entities.append( gen_fn( entity, diff --git a/python/tests/test_algorithms.py b/python/tests/test_algorithms.py index aec08c906..47112b169 100644 --- a/python/tests/test_algorithms.py +++ b/python/tests/test_algorithms.py @@ -1,6 +1,7 @@ import pytest import pandas as pd import pandas.core.frame + from raphtory import Graph, PersistentGraph from raphtory import algorithms from raphtory import graph_loader @@ -518,3 +519,29 @@ def test_temporal_SEIR(): for i, (n, v) in enumerate(res): assert n == g.node(i + 1) assert v.infected == i + + +def test_max_weight_matching(): + g = Graph() + g.add_edge(0, 1, 2, {"weight": 5}) + g.add_edge(0, 2, 3, {"weight": 11}) + g.add_edge(0, 3, 4, {"weight": 5}) + + # Run max weight matching with max cardinality set to false + max_weight = algorithms.max_weight_matching(g, "weight", False) + max_cardinality = algorithms.max_weight_matching(g, "weight") + + assert len(max_weight) == 1 + assert len(max_cardinality) == 2 + assert (2, 3) in max_weight + assert (1, 2) in max_cardinality + assert (3, 4) in max_cardinality + + assert max_weight.edges().id == [(2, 3)] + assert sorted(max_cardinality.edges().id) == [(1, 2), (3, 4)] + + assert max_weight.src(3).id == 2 + assert max_weight.src(2) is None + + assert max_weight.dst(2).id == 3 + assert max_weight.dst(3) is None diff --git a/raphtory-api/src/core/entities/mod.rs b/raphtory-api/src/core/entities/mod.rs index b334b31d6..068a8320a 100644 --- a/raphtory-api/src/core/entities/mod.rs +++ b/raphtory-api/src/core/entities/mod.rs @@ -86,6 +86,32 @@ pub enum GID { U64(u64), Str(String), } +impl PartialEq for GID { + fn eq(&self, other: &str) -> bool { + match self { + GID::U64(_) => false, + GID::Str(id) => id == other, + } + } +} + +impl PartialEq for GID { + fn eq(&self, other: &String) -> bool { + match self { + GID::U64(_) => false, + GID::Str(id) => id == other, + } + } +} + +impl PartialEq for GID { + fn eq(&self, other: &u64) -> bool { + match self { + GID::Str(_) => false, + GID::U64(id) => id == other, + } + } +} impl Default for GID { fn default() -> Self { diff --git a/raphtory/src/algorithms/bipartite/max_weight_matching.rs b/raphtory/src/algorithms/bipartite/max_weight_matching.rs new file mode 100644 index 000000000..0e736f689 --- /dev/null +++ b/raphtory/src/algorithms/bipartite/max_weight_matching.rs @@ -0,0 +1,1588 @@ +// Licensed under the Apache License, Version 2.0 (the "License"); you may +// not use this file except in compliance with the License. You may obtain +// a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +// License for the specific language governing permissions and limitations +// under the License. + +// Needed to pass shared state between functions +// closures don't work because of recurssion +#![allow(clippy::too_many_arguments)] +// Allow single character names to match naming convention from +// paper +#![allow(clippy::many_single_char_names)] + +use crate::{ + core::{entities::nodes::node_ref::AsNodeRef, utils::iter::GenLockedIter}, + db::{ + api::{ + storage::graph::edges::edge_storage_ops::EdgeStorageOps, + view::{DynamicGraph, IntoDynBoxed, IntoDynamic, StaticGraphViewOps}, + }, + graph::{edge::EdgeView, edges::Edges, node::NodeView}, + }, + prelude::{EdgeViewOps, GraphViewOps, Prop, PropUnwrap}, +}; +use hashbrown::HashMap; +use raphtory_api::core::entities::{EID, VID}; +use std::{ + cmp::max, + fmt::{Debug, Display, Formatter}, + mem, + sync::Arc, +}; + +/// Return 2 * slack of edge k (does not work inside blossoms). +fn slack(edge_index: usize, dual_var: &[i64], edges: &[(usize, usize, i64)]) -> i64 { + let (source_index, target_index, weight) = edges[edge_index]; + dual_var[source_index] + dual_var[target_index] - 2 * weight +} + +/// Generate the leaf vertices of a blossom. +fn blossom_leaves(blossom: usize, num_nodes: usize, blossom_children: &[Vec]) -> Vec { + let mut out_vec: Vec = Vec::new(); + if blossom < num_nodes { + out_vec.push(blossom); + } else { + let child_blossom = &blossom_children[blossom]; + for c in child_blossom { + let child = *c; + if child < num_nodes { + out_vec.push(child); + } else { + for v in blossom_leaves(child, num_nodes, blossom_children) { + out_vec.push(v); + } + } + } + } + out_vec +} + +/// Assign label t to the top-level blossom containing vertex w +/// and record the fact that w was reached through the edge with +/// remote endpoint p. +fn assign_label( + w: usize, + t: usize, + p: Option, + num_nodes: usize, + in_blossoms: &[usize], + labels: &mut Vec>, + label_ends: &mut Vec>, + best_edge: &mut Vec>, + queue: &mut Vec, + blossom_children: &[Vec], + blossom_base: &[Option], + endpoints: &[usize], + mate: &HashMap, +) -> () { + let b = in_blossoms[w]; + assert!(labels[w] == Some(0) && labels[b] == Some(0)); + labels[w] = Some(t); + labels[b] = Some(t); + label_ends[b] = p; + label_ends[w] = p; + best_edge[w] = None; + best_edge[b] = None; + if t == 1 { + // b became an S-vertex/blossom; add it(s verticies) to the queue + queue.append(&mut blossom_leaves(b, num_nodes, blossom_children)); + } else if t == 2 { + // b became a T-vertex/blossom; assign label S to its mate. + // (If b is a non-trivial blossom, its base is the only vertex + // with an external mate.) + let blossom_index: usize = b; + let base: usize = blossom_base[blossom_index].unwrap(); + assert!(mate.get(&base).is_some()); + assign_label( + endpoints[mate[&base]], + 1, + mate.get(&base).map(|p| p ^ 1), + num_nodes, + in_blossoms, + labels, + label_ends, + best_edge, + queue, + blossom_children, + blossom_base, + endpoints, + mate, + ); + } + () +} + +/// Trace back from vertices v and w to discover either a new blossom +/// or an augmenting path. Return the base vertex of the new blossom or None. +fn scan_blossom( + node_a: usize, + node_b: usize, + in_blossoms: &[usize], + blossom_base: &[Option], + endpoints: &[usize], + labels: &mut [Option], + label_ends: &[Option], + mate: &HashMap, +) -> Option { + let mut v: Option = Some(node_a); + let mut w: Option = Some(node_b); + // Trace back from v and w, placing breadcrumbs as we go + let mut path: Vec = Vec::new(); + let mut base: Option = None; + while v.is_some() || w.is_some() { + // Look for a breadcrumb in v's blossom or put a new breadcrump + let mut blossom = in_blossoms[v.unwrap()]; + if labels[blossom].is_none() || labels[blossom].unwrap() & 4 > 0 { + base = blossom_base[blossom]; + break; + } + assert_eq!(labels[blossom], Some(1)); + path.push(blossom); + labels[blossom] = Some(5); + // Trace one step bacl. + assert_eq!( + label_ends[blossom], + mate.get(&blossom_base[blossom].unwrap()).copied() + ); + if label_ends[blossom].is_none() { + // The base of blossom is single; stop tracing this path + v = None; + } else { + let tmp = endpoints[label_ends[blossom].unwrap()]; + blossom = in_blossoms[tmp]; + assert_eq!(labels[blossom], Some(2)); + // blossom is a T-blossom; trace one more step back. + assert!(label_ends[blossom].is_some()); + v = Some(endpoints[label_ends[blossom].unwrap()]); + } + // Swap v and w so that we alternate between both paths. + if w.is_some() { + mem::swap(&mut v, &mut w); + } + } + // Remvoe breadcrumbs. + for blossom in path { + labels[blossom] = Some(1); + } + // Return base vertex, if we found one. + base +} + +/// Construct a new blossom with given base, containing edge k which +/// connects a pair of S vertices. Label the new blossom as S; set its dual +/// variable to zero; relabel its T-vertices to S and add them to the queue. +fn add_blossom( + base: usize, + edge: usize, + blossom_children: &mut [Vec], + num_nodes: usize, + edges: &[(usize, usize, i64)], + in_blossoms: &mut [usize], + dual_var: &mut [i64], + labels: &mut [Option], + label_ends: &mut [Option], + best_edge: &mut [Option], + queue: &mut Vec, + blossom_base: &mut [Option], + endpoints: &[usize], + blossom_endpoints: &mut [Vec], + unused_blossoms: &mut Vec, + blossom_best_edges: &mut [Vec], + blossom_parents: &mut [Option], + neighbor_endpoints: &[Vec], + mate: &HashMap, +) -> () { + let (mut v, mut w, _weight) = edges[edge]; + let blossom_b = in_blossoms[base]; + let mut blossom_v = in_blossoms[v]; + let mut blossom_w = in_blossoms[w]; + // Create blossom + let blossom = unused_blossoms.pop().unwrap(); + blossom_base[blossom] = Some(base); + blossom_parents[blossom] = None; + blossom_parents[blossom_b] = Some(blossom); + // Make list of sub-blossoms and their interconnecting edge endpoints. + blossom_children[blossom].clear(); + blossom_endpoints[blossom].clear(); + // Trace back from blossom_v to base. + while blossom_v != blossom_b { + // Add blossom_v to the new blossom + blossom_parents[blossom_v] = Some(blossom); + blossom_children[blossom].push(blossom_v); + let blossom_v_endpoint_label = label_ends[blossom_v].unwrap(); + blossom_endpoints[blossom].push(blossom_v_endpoint_label); + assert!( + labels[blossom_v] == Some(2) + || (labels[blossom_v] == Some(1) + && label_ends[blossom_v] + == mate.get(&blossom_base[blossom_v].unwrap()).copied()) + ); + // Trace one step back. + assert!(label_ends[blossom_v].is_some()); + v = endpoints[blossom_v_endpoint_label]; + blossom_v = in_blossoms[v]; + } + // Reverse lists, add endpoint that connects the pair of S vertices. + blossom_children[blossom].push(blossom_b); + blossom_children[blossom].reverse(); + blossom_endpoints[blossom].reverse(); + blossom_endpoints[blossom].push(2 * edge); + // Trace back from w to base. + while blossom_w != blossom_b { + // Add blossom_w to the new blossom + blossom_parents[blossom_w] = Some(blossom); + blossom_children[blossom].push(blossom_w); + let blossom_w_endpoint_label = label_ends[blossom_w].unwrap(); + blossom_endpoints[blossom].push(blossom_w_endpoint_label ^ 1); + assert!( + labels[blossom_w] == Some(2) + || (labels[blossom_w] == Some(1) + && label_ends[blossom_w] + == mate.get(&blossom_base[blossom_w].unwrap()).copied()) + ); + // Trace one step back + assert!(label_ends[blossom_w].is_some()); + w = endpoints[blossom_w_endpoint_label]; + blossom_w = in_blossoms[w]; + } + // Set label to S. + assert_eq!(labels[blossom_b], Some(1)); + labels[blossom] = Some(1); + label_ends[blossom] = label_ends[blossom_b]; + // Set dual variable to 0 + dual_var[blossom] = 0; + // Relabel vertices + for node in blossom_leaves(blossom, num_nodes, blossom_children) { + if labels[in_blossoms[node]] == Some(2) { + // This T-vertex now turns into an S-vertex because it becomes + // part of an S-blossom; add it to the queue + queue.push(node); + } + in_blossoms[node] = blossom; + } + // Compute blossom_best_edges[blossom] + let mut best_edge_to: HashMap = HashMap::with_capacity(2 * num_nodes); + for bv_ref in &blossom_children[blossom] { + let bv = *bv_ref; + // This sub-blossom does not have a list of least-slack edges; + // get the information from the vertices. + let nblists: Vec> = if blossom_best_edges[bv].is_empty() { + let mut tmp: Vec> = Vec::new(); + for node in blossom_leaves(bv, num_nodes, blossom_children) { + tmp.push(neighbor_endpoints[node].iter().map(|p| p / 2).collect()); + } + tmp + } else { + // Walk this sub-blossom's least-slack edges. + vec![blossom_best_edges[bv].clone()] + }; + for nblist in nblists { + for edge_index in nblist { + let (mut i, mut j, _edge_weight) = edges[edge_index]; + if in_blossoms[j] == blossom { + mem::swap(&mut i, &mut j); + } + let blossom_j = in_blossoms[j]; + if blossom_j != blossom + && labels[blossom_j] == Some(1) + && (best_edge_to.get(&blossom_j).is_none() + || slack(edge_index, dual_var, edges) + < slack(best_edge_to[&blossom_j], dual_var, edges)) + { + best_edge_to.insert(blossom_j, edge_index); + } + } + } + // Forget about least-slack edges of the sub-blossom. + blossom_best_edges[bv].clear(); + best_edge[bv] = None; + } + blossom_best_edges[blossom] = best_edge_to.values().copied().collect(); + //select best_edge[blossom] + best_edge[blossom] = None; + for edge_index in &blossom_best_edges[blossom] { + if best_edge[blossom].is_none() + || slack(*edge_index, dual_var, edges) + < slack(best_edge[blossom].unwrap(), dual_var, edges) + { + best_edge[blossom] = Some(*edge_index); + } + } + () +} + +/// Expand the given top level blossom +fn expand_blossom( + blossom: usize, + end_stage: bool, + num_nodes: usize, + blossom_children: &mut Vec>, + blossom_parents: &mut Vec>, + in_blossoms: &mut Vec, + dual_var: &[i64], + labels: &mut Vec>, + label_ends: &mut Vec>, + best_edge: &mut Vec>, + queue: &mut Vec, + blossom_base: &mut Vec>, + endpoints: &[usize], + mate: &HashMap, + blossom_endpoints: &mut Vec>, + allowed_edge: &mut Vec, + unused_blossoms: &mut Vec, +) -> () { + // Convert sub-blossoms into top-level blossoms. + for s in blossom_children[blossom].clone() { + blossom_parents[s] = None; + if s < num_nodes { + in_blossoms[s] = s + } else if end_stage && dual_var[s] == 0 { + // Recursively expand this sub-blossom + expand_blossom( + s, + end_stage, + num_nodes, + blossom_children, + blossom_parents, + in_blossoms, + dual_var, + labels, + label_ends, + best_edge, + queue, + blossom_base, + endpoints, + mate, + blossom_endpoints, + allowed_edge, + unused_blossoms, + ); + } else { + for v in blossom_leaves(s, num_nodes, blossom_children) { + in_blossoms[v] = s; + } + } + } + // if we expand a T-blossom during a stage, its a sub-blossoms must be + // relabeled + if !end_stage && labels[blossom] == Some(2) { + // start at the sub-blossom through which the expanding blossom + // obtained its label, and relabel sub-blossoms until we reach the + // base. + assert!(label_ends[blossom].is_some()); + let entry_child = in_blossoms[endpoints[label_ends[blossom].unwrap() ^ 1]]; + // Decied in which direction we will go around the blossom. + let i = blossom_children[blossom] + .iter() + .position(|x| *x == entry_child) + .unwrap(); + let mut j = i as i64; + let j_step: i64; + let endpoint_trick: usize = if i & 1 != 0 { + // Start index is odd; go forward and wrap + j -= blossom_children[blossom].len() as i64; + j_step = 1; + 0 + } else { + // start index is even; go backward + j_step = -1; + 1 + }; + // Move along the blossom until we get to the base + let mut p = label_ends[blossom].unwrap(); + while j != 0 { + // Relabel the T-sub-blossom. + labels[endpoints[p ^ 1]] = Some(0); + if j < 0 { + let length = blossom_endpoints[blossom].len(); + let index = length - j.unsigned_abs() as usize; + labels[endpoints + [blossom_endpoints[blossom][index - endpoint_trick] ^ endpoint_trick ^ 1]] = + Some(0); + } else { + labels[endpoints[blossom_endpoints[blossom][j as usize - endpoint_trick] + ^ endpoint_trick + ^ 1]] = Some(0); + } + assign_label( + endpoints[p ^ 1], + 2, + Some(p), + num_nodes, + in_blossoms, + labels, + label_ends, + best_edge, + queue, + blossom_children, + blossom_base, + endpoints, + mate, + ); + // Step to the next S-sub-blossom and note it's forwward endpoint. + let endpoint_index = if j < 0 { + let tmp = j - endpoint_trick as i64; + let length = blossom_endpoints[blossom].len(); + let index = length - tmp.unsigned_abs() as usize; + blossom_endpoints[blossom][index] + } else { + blossom_endpoints[blossom][j as usize - endpoint_trick] + }; + allowed_edge[endpoint_index / 2] = true; + j += j_step; + p = if j < 0 { + let tmp = j - endpoint_trick as i64; + let length = blossom_endpoints[blossom].len(); + let index = length - tmp.unsigned_abs() as usize; + blossom_endpoints[blossom][index] ^ endpoint_trick + } else { + blossom_endpoints[blossom][j as usize - endpoint_trick] ^ endpoint_trick + }; + // Step to the next T-sub-blossom. + allowed_edge[p / 2] = true; + j += j_step; + } + // Relabel the base T-sub-blossom WITHOUT stepping through to + // its mate (so don't call assign_label()) + let blossom_v = if j < 0 { + let length = blossom_children[blossom].len(); + let index = length - j.unsigned_abs() as usize; + blossom_children[blossom][index] + } else { + blossom_children[blossom][j as usize] + }; + labels[endpoints[p ^ 1]] = Some(2); + labels[blossom_v] = Some(2); + label_ends[endpoints[p ^ 1]] = Some(p); + label_ends[blossom_v] = Some(p); + best_edge[blossom_v] = None; + // Continue along the blossom until we get back to entry_child + j += j_step; + let mut j_index = if j < 0 { + let length = blossom_children[blossom].len(); + length - j.unsigned_abs() as usize + } else { + j as usize + }; + while blossom_children[blossom][j_index] != entry_child { + // Examine the vertices of the sub-blossom to see whether + // it is reachable from a neighboring S-vertex outside the + // expanding blososm. + let bv = blossom_children[blossom][j_index]; + if labels[bv] == Some(1) { + // This sub-blossom just got label S through one of its + // neighbors; leave it + j += j_step; + if j < 0 { + let length = blossom_children[blossom].len(); + j_index = length - j.unsigned_abs() as usize; + } else { + j_index = j as usize; + } + continue; + } + let mut v: usize = 0; + for tmp in blossom_leaves(bv, num_nodes, blossom_children) { + v = tmp; + if labels[v].unwrap() != 0 { + break; + } + } + // If the sub-blossom contains a reachable vertex, assign label T + // to the sub-blossom + if labels[v] != Some(0) { + assert_eq!(labels[v], Some(2)); + assert_eq!(in_blossoms[v], bv); + labels[v] = Some(0); + labels[endpoints[mate[&blossom_base[bv].unwrap()]]] = Some(0); + assign_label( + v, + 2, + label_ends[v], + num_nodes, + in_blossoms, + labels, + label_ends, + best_edge, + queue, + blossom_children, + blossom_base, + endpoints, + mate, + ); + } + j += j_step; + if j < 0 { + let length = blossom_children[blossom].len(); + j_index = length - j.unsigned_abs() as usize; + } else { + j_index = j as usize; + } + } + } + // Recycle the blossom number. + labels[blossom] = None; + label_ends[blossom] = None; + blossom_children[blossom].clear(); + blossom_endpoints[blossom].clear(); + blossom_base[blossom] = None; + best_edge[blossom] = None; + unused_blossoms.push(blossom); + () +} + +/// Swap matched/unmatched edges over an alternating path through blossom b +/// between vertex v and the base vertex. Keep blossom bookkeeping consistent. +fn augment_blossom( + blossom: usize, + node: usize, + num_nodes: usize, + blossom_parents: &[Option], + endpoints: &[usize], + blossom_children: &mut Vec>, + blossom_endpoints: &mut Vec>, + blossom_base: &mut Vec>, + mate: &mut HashMap, +) { + // Bubble up through the blossom tree from vertex v to an immediate + // sub-blossom of b. + let mut tmp = node; + while blossom_parents[tmp] != Some(blossom) { + tmp = blossom_parents[tmp].unwrap(); + } + // Recursively deal with the first sub-blossom. + if tmp >= num_nodes { + augment_blossom( + tmp, + node, + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Decide in which direction we will go around the blossom. + let i = blossom_children[blossom] + .iter() + .position(|x| *x == tmp) + .unwrap(); + let mut j: i64 = i as i64; + let j_step: i64; + let endpoint_trick: usize = if i & 1 != 0 { + // start index is odd; go forward and wrap. + j -= blossom_children[blossom].len() as i64; + j_step = 1; + 0 + } else { + // Start index is even; go backward. + j_step = -1; + 1 + }; + // Move along the blossom until we get to the base. + while j != 0 { + // Step to the next sub-blossom and augment it recursively. + j += j_step; + + tmp = if j < 0 { + let length = blossom_children[blossom].len(); + let index = length - j.unsigned_abs() as usize; + blossom_children[blossom][index] + } else { + blossom_children[blossom][j as usize] + }; + let p = if j < 0 { + let length = blossom_endpoints[blossom].len(); + let index = length - j.unsigned_abs() as usize - endpoint_trick; + blossom_endpoints[blossom][index] ^ endpoint_trick + } else { + blossom_endpoints[blossom][j as usize - endpoint_trick] ^ endpoint_trick + }; + if tmp > num_nodes { + augment_blossom( + tmp, + endpoints[p], + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + j += j_step; + if j < 0 { + let length = blossom_children[blossom].len(); + let index = length - j.unsigned_abs() as usize; + tmp = blossom_children[blossom][index]; + } else { + tmp = blossom_children[blossom][j as usize]; + } + if tmp > num_nodes { + augment_blossom( + tmp, + endpoints[p ^ 1], + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Match the edge connecting those sub-blossoms. + mate.insert(endpoints[p], p ^ 1); + mate.insert(endpoints[p ^ 1], p); + } + // Rotate the list of sub-blossoms to put the new base at the front. + let mut children: Vec = blossom_children[blossom][i..].to_vec(); + children.extend_from_slice(&blossom_children[blossom][..i]); + blossom_children[blossom] = children; + let mut endpoint: Vec = blossom_endpoints[blossom][i..].to_vec(); + endpoint.extend_from_slice(&blossom_endpoints[blossom][..i]); + blossom_endpoints[blossom] = endpoint; + blossom_base[blossom] = blossom_base[blossom_children[blossom][0]]; + assert_eq!(blossom_base[blossom], Some(node)); +} + +/// Swap matched/unmatched edges over an alternating path between two +/// single vertices. The augmenting path runs through edge k, which +/// connects a pair of S vertices. +fn augment_matching( + edge: usize, + num_nodes: usize, + edges: &[(usize, usize, i64)], + in_blossoms: &[usize], + labels: &[Option], + label_ends: &[Option], + blossom_parents: &[Option], + endpoints: &[usize], + blossom_children: &mut Vec>, + blossom_endpoints: &mut Vec>, + blossom_base: &mut Vec>, + mate: &mut HashMap, +) { + let (v, w, _weight) = edges[edge]; + for (s_ref, p_ref) in [(v, 2 * edge + 1), (w, 2 * edge)].iter() { + // Match vertex s to remote endpoint p. Then trace back from s + // until we find a single vertex, swapping matched and unmatched + // edges as we go. + let mut s: usize = *s_ref; + let mut p: usize = *p_ref; + loop { + let blossom_s = in_blossoms[s]; + assert_eq!(labels[blossom_s], Some(1)); + assert_eq!( + label_ends[blossom_s], + mate.get(&blossom_base[blossom_s].unwrap()).copied() + ); + // Augment through the S-blossom from s to base. + if blossom_s >= num_nodes { + augment_blossom( + blossom_s, + s, + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Update mate[s] + mate.insert(s, p); + // Trace one step back. + if label_ends[blossom_s].is_none() { + // Reached a single vertex; stop + break; + } + let t = endpoints[label_ends[blossom_s].unwrap()]; + let blossom_t = in_blossoms[t]; + assert_eq!(labels[blossom_t], Some(2)); + // Trace one step back + assert!(label_ends[blossom_t].is_some()); + s = endpoints[label_ends[blossom_t].unwrap()]; + let j = endpoints[label_ends[blossom_t].unwrap() ^ 1]; + // Augment through the T-blossom from j to base. + assert_eq!(blossom_base[blossom_t], Some(t)); + if blossom_t >= num_nodes { + augment_blossom( + blossom_t, + j, + num_nodes, + blossom_parents, + endpoints, + blossom_children, + blossom_endpoints, + blossom_base, + mate, + ); + } + // Update mate[j] + mate.insert(j, label_ends[blossom_t].unwrap()); + // Keep the opposite endpoint; + // it will be assigned to mate[s] in the next step. + p = label_ends[blossom_t].unwrap() ^ 1; + } + } +} + +/// Swap matched/unmatched edges over an alternating path between two +/// single vertices. The augmenting path runs through the edge, which +/// connects a pair of S vertices. +fn verify_optimum( + max_cardinality: bool, + num_nodes: usize, + num_edges: usize, + edges: &[(usize, usize, i64)], + endpoints: &[usize], + dual_var: &[i64], + blossom_parents: &[Option], + blossom_endpoints: &[Vec], + blossom_base: &[Option], + mate: &HashMap, +) { + let dual_var_node_min: i64 = *dual_var[..num_nodes].iter().min().unwrap(); + let node_dual_offset: i64 = if max_cardinality { + // Vertices may have negative dual; + // find a constant non-negative number to add to all vertex duals. + max(0, -dual_var_node_min) + } else { + 0 + }; + assert!(dual_var_node_min + node_dual_offset >= 0); + assert!(*dual_var[num_nodes..].iter().min().unwrap() >= 0); + // 0. all edges have non-negative slack and + // 1. all matched edges have zero slack; + for (edge, (i, j, weight)) in edges.iter().enumerate().take(num_edges) { + let mut s = dual_var[*i] + dual_var[*j] - 2 * weight; + let mut i_blossoms: Vec = vec![*i]; + let mut j_blossoms: Vec = vec![*j]; + while blossom_parents[*i_blossoms.last().unwrap()].is_some() { + i_blossoms.push(blossom_parents[*i_blossoms.last().unwrap()].unwrap()); + } + while blossom_parents[*j_blossoms.last().unwrap()].is_some() { + j_blossoms.push(blossom_parents[*j_blossoms.last().unwrap()].unwrap()); + } + i_blossoms.reverse(); + j_blossoms.reverse(); + for (blossom_i, blossom_j) in i_blossoms.iter().zip(j_blossoms.iter()) { + if blossom_i != blossom_j { + break; + } + s += 2 * dual_var[*blossom_i]; + } + assert!(s >= 0); + + if (mate.get(i).is_some() && mate.get(i).unwrap() / 2 == edge) + || (mate.get(j).is_some() && mate.get(j).unwrap() / 2 == edge) + { + assert!(mate[i] / 2 == edge && mate[j] / 2 == edge); + assert_eq!(s, 0); + } + } + // 2. all single vertices have zero dual value; + for (node, dual_var_node) in dual_var.iter().enumerate().take(num_nodes) { + assert!(mate.get(&node).is_some() || dual_var_node + node_dual_offset == 0); + } + // 3. all blossoms with positive dual value are full. + for blossom in num_nodes..2 * num_nodes { + if blossom_base[blossom].is_some() && dual_var[blossom] > 0 { + assert_eq!(blossom_endpoints[blossom].len() % 2, 1); + for p in blossom_endpoints[blossom].iter().skip(1).step_by(2) { + assert_eq!(mate.get(&endpoints[*p]).copied(), Some(p ^ 1)); + assert_eq!(mate.get(&endpoints[*p ^ 1]), Some(p)); + } + } + } +} + +/// Compute a maximum-weighted matching in the general undirected weighted +/// graph given by "edges". If `max_cardinality` is true, only +/// maximum-cardinality matchings are considered as solutions. +/// +/// The algorithm is based on "Efficient Algorithms for Finding Maximum +/// Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986. +/// +/// Based on networkx implementation +/// +/// +/// With reference to the standalone protoype implementation from: +/// +/// +/// +/// +/// The function takes time O(n**3) +/// +/// Arguments: +/// +/// * `graph` - The graph to compute the maximum weight matching for +/// * `max_cardinality` - If set to true compute the maximum-cardinality matching +/// with maximum weight among all maximum-cardinality matchings +/// * `verify_optimum_flag`: If true prior to returning an additional routine +/// to verify the optimal solution was found will be run after computing +/// the maximum weight matching. If it's true and the found matching is not +/// an optimal solution this function will panic. This option should +/// normally be only set true during testing. +/// +/// # Example +/// ```rust +/// +/// use hashbrown::HashSet; +/// use raphtory::core::entities::properties::props::Props; +/// use raphtory::prelude::{AdditionOps, Prop}; +/// use raphtory::algorithms::bipartite::max_weight_matching::max_weight_matching; +/// use raphtory_api::core::entities::GID; +/// +/// // Create a path graph +/// let g = raphtory::prelude::Graph::new(); +/// let vs = vec![ +/// (1, 2, 5), +/// (2, 3, 11), +/// (3, 4, 5)]; +/// for (src, dst,weight) in &vs { +/// g.add_edge(0, *src, *dst, [("weight",Prop::I64(*weight),)], None).unwrap(); +/// } +/// +/// // Run max weight matching with max cardinality set to false +/// let res = max_weight_matching( +/// &g, Some("weight"), false, true +/// ); +/// // Run max weight matching with max cardinality set to true +/// let maxc_res = max_weight_matching( +/// &g, Some("weight"), true, true +/// ); +/// +/// let matching = res; +/// let maxc_matching = maxc_res; +/// // Check output +/// assert_eq!(matching.len(), 1); +/// assert!(matching.contains(2, 3)); +/// assert_eq!(maxc_matching.len(), 2); +/// assert!(maxc_matching.contains(1, 2)); +/// assert!(maxc_matching.contains(3, 4)); +/// ``` +pub fn max_weight_matching<'graph, G: GraphViewOps<'graph>>( + g: &'graph G, + weight_prop: Option<&str>, + max_cardinality: bool, + verify_optimum_flag: bool, +) -> Matching { + let num_edges = g.count_edges(); + let num_nodes = g.count_nodes(); + // Exit fast for graph without edges + if num_edges == 0 { + return Matching::empty(g.clone()); + } + + let node_map: HashMap = g + .nodes() + .into_iter() + .enumerate() + .map(|(index, node)| (node.node, index)) + .collect(); + let mut edges: Vec<(usize, usize, i64)> = Vec::with_capacity(num_edges); + let mut max_weight: i64 = 0; + + g.edges().iter().for_each(|edge| { + let weight = match weight_prop { + None => 1, + Some(weight_prop) => edge + .properties() + .get(weight_prop) + .unwrap_or(Prop::I64(1)) + .into_i64() + .unwrap_or(1), + }; + if weight > max_weight { + max_weight = weight; + }; + edges.push(( + node_map[&edge.src().node], + node_map[&edge.dst().node], + weight, + )); + }); + // If p is an edge endpoint + // endpoints[p] is the node index to which endpoint p is attached + let endpoints: Vec = (0..2 * num_edges) + .map(|endpoint| { + let edge_tuple = edges[endpoint / 2]; + let out_value: usize = if endpoint % 2 == 0 { + edge_tuple.0 + } else { + edge_tuple.1 + }; + out_value + }) + .collect(); + // If v is a node/vertex + // neighbor_endpoints[v] is the list of remote endpoints of the edges + // attached to v. + // Not modified by algorithm (only mut to initially construct contents). + let mut neighbor_endpoints: Vec> = (0..num_nodes).map(|_| Vec::new()).collect(); + for edge in 0..num_edges { + neighbor_endpoints[edges[edge].0].push(2 * edge + 1); + neighbor_endpoints[edges[edge].1].push(2 * edge); + } + // If v is a vertex, + // mate[v] is the remote endpoint of its matched edge, or None if it is + // single (i.e. endpoint[mate[v]] is v's partner vertex). + // Initially all vertices are single; updated during augmentation. + let mut mate: HashMap = HashMap::with_capacity(num_nodes); + // If b is a top-level blossom + // label[b] is 0 if b is unlabeled (free); + // 1 if b is a S-vertex/blossom; + // 2 if b is a T-vertex/blossom; + // The label of a vertex/node is found by looking at the label of its + // top-level containing blossom + // If v is a node/vertex inside a T-Blossom, + // label[v] is 2 if and only if v is reachable from an S_Vertex outside + // the blossom. + let mut labels: Vec>; + // If b is a labeled top-level blossom, + // label_ends[b] is the remote endpoint of the edge through which b + // obtained its label, or None if b's base vertex is single. + // If v is a vertex inside a T-blossom and label[v] == 2, + // label_ends[v] is the remote endpoint of the edge through which v is + // reachable from outside the blossom + let mut label_ends: Vec>; + // If v is a vertex/node + // in_blossoms[v] is the top-level blossom to which v belongs. + // If v is a top-level vertex, v is itself a blossom (a trivial blossom) + // and in_blossoms[v] == v. + // Initially all nodes are top-level trivial blossoms. + let mut in_blossoms: Vec = (0..num_nodes).collect(); + // if b is a sub-blossom + // blossom_parents[b] is its immediate parent + // If b is a top-level blossom, blossom_parents[b] is None + let mut blossom_parents: Vec> = vec![None; 2 * num_nodes]; + // If b is a non-trivial (sub-)blossom, + // blossom_children[b] is an ordered list of its sub-blossoms, starting with + // the base and going round the blossom. + let mut blossom_children: Vec> = (0..2 * num_nodes).map(|_| Vec::new()).collect(); + // If b is a (sub-)blossom, + // blossombase[b] is its base VERTEX (i.e. recursive sub-blossom). + let mut blossom_base: Vec> = (0..num_nodes).map(Some).collect(); + blossom_base.append(&mut vec![None; num_nodes]); + // If b is a non-trivial (sub-)blossom, + // blossom_endpoints[b] is a list of endpoints on its connecting edges, + // such that blossom_endpoints[b][i] is the local endpoint of + // blossom_children[b][i] on the edge that connects it to + // blossom_children[b][wrap(i+1)]. + let mut blossom_endpoints: Vec> = (0..2 * num_nodes).map(|_| Vec::new()).collect(); + // If v is a free vertex (or an unreached vertex inside a T-blossom), + // best_edge[v] is the edge to an S-vertex with least slack, + // or None if there is no such edge. If b is a (possibly trivial) + // top-level S-blossom, + // best_edge[b] is the least-slack edge to a different S-blossom, + // or None if there is no such edge. + let mut best_edge: Vec>; + // If b is a non-trivial top-level S-blossom, + // blossom_best_edges[b] is a list of least-slack edges to neighboring + // S-blossoms, or None if no such list has been computed yet. + // This is used for efficient computation of delta3. + let mut blossom_best_edges: Vec> = (0..2 * num_nodes).map(|_| Vec::new()).collect(); + let mut unused_blossoms: Vec = (num_nodes..2 * num_nodes).collect(); + // If v is a vertex, + // dual_var[v] = 2 * u(v) where u(v) is the v's variable in the dual + // optimization problem (multiplication by two ensures integer values + // throughout the algorithm if all edge weights are integers). + // If b is a non-trivial blossom, + // dual_var[b] = z(b) where z(b) is b's variable in the dual optimization + // problem. + // dual_var is for vertex v in 0..num_nodes and blossom b in + // num_nodes..2* num nodes + let mut dual_var: Vec = vec![max_weight; num_nodes]; + dual_var.append(&mut vec![0; num_nodes]); + // If allowed_edge[k] is true, edge k has zero slack in the optimization + // problem; if allowed_edge[k] is false, the edge's slack may or may not + // be zero. + let mut allowed_edge: Vec; + // Queue of newly discovered S-vertices + let mut queue: Vec = Vec::with_capacity(num_nodes); + + // Main loop: continue until no further improvement is possible + for _ in 0..num_nodes { + // Each iteration of this loop is a "stage". + // A stage finds an augmenting path and uses that to improve + // the matching. + + // Removal labels from top-level blossoms/vertices + labels = vec![Some(0); 2 * num_nodes]; + label_ends = vec![None; 2 * num_nodes]; + // Forget all about least-slack edges. + best_edge = vec![None; 2 * num_nodes]; + blossom_best_edges.splice(num_nodes.., (0..num_nodes).map(|_| Vec::new())); + // Loss of labeling means that we can not be sure that currently + // allowable edges remain allowable througout this stage. + allowed_edge = vec![false; num_edges]; + // Make queue empty + queue.clear(); + // Label single blossom/vertices with S and put them in queue. + for v in 0..num_nodes { + if mate.get(&v).is_none() && labels[in_blossoms[v]] == Some(0) { + assign_label( + v, + 1, + None, + num_nodes, + &in_blossoms, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &blossom_children, + &blossom_base, + &endpoints, + &mate, + ); + } + } + // Loop until we succeed in augmenting the matching. + let mut augmented = false; + loop { + // Each iteration of this loop is a "substage". + // A substage tries to find an augmenting path; + // if found, the path is used to improve the matching and + // the stage ends. If there is no augmenting path, the + // primal-dual method is used to find some slack out of + // the dual variables. + + // Continue labeling until all vertices which are reachable + // through an alternating path have got a label. + while !queue.is_empty() && !augmented { + // Take an S vertex from the queue + let v = queue.pop().unwrap(); + assert_eq!(labels[in_blossoms[v]], Some(1)); + + // Scan its neighbors + for p in &neighbor_endpoints[v] { + let k = *p / 2; + let mut kslack = 0; + let w = endpoints[*p]; + // w is a neighbor of v + if in_blossoms[v] == in_blossoms[w] { + // this edge is internal to a blossom; ignore it + continue; + } + if !allowed_edge[k] { + kslack = slack(k, &dual_var, &edges); + if kslack <= 0 { + // edge k has zero slack -> it is allowable + allowed_edge[k] = true; + } + } + if allowed_edge[k] { + if labels[in_blossoms[w]] == Some(0) { + // (C1) w is a free vertex; + // label w with T and label its mate with S (R12). + assign_label( + w, + 2, + Some(*p ^ 1), + num_nodes, + &in_blossoms, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &blossom_children, + &blossom_base, + &endpoints, + &mate, + ); + } else if labels[in_blossoms[w]] == Some(1) { + // (C2) w is an S-vertex (not in the same blossom); + // follow back-links to discover either an + // augmenting path or a new blossom. + let base = scan_blossom( + v, + w, + &in_blossoms, + &blossom_base, + &endpoints, + &mut labels, + &label_ends, + &mate, + ); + match base { + // Found a new blossom; add it to the blossom + // bookkeeping and turn it into an S-blossom. + Some(base) => add_blossom( + base, + k, + &mut blossom_children, + num_nodes, + &edges, + &mut in_blossoms, + &mut dual_var, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &mut blossom_base, + &endpoints, + &mut blossom_endpoints, + &mut unused_blossoms, + &mut blossom_best_edges, + &mut blossom_parents, + &neighbor_endpoints, + &mate, + ), + // Found an augmenting path; augment the + // matching and end this stage. + None => { + augment_matching( + k, + num_nodes, + &edges, + &in_blossoms, + &labels, + &label_ends, + &blossom_parents, + &endpoints, + &mut blossom_children, + &mut blossom_endpoints, + &mut blossom_base, + &mut mate, + ); + augmented = true; + break; + } + }; + } else if labels[w] == Some(0) { + // w is inside a T-blossom, but w itself has not + // yet been reached from outside the blossom; + // mark it as reached (we need this to relabel + // during T-blossom expansion). + assert_eq!(labels[in_blossoms[w]], Some(2)); + labels[w] = Some(2); + label_ends[w] = Some(*p ^ 1); + } + } else if labels[in_blossoms[w]] == Some(1) { + // Keep track of the least-slack non-allowable edge to + // a different S-blossom + let blossom = in_blossoms[v]; + if best_edge[blossom].is_none() + || kslack < slack(best_edge[blossom].unwrap(), &dual_var, &edges) + { + best_edge[blossom] = Some(k); + } + } else if labels[w] == Some(0) { + // w is a free vertex (or an unreached vertex inside + // a T-blossom) but we can not reach it yet; + // keep track of the least-slack edge that reaches w. + if best_edge[w].is_none() + || kslack < slack(best_edge[w].unwrap(), &dual_var, &edges) + { + best_edge[w] = Some(k) + } + } + } + } + if augmented { + break; + } + // There is no augmenting path under these constraints; + // compute delta and reduce slack in the optimization problem. + // (Note that our vertex dual variables, edge slacks and delta's + // are pre-multiplied by two.) + let mut delta_type = -1; + let mut delta: Option = None; + let mut delta_edge: Option = None; + let mut delta_blossom: Option = None; + + // Compute delta1: the minimum value of any vertex dual. + if !max_cardinality { + delta_type = 1; + delta = Some(*dual_var[..num_nodes].iter().min().unwrap()); + } + + // Compute delta2: the minimum slack on any edge between + // an S-vertex and a free vertex. + for v in 0..num_nodes { + if labels[in_blossoms[v]] == Some(0) && best_edge[v].is_some() { + let d = slack(best_edge[v].unwrap(), &dual_var, &edges); + if delta_type == -1 || Some(d) < delta { + delta = Some(d); + delta_type = 2; + delta_edge = best_edge[v]; + } + } + } + + // Compute delta3: half the minimum slack on any edge between a + // pair of S-blossoms + for blossom in 0..2 * num_nodes { + if blossom_parents[blossom].is_none() + && labels[blossom] == Some(1) + && best_edge[blossom].is_some() + { + let kslack = slack(best_edge[blossom].unwrap(), &dual_var, &edges); + assert_eq!(kslack % 2, 0); + let d = Some(kslack / 2); + if delta_type == -1 || d < delta { + delta = d; + delta_type = 3; + delta_edge = best_edge[blossom]; + } + } + } + + // Compute delta4: minimum z variable of any T-blossom + for blossom in num_nodes..2 * num_nodes { + if blossom_base[blossom].is_some() + && blossom_parents[blossom].is_none() + && labels[blossom] == Some(2) + && (delta_type == -1 || dual_var[blossom] < delta.unwrap()) + { + delta = Some(dual_var[blossom]); + delta_type = 4; + delta_blossom = Some(blossom); + } + } + if delta_type == -1 { + // No further improvement possible; max-cardinality optimum + // reached. Do a final delta update to make the optimum + // verifyable + assert!(max_cardinality); + delta_type = 1; + delta = Some(max(0, *dual_var[..num_nodes].iter().min().unwrap())); + } + + // Update dual variables according to delta. + for v in 0..num_nodes { + if labels[in_blossoms[v]] == Some(1) { + // S-vertex: 2*u = 2*u - 2*delta + dual_var[v] -= delta.unwrap(); + } else if labels[in_blossoms[v]] == Some(2) { + // T-vertex: 2*u = 2*u + 2*delta + dual_var[v] += delta.unwrap(); + } + } + for b in num_nodes..2 * num_nodes { + if blossom_base[b].is_some() && blossom_parents[b].is_none() { + if labels[b] == Some(1) { + // top-level S-blossom: z = z + 2*delta + dual_var[b] += delta.unwrap(); + } else if labels[b] == Some(2) { + // top-level T-blossom: z = z - 2*delta + dual_var[b] -= delta.unwrap(); + } + } + } + // Take action at the point where minimum delta occured. + if delta_type == 1 { + // No further improvement possible; optimum reached + break; + } else if delta_type == 2 { + // Use the least-slack edge to continue the search. + allowed_edge[delta_edge.unwrap()] = true; + let (mut i, mut j, _weight) = edges[delta_edge.unwrap()]; + if labels[in_blossoms[i]] == Some(0) { + mem::swap(&mut i, &mut j); + } + assert_eq!(labels[in_blossoms[i]], Some(1)); + queue.push(i); + } else if delta_type == 3 { + // Use the least-slack edge to continue the search. + allowed_edge[delta_edge.unwrap()] = true; + let (i, _j, _weight) = edges[delta_edge.unwrap()]; + assert_eq!(labels[in_blossoms[i]], Some(1)); + queue.push(i); + } else if delta_type == 4 { + // Expand the least-z blossom + expand_blossom( + delta_blossom.unwrap(), + false, + num_nodes, + &mut blossom_children, + &mut blossom_parents, + &mut in_blossoms, + &dual_var, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &mut blossom_base, + &endpoints, + &mate, + &mut blossom_endpoints, + &mut allowed_edge, + &mut unused_blossoms, + ); + } + // end of this substage + } + // Stop when no more augment paths can be found + if !augmented { + break; + } + + // End of a stage; expand all S-blossoms which have a dual_var == 0 + for blossom in num_nodes..2 * num_nodes { + if blossom_parents[blossom].is_none() + && blossom_base[blossom].is_some() + && labels[blossom] == Some(1) + && dual_var[blossom] == 0 + { + expand_blossom( + blossom, + true, + num_nodes, + &mut blossom_children, + &mut blossom_parents, + &mut in_blossoms, + &dual_var, + &mut labels, + &mut label_ends, + &mut best_edge, + &mut queue, + &mut blossom_base, + &endpoints, + &mate, + &mut blossom_endpoints, + &mut allowed_edge, + &mut unused_blossoms, + ); + } + } + } + if verify_optimum_flag { + verify_optimum( + max_cardinality, + num_nodes, + num_edges, + &edges, + &endpoints, + &dual_var, + &blossom_parents, + &blossom_endpoints, + &blossom_base, + &mate, + ); + } + + // Transform mate[] such that mate[v] is the vertex to which v is paired + // Also handle holes in node indices from PyGraph node removals by mapping + // linear index to node index. + Matching::from_mates(g.clone(), mate, endpoints) +} + +#[derive(Clone)] +pub struct Matching { + graph: G, + forward_map: Arc>, + reverse_map: Arc>, +} + +impl Matching { + pub(crate) fn into_dyn(self) -> Matching { + Matching { + graph: self.graph.into_dynamic(), + forward_map: self.forward_map, + reverse_map: self.reverse_map, + } + } +} + +impl<'graph, G: GraphViewOps<'graph>> Debug for Matching { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + f.debug_struct("Matching") + .field("forward_map", &self.forward_map) + .field("reverse_map", &self.reverse_map) + .finish() + } +} + +impl<'graph, G: GraphViewOps<'graph>> Display for Matching { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + write!(f, "Matching([")?; + for (src, dst) in self.edges().id() { + write!(f, "({src}, {dst})")?; + } + write!(f, "])") + } +} + +impl<'graph, G: GraphViewOps<'graph>> Matching { + fn empty(graph: G) -> Self { + Matching { + graph, + forward_map: Default::default(), + reverse_map: Default::default(), + } + } + fn from_mates(graph: G, mates: HashMap, endpoints: Vec) -> Self { + let mut forward_map = HashMap::with_capacity(mates.len() / 2); + let mut reverse_map = HashMap::with_capacity(mates.len() / 2); + let node_map: Vec<_> = graph.nodes().iter().map(|node| node.node).collect(); + let edge_map: Vec<_> = graph.edges().iter().map(|edge| edge.edge.pid()).collect(); + for (node, edge) in mates.iter() { + let eid = edge_map[*edge / 2]; + if edge % 2 == 0 { + reverse_map.insert(node_map[*node], (node_map[endpoints[*edge]], eid)); + } else { + forward_map.insert(node_map[*node], (node_map[endpoints[*edge]], eid)); + } + } + Self { + graph, + forward_map: forward_map.into(), + reverse_map: reverse_map.into(), + } + } + + pub fn len(&self) -> usize { + self.forward_map.len() + } + + pub fn is_empty(&self) -> bool { + self.forward_map.is_empty() + } + + pub fn edges_iter<'a>(&'a self) -> impl Iterator> + where + 'graph: 'a, + { + let storage = self.graph.core_edges(); + self.forward_map + .values() + .map(move |(_, eid)| EdgeView::new(&self.graph, storage.as_ref().edge(*eid).out_ref())) + } + + pub fn edges(&self) -> Edges<'graph, G> { + let storage = self.graph.core_graph().clone(); + let forward_map = self.forward_map.clone(); + let edges_iter = Arc::new(move || { + let storage = storage.clone(); + let forward_map = forward_map.clone(); + GenLockedIter::from(forward_map, move |forward_map| { + forward_map + .values() + .map(move |(_, eid)| storage.edge_entry(*eid).out_ref()) + .into_dyn_boxed() + }) + .into_dyn_boxed() + }); + Edges { + base_graph: self.graph.clone(), + graph: self.graph.clone(), + edges: edges_iter, + } + } + + pub fn contains(&self, src: N, dst: N) -> bool { + if let Some(src) = self.graph.internalise_node(src.as_node_ref()) { + if let Some(dst) = self.graph.internalise_node(dst.as_node_ref()) { + if let Some((nbr, _)) = self.forward_map.get(&src) { + return nbr == &dst; + } + } + } + false + } + + pub fn src<'a>(&'a self, dst: impl AsNodeRef) -> Option> + where + 'graph: 'a, + { + let dst = (&&self.graph).node(dst)?.node; + let (src, _) = self.reverse_map.get(&dst)?; + Some(NodeView::new_internal(&self.graph, *src)) + } + + pub fn edge_for_src<'a>(&'a self, src: impl AsNodeRef) -> Option> + where + 'graph: 'a, + { + let src = (&&self.graph).node(src)?.node; + let (_, eid) = self.forward_map.get(&src)?; + Some(EdgeView::new( + &self.graph, + self.graph.core_edge(*eid).out_ref(), + )) + } + pub fn dst<'a>(&'a self, src: impl AsNodeRef) -> Option> + where + 'graph: 'a, + { + let src = (&&self.graph).node(src)?.node; + let (dst, _) = self.forward_map.get(&src)?; + Some(NodeView::new_internal(&self.graph, *dst)) + } + + pub fn edge_for_dst<'a>(&'a self, dst: impl AsNodeRef) -> Option> + where + 'graph: 'a, + { + let dst = (&&self.graph).node(dst)?.node; + let (_, eid) = self.reverse_map.get(&dst)?; + Some(EdgeView::new( + &self.graph, + self.graph.core_edge(*eid).out_ref(), + )) + } +} + +#[cfg(test)] +mod test { + use crate::{algorithms::bipartite::max_weight_matching::max_weight_matching, prelude::*}; + use itertools::Itertools; + + #[test] + fn test() { + let g = Graph::new(); + let vs = vec![(1, 2, 5), (2, 3, 11), (3, 4, 5)]; + for (src, dst, weight) in &vs { + g.add_edge(0, *src, *dst, [("weight", Prop::I64(*weight))], None) + .unwrap(); + } + + // Run max weight matching with max cardinality set to false + let res = max_weight_matching(&g, Some("weight"), false, true); + // Run max weight matching with max cardinality set to true + let maxc_res = max_weight_matching(&g, Some("weight"), true, true); + + let matching = res; + println!("{}", matching); + let maxc_matching = maxc_res; + // Check output + assert_eq!(matching.len(), 1); + assert!(matching.contains(2, 3)); + assert_eq!(maxc_matching.len(), 2); + assert!(maxc_matching.contains(1, 2)); + assert!(maxc_matching.contains(3, 4)); + + assert_eq!(matching.src(3).unwrap().id(), 2); + assert_eq!(matching.src(2), None); + + assert_eq!(matching.dst(2).unwrap().id(), 3); + assert_eq!(matching.dst(3), None); + + assert_eq!(matching.edge_for_src(2).unwrap(), g.edge(2, 3).unwrap()); + assert_eq!(matching.edge_for_src(1), None); + + assert_eq!(matching.edge_for_dst(3).unwrap(), g.edge(2, 3).unwrap()); + assert_eq!(matching.edge_for_dst(2), None); + + assert_eq!(matching.edges().collect(), vec![g.edge(2, 3).unwrap()]); + assert_eq!( + matching.edges_iter().collect_vec(), + vec![g.edge(2, 3).unwrap()] + ); + } +} diff --git a/raphtory/src/algorithms/bipartite/mod.rs b/raphtory/src/algorithms/bipartite/mod.rs new file mode 100644 index 000000000..0281237db --- /dev/null +++ b/raphtory/src/algorithms/bipartite/mod.rs @@ -0,0 +1 @@ +pub mod max_weight_matching; diff --git a/raphtory/src/algorithms/community_detection/louvain.rs b/raphtory/src/algorithms/community_detection/louvain.rs index 6b2a4cfed..808ea81cc 100644 --- a/raphtory/src/algorithms/community_detection/louvain.rs +++ b/raphtory/src/algorithms/community_detection/louvain.rs @@ -70,9 +70,11 @@ mod test { test_storage, }; use proptest::prelude::*; - use raphtory_api::core::utils::logging::global_info_logger; use tracing::info; + #[cfg(feature = "io")] + use raphtory_api::core::utils::logging::global_info_logger; + #[test] fn test_louvain() { let edges = vec![ diff --git a/raphtory/src/algorithms/mod.rs b/raphtory/src/algorithms/mod.rs index 7003b8361..8f21dd870 100644 --- a/raphtory/src/algorithms/mod.rs +++ b/raphtory/src/algorithms/mod.rs @@ -30,6 +30,7 @@ pub mod algorithm_result; pub mod centrality; pub mod community_detection; +pub mod bipartite; pub mod components; pub mod cores; pub mod dynamics; diff --git a/raphtory/src/db/graph/graph.rs b/raphtory/src/db/graph/graph.rs index ac763a756..0b07cd0a1 100644 --- a/raphtory/src/db/graph/graph.rs +++ b/raphtory/src/db/graph/graph.rs @@ -1941,14 +1941,11 @@ mod db_tests { test_graph(&graph, |graph| { assert_eq!( graph.nodes().id().collect::>(), - vec![1u64.into(), 2u64.into(), 3u64.into()] + vec![1u64, 2u64, 3u64] ); let g_at = graph.at(1); - assert_eq!( - g_at.nodes().id().collect::>(), - vec![1u64.into(), 2u64.into()] - ); + assert_eq!(g_at.nodes().id().collect::>(), vec![1u64, 2u64]); }); } diff --git a/raphtory/src/python/algorithm/max_weight_matching.rs b/raphtory/src/python/algorithm/max_weight_matching.rs new file mode 100644 index 000000000..ef7d8f788 --- /dev/null +++ b/raphtory/src/python/algorithm/max_weight_matching.rs @@ -0,0 +1,121 @@ +use crate::{ + algorithms::bipartite::max_weight_matching::Matching, + db::{ + api::view::{DynamicGraph, IntoDynamic, StaticGraphViewOps}, + graph::{edge::EdgeView, edges::Edges, node::NodeView}, + }, + prelude::GraphViewOps, + python::{ + types::{repr::Repr, wrappers::iterators::PyBorrowingIterator}, + utils::PyNodeRef, + }, +}; +use pyo3::prelude::*; + +/// A Matching (i.e., a set of edges that do not share any nodes) +#[pyclass(frozen, name = "Matching")] +pub struct PyMatching { + inner: Matching, +} + +impl IntoPy for Matching { + fn into_py(self, py: Python) -> PyObject { + PyMatching { + inner: self.into_dyn(), + } + .into_py(py) + } +} + +#[pymethods] +impl PyMatching { + /// The number of edges in the matching + fn __len__(&self) -> usize { + self.inner.len() + } + + /// Returns True if the matching is not empty + fn __bool__(&self) -> bool { + !self.inner.is_empty() + } + + fn __iter__(&self) -> PyBorrowingIterator { + py_borrowing_iter!(self.inner.clone(), Matching, |inner| inner + .edges_iter()) + } + + /// Get the matched source node for a destination node + /// + /// Arguments: + /// dst (InputNode): The destination node + /// + /// Returns: + /// Optional[Node]: The matched source node if it exists + /// + fn src(&self, dst: PyNodeRef) -> Option> { + self.inner.src(dst).map(|n| n.cloned()) + } + + /// Get the matched destination node for a source node + /// + /// Arguments: + /// src (InputNode): The source node + /// + /// Returns: + /// Optional[Node]: The matched destination node if it exists + /// + fn dst(&self, src: PyNodeRef) -> Option> { + self.inner.dst(src).map(|n| n.cloned()) + } + + /// Get a view of the matched edges + /// + /// Returns: + /// Edges: The edges in the matching + fn edges(&self) -> Edges<'static, DynamicGraph> { + self.inner.edges() + } + + /// Get the matched edge for a source node + /// + /// Arguments: + /// src (InputNode): The source node + /// + /// Returns: + /// Optional[Edge]: The matched edge if it exists + fn edge_for_src(&self, src: PyNodeRef) -> Option> { + self.inner.edge_for_src(src).map(|e| e.cloned()) + } + + /// Get the matched edge for a destination node + /// + /// Arguments: + /// dst (InputNode): The source node + /// + /// Returns: + /// Optional[Edge]: The matched edge if it exists + fn edge_for_dst(&self, dst: PyNodeRef) -> Option> { + self.inner.edge_for_dst(dst).map(|e| e.cloned()) + } + + /// Check if an edge is part of the matching + /// + /// Arguments: + /// edge (Tuple[InputNode, InputNode]): The edge to check + /// + /// Returns: + /// bool: Returns True if the edge is part of the matching, False otherwise + fn __contains__(&self, edge: (PyNodeRef, PyNodeRef)) -> bool { + self.inner.contains(edge.0, edge.1) + } + + fn __repr__(&self) -> String { + self.inner.repr() + } +} + +impl<'graph, G: GraphViewOps<'graph>> Repr for Matching { + fn repr(&self) -> String { + format!("{self}") + } +} diff --git a/raphtory/src/python/algorithm/mod.rs b/raphtory/src/python/algorithm/mod.rs index a10f38cf7..078eb0e19 100644 --- a/raphtory/src/python/algorithm/mod.rs +++ b/raphtory/src/python/algorithm/mod.rs @@ -1 +1,2 @@ pub(crate) mod epidemics; +pub(crate) mod max_weight_matching; diff --git a/raphtory/src/python/packages/algorithms.rs b/raphtory/src/python/packages/algorithms.rs index 19633452c..7aef68bd3 100644 --- a/raphtory/src/python/packages/algorithms.rs +++ b/raphtory/src/python/packages/algorithms.rs @@ -3,6 +3,7 @@ use crate::{ algorithms::{ algorithm_result::AlgorithmResult, + bipartite::max_weight_matching::max_weight_matching as mwm, centrality::{ betweenness::betweenness_centrality as betweenness_rs, degree_centrality::degree_centrality as degree_centrality_rs, hits::hits as hits_rs, @@ -62,6 +63,7 @@ use rand::{prelude::StdRng, SeedableRng}; use raphtory_api::core::{entities::GID, Direction}; use std::collections::{HashMap, HashSet}; +use crate::algorithms::bipartite::max_weight_matching::Matching; #[cfg(feature = "storage")] use crate::python::graph::disk_graph::PyDiskGraph; #[cfg(feature = "storage")] @@ -862,3 +864,50 @@ pub fn temporal_rich_club_coefficient( .collect::>>()?; Ok(temporal_rich_club_rs(graph.graph, views, k, delta)) } + +/// Compute a maximum-weighted matching in the general undirected weighted +/// graph given by "edges". If `max_cardinality` is true, only +/// maximum-cardinality matchings are considered as solutions. +/// +/// The algorithm is based on "Efficient Algorithms for Finding Maximum +/// Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986. +/// +/// Based on networkx implementation +/// +/// +/// With reference to the standalone protoype implementation from: +/// +/// +/// +/// +/// The function takes time O(n**3) +/// +/// Arguments: +/// graph (GraphView): The graph to compute the maximum weight matching for +/// weight_prop (str, optional): The property on the edge to use for the weight. If not +/// provided, +/// max_cardinality (bool): If set to true compute the maximum-cardinality matching +/// with maximum weight among all maximum-cardinality matchings. Defaults to True. +/// verify_optimum_flag (bool): If true prior to returning an additional routine +/// to verify the optimal solution was found will be run after computing +/// the maximum weight matching. If it's true and the found matching is not +/// an optimal solution this function will panic. This option should +/// normally be only set true during testing. Defaults to False. +/// +/// Returns: +/// Matching: The matching +#[pyfunction] +#[pyo3(signature = (graph, weight_prop=None, max_cardinality=true, verify_optimum_flag=false))] +pub fn max_weight_matching( + graph: &PyGraphView, + weight_prop: Option<&str>, + max_cardinality: bool, + verify_optimum_flag: bool, +) -> Matching { + mwm( + &graph.graph, + weight_prop, + max_cardinality, + verify_optimum_flag, + ) +} diff --git a/raphtory/src/python/packages/base_modules.rs b/raphtory/src/python/packages/base_modules.rs index 9304345a6..46f8814a0 100644 --- a/raphtory/src/python/packages/base_modules.rs +++ b/raphtory/src/python/packages/base_modules.rs @@ -5,6 +5,7 @@ use crate::python::graph::disk_graph::PyDiskGraph; use crate::{ add_classes, add_functions, python::{ + algorithm::max_weight_matching::PyMatching, graph::{ algorithm_result::AlgorithmResult, edge::{PyEdge, PyMutableEdge}, @@ -58,7 +59,7 @@ pub fn add_raphtory_classes(m: &Bound) -> PyResult<()> { #[cfg(feature = "storage")] add_classes!(m, PyDiskGraph); - return Ok(()); + Ok(()) } pub fn base_algorithm_module(py: Python<'_>) -> Result, PyErr> { @@ -102,8 +103,10 @@ pub fn base_algorithm_module(py: Python<'_>) -> Result, PyErr> { louvain, fruchterman_reingold, cohesive_fruchterman_reingold, + max_weight_matching ); + add_classes!(&algorithm_module, PyMatching); #[cfg(feature = "storage")] add_functions!(&algorithm_module, connected_components); Ok(algorithm_module) @@ -120,7 +123,7 @@ pub fn base_graph_loader_module(py: Python<'_>) -> Result, PyErr reddit_hyperlink_graph_local, karate_club_graph, ); - return Ok(graph_loader_module); + Ok(graph_loader_module) } pub fn base_graph_gen_module(py: Python<'_>) -> Result, PyErr> { @@ -130,7 +133,7 @@ pub fn base_graph_gen_module(py: Python<'_>) -> Result, PyErr> { random_attachment, ba_preferential_attachment, ); - return Ok(graph_gen_module); + Ok(graph_gen_module) } pub fn base_vectors_module(py: Python<'_>) -> Result, PyErr> { @@ -138,5 +141,5 @@ pub fn base_vectors_module(py: Python<'_>) -> Result, PyErr> { vectors_module.add_class::()?; vectors_module.add_class::()?; vectors_module.add_class::()?; - return Ok(vectors_module); + Ok(vectors_module) }