diff --git a/examples/benchmark_parallel_core_decomposition.py b/examples/benchmark_parallel_core_decomposition.py new file mode 100644 index 000000000..2c700811f --- /dev/null +++ b/examples/benchmark_parallel_core_decomposition.py @@ -0,0 +1,257 @@ +#!/usr/bin/env python3 + +"""Compare ParK with adaptive parallel core decomposition (including its ParK sparse path).""" + +import argparse +import math +import statistics +import time + +import networkit as nk + + +def validate(algorithm, expected_scores, expected_max_core, name): + actual_scores = algorithm.scores() + actual_max_core = algorithm.maxCoreNumber() + if actual_max_core != expected_max_core: + raise RuntimeError( + f"{name} returned max core {actual_max_core}; expected {expected_max_core}") + if len(actual_scores) != len(expected_scores): + raise RuntimeError( + f"{name} returned {len(actual_scores)} scores; expected {len(expected_scores)}") + if actual_scores != expected_scores: + mismatch = next( + node for node, (expected, actual) in enumerate(zip(expected_scores, actual_scores)) + if expected != actual) + raise RuntimeError( + f"{name} returned core {actual_scores[mismatch]} for node {mismatch}; " + f"expected {expected_scores[mismatch]}") + + +def measure(name, factory, expected_scores, expected_max_core, repeats): + samples = [] + for _ in range(repeats): + algorithm = factory() + started = time.perf_counter() + algorithm.run() + samples.append(time.perf_counter() - started) + validate(algorithm, expected_scores, expected_max_core, name) + return statistics.median(samples), samples + + +def benchmark(name, graph, repeats): + oracle = nk.centrality.CoreDecomposition( + graph, normalized=False, enforceBucketQueueAlgorithm=True) + oracle.run() + expected_scores = oracle.scores() + expected_max_core = oracle.maxCoreNumber() + + park_median, park_samples = measure( + "ParK", lambda: nk.centrality.CoreDecomposition(graph), expected_scores, + expected_max_core, repeats) + parallel_median, parallel_samples = measure( + "ParallelCore", lambda: nk.centrality.ParallelCoreDecomposition(graph), expected_scores, + expected_max_core, repeats) + + print(f"\n{name}: n={graph.numberOfNodes()} m={graph.numberOfEdges()}") + print( + f"validated={2 * repeats} timed runs against exact per-node scores; " + f"max_core={expected_max_core}") + print(f"ParK: median={park_median:.6f}s samples={park_samples}") + print(f"ParallelCore: median={parallel_median:.6f}s samples={parallel_samples}") + ratio = park_median / parallel_median if parallel_median > 0 else math.inf + print(f"ParK/ParallelCore={ratio:.3f}x") + + +def sparse_er(): + nk.setSeed(42, False) + return nk.generators.ErdosRenyiGenerator(200_000, 0.00005, False).generate() + + +def vgc_path(): + graph = nk.Graph(1_000_000) + for node in range(1, graph.numberOfNodes()): + graph.addEdge(node - 1, node) + return graph + + +def sampling_star(): + graph = nk.Graph(1_000_001) + for leaf in range(1, graph.numberOfNodes()): + graph.addEdge(0, leaf) + return graph + + +def sampling_bipartite(): + graph = nk.Graph(20_020) + for high_degree_node in range(20): + for low_degree_node in range(20, graph.numberOfNodes()): + graph.addEdge(high_degree_node, low_degree_node) + return graph + + +def hbs_clique(): + return nk.generators.ErdosRenyiGenerator(2_000, 1.0, False).generate() + + +def graph_from_symmetric_csr(indices, indptr): + """Build an undirected zero-copy GraphR from NumPy uint64 CSR storage.""" + import pyarrow as pa + + return nk.Graph.fromCSR( + len(indptr) - 1, + False, + pa.array(indices, type=pa.uint64()), + pa.array(indptr, type=pa.uint64()), + ) + + +def large_bipartite(target_edges): + """Build K(left, right) with exactly ``target_edges`` when divisible by left.""" + import numpy as np + + left = 1_000 + if target_edges < left or target_edges % left: + raise ValueError(f"large bipartite edge count must be divisible by {left}") + right = target_edges // left + nodes = left + right + adjacency_entries = 2 * target_edges + + indices = np.empty(adjacency_entries, dtype=np.uint64) + right_nodes = np.arange(left, nodes, dtype=np.uint64) + for node in range(left): + begin = node * right + indices[begin:begin + right] = right_nodes + + left_nodes = np.arange(left, dtype=np.uint64) + reverse_begin = target_edges + rows_per_chunk = 16_384 + for row_begin in range(0, right, rows_per_chunk): + row_end = min(row_begin + rows_per_chunk, right) + begin = reverse_begin + row_begin * left + end = reverse_begin + row_end * left + indices[begin:end].reshape(row_end - row_begin, left)[:] = left_nodes + + indptr = np.empty(nodes + 1, dtype=np.uint64) + indptr[:left + 1] = np.arange(left + 1, dtype=np.uint64) * right + indptr[left + 1:] = target_edges + np.arange(1, right + 1, dtype=np.uint64) * left + return graph_from_symmetric_csr(indices, indptr) + + +def large_clique(target_edges): + """Build the largest complete graph with no more than ``target_edges`` edges.""" + import numpy as np + + nodes = (1 + math.isqrt(1 + 8 * target_edges)) // 2 + edges = nodes * (nodes - 1) // 2 + adjacency_entries = 2 * edges + indices = np.empty(adjacency_entries, dtype=np.uint64) + all_nodes = np.arange(nodes, dtype=np.uint64) + row_width = nodes - 1 + for node in range(nodes): + begin = node * row_width + indices[begin:begin + node] = all_nodes[:node] + indices[begin + node:begin + row_width] = all_nodes[node + 1:] + + indptr = np.arange(nodes + 1, dtype=np.uint64) * row_width + return graph_from_symmetric_csr(indices, indptr) + + +def large_sparse_er(target_edges): + """Build a sparse Erdos-Renyi graph with about ``target_edges`` edges.""" + nodes = max(1_000, target_edges // 20) + probability = (2.0 * target_edges) / (nodes * (nodes - 1)) + if probability > 1.0: + raise ValueError("large sparse ER target is too dense for the selected node count") + nk.setSeed(42, False) + return nk.generators.ErdosRenyiGenerator(nodes, probability, False).generate() + + +def large_star(target_edges): + """Build a zero-copy star with exactly ``target_edges`` edges.""" + import numpy as np + + nodes = target_edges + 1 + indices = np.empty(2 * target_edges, dtype=np.uint64) + indices[:target_edges] = np.arange(1, nodes, dtype=np.uint64) + indices[target_edges:] = 0 + + indptr = np.empty(nodes + 1, dtype=np.uint64) + indptr[0] = 0 + indptr[1:] = target_edges + np.arange(nodes, dtype=np.uint64) + return graph_from_symmetric_csr(indices, indptr) + + +def large_path(target_edges): + """Build a zero-copy path with exactly ``target_edges`` edges.""" + import numpy as np + + nodes = target_edges + 1 + indices = np.empty(2 * target_edges, dtype=np.uint64) + indices[0] = 1 + if nodes > 2: + internal_nodes = np.arange(1, nodes - 1, dtype=np.uint64) + indices[1:-1:2] = internal_nodes - 1 + indices[2:-1:2] = internal_nodes + 1 + indices[-1] = nodes - 2 + + indptr = np.empty(nodes + 1, dtype=np.uint64) + indptr[0] = 0 + indptr[1:nodes] = 2 * np.arange(1, nodes, dtype=np.uint64) - 1 + indptr[nodes] = 2 * target_edges + return graph_from_symmetric_csr(indices, indptr) + + +CASES = { + "sparse-er": sparse_er, + "vgc-path": vgc_path, + "sampling-star": sampling_star, + "sampling-bipartite": sampling_bipartite, + "hbs-clique": hbs_clique, +} + +LARGE_CASES = { + "large-bipartite": large_bipartite, + "large-clique": large_clique, + "large-sparse-er": large_sparse_er, + "large-star": large_star, + "large-path": large_path, +} + +LARGE_ALL_CASES = ("large-bipartite", "large-clique") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--case", choices=["all", "large-all", *CASES, *LARGE_CASES], default="all") + parser.add_argument("--repeats", type=int, default=5) + parser.add_argument("--threads", type=int, default=15) + parser.add_argument( + "--large-edges", type=int, default=100_000_000, + help="target undirected edge count for large CSR cases (default: 100 million)") + args = parser.parse_args() + + if args.repeats < 1: + parser.error("--repeats must be positive") + if args.threads < 1: + parser.error("--threads must be positive") + if args.large_edges < 1: + parser.error("--large-edges must be positive") + + nk.setNumberOfThreads(args.threads) + if args.case == "all": + selected = [(name, factory) for name, factory in CASES.items()] + elif args.case == "large-all": + selected = [(name, LARGE_CASES[name]) for name in LARGE_ALL_CASES] + elif args.case in LARGE_CASES: + selected = [(args.case, LARGE_CASES[args.case])] + else: + selected = [(args.case, CASES[args.case])] + for name, factory in selected: + graph = factory(args.large_edges) if name in LARGE_CASES else factory() + benchmark(name, graph, args.repeats) + + +if __name__ == "__main__": + main() diff --git a/include/networkit/auxiliary/HierarchicalCoreBuckets.hpp b/include/networkit/auxiliary/HierarchicalCoreBuckets.hpp new file mode 100644 index 000000000..dc9f10aa2 --- /dev/null +++ b/include/networkit/auxiliary/HierarchicalCoreBuckets.hpp @@ -0,0 +1,72 @@ +#ifndef NETWORKIT_AUXILIARY_HIERARCHICAL_CORE_BUCKETS_HPP_ +#define NETWORKIT_AUXILIARY_HIERARCHICAL_CORE_BUCKETS_HPP_ + +#include +#include +#include + +namespace Aux { +namespace CoreBuckets { +/** + * @ingroup auxiliary + * Bucket routing for online core peeling: eight exact levels and logarithmic intervals. + * Callers own synchronization, stale-entry filtering and intermediate redistribution. + * Base levels are multiples of eight; routed degrees must belong to the same aligned + * 512-level window as the base. The peeling engine advances that window separately. + */ +constexpr NetworKit::count kLog2SingleBuckets = 3; +constexpr NetworKit::count kNumSingleBuckets = NetworKit::count{1} << kLog2SingleBuckets; +constexpr NetworKit::count kNumIntermediateBuckets = 6; +constexpr NetworKit::count kNumBuckets = kNumSingleBuckets + kNumIntermediateBuckets; +constexpr NetworKit::count kBucketMask = kNumSingleBuckets - 1; +constexpr NetworKit::count kBucketStride = kNumSingleBuckets << kNumIntermediateBuckets; + +using Storage = std::array, kNumBuckets>; +inline NetworKit::count mostSignificantBit(NetworKit::count value) { + NetworKit::count result = 0; + while (value >>= 1) { + ++result; + } + return result; +} + +template +/** Insert a candidate at its current degree; candidates outside the window are ignored. */ +void add(Buckets &target, NetworKit::node u, NetworKit::count degree, NetworKit::count baseLevel) { + if (degree < baseLevel || degree >= baseLevel + kBucketStride) { + return; + } + + if (degree < baseLevel + kNumSingleBuckets) { + target[degree & kBucketMask].push_back(u); + return; + } + + const NetworKit::count differingBit = mostSignificantBit(degree ^ baseLevel); + const NetworKit::count bucket = kNumSingleBuckets + differingBit - kLog2SingleBuckets; + target[bucket].push_back(u); +} + +template +/** Record a decrement only when it crosses a bucket boundary (or reaches an exact bucket). */ +void move(Buckets &target, NetworKit::node u, NetworKit::count degree, NetworKit::count baseLevel) { + if (degree < baseLevel || degree >= baseLevel + kBucketStride) { + return; + } + + if (degree < baseLevel + kNumSingleBuckets) { + target[degree & kBucketMask].push_back(u); + return; + } + + const NetworKit::count differingBit = mostSignificantBit(degree ^ baseLevel); + const NetworKit::count previousDifferingBit = mostSignificantBit((degree + 1) ^ baseLevel); + if (differingBit != previousDifferingBit) { + const NetworKit::count bucket = kNumSingleBuckets + differingBit - kLog2SingleBuckets; + target[bucket].push_back(u); + } +} + +} // namespace CoreBuckets +} // namespace Aux +#endif // NETWORKIT_AUXILIARY_HIERARCHICAL_CORE_BUCKETS_HPP_ diff --git a/include/networkit/centrality/ParallelCoreDecomposition.hpp b/include/networkit/centrality/ParallelCoreDecomposition.hpp new file mode 100644 index 000000000..736fd747e --- /dev/null +++ b/include/networkit/centrality/ParallelCoreDecomposition.hpp @@ -0,0 +1,92 @@ +/* + * ParallelCoreDecomposition.hpp + * + * Exact shared-memory k-core decomposition based on the SIGMOD 2025 algorithm. + */ + +#ifndef NETWORKIT_CENTRALITY_PARALLEL_CORE_DECOMPOSITION_HPP_ +#define NETWORKIT_CENTRALITY_PARALLEL_CORE_DECOMPOSITION_HPP_ + +#include +#include +#include +#include +#include + +namespace NetworKit { + +/** + * @ingroup centrality + * Computes an exact k-core decomposition in parallel. + * + * The implementation is based on the online peeling framework, hierarchical bucketing structure, + * probabilistic high-degree sampling, and vertical granularity control described by Liu et al. in + * "Parallel k-Core Decomposition: Theory and Practice" (Proc. ACM Manag. Data 3, 3, SIGMOD + * 2025). Sampling errors are detected and cause a restart without sampling, making the + * implementation Las Vegas rather than Monte Carlo. + * Low-average-degree graphs with maximum degree below the edge-parallelism threshold use + * the existing ParK peeling implementation to avoid sampling/bucketing overhead. Selection + * is automatic and preserves exact core numbers and normalization semantics. + * + * The input graph must be undirected, must not contain self-loops, and must have continuous node + * ids. Both GraphW and zero-copy CSR-backed GraphR inputs are supported. + */ +class ParallelCoreDecomposition final : public Centrality { +public: + /** + * Create a parallel core decomposition for @a G. + * + * @param G The undirected input graph. + * @param normalized If @c true, divide scores by the maximum degree. Partitions and + * covers always use the original integer core numbers. + * @param useSampling If @c true, use the paper's high-degree sampling optimization. + */ + ParallelCoreDecomposition(const Graph &G, bool normalized = false, bool useSampling = true); + + /** Compute the core number of every node. */ + void run() override; + + /** + * Get the nested k-cores as a cover. + * + * @return The k-cores as a Cover. + */ + Cover getCover() const; + + /** + * Get the k-shells as a partition. + * + * @return The k-shells as a Partition. + */ + Partition getPartition() const; + + /** @return The maximum core number. */ + index maxCoreNumber() const; + + /** + * Get the number of sampling-error recovery restarts performed by the last run. + * + * @return Zero unless a sampling validation failed and the algorithm restarted without + * sampling. + */ + count numberOfRestarts() const; + + /** + * @return The theoretical maximum score: 1 for normalized scores, n - 1 otherwise. + * Empty and single-node graphs return 0. + */ + double maximum() override; + +private: + index coreNumber(node u) const; + + // Preserve integer membership only when scoreData is normalized. + std::vector coreNumbers_; + index maxCore{0}; + bool useSampling; + count restarts{0}; +}; + +} // namespace NetworKit + +#endif // NETWORKIT_CENTRALITY_PARALLEL_CORE_DECOMPOSITION_HPP_ diff --git a/networkit/centrality.pyx b/networkit/centrality.pyx index e3c94bc4d..89f3d1b70 100644 --- a/networkit/centrality.pyx +++ b/networkit/centrality.pyx @@ -2276,6 +2276,58 @@ cdef class CoreDecomposition(Centrality): """ return (<_CoreDecomposition*>(self._this)).getNodeOrder() +cdef extern from "": + + cdef cppclass _ParallelCoreDecomposition "NetworKit::ParallelCoreDecomposition" (_Centrality): + _ParallelCoreDecomposition(_Graph, bool_t, bool_t) except + + _Cover getCover() except + + _Partition getPartition() except + + index maxCoreNumber() except + + count numberOfRestarts() except + + +cdef class ParallelCoreDecomposition(Centrality): + """ + ParallelCoreDecomposition(G, normalized=False, useSampling=True) + + Computes an exact parallel k-core decomposition of an undirected graph using the SIGMOD 2025 + online peeling algorithm with hierarchical buckets, high-degree sampling, and vertical + granularity control. If sampling validation detects an error, the algorithm restarts without + sampling to preserve exact results. + Graphs with low average and maximum degree automatically use the existing ParK peeling + implementation to avoid bucketing and sampling overhead. + + Parameters + ---------- + G : networkit.Graph + The undirected input graph. Node ids must be continuous and the graph may not contain + self-loops. + normalized : bool, optional + Divide each core number by the maximum degree. Default: False + useSampling : bool, optional + Use high-degree sampling to reduce atomic contention. Default: True + """ + + def __cinit__(self, Graph G, bool_t normalized=False, bool_t useSampling=True): + self._G = G + self._this = new _ParallelCoreDecomposition( + dereference(G._view()), normalized, useSampling) + + def maxCoreNumber(self): + """Return the maximum core number.""" + return (<_ParallelCoreDecomposition*>(self._this)).maxCoreNumber() + + def numberOfRestarts(self): + """Return the number of sampling-error recovery restarts in the last run.""" + return (<_ParallelCoreDecomposition*>(self._this)).numberOfRestarts() + + def getCover(self): + """Return the nested k-cores as a cover.""" + return Cover().setThis((<_ParallelCoreDecomposition*>(self._this)).getCover()) + + def getPartition(self): + """Return the k-shells as a partition.""" + return Partition().setThis((<_ParallelCoreDecomposition*>(self._this)).getPartition()) + cdef extern from "": cdef cppclass _EigenvectorCentrality "NetworKit::EigenvectorCentrality" (_Centrality): diff --git a/networkit/cpp/centrality/CMakeLists.txt b/networkit/cpp/centrality/CMakeLists.txt index aa4622f07..00e8d365b 100644 --- a/networkit/cpp/centrality/CMakeLists.txt +++ b/networkit/cpp/centrality/CMakeLists.txt @@ -35,6 +35,7 @@ networkit_add_module(centrality LocalPartitionCoverage.cpp LocalSquareClusteringCoefficient.cpp PageRank.cpp + ParallelCoreDecomposition.cpp PermanenceCentrality.cpp Sfigality.cpp SpanningEdgeCentrality.cpp diff --git a/networkit/cpp/centrality/CoreSampler.hpp b/networkit/cpp/centrality/CoreSampler.hpp new file mode 100644 index 000000000..72be93b7f --- /dev/null +++ b/networkit/cpp/centrality/CoreSampler.hpp @@ -0,0 +1,47 @@ +#ifndef NETWORKIT_CENTRALITY_CORE_SAMPLER_HPP_ +#define NETWORKIT_CENTRALITY_CORE_SAMPLER_HPP_ +#include +#include +#include +#include +#include + +namespace NetworKit { +namespace CoreDecompositionDetail { +// reset() is called only between parallel sampling phases. +class Sampler final { +public: + void reset(count newExpectedHits, double newRate) { + expectedHits = newExpectedHits; + rate = std::min(1.0, newRate); + threshold = + static_cast(rate * static_cast(std::numeric_limits::max())); + hits.store(0, std::memory_order_relaxed); + } + + bool sample(uint32_t randomNumber) { + count current = hits.load(std::memory_order_relaxed); + if (current >= expectedHits || randomNumber >= threshold) { + return false; + } + + current = hits.fetch_add(1, std::memory_order_relaxed); + return current + 1 == expectedHits; + } + + count numberOfHits() const { return hits.load(std::memory_order_relaxed); } + + count numberOfExpectedHits() const { return expectedHits; } + + double sampleRate() const { return rate; } + +private: + std::atomic hits{0}; + count expectedHits{0}; + uint32_t threshold{0}; + double rate{0.0}; +}; + +} // namespace CoreDecompositionDetail +} // namespace NetworKit +#endif diff --git a/networkit/cpp/centrality/ParallelCoreDecomposition.cpp b/networkit/cpp/centrality/ParallelCoreDecomposition.cpp new file mode 100644 index 000000000..a5cca7f97 --- /dev/null +++ b/networkit/cpp/centrality/ParallelCoreDecomposition.cpp @@ -0,0 +1,87 @@ +/* + * ParallelCoreDecomposition.cpp + * + * Exact shared-memory k-core decomposition based on Liu et al., SIGMOD 2025. + */ + +#include "ParallelCoreDecompositionInternal.hpp" + +namespace NetworKit { + +ParallelCoreDecomposition::ParallelCoreDecomposition(const Graph &G, bool normalized, + bool useSampling) + : Centrality(G, normalized), useSampling(useSampling) { + if (G.isDirected()) { + throw std::runtime_error("ParallelCoreDecomposition requires an undirected graph"); + } + if (G.numberOfSelfLoops()) { + throw std::runtime_error("ParallelCoreDecomposition does not support self-loops; call " + "Graph.removeSelfLoops() first"); + } + if (G.numberOfNodes() != G.upperNodeIdBound()) { + throw std::runtime_error("ParallelCoreDecomposition requires continuous node ids"); + } +} + +void ParallelCoreDecomposition::run() { + const auto result = CoreDecompositionDetail::runWithRecovery(G, scoreData, useSampling); + restarts = result.restarts; + maxCore = result.maxCore; + + if (normalized) { + coreNumbers_.resize(scoreData.size()); + count maximumDegree = 0; + G.forNodes([&](node u) { maximumDegree = std::max(maximumDegree, G.degree(u)); }); + G.parallelForNodes([&](node u) { + coreNumbers_[u] = static_cast(scoreData[u]); + if (maximumDegree > 0) + scoreData[u] /= maximumDegree; + }); + } + + hasRun = true; +} + +Cover ParallelCoreDecomposition::getCover() const { + assureFinished(); + Cover result(G.upperNodeIdBound()); + result.setUpperBound(G.upperNodeIdBound()); + G.forNodes([&](node u) { + for (index core = 0; core <= coreNumber(u); ++core) { + result.addToSubset(core, u); + } + }); + return result; +} + +Partition ParallelCoreDecomposition::getPartition() const { + assureFinished(); + Partition result(G.upperNodeIdBound()); + result.allToSingletons(); + G.forNodes([&](node u) { result.moveToSubset(coreNumber(u), u); }); + return result; +} + +index ParallelCoreDecomposition::coreNumber(node u) const { + return normalized ? coreNumbers_[u] : static_cast(scoreData[u]); +} + +index ParallelCoreDecomposition::maxCoreNumber() const { + assureFinished(); + return maxCore; +} + +count ParallelCoreDecomposition::numberOfRestarts() const { + assureFinished(); + return restarts; +} + +double ParallelCoreDecomposition::maximum() { + const count n = G.numberOfNodes(); + if (n <= 1) { + return 0.0; + } + return normalized ? 1.0 : static_cast(n - 1); +} + +} // namespace NetworKit diff --git a/networkit/cpp/centrality/ParallelCoreDecompositionInternal.hpp b/networkit/cpp/centrality/ParallelCoreDecompositionInternal.hpp new file mode 100644 index 000000000..2b8f1c7b0 --- /dev/null +++ b/networkit/cpp/centrality/ParallelCoreDecompositionInternal.hpp @@ -0,0 +1,728 @@ +#ifndef NETWORKIT_CENTRALITY_PARALLEL_CORE_INTERNAL_HPP_ +#define NETWORKIT_CENTRALITY_PARALLEL_CORE_INTERNAL_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "CoreSampler.hpp" +#include +#include +#include + +namespace NetworKit { +namespace CoreDecompositionDetail { + +constexpr count kSampleThreshold = 2000; +constexpr double kInitialReduceRatio = 0.1; +constexpr count kLog2ErrorFactor = 32; +constexpr count kBaselineExpectedHits = + static_cast(kLog2ErrorFactor / (kInitialReduceRatio * kInitialReduceRatio)); +constexpr double kBiasFactor = 0.5; +constexpr double kErrorRateTolerance = 1e-10; + +constexpr count kLocalQueueSize = 128; +constexpr count kNeighborParallelismThreshold = 128; +constexpr count kParallelFrontierThreshold = 256; +constexpr count kHbsSwitchCore = 16; + +using BucketArray = Aux::CoreBuckets::Storage; +using namespace Aux::CoreBuckets; + +struct ThreadScratch final { + BucketArray buckets; + std::vector recounts; +}; + +struct EdgeBlock final { + node source; + index begin; + index end; +}; + +inline uint32_t hash32(uint32_t value) { + value = (value + 0x7ed55d16U) + (value << 12); + value = (value ^ 0xc761c23cU) ^ (value >> 19); + value = (value + 0x165667b1U) + (value << 5); + value = (value + 0xd3a2646cU) ^ (value << 9); + value = (value + 0xfd7046c5U) + (value << 3); + value = (value ^ 0xb55a4f09U) ^ (value >> 16); + return value; +} + +inline count log2Up(count value) { + if (value <= 1) { + return 0; + } + + count result = 0; + --value; + while (value > 0) { + ++result; + value >>= 1; + } + return result; +} + +class ParallelKCoreSolver final { +public: + ParallelKCoreSolver(const Graph &G, std::vector &scores, bool samplingEnabled, + bool (*acceptRecount)(node, count) = nullptr) + : G(G), scores(scores), samplingEnabled(samplingEnabled), acceptRecount(acceptRecount), + z(G.upperNodeIdBound()), degrees(new std::atomic[z]), + active(new std::atomic[z]), threadScratch(omp_get_max_threads()), + filterCounts(omp_get_max_threads()), filterOffsets(omp_get_max_threads()), + seenAtEpoch(z, 0) {} + + bool run() { + const count n = G.numberOfNodes(); + if (n == 0) { + return true; + } + + initialize(); + + std::vector remaining(z); +#pragma omp parallel for schedule(static) + for (omp_index u = 0; u < static_cast(z); ++u) { + remaining[u] = static_cast(u); + } + + const count averageDegree = G.numberOfEdges() * 2 / n; + count firstHbsLevel = 0; + if (averageDegree < kHbsSwitchCore) { + for (count level = 0; level < kHbsSwitchCore && !remaining.empty(); ++level) { + buckets[0].clear(); + if (!buildSparseBucket(remaining, level)) { + return false; + } + if (!processLevel(level, level, false)) { + return false; + } + firstHbsLevel = level + 1; + } + filterRemaining(remaining); + } + + for (count baseLevel = 0; !remaining.empty(); baseLevel += kBucketStride) { + if (baseLevel > maximumInitialDegree) { + return false; + } + + clearBuckets(); + const count levelBegin = std::max(baseLevel, firstHbsLevel); + if (!buildBuckets(remaining, levelBegin, baseLevel + kBucketStride, baseLevel, true)) { + return false; + } + + for (count level = levelBegin; level < baseLevel + kBucketStride && !remaining.empty(); + ++level) { + if (level != baseLevel && (level & kBucketMask) == 0) { + redistributeIntermediateBucket(level); + } + + const count bucketBase = level & ~kBucketMask; + if (!processLevel(level, bucketBase, true)) { + return false; + } + } + + filterRemaining(remaining); + } + + return true; + } + + index maxCoreNumber() const { return maxCore; } + +private: + const Graph &G; + std::vector &scores; + bool samplingEnabled; + // Per-run internal validation hook; never shared globally or exposed through the public API. + bool (*acceptRecount)(node, count); + count z; + std::unique_ptr[]> degrees; + std::unique_ptr samplers; + std::unique_ptr[]> recountQueued; + // Atomics protect accesses within a peeling round. Relaxed ordering is sufficient: + // OpenMP region/barrier boundaries publish sampler configuration and scratch buffers + // between initialization, peeling and recount phases. No flag publishes other data. + std::unique_ptr[]> active; + std::unique_ptr[]> sampleMode; + BucketArray buckets; + std::vector threadScratch; + std::vector filterCounts; + std::vector filterOffsets; + std::vector filterBuffer; + std::vector seenAtEpoch; + count scratchInUse{1}; + bool samplingPossible{false}; + count currentEpoch{0}; + index maxCore{0}; + count maximumInitialDegree{0}; + + void initialize() { + count maxDegree = 0; +// MSVC classic OpenMP (NETWORKIT_OMP2) lacks reduction(max); fall back to serial. +#ifndef NETWORKIT_OMP2 +#pragma omp parallel for reduction(max : maxDegree) schedule(static) +#endif + for (omp_index u = 0; u < static_cast(z); ++u) { + const count degree = G.degree(static_cast(u)); + degrees[u].store(degree, std::memory_order_relaxed); + active[u].store(1, std::memory_order_relaxed); + maxDegree = std::max(maxDegree, degree); + } + maximumInitialDegree = maxDegree; + samplingPossible = samplingEnabled && maxDegree * kInitialReduceRatio >= kSampleThreshold; + + if (samplingPossible) { + // Most sparse inputs never sample. Avoid O(n) sampler construction and storage + // until the degree scan establishes that sampling can actually be used. + samplers.reset(new Sampler[z]); + recountQueued.reset(new std::atomic[z]); + sampleMode.reset(new std::atomic[z]); +#pragma omp parallel for schedule(static) + for (omp_index u = 0; u < static_cast(z); ++u) { + recountQueued[u].store(false, std::memory_order_relaxed); + setSampler(static_cast(u), 0); + } + } + } + + void setSampler(node u, count level) { + if (!samplingPossible) { + return; + } + + const count degree = degrees[u].load(std::memory_order_relaxed); + const double reducedDegree = degree * kInitialReduceRatio; + const double remainingAfterReduction = + (degree - std::min(degree, level)) * (1.0 - kInitialReduceRatio); + + if (reducedDegree >= kSampleThreshold && level < reducedDegree * kBiasFactor + && kBaselineExpectedHits < remainingAfterReduction) { + sampleMode[u].store(1, std::memory_order_relaxed); + const count distance = degree - level; + const count logarithm = log2Up(distance); + const count expectedHits = kLog2ErrorFactor * logarithm * logarithm / 2; + const double rate = expectedHits / ((1.0 - kInitialReduceRatio) * degree); + samplers[u].reset(expectedHits, rate); + } else { + sampleMode[u].store(0, std::memory_order_relaxed); + } + } + + bool sampleIsSafe(node u, count horizon) const { + const count degree = degrees[u].load(std::memory_order_relaxed); + if (degree * kInitialReduceRatio * kBiasFactor < horizon || degree <= horizon) { + return false; + } + + const double rate = samplers[u].sampleRate(); + const double distance = static_cast(degree - horizon); + const double expected = distance * rate; + if (expected <= 0.0) { + return false; + } + + const double hits = static_cast(std::max(1, samplers[u].numberOfHits())); + const double exponent = -expected + 2.0 * hits - hits * hits / expected; + return std::exp(exponent) < kErrorRateTolerance; + } + + void resetScratch(count numberOfThreads) { + for (count thread = 0; thread < numberOfThreads; ++thread) { + auto &scratch = threadScratch[thread]; + scratch.recounts.clear(); + for (auto &bucket : scratch.buckets) { + bucket.clear(); + } + } + } + + void mergeScratchBuckets(count numberOfThreads) { + for (count thread = 0; thread < numberOfThreads; ++thread) { + auto &scratch = threadScratch[thread]; + for (count bucket = 0; bucket < kNumBuckets; ++bucket) { + auto &source = scratch.buckets[bucket]; + buckets[bucket].insert(buckets[bucket].end(), source.begin(), source.end()); + } + } + } + + std::vector collectScratchRecounts(count numberOfThreads) { + count size = 0; + for (count thread = 0; thread < numberOfThreads; ++thread) { + size += threadScratch[thread].recounts.size(); + } + + std::vector recounts; + recounts.reserve(size); + for (count thread = 0; thread < numberOfThreads; ++thread) { + const auto &source = threadScratch[thread].recounts; + recounts.insert(recounts.end(), source.begin(), source.end()); + } + return recounts; + } + + void clearBuckets() { + for (auto &bucket : buckets) { + bucket.clear(); + } + } + + bool buildBuckets(const std::vector &remaining, count level, count horizon, + count baseLevel, bool useHbs) { + resetScratch(threadScratch.size()); +#pragma omp parallel for schedule(static) + for (omp_index i = 0; i < static_cast(remaining.size()); ++i) { + const node u = remaining[i]; + const count degree = degrees[u].load(std::memory_order_relaxed); + auto &scratch = threadScratch[omp_get_thread_num()]; + Aux::CoreBuckets::add(scratch.buckets, u, degree, baseLevel); + if (samplingPossible && sampleMode[u].load(std::memory_order_relaxed) + && !sampleIsSafe(u, horizon)) { + bool expected = false; + if (recountQueued[u].compare_exchange_strong(expected, true, + std::memory_order_relaxed)) { + scratch.recounts.push_back(u); + } + } + } + mergeScratchBuckets(threadScratch.size()); + return recountVertices(collectScratchRecounts(threadScratch.size()), level, baseLevel, + useHbs); + } + + bool buildSparseBucket(std::vector &remaining, count level) { + resetScratch(threadScratch.size()); + // Compact the previous level and build the next frontier in the same scan. + filterRemaining(remaining, true, level); + mergeScratchBuckets(threadScratch.size()); + return recountVertices(collectScratchRecounts(threadScratch.size()), level, level, false); + } + + std::vector uniqueActive(std::vector candidates) { + ++currentEpoch; + if (currentEpoch == 0) { + std::fill(seenAtEpoch.begin(), seenAtEpoch.end(), 0); + ++currentEpoch; + } + + size_t output = 0; + for (node u : candidates) { + if (active[u].load(std::memory_order_relaxed) && seenAtEpoch[u] != currentEpoch) { + seenAtEpoch[u] = currentEpoch; + candidates[output++] = u; + } + } + candidates.resize(output); + return candidates; + } + + std::vector takeExactBucket(count level, bool useHbs) { + const count bucket = useHbs ? (level & kBucketMask) : 0; + std::vector candidates; + candidates.swap(buckets[bucket]); + candidates = uniqueActive(std::move(candidates)); + + candidates.erase(std::remove_if(candidates.begin(), candidates.end(), + [&](node u) { + return degrees[u].load(std::memory_order_relaxed) + != level; + }), + candidates.end()); + return candidates; + } + + void redistributeIntermediateBucket(count level) { + count intermediate = kNumIntermediateBuckets; + for (count i = kNumIntermediateBuckets; i-- > 0;) { + const count mask = (kNumSingleBuckets << i) - 1; + if ((level & mask) == 0) { + intermediate = i; + break; + } + } + if (intermediate == kNumIntermediateBuckets) { + return; + } + + const count bucket = kNumSingleBuckets + intermediate; + std::vector candidates; + candidates.swap(buckets[bucket]); + candidates = uniqueActive(std::move(candidates)); + + resetScratch(threadScratch.size()); +#pragma omp parallel for schedule(static) + for (omp_index i = 0; i < static_cast(candidates.size()); ++i) { + const node u = candidates[i]; + const count degree = degrees[u].load(std::memory_order_relaxed); + Aux::CoreBuckets::add(threadScratch[omp_get_thread_num()].buckets, u, degree, level); + } + mergeScratchBuckets(threadScratch.size()); + } + + void decrementNeighbor(node source, node target, count level, count bucketBase, bool useHbs, + bool allowLocalQueue, std::array &localQueue, + count &rear, BucketArray &localBuckets, + std::vector &localRecounts) { + count oldDegree = degrees[target].load(std::memory_order_relaxed); + if (oldDegree <= level) { + return; + } + + if (samplingPossible && sampleMode[target].load(std::memory_order_relaxed)) { + // Unsigned multiplication wraps and the 32-bit key intentionally collides. + // Exact recounts and recovery protect correctness from sampling collisions. + const uint32_t key = static_cast(source * z + target); + if (samplers[target].sample(hash32(key))) { + bool expected = false; + if (recountQueued[target].compare_exchange_strong(expected, true, + std::memory_order_relaxed)) { + localRecounts.push_back(target); + } + } + return; + } + + while (oldDegree > level) { + if (degrees[target].compare_exchange_weak(oldDegree, oldDegree - 1, + std::memory_order_relaxed)) { + const count newDegree = oldDegree - 1; + if (newDegree == level && allowLocalQueue && rear < kLocalQueueSize) { + localQueue[rear++] = target; + } else if (useHbs) { + Aux::CoreBuckets::move(localBuckets, target, newDegree, bucketBase); + } else if (newDegree == level) { + localBuckets[0].push_back(target); + } + return; + } + } + } + + count processRootSequential(node root, count level, count bucketBase, bool useHbs, + bool deferHighDegree, BucketArray &localBuckets, + std::vector &localRecounts) { + std::array localQueue; + count front = 0; + count rear = 1; + count processed = 0; + localQueue[0] = root; + + while (front < rear) { + const node u = localQueue[front++]; + // In the serial small-frontier path, continue a chain without another round. + // Parallel roots retain a bounded work budget so one worker cannot monopolize + // a long chain. High-degree vertices are still deferred to edge blocks. + if (deferHighDegree && front == rear) { + front = 0; + rear = 0; + } + if (!active[u].load(std::memory_order_relaxed)) { + continue; + } + + // Use original degree: even a low residual degree still scans the full CSR row. + if (deferHighDegree && G.degree(u) >= kNeighborParallelismThreshold) { + const count bucket = useHbs ? (level & kBucketMask) : 0; + localBuckets[bucket].push_back(u); + continue; + } + + active[u].store(0, std::memory_order_relaxed); + scores[u] = level; + ++processed; + G.forNeighborsOf(u, [&](node v) { + decrementNeighbor(u, v, level, bucketBase, useHbs, true, localQueue, rear, + localBuckets, localRecounts); + }); + } + return processed; + } + + count processLargeFrontier(const std::vector &frontier, count level, count bucketBase, + bool useHbs) { + scratchInUse = threadScratch.size(); + resetScratch(scratchInUse); + count processed = 0; +#pragma omp parallel for schedule(dynamic, 64) reduction(+ : processed) + for (omp_index i = 0; i < static_cast(frontier.size()); ++i) { + auto &scratch = threadScratch[omp_get_thread_num()]; + processed += processRootSequential(frontier[i], level, bucketBase, useHbs, false, + scratch.buckets, scratch.recounts); + } + mergeScratchBuckets(scratchInUse); + return processed; + } + + count processSmallFrontier(const std::vector &frontier, count level, count bucketBase, + bool useHbs) { + std::vector highDegreeRoots; + for (node root : frontier) { + if (G.degree(root) >= kNeighborParallelismThreshold) { + highDegreeRoots.push_back(root); + } + } + + scratchInUse = highDegreeRoots.empty() ? 1 : threadScratch.size(); + const bool directBuckets = highDegreeRoots.empty(); + if (directBuckets) { + threadScratch[0].recounts.clear(); + } else { + resetScratch(scratchInUse); + } + count processed = 0; + for (node root : frontier) { + if (G.degree(root) < kNeighborParallelismThreshold) { + processed += processRootSequential( + root, level, bucketBase, useHbs, true, + directBuckets ? buckets : threadScratch[0].buckets, threadScratch[0].recounts); + } + } + + std::vector edgeBlocks; + for (node root : highDegreeRoots) { + if (!active[root].load(std::memory_order_relaxed)) { + continue; + } + + active[root].store(0, std::memory_order_relaxed); + scores[root] = level; + ++processed; + const count degree = G.degree(root); + for (index begin = 0; begin < degree; begin += kNeighborParallelismThreshold) { + edgeBlocks.push_back( + {root, begin, std::min(degree, begin + kNeighborParallelismThreshold)}); + } + } + + if (!edgeBlocks.empty()) { +#pragma omp parallel + { + auto &scratch = threadScratch[omp_get_thread_num()]; + std::array unusedQueue; + count unusedRear = 0; +#pragma omp for schedule(dynamic, 8) + for (omp_index blockIndex = 0; + blockIndex < static_cast(edgeBlocks.size()); ++blockIndex) { + const auto &block = edgeBlocks[blockIndex]; + for (index neighborIndex = block.begin; neighborIndex < block.end; + ++neighborIndex) { + const node v = G.getIthNeighbor(Unsafe{}, block.source, neighborIndex); + decrementNeighbor(block.source, v, level, bucketBase, useHbs, false, + unusedQueue, unusedRear, scratch.buckets, + scratch.recounts); + } + } + } + } + + if (!directBuckets) + mergeScratchBuckets(scratchInUse); + return processed; + } + + bool processLevel(count level, count bucketBase, bool useHbs) { + std::vector frontier = takeExactBucket(level, useHbs); + while (!frontier.empty()) { + const count processed = frontier.size() > kParallelFrontierThreshold + ? processLargeFrontier(frontier, level, bucketBase, useHbs) + : processSmallFrontier(frontier, level, bucketBase, useHbs); + if (processed > 0) { + maxCore = std::max(maxCore, level); + } + + const std::vector recounts = collectScratchRecounts(scratchInUse); + if (!recounts.empty() && !recountVertices(recounts, level, bucketBase, useHbs)) { + return false; + } + frontier = takeExactBucket(level, useHbs); + } + return true; + } + + bool recountVertices(const std::vector &vertices, count level, count bucketBase, + bool useHbs) { + if (vertices.empty()) { + return true; + } + + std::atomic samplingError{false}; + resetScratch(threadScratch.size()); + +#pragma omp parallel for schedule(dynamic, 1) if (vertices.size() > 1) + for (omp_index i = 0; i < static_cast(vertices.size()); ++i) { + const node u = vertices[i]; + recountQueued[u].store(false, std::memory_order_relaxed); + if (!active[u].load(std::memory_order_relaxed)) { + continue; + } + + if (samplingEnabled && acceptRecount && !acceptRecount(u, level)) { + samplingError.store(true, std::memory_order_relaxed); + continue; + } + + count exactDegree = 0; + G.forNeighborsOf( + u, [&](node v) { exactDegree += active[v].load(std::memory_order_relaxed) != 0; }); + + if (exactDegree < level) { + count previousRoundDegree = 0; + G.forNeighborsOf(u, [&](node v) { + previousRoundDegree += active[v].load(std::memory_order_relaxed) + || degrees[v].load(std::memory_order_relaxed) == level; + }); + + if (previousRoundDegree < level) { + samplingError.store(true, std::memory_order_relaxed); + continue; + } + exactDegree = level; + } + + degrees[u].store(exactDegree, std::memory_order_relaxed); + auto &scratch = threadScratch[omp_get_thread_num()]; + if (exactDegree == level) { + const count bucket = useHbs ? (level & kBucketMask) : 0; + scratch.buckets[bucket].push_back(u); + } else if (useHbs) { + Aux::CoreBuckets::add(scratch.buckets, u, exactDegree, bucketBase); + } + setSampler(u, level); + } + + mergeScratchBuckets(threadScratch.size()); + return !samplingError.load(std::memory_order_relaxed); + } + + void filterRemaining(std::vector &remaining, bool buildSparse = false, count level = 0) { + if (remaining.empty()) + return; + const count inputSize = remaining.size(); + filterBuffer.resize(inputSize); + count outputSize = 0; + +#pragma omp parallel shared(outputSize) + { + const int thread = omp_get_thread_num(); + const int threadCount = omp_get_num_threads(); + const count begin = inputSize * thread / threadCount; + const count end = inputSize * (thread + 1) / threadCount; + count localCount = 0; + for (count i = begin; i < end; ++i) { + const node u = remaining[i]; + if (!active[u].load(std::memory_order_relaxed)) + continue; + ++localCount; + if (buildSparse) { + auto &scratch = threadScratch[thread]; + if (degrees[u].load(std::memory_order_relaxed) == level) { + scratch.buckets[0].push_back(u); + } + if (samplingPossible && sampleMode[u].load(std::memory_order_relaxed) + && !sampleIsSafe(u, level + kBucketStride)) { + bool expected = false; + if (recountQueued[u].compare_exchange_strong(expected, true, + std::memory_order_relaxed)) { + scratch.recounts.push_back(u); + } + } + } + } + filterCounts[thread] = localCount; + +#pragma omp barrier +#pragma omp single + { + outputSize = 0; + for (int owner = 0; owner < threadCount; ++owner) { + filterOffsets[owner] = outputSize; + outputSize += filterCounts[owner]; + } + } +#pragma omp barrier + + count output = filterOffsets[thread]; + for (count i = begin; i < end; ++i) { + const node u = remaining[i]; + if (active[u].load(std::memory_order_relaxed)) { + filterBuffer[output++] = u; + } + } + } + + filterBuffer.resize(outputSize); + remaining.swap(filterBuffer); + } +}; + +struct CoreResult { + index maxCore; + count restarts; +}; + +// Reuse the established sparse peeling implementation without copying its score vector. +class SparseCoreSolver final : public CoreDecomposition { +public: + explicit SparseCoreSolver(const Graph &graph) : CoreDecomposition(graph, false, false) {} + + void runInto(std::vector &scores) { + scoreData.swap(scores); + run(); + scoreData.swap(scores); + } +}; + +inline bool preferSparsePeeling(const Graph &graph) { + const count n = graph.numberOfNodes(); + if (n == 0 || graph.numberOfEdges() / n >= kHbsSwitchCore / 2) + return false; + count maximumDegree = 0; +// MSVC classic OpenMP (NETWORKIT_OMP2) lacks reduction(max); fall back to serial. +#ifndef NETWORKIT_OMP2 +#pragma omp parallel for reduction(max : maximumDegree) schedule(static) if (n > 256) +#endif + for (omp_index u = 0; u < static_cast(n); ++u) { + maximumDegree = std::max(maximumDegree, graph.degree(static_cast(u))); + } + // Preserve sampling/edge-block peeling for stars and dense components in sparse graphs. + return maximumDegree < kNeighborParallelismThreshold; +} + +// Shared by the public algorithm and white-box recovery tests. +inline CoreResult runWithRecovery(const Graph &G, std::vector &scores, bool sampling, + bool (*acceptRecount)(node, count) = nullptr) { + if (!acceptRecount && preferSparsePeeling(G)) { + SparseCoreSolver sparse(G); + sparse.runInto(scores); + return {sparse.maxCoreNumber(), 0}; + } + scores.assign(G.upperNodeIdBound(), 0.0); + ParallelKCoreSolver solver(G, scores, sampling, acceptRecount); + if (solver.run()) { + return {solver.maxCoreNumber(), 0}; + } + scores.assign(G.upperNodeIdBound(), 0.0); + ParallelKCoreSolver exactSolver(G, scores, false); + if (!exactSolver.run()) { + throw std::runtime_error("ParallelCoreDecomposition exact recovery run failed"); + } + return {exactSolver.maxCoreNumber(), 1}; +} + +} // namespace CoreDecompositionDetail +} // namespace NetworKit +#endif diff --git a/networkit/cpp/centrality/test/CentralityGTest.cpp b/networkit/cpp/centrality/test/CentralityGTest.cpp index 560fd2131..4b2df8e64 100644 --- a/networkit/cpp/centrality/test/CentralityGTest.cpp +++ b/networkit/cpp/centrality/test/CentralityGTest.cpp @@ -10,7 +10,11 @@ #include #include +#include +#include "../ParallelCoreDecompositionInternal.hpp" +#include #include +#include #include #include @@ -46,6 +50,7 @@ #include #include #include +#include #include #include #include @@ -1523,6 +1528,222 @@ TEST_F(CentralityGTest, testCoreDecomposition) { EXPECT_ANY_THROW(CoreDecomposition CoreDec(H)); } +TEST_F(CentralityGTest, testParallelCoreDecompositionMatchesSequentialOracle) { + Aux::Random::setSeed(42, false); + GraphW G = ErdosRenyiGenerator(1000, 0.02, false).generate(); + + CoreDecomposition reference(G, false, true); + reference.run(); + + ParallelCoreDecomposition parallel(G); + parallel.run(); + + EXPECT_EQ(reference.scores(), parallel.scores()); + EXPECT_EQ(reference.maxCoreNumber(), parallel.maxCoreNumber()); + EXPECT_LE(parallel.numberOfRestarts(), 1u); +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionHierarchicalBuckets) { + GraphW G = ErdosRenyiGenerator(520, 1.0, false).generate(); + ParallelCoreDecomposition parallel(G); + parallel.run(); + + EXPECT_EQ(519u, parallel.maxCoreNumber()); + G.forNodes([&](node u) { EXPECT_EQ(519, parallel.score(u)); }); +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionSampling) { + GraphW G(25001); + for (node leaf = 1; leaf < G.numberOfNodes(); ++leaf) { + G.addEdge(0, leaf); + } + + ParallelCoreDecomposition parallel(G); + parallel.run(); + + EXPECT_EQ(1u, parallel.maxCoreNumber()); + G.forNodes([&](node u) { EXPECT_EQ(1, parallel.score(u)); }); + EXPECT_LE(parallel.numberOfRestarts(), 1u); +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionThreadMatrixAndGraphR) { + struct RestoreThreads { + int threads = omp_get_max_threads(); + ~RestoreThreads() { omp_set_num_threads(threads); } + } restore; + + // A high-degree center exercises sampling and the GraphR indexed-neighbor path. + GraphW writable(25001); + for (node u = 1; u < writable.numberOfNodes(); ++u) { + writable.addEdge(0, u); + } + arrow::UInt64Builder indicesBuilder, offsetsBuilder; + uint64_t offset = 0; + ASSERT_TRUE(offsetsBuilder.Append(offset).ok()); + writable.forNodes([&](node u) { + writable.forNeighborsOf(u, [&](node v) { + ASSERT_TRUE(indicesBuilder.Append(v).ok()); + ++offset; + }); + ASSERT_TRUE(offsetsBuilder.Append(offset).ok()); + }); + std::shared_ptr indices, offsets; + ASSERT_TRUE(indicesBuilder.Finish(&indices).ok()); + ASSERT_TRUE(offsetsBuilder.Finish(&offsets).ok()); + GraphR readonly(writable.numberOfNodes(), false, indices, offsets); + GraphW clique = ErdosRenyiGenerator(520, 1.0, false).generate(); + + for (int threads : {1, 2, 4, 15}) { + omp_set_num_threads(threads); + for (const Graph *graph : + {&static_cast(writable), &static_cast(readonly), + &static_cast(clique)}) { + for (bool sampling : {false, true}) { + SCOPED_TRACE(threads); + CoreDecomposition reference(*graph, false, true); + reference.run(); + ParallelCoreDecomposition parallel(*graph, false, sampling); + parallel.run(); + EXPECT_EQ(reference.scores(), parallel.scores()); + EXPECT_EQ(reference.maxCoreNumber(), parallel.maxCoreNumber()); + if (!sampling) + EXPECT_EQ(0u, parallel.numberOfRestarts()); + } + } + } +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionSamplingRecovery) { + GraphW graph(25001); + for (node u = 1; u < graph.numberOfNodes(); ++u) + graph.addEdge(0, u); + CoreDecomposition reference(graph, false, true); + reference.run(); + std::vector scores(graph.numberOfNodes(), -1.0); + // Accept initial recounts, then reject a sampled recount after level-one peeling. + const auto rejectAfterPeeling = +[](node, count level) { return level == 0; }; + const auto recovered = + CoreDecompositionDetail::runWithRecovery(graph, scores, true, rejectAfterPeeling); + EXPECT_EQ(1u, recovered.restarts); + EXPECT_EQ(reference.scores(), scores); + EXPECT_EQ(reference.maxCoreNumber(), recovered.maxCore); + // Exact execution must ignore sampling fault injection and start from fresh state. + const auto exact = + CoreDecompositionDetail::runWithRecovery(graph, scores, false, rejectAfterPeeling); + EXPECT_EQ(0u, exact.restarts); + EXPECT_EQ(reference.scores(), scores); +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionInternalComponents) { + Aux::CoreBuckets::Storage buckets; + for (count degree : {0u, 7u, 8u, 15u, 16u, 511u, 512u}) { + Aux::CoreBuckets::add(buckets, degree, degree, 0); + } + EXPECT_EQ((std::vector{0}), buckets[0]); + EXPECT_EQ((std::vector{7}), buckets[7]); + EXPECT_EQ((std::vector{8, 15}), buckets[8]); + EXPECT_EQ((std::vector{16}), buckets[9]); + EXPECT_EQ((std::vector{511}), buckets[13]); + Aux::CoreBuckets::Storage moved; + Aux::CoreBuckets::move(moved, 1, 15, 0); // Cross from [16, 32) to [8, 16). + Aux::CoreBuckets::move(moved, 2, 14, 0); // Stay in the same interval. + EXPECT_EQ((std::vector{1}), moved[8]); + Aux::CoreBuckets::add(moved, 3, 519, 512); + EXPECT_EQ((std::vector{3}), moved[7]); + + CoreDecompositionDetail::Sampler sampler; + sampler.reset(100, 1.0); + count triggers = 0; +#pragma omp parallel for reduction(+ : triggers) + for (int i = 0; i < 1000; ++i) + triggers += sampler.sample(0); + EXPECT_EQ(1u, triggers); + EXPECT_GE(sampler.numberOfHits(), 100u); + sampler.reset(10, 0.0); + EXPECT_FALSE(sampler.sample(0)); + EXPECT_EQ(0u, sampler.numberOfHits()); +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionSparseSelectionAndLongChains) { + struct RestoreThreads { + int threads = omp_get_max_threads(); + ~RestoreThreads() { omp_set_num_threads(threads); } + } restore; + // Disconnected long path, a cycle and isolates distinguish core numbers 0, 1 and 2. + GraphW graph(10010); + for (node u = 1; u < 10000; ++u) + graph.addEdge(u - 1, u); + graph.addEdge(10000, 10001); + graph.addEdge(10001, 10002); + graph.addEdge(10002, 10000); + EXPECT_TRUE(CoreDecompositionDetail::preferSparsePeeling(graph)); + CoreDecomposition reference(graph, false, true); + reference.run(); + for (int threads : {1, 2, 4, 15}) { + omp_set_num_threads(threads); + for (bool sampling : {false, true}) { + ParallelCoreDecomposition algorithm(graph, false, sampling); + algorithm.run(); + EXPECT_EQ(reference.scores(), algorithm.scores()); + EXPECT_EQ(2u, algorithm.maxCoreNumber()); + EXPECT_EQ(0u, algorithm.numberOfRestarts()); + // Also exercise the engine directly: serial queue reuse must terminate and agree + // with the oracle even when adaptive selection would normally use sparse peeling. + std::vector scores(graph.numberOfNodes(), 0.0); + CoreDecompositionDetail::ParallelKCoreSolver solver(graph, scores, sampling); + ASSERT_TRUE(solver.run()); + EXPECT_EQ(reference.scores(), scores); + } + } + GraphW star(1000); + for (node u = 1; u < star.numberOfNodes(); ++u) + star.addEdge(0, u); + EXPECT_FALSE(CoreDecompositionDetail::preferSparsePeeling(star)); + GraphW clique = ErdosRenyiGenerator(520, 1.0, false).generate(); + EXPECT_FALSE(CoreDecompositionDetail::preferSparsePeeling(clique)); +} + +TEST_F(CentralityGTest, testParallelCoreDecompositionNormalizedAndEmpty) { + for (count n : {0u, 1u, 4u}) { + GraphW graph(n); + if (n == 4) { + graph.addEdge(0, 1); + graph.addEdge(1, 2); + graph.addEdge(2, 0); + graph.addEdge(2, 3); + } + for (bool normalized : {false, true}) { + ParallelCoreDecomposition algorithm(graph, normalized); + EXPECT_EQ(n == 4 ? (normalized ? 1.0 : 3.0) : 0.0, algorithm.maximum()); + algorithm.run(); + const auto firstScores = algorithm.scores(); + algorithm.run(); + EXPECT_EQ(firstScores, algorithm.scores()); + EXPECT_EQ(n == 4 ? 2u : 0u, algorithm.maxCoreNumber()); + if (n == 4) { + EXPECT_NEAR(0.2, algorithm.centralization(), 1e-12); + const double divisor = normalized ? 3.0 : 1.0; + EXPECT_EQ((std::vector{2 / divisor, 2 / divisor, 2 / divisor, 1 / divisor}), + algorithm.scores()); + } else { + EXPECT_EQ(std::vector(n, 0.0), algorithm.scores()); + } + const auto partition = algorithm.getPartition(); + const auto cover = algorithm.getCover(); + EXPECT_EQ(n, partition.numberOfElements()); + EXPECT_EQ(n, cover.numberOfElements()); + for (node u = 0; u < n; ++u) { + const index expected = n == 4 ? (u == 3 ? 1 : 2) : 0; + EXPECT_EQ(expected, partition.subsetOf(u)); + std::set memberships; + for (index core = 0; core <= expected; ++core) + memberships.insert(core); + EXPECT_EQ(memberships, cover.subsetsOf(u)); + } + } + } +} + TEST_F(CentralityGTest, benchCoreDecompositionLocal) { METISGraphReader reader; std::vector filenames = {"caidaRouterLevel", "wing", "astro-ph", "PGPgiantcompo"}; diff --git a/networkit/test/test_centrality.py b/networkit/test/test_centrality.py index ecdf9bce1..753c3b29b 100755 --- a/networkit/test/test_centrality.py +++ b/networkit/test/test_centrality.py @@ -105,6 +105,128 @@ def testCoreDecomposition(self): CLL = nk.centrality.CoreDecomposition(tmp) CLL.run() self.assertTrue(self.checkCovers(CL.getCover(),CLL.getCover())) + + def assertParallelCoreMatchesReference(self, graph, use_sampling=True, normalized=False): + reference = nk.centrality.CoreDecomposition( + graph, normalized=normalized, enforceBucketQueueAlgorithm=True) + reference.run() + parallel = nk.centrality.ParallelCoreDecomposition( + graph, normalized=normalized, useSampling=use_sampling) + parallel.run() + self.assertListEqual(reference.scores(), parallel.scores()) + self.assertEqual(reference.maxCoreNumber(), parallel.maxCoreNumber()) + return parallel + + def testParallelCoreDecompositionRandomGraphs(self): + for seed, probability in [(1, 0.01), (2, 0.05), (3, 0.2), (4, 0.6)]: + nk.setSeed(seed, False) + graph = nk.generators.ErdosRenyiGenerator(300, probability, False).generate() + self.assertParallelCoreMatchesReference(graph) + + def testParallelCoreDecompositionNormalizedAndEmpty(self): + for size in [0, 1, 4]: + graph = nk.Graph(size) + if size == 4: + for source, target in [(0, 1), (1, 2), (2, 0), (2, 3)]: + graph.addEdge(source, target) + for normalized in [False, True]: + with self.subTest(size=size, normalized=normalized): + # The legacy normalized oracle divides by maximum degree, which is + # undefined for edgeless inputs. Use explicit expected scores there. + algorithm = nk.centrality.ParallelCoreDecomposition( + graph, normalized=normalized).run() + self.assertEqual(algorithm.maximum(), + (1.0 if normalized else 3.0) if size == 4 else 0.0) + if size == 4: + self.assertAlmostEqual(algorithm.centralization(), 0.2) + self.assertParallelCoreMatchesReference(graph, normalized=normalized) + divisor = 3.0 if normalized else 1.0 + self.assertEqual(algorithm.scores(), + [2 / divisor, 2 / divisor, 2 / divisor, 1 / divisor]) + else: + self.assertEqual(algorithm.scores(), [0.0] * size) + expected_cores = [2, 2, 2, 1] if size == 4 else [0] * size + for run in range(2): + partition = algorithm.getPartition() + cover = algorithm.getCover() + self.assertEqual([partition.subsetOf(u) for u in range(size)], expected_cores) + for u, core in enumerate(expected_cores): + self.assertEqual(sorted(cover.subsetsOf(u)), list(range(core + 1))) + if run == 0: + algorithm.run() + + def testParallelCoreDecompositionHierarchicalBuckets(self): + graph = nk.generators.ErdosRenyiGenerator(520, 1.0, False).generate() + parallel = self.assertParallelCoreMatchesReference(graph) + self.assertEqual(parallel.maxCoreNumber(), 519) + + def testParallelCoreDecompositionSampling(self): + graph = nk.Graph(25001) + for leaf in range(1, graph.numberOfNodes()): + graph.addEdge(0, leaf) + parallel = self.assertParallelCoreMatchesReference(graph) + self.assertLessEqual(parallel.numberOfRestarts(), 1) + + def testParallelCoreDecompositionSamplingDenseBipartite(self): + high_degree_nodes = 20 + low_degree_nodes = 20000 + graph = nk.Graph(high_degree_nodes + low_degree_nodes) + for high_degree_node in range(high_degree_nodes): + for low_degree_node in range(high_degree_nodes, graph.numberOfNodes()): + graph.addEdge(high_degree_node, low_degree_node) + parallel = self.assertParallelCoreMatchesReference(graph) + self.assertEqual(parallel.maxCoreNumber(), high_degree_nodes) + self.assertLessEqual(parallel.numberOfRestarts(), 1) + + def testParallelCoreDecompositionWithoutSampling(self): + nk.setSeed(5, False) + graph = nk.generators.ErdosRenyiGenerator(500, 0.15, False).generate() + parallel = self.assertParallelCoreMatchesReference(graph, use_sampling=False) + self.assertEqual(parallel.numberOfRestarts(), 0) + + def testParallelCoreDecompositionGraphR(self): + import pyarrow as pa + + indices = pa.array([1, 2, 0, 2, 0, 1, 3, 2], type=pa.uint64()) + indptr = pa.array([0, 2, 4, 7, 8], type=pa.uint64()) + graph = nk.Graph.fromCSR(4, False, indices, indptr) + parallel = self.assertParallelCoreMatchesReference(graph) + self.assertListEqual(parallel.scores(), [2.0, 2.0, 2.0, 1.0]) + normalized = nk.centrality.ParallelCoreDecomposition(graph, normalized=True).run() + self.assertEqual([normalized.getPartition().subsetOf(u) for u in range(4)], [2, 2, 2, 1]) + for u in range(4): + self.assertEqual(normalized.getCover().subsetsOf(u), parallel.getCover().subsetsOf(u)) + + def testParallelCoreDecompositionGraphRSampling(self): + import pyarrow as pa + + number_of_leaves = 25000 + indices = pa.array( + list(range(1, number_of_leaves + 1)) + [0] * number_of_leaves, + type=pa.uint64()) + indptr = pa.array( + [0, number_of_leaves] + + list(range(number_of_leaves + 1, 2 * number_of_leaves + 1)), + type=pa.uint64()) + graph = nk.Graph.fromCSR(number_of_leaves + 1, False, indices, indptr) + parallel = self.assertParallelCoreMatchesReference(graph) + self.assertEqual(parallel.maxCoreNumber(), 1) + self.assertLessEqual(parallel.numberOfRestarts(), 1) + + def testParallelCoreDecompositionRejectsUnsupportedGraphs(self): + directed = nk.Graph(3, False, True) + with self.assertRaises(RuntimeError): + nk.centrality.ParallelCoreDecomposition(directed) + + self_loop = nk.Graph(3) + self_loop.addEdge(0, 0) + with self.assertRaises(RuntimeError): + nk.centrality.ParallelCoreDecomposition(self_loop) + + non_continuous = nk.Graph(3) + non_continuous.removeNode(1) + with self.assertRaises(RuntimeError): + nk.centrality.ParallelCoreDecomposition(non_continuous) def testComplexPathsAllNodes(self): CP = nk.centrality.ComplexPaths(self.L, 3)