diff --git a/setup.py b/setup.py index 3ba5eaa3..2a9626ab 100644 --- a/setup.py +++ b/setup.py @@ -358,7 +358,8 @@ def format_documentation_entry( docs[f"{name}::{member_name}"] = format_documentation_entry( f"{name}::{member_name}", member_brief, member_detailed ) - with open("./src/dsf/.docstrings.hpp", "w") as f: + + with open("./src/dsf/.docstrings.hpp", "w", encoding="UTF-8") as f: f.write("#pragma once\n\n#include \n#include \n\n") f.write("namespace dsf {\n") f.write( diff --git a/src/dsf/base/Network.hpp b/src/dsf/base/Network.hpp index 0dd3a02b..aad12b45 100644 --- a/src/dsf/base/Network.hpp +++ b/src/dsf/base/Network.hpp @@ -6,17 +6,27 @@ #include "../utility/Typedef.hpp" #include +#include #include #include #include +#include +#include #include +#include +#include #include -#include +#include +#include #include #include #include #include +#include +#include +#include + namespace dsf { template requires(std::is_base_of_v && std::is_base_of_v) @@ -136,6 +146,24 @@ namespace dsf { /// @brief Compute edge weighted betweenness centralities using Brandes' algorithm void computeEdgeBetweennessCentralities(); + + /// @brief Compute edge betweenness centralities using Yen's K-shortest paths. + /// + /// @details For every ordered source–target pair (s, t) with s ≠ t the method + /// finds up to K loopless shortest paths using Yen's algorithm. + /// Each edge e that appears in at least one of those paths receives + /// an increment equal to the number of those K paths that traverse it. + /// After all pairs are processed the raw counts are normalised by + /// (N−1)·(N−2), the number of ordered pairs in a directed graph of N + /// nodes (same denominator used by the standard Brandes formulation). + /// + /// The computation is parallelised over source nodes with Intel TBB: + /// each thread maintains its own accumulator map, which are merged + /// sequentially once all threads finish. + /// + /// @param K Maximum number of shortest paths per (s, t) pair. + /// K = 1 reproduces single-shortest-path betweenness. + void computeEdgeKBetweennessCentralities(std::size_t const K); }; template @@ -219,8 +247,6 @@ namespace dsf { return *it->second; } - constexpr double kPathBudgetEpsilon = 1e-9; - template requires(std::is_base_of_v && std::is_base_of_v) inline std::unordered_map @@ -328,7 +354,7 @@ namespace dsf { double nodeBudget = nodeDistToTarget; if (m_weightThreshold.has_value()) { - nodeBudget *= (1. + *m_weightThreshold); + nodeBudget *= (1.0 + *m_weightThreshold); } std::vector hops; hops.reserve(pNode->outgoingEdges().size()); @@ -369,20 +395,18 @@ namespace dsf { requires(std::is_base_of_v && std::is_base_of_v) inline PathCollection Network::shortestPath(Id const sourceId, Id const targetId) const { - // If source equals target, return empty map (no intermediate hops needed) if (sourceId == targetId) { return PathCollection{}; } - // Check if source node exists if (!this->nodes().contains(sourceId)) { throw std::out_of_range( std::format("Source node with id {} does not exist in the graph", sourceId)); } - // Check if target node exists if (!this->nodes().contains(targetId)) { throw std::out_of_range( std::format("Target node with id {} does not exist in the graph", targetId)); } + auto const distToTarget = m_computeDistancesToTarget(targetId); auto const sourceBestDistance = distToTarget.at(sourceId); @@ -392,7 +416,7 @@ namespace dsf { double sourceBudget = sourceBestDistance; if (m_weightThreshold.has_value()) { - sourceBudget *= (1. + *m_weightThreshold); + sourceBudget *= (1.0 + *m_weightThreshold); } auto const distFromSource = m_computeDistancesFromSource(sourceId); @@ -429,8 +453,7 @@ namespace dsf { auto const projectedDistFromSource = nodeDistFromSource + edgeWeight; // Keep intermediate transitions source-distance-consistent so all - // prefixes to a node share the same cost label. For target hops this - // constraint is unnecessary because no further expansion occurs. + // prefixes to a node share the same cost label. if (nextNodeId != targetId && (projectedDistFromSource > nextDistFromSource + 1e-12 || projectedDistFromSource + 1e-12 < nextDistFromSource)) { @@ -535,25 +558,39 @@ namespace dsf { pNode->setAttribute("betweennessCentrality", 0.0); } - struct PathDataHelper { - std::vector P; - double sigma{0.0}; - double dist{std::numeric_limits::infinity()}; - double delta{0.0}; - }; + auto const N_NODES = m_nodes.size(); - for (auto const& [sourceId, sourceNode] : m_nodes) { - std::unordered_map pathData; - pathData.reserve(this->nNodes()); + // Build an indexed list of source node ids so TBB can address them by integer. + std::vector sourceIds; + sourceIds.reserve(N_NODES); + for (auto const& [id, _] : m_nodes) { + sourceIds.push_back(id); + } + + // Each TBB thread accumulates its own BC increments; we merge afterwards. + tbb::combinable> localBC; + + tbb::parallel_for(std::size_t(0), sourceIds.size(), [&](std::size_t idx) { + auto& local = localBC.local(); + auto const sourceId = sourceIds[idx]; + + struct PathDataHelper { + std::vector P; // predecessor node ids on shortest paths + double sigma{0.0}; + double dist{std::numeric_limits::infinity()}; + double delta{0.0}; + }; + std::unordered_map pathData; + pathData.reserve(N_NODES); for (auto const& [nId, _] : m_nodes) { - pathData.emplace(nId, PathDataHelper()); + pathData.emplace(nId, PathDataHelper{}); } - // Initialize source node data + // Initialise source node data. { - auto& sourceData = pathData[sourceId]; - sourceData.sigma = 1.0; - sourceData.dist = 0.0; + auto& src = pathData.at(sourceId); + src.sigma = 1.0; + src.dist = 0.0; } std::stack S; @@ -574,10 +611,13 @@ namespace dsf { } visited.insert(v); S.push(v); - auto const& vData = pathData[v]; + auto const& vData = pathData.at(v); for (auto const& edgeId : m_nodes.at(v)->outgoingEdges()) { auto const& edgeObj = *m_edges.at(edgeId); + if (!edgeObj.isActive()) { + continue; + } auto const w = edgeObj.target(); auto& wData = pathData.at(w); if (visited.contains(w)) { @@ -598,51 +638,74 @@ namespace dsf { } } + // Back-propagate dependency scores. while (!S.empty()) { auto const w = S.top(); - auto const& wData = pathData[w]; + auto& wData = pathData.at(w); S.pop(); for (auto const v : wData.P) { auto& vData = pathData.at(v); vData.delta += (vData.sigma / wData.sigma) * (1.0 + wData.delta); } if (w != sourceId) { - auto& currentNode = this->node(w); - auto const optBC = - currentNode.template getAttribute("betweennessCentrality"); - // BC must exist since we initialized it for all nodes at the start - currentNode.setAttribute("betweennessCentrality", *optBC + wData.delta); + local[w] += wData.delta; } } - } + }); + + // Merge thread-local accumulations into node attributes (sequential). + localBC.combine_each([&](auto const& localMap) { + for (auto const& [nodeId, value] : localMap) { + auto& nd = this->node(nodeId); + auto const bc = nd.template getAttribute("betweennessCentrality"); + nd.setAttribute("betweennessCentrality", *bc + value); + } + }); } + // ───────────────────────────────────────────────────────────────────────── + // Brandes edge betweenness (parallelised over sources with TBB) + // ───────────────────────────────────────────────────────────────────────── template requires(std::is_base_of_v && std::is_base_of_v) void Network::computeEdgeBetweennessCentralities() { + // Initialise all edge BC values to 0. for (auto& [edgeId, pEdge] : m_edges) { pEdge->setAttribute("betweennessCentrality", 0.0); } - struct PathDataHelper { - std::vector P; // predecessor edge ids on shortest paths - double sigma{0.0}; - double dist{std::numeric_limits::infinity()}; - double delta{0.0}; - }; + auto const N_NODES = m_nodes.size(); - for (auto const& [sourceId, sourceNode] : m_nodes) { - std::unordered_map pathData; - pathData.reserve(this->nNodes()); + std::vector sourceIds; + sourceIds.reserve(N_NODES); + for (auto const& [id, _] : m_nodes) { + sourceIds.push_back(id); + } + // Each TBB thread accumulates its own edge BC increments. + tbb::combinable> localBC; + + tbb::parallel_for(size_t(0), sourceIds.size(), [&](size_t idx) { + auto& local = localBC.local(); + auto const sourceId = sourceIds[idx]; + + struct PathDataHelper { + std::vector P; // predecessor edge ids on shortest paths + double sigma{0.0}; + double dist{std::numeric_limits::infinity()}; + double delta{0.0}; + }; + + std::unordered_map pathData; + pathData.reserve(N_NODES); for (auto const& [nId, _] : m_nodes) { - pathData.emplace(nId, PathDataHelper()); + pathData.emplace(nId, PathDataHelper{}); } - // Initialize source node data + // Initialise source node data. { - auto& sourceData = pathData.at(sourceId); - sourceData.sigma = 1.0; - sourceData.dist = 0.0; + auto& src = pathData.at(sourceId); + src.sigma = 1.0; + src.dist = 0.0; } std::stack S; @@ -658,18 +721,23 @@ namespace dsf { auto [d, v] = pq.top(); pq.pop(); - if (visited.contains(v)) + if (visited.contains(v)) { continue; + } visited.insert(v); S.push(v); auto const& vData = pathData.at(v); - for (auto const& eId : this->node(v).outgoingEdges()) { - auto const& edgeObj = this->edge(eId); + for (auto const& eId : m_nodes.at(v)->outgoingEdges()) { + auto const& edgeObj = *m_edges.at(eId); + if (!edgeObj.isActive()) { + continue; + } auto const w = edgeObj.target(); auto& wData = pathData.at(w); - if (visited.contains(w)) + if (visited.contains(w)) { continue; + } double const edgeWeight = m_weightFunction(edgeObj); double const newDist = vData.dist + edgeWeight; @@ -685,21 +753,344 @@ namespace dsf { } } + // Back-propagate dependency scores; accumulate per-edge contributions. while (!S.empty()) { auto const w = S.top(); auto& wData = pathData.at(w); S.pop(); for (auto const eId : wData.P) { - auto& vData = pathData.at(this->edge(eId).source()); + auto& vData = pathData.at(m_edges.at(eId)->source()); double const contrib = (vData.sigma / wData.sigma) * (1.0 + wData.delta); vData.delta += contrib; - auto& currentEdge = this->edge(eId); - auto const optBC = - currentEdge.template getAttribute("betweennessCentrality"); - // BC must exist since we initialized it for all edges at the start - currentEdge.setAttribute("betweennessCentrality", *optBC + contrib); + local[eId] += contrib; + } + } + }); + + // Merge thread-local accumulations into edge attributes (sequential). + localBC.combine_each([&](auto const& localMap) { + for (auto const& [edgeId, value] : localMap) { + auto& ed = this->edge(edgeId); + auto const bc = + ed.template getAttribute("betweennessCentrality").value_or(0.0); + ed.setAttribute("betweennessCentrality", bc + value); + } + }); + } + + // ───────────────────────────────────────────────────────────────────────── + // Yen K-shortest-paths edge betweenness (parallelised over sources, TBB) + // ───────────────────────────────────────────────────────────────────────── + template + requires(std::is_base_of_v && std::is_base_of_v) + void Network::computeEdgeKBetweennessCentralities(std::size_t const K) { + for (auto& [eId, pEdge] : m_edges) { + pEdge->setAttribute("betweennessCentrality", 0.0); + } + + auto const N_NODES = m_nodes.size(); + if (N_NODES < 2 || K == 0) { + return; + } + + // K == 1 reduces exactly to standard Brandes edge betweenness. + // K=1 delegates to the Brandes algorithm (computeEdgeBetweennessCentralities), + // which accounts for tied shortest paths via σ fractions. Both K=1 and K>1 + // paths use the same (N-1)*(N-2) normalization, making K=1 equivalent to + // single-source Brandes betweenness. Results are normalized by (N-1)*(N-2). + // K>1 uses Yen's algorithm to find K shortest loopless paths and counts + // edge contributions across all K paths. + if (K == 1) { + computeEdgeBetweennessCentralities(); + // Re-normalise: computeEdgeBetweennessCentralities leaves raw Brandes + // values; apply the same norm as the K>1 path below. + auto const norm = static_cast((N_NODES - 1) * (N_NODES - 2)); + if (norm > 0.0) { + for (auto& [_, pEdge] : m_edges) { + auto const bc = + pEdge->template getAttribute("betweennessCentrality").value_or(0.0); + pEdge->setAttribute("betweennessCentrality", bc / norm); + } + } + return; + } + + // K > 1: Yen's algorithm, parallelised over source nodes. + std::vector nodeIds; + nodeIds.reserve(N_NODES); + for (auto const& [id, _] : m_nodes) { + nodeIds.push_back(id); + } + + std::unordered_map nodeIndex; + nodeIndex.reserve(N_NODES); + for (std::size_t i = 0; i < N_NODES; ++i) { + nodeIndex.emplace(nodeIds[i], i); + } + + tbb::combinable> localAccum; + + tbb::parallel_for(std::size_t(0), N_NODES, [&](std::size_t const sourceIndex) { + auto& local = localAccum.local(); + + struct PathData { + std::vector edges; + double cost{0.0}; + }; + + for (std::size_t targetIndex = 0; targetIndex < N_NODES; ++targetIndex) { + if (sourceIndex == targetIndex) { + continue; + } + + std::vector acceptedPaths; + + // ── Global candidate set for Yen's algorithm across all ranks ───── + // Standard Yen's algorithm maintains a single priority queue of candidate + // paths across all K rank iterations. Candidates generated but not selected + // in earlier ranks remain available for later ranks (B-heap approach). + std::priority_queue>, + std::vector>>, + std::greater<>> + globalCandidates; + + // ── Initial shortest path (Dijkstra) ────────────────────────────── + std::vector dist(N_NODES, std::numeric_limits::infinity()); + std::vector parentEdge(N_NODES, Id{}); + std::priority_queue, + std::vector>, + std::greater<>> + pq; + dist[sourceIndex] = 0.0; + pq.push({0.0, sourceIndex}); + + while (!pq.empty()) { + auto const [currentDist, currentIndex] = pq.top(); + pq.pop(); + if (currentDist > dist[currentIndex]) { + continue; + } + if (currentIndex == targetIndex) { + break; + } + + auto const currentNodeId = nodeIds[currentIndex]; + for (auto const edgeId : m_nodes.at(currentNodeId)->outgoingEdges()) { + auto const& edgeObj = this->edge(edgeId); + if (!edgeObj.isActive()) { + continue; + } + auto const nextIdx = nodeIndex.at(edgeObj.target()); + auto const nd = currentDist + m_weightFunction(edgeObj); + if (nd < dist[nextIdx]) { + dist[nextIdx] = nd; + parentEdge[nextIdx] = edgeId; + pq.push({nd, nextIdx}); + } + } + } + + if (!std::isfinite(dist[targetIndex])) { + continue; + } + + // Reconstruct first shortest path. + PathData firstPath; + firstPath.cost = dist[targetIndex]; + for (size_t cursor = targetIndex; cursor != sourceIndex;) { + Id const eId = parentEdge[cursor]; + firstPath.edges.push_back(eId); + cursor = nodeIndex.at(m_edges.at(eId)->source()); } + std::reverse(firstPath.edges.begin(), firstPath.edges.end()); + acceptedPaths.push_back(std::move(firstPath)); + + // ── Yen's spur loop ─────────────────────────────────────────────── + for (size_t pathRank = 1; pathRank < K; ++pathRank) { + auto const& lastPath = acceptedPaths.back().edges; + + for (size_t prefixSize = 0; prefixSize < lastPath.size(); ++prefixSize) { + std::vector rootPath( + lastPath.begin(), + lastPath.begin() + static_cast(prefixSize)); + + // Determine spur node. + Id spurNodeId = nodeIds[sourceIndex]; + for (auto const eId : rootPath) { + spurNodeId = m_edges.at(eId)->target(); + } + size_t const spurIndex = nodeIndex.at(spurNodeId); + + // Build banned-node set: all root-path nodes except the spur node. + std::unordered_set bannedNodes; + Id visitedNodeId = nodeIds[sourceIndex]; + for (auto const eId : rootPath) { + Id const srcNodeId = visitedNodeId; + visitedNodeId = m_edges.at(eId)->target(); + if (srcNodeId != spurNodeId) { + bannedNodes.insert(srcNodeId); + } + } + + // Ban the outgoing edge at the spur position for every accepted + // path that shares the same root prefix. + std::unordered_set bannedEdges; + for (auto const& path : acceptedPaths) { + if (path.edges.size() <= prefixSize) { + continue; + } + bool samePrefix = true; + for (size_t i = 0; i < prefixSize; ++i) { + if (path.edges[i] != rootPath[i]) { + samePrefix = false; + break; + } + } + if (samePrefix) { + bannedEdges.insert(path.edges[prefixSize]); + } + } + + if (bannedNodes.contains(spurNodeId)) { + continue; + } + + // Dijkstra from the spur node with banned nodes/edges. + std::vector spurDist(N_NODES, + std::numeric_limits::infinity()); + std::vector spurParentEdge(N_NODES, Id{}); + std::priority_queue, + std::vector>, + std::greater<>> + spurQueue; + spurDist[spurIndex] = 0.0; + spurQueue.push({0.0, spurIndex}); + + while (!spurQueue.empty()) { + auto const [currentDist, currentIndex] = spurQueue.top(); + spurQueue.pop(); + if (currentDist > spurDist[currentIndex]) { + continue; + } + if (currentIndex == targetIndex) { + break; + } + + auto const curNodeId = nodeIds[currentIndex]; + if (bannedNodes.contains(curNodeId)) { + continue; + } + + for (auto const eId : m_nodes.at(curNodeId)->outgoingEdges()) { + if (bannedEdges.contains(eId)) { + continue; + } + auto const& edgeObj = *m_edges.at(eId); + if (!edgeObj.isActive()) { + continue; + } + Id const nextNode = edgeObj.target(); + if (bannedNodes.contains(nextNode)) { + continue; + } + auto const nextIdx = nodeIndex.at(nextNode); + auto const nd = currentDist + m_weightFunction(edgeObj); + if (nd < spurDist[nextIdx]) { + spurDist[nextIdx] = nd; + spurParentEdge[nextIdx] = eId; + spurQueue.push({nd, nextIdx}); + } + } + } + + if (!std::isfinite(spurDist[targetIndex])) { + continue; + } + + // Build candidate = rootPath + spurPath. + PathData candidatePath; + candidatePath.edges = rootPath; + candidatePath.cost = 0.0; + for (auto const eId : rootPath) { + candidatePath.cost += m_weightFunction(*m_edges.at(eId)); + } + + std::vector spurPath; + for (std::size_t cursor = targetIndex; cursor != spurIndex;) { + Id const eId = spurParentEdge[cursor]; + spurPath.push_back(eId); + cursor = nodeIndex.at(m_edges.at(eId)->source()); + } + std::reverse(spurPath.begin(), spurPath.end()); + + candidatePath.cost += spurDist[targetIndex]; + candidatePath.edges.insert( + candidatePath.edges.end(), spurPath.begin(), spurPath.end()); + + bool duplicate = false; + for (auto const& path : acceptedPaths) { + if (path.edges == candidatePath.edges) { + duplicate = true; + break; + } + } + if (!duplicate) { + globalCandidates.push({candidatePath.cost, std::move(candidatePath.edges)}); + } + } // end prefix loop + + bool foundPath = false; + while (!globalCandidates.empty() && !foundPath) { + auto [candidateCost, candidateEdges] = globalCandidates.top(); + globalCandidates.pop(); + + bool alreadyAccepted = false; + for (auto const& path : acceptedPaths) { + if (path.edges == candidateEdges) { + alreadyAccepted = true; + break; + } + } + + if (!alreadyAccepted) { + acceptedPaths.push_back(PathData{std::move(candidateEdges), candidateCost}); + foundPath = true; + } + } + + if (!foundPath && globalCandidates.empty()) { + break; // No more candidates; K-path generation complete + } + } // end Yen rank loop + + // Accumulate edge counts. + for (auto const& path : acceptedPaths) { + for (auto const eId : path.edges) { + local[eId] += 1.0; + } + } + } // end target loop + }); // end parallel_for + + // Merge thread-local accumulations into edge attributes (sequential). + localAccum.combine_each([&](auto const& localMap) { + for (auto const& [edgeId, value] : localMap) { + auto& ed = this->edge(edgeId); + auto const bc = + ed.template getAttribute("betweennessCentrality").value_or(0.0); + ed.setAttribute("betweennessCentrality", bc + value); + } + }); + + // Normalise by (N-1)·(N-2). + auto const norm = static_cast((N_NODES - 1) * (N_NODES - 2)); + if (norm > 0.0) { + for (auto& [_, pEdge] : m_edges) { + auto const bc = + pEdge->template getAttribute("betweennessCentrality").value_or(0.0); + pEdge->setAttribute("betweennessCentrality", bc / norm); } } } + } // namespace dsf \ No newline at end of file diff --git a/src/dsf/bindings.cpp b/src/dsf/bindings.cpp index 3e4a558a..27d65535 100644 --- a/src/dsf/bindings.cpp +++ b/src/dsf/bindings.cpp @@ -329,6 +329,11 @@ PYBIND11_MODULE(dsf_cpp, m) { "computeEdgeBetweennessCentralities", &dsf::mobility::RoadNetwork::computeEdgeBetweennessCentralities, dsf::g_docstrings.at("dsf::Network::computeEdgeBetweennessCentralities").c_str()) + .def("computeEdgeKBetweennessCentralities", + &dsf::mobility::RoadNetwork::computeEdgeKBetweennessCentralities, + pybind11::arg("k"), + dsf::g_docstrings.at("dsf::Network::computeEdgeKBetweennessCentralities") + .c_str()) .def( "nodeBetweennessCentralities", [](const dsf::mobility::RoadNetwork& self) { diff --git a/src/dsf/mobility/RoadNetwork.cpp b/src/dsf/mobility/RoadNetwork.cpp index 39dca561..d75ff4bd 100644 --- a/src/dsf/mobility/RoadNetwork.cpp +++ b/src/dsf/mobility/RoadNetwork.cpp @@ -3,7 +3,9 @@ #include "../geometry/PolyLine.hpp" #include +#include #include +#include #include #include @@ -22,6 +24,8 @@ static constexpr auto EDGE_DEFAULT_ATTRIBUTES = "nlanes", "name", "type", + "capacity", + "status", "coilcode", "priority", "geometry"}); @@ -63,6 +67,10 @@ namespace dsf::mobility { (std::find(colNames.begin(), colNames.end(), "coilcode") != colNames.end()); bool const bHasPriority = (std::find(colNames.begin(), colNames.end(), "priority") != colNames.end()); + bool const bHasCapacity = + (std::find(colNames.begin(), colNames.end(), "capacity") != colNames.end()); + bool const bHasStatus = + (std::find(colNames.begin(), colNames.end(), "status") != colNames.end()); for (auto& row : reader) { auto const sourceId = row["source"].get(); @@ -150,6 +158,43 @@ namespace dsf::mobility { addCoil(streetId, strCoilCode); } } + + // Parse capacity field if present + if (bHasCapacity) { + try { + int capacityValue = row["capacity"].get(); + edge(streetId).setCapacity(capacityValue); + } catch (...) { + spdlog::warn("Invalid capacity for edge {}. Using default.", streetId); + } + } + + // Parse status field if present + if (bHasStatus) { + try { + auto statusStr = row["status"].get(); + std::transform( + statusStr.begin(), statusStr.end(), statusStr.begin(), [](unsigned char c) { + return std::tolower(c); + }); + if (statusStr == "closed") { + edge(streetId).setStatus(RoadStatus::CLOSED); + } else if (statusStr == "open" || statusStr.empty()) { + edge(streetId).setStatus(RoadStatus::OPEN); + } else { + spdlog::warn( + "Unknown status '{}' for edge {}. Valid values: 'open', 'closed'. Using " + "'open'.", + statusStr, + streetId); + edge(streetId).setStatus(RoadStatus::OPEN); + } + } catch (...) { + spdlog::warn("Invalid status for edge {}. Using default (OPEN).", streetId); + edge(streetId).setStatus(RoadStatus::OPEN); + } + } + // Handle additional attributes for (auto const attrName : additionalAttributes) { auto const strAttrName{std::string(attrName)}; diff --git a/test/mobility/Test_graph.cpp b/test/mobility/Test_graph.cpp index ff513ee6..f246b42d 100644 --- a/test/mobility/Test_graph.cpp +++ b/test/mobility/Test_graph.cpp @@ -2218,3 +2218,225 @@ TEST_CASE("EdgeBetweennessCentrality") { doctest::Approx(1.0)); } } + +TEST_CASE("computeEdgeKBetweennessCentralities") { + Road::setMeanVehicleLength(5.); + + SUBCASE("K=1 matches normalized standard edge BC on linear chain") { + RoadNetwork standardGraph{}; + Street stdS01(0, std::make_pair(0, 1), 10.0); + Street stdS12(1, std::make_pair(1, 2), 10.0); + Street stdS23(2, std::make_pair(2, 3), 10.0); + standardGraph.addStreets(stdS01, stdS12, stdS23); + standardGraph.setEdgeWeight("uniform"); + standardGraph.computeEdgeBetweennessCentralities(); + + size_t const nNodes = standardGraph.nNodes(); + double const norm = static_cast((nNodes - 1) * (nNodes - 2)); + REQUIRE(norm > 0.0); + + auto const stdBc0 = standardGraph.edge(static_cast(0)) + .getAttribute("betweennessCentrality"); + auto const stdBc1 = standardGraph.edge(static_cast(1)) + .getAttribute("betweennessCentrality"); + auto const stdBc2 = standardGraph.edge(static_cast(2)) + .getAttribute("betweennessCentrality"); + REQUIRE(stdBc0.has_value()); + REQUIRE(stdBc1.has_value()); + REQUIRE(stdBc2.has_value()); + + double const expected0 = *stdBc0 / norm; + double const expected1 = *stdBc1 / norm; + double const expected2 = *stdBc2 / norm; + + RoadNetwork yenGraph{}; + Street yenS01(0, std::make_pair(0, 1), 10.0); + Street yenS12(1, std::make_pair(1, 2), 10.0); + Street yenS23(2, std::make_pair(2, 3), 10.0); + yenGraph.addStreets(yenS01, yenS12, yenS23); + yenGraph.setEdgeWeight("uniform"); + yenGraph.computeEdgeKBetweennessCentralities(1); + + auto const yenBc0 = + yenGraph.edge(static_cast(0)).getAttribute("betweennessCentrality"); + auto const yenBc1 = + yenGraph.edge(static_cast(1)).getAttribute("betweennessCentrality"); + auto const yenBc2 = + yenGraph.edge(static_cast(2)).getAttribute("betweennessCentrality"); + REQUIRE(yenBc0.has_value()); + REQUIRE(yenBc1.has_value()); + REQUIRE(yenBc2.has_value()); + + CHECK_EQ(*yenBc0, doctest::Approx(expected0)); + CHECK_EQ(*yenBc1, doctest::Approx(expected1)); + CHECK_EQ(*yenBc2, doctest::Approx(expected2)); + } + + SUBCASE("K=2 captures additional shortest-path alternatives") { + RoadNetwork graphK1{}; + Street k1S01(0, std::make_pair(0, 1), 10.0); + Street k1S02(1, std::make_pair(0, 2), 10.0); + Street k1S13(2, std::make_pair(1, 3), 10.0); + Street k1S23(3, std::make_pair(2, 3), 10.0); + graphK1.addStreets(k1S01, k1S02, k1S13, k1S23); + graphK1.setEdgeWeight("uniform"); + graphK1.computeEdgeKBetweennessCentralities(1); + + RoadNetwork graphK2{}; + Street k2S01(0, std::make_pair(0, 1), 10.0); + Street k2S02(1, std::make_pair(0, 2), 10.0); + Street k2S13(2, std::make_pair(1, 3), 10.0); + Street k2S23(3, std::make_pair(2, 3), 10.0); + graphK2.addStreets(k2S01, k2S02, k2S13, k2S23); + graphK2.setEdgeWeight("uniform"); + graphK2.computeEdgeKBetweennessCentralities(2); + + bool hasStrictIncrease = false; + Id const edgeIds[] = { + static_cast(0), static_cast(1), static_cast(2), static_cast(3)}; + + for (Id const edgeId : edgeIds) { + auto const k1Bc = + graphK1.edge(edgeId).getAttribute("betweennessCentrality"); + auto const k2Bc = + graphK2.edge(edgeId).getAttribute("betweennessCentrality"); + REQUIRE(k1Bc.has_value()); + REQUIRE(k2Bc.has_value()); + + CHECK(*k2Bc + 1e-12 >= *k1Bc); + if (*k2Bc > *k1Bc + 1e-12) { + hasStrictIncrease = true; + } + } + + CHECK(hasStrictIncrease); + } + + SUBCASE("K=1 matches Brandes on a diamond graph with tied shortest paths") { + RoadNetwork standardGraph{}; + Street stdS01(0, std::make_pair(0, 1), 10.0); + Street stdS02(1, std::make_pair(0, 2), 10.0); + Street stdS13(2, std::make_pair(1, 3), 10.0); + Street stdS23(3, std::make_pair(2, 3), 10.0); + standardGraph.addStreets(stdS01, stdS02, stdS13, stdS23); + standardGraph.setEdgeWeight("uniform"); + standardGraph.computeEdgeBetweennessCentralities(); + + std::size_t const nNodes = standardGraph.nNodes(); + double const norm = static_cast((nNodes - 1) * (nNodes - 2)); + REQUIRE(norm > 0.0); + + RoadNetwork yenGraph{}; + Street yenS01(0, std::make_pair(0, 1), 10.0); + Street yenS02(1, std::make_pair(0, 2), 10.0); + Street yenS13(2, std::make_pair(1, 3), 10.0); + Street yenS23(3, std::make_pair(2, 3), 10.0); + yenGraph.addStreets(yenS01, yenS02, yenS13, yenS23); + yenGraph.setEdgeWeight("uniform"); + yenGraph.computeEdgeKBetweennessCentralities(1); + + Id const edgeIds[] = { + static_cast(0), static_cast(1), static_cast(2), static_cast(3)}; + + for (Id const edgeId : edgeIds) { + auto const standardBc = + standardGraph.edge(edgeId).getAttribute("betweennessCentrality"); + auto const yenBc = + yenGraph.edge(edgeId).getAttribute("betweennessCentrality"); + REQUIRE(standardBc.has_value()); + REQUIRE(yenBc.has_value()); + CHECK(*yenBc == doctest::Approx(*standardBc / norm)); + } + } + + SUBCASE("K=2 counts both equal shortest paths exactly on a diamond graph") { + RoadNetwork graph{}; + Street s01(0, std::make_pair(0, 1), 10.0); + Street s02(1, std::make_pair(0, 2), 10.0); + Street s13(2, std::make_pair(1, 3), 10.0); + Street s23(3, std::make_pair(2, 3), 10.0); + graph.addStreets(s01, s02, s13, s23); + graph.setEdgeWeight("uniform"); + graph.computeEdgeKBetweennessCentralities(2); + + Id const edgeIds[] = { + static_cast(0), static_cast(1), static_cast(2), static_cast(3)}; + + for (Id const edgeId : edgeIds) { + auto const bc = graph.edge(edgeId).getAttribute("betweennessCentrality"); + REQUIRE(bc.has_value()); + CHECK(*bc == doctest::Approx(1.0 / 3.0)); + } + } + + SUBCASE("K=3 finds three equal shortest paths with persistent candidate heap") { + // Graph structure: 3 edge-disjoint paths from source 0 to target 4 + // Each path has 2 edges with uniform cost + // Path 1: 0 -> 1 -> 4 (edges 0, 1) + // Path 2: 0 -> 2 -> 4 (edges 2, 3) + // Path 3: 0 -> 3 -> 4 (edges 4, 5) + // + // For (s, t) = (0, 4), all 3 paths have cost 20, so K=3 should find all of them. + // This tests that the global candidate heap persists across rank iterations + // and correctly accumulates spur candidates from earlier ranks. + RoadNetwork graph{}; + Street s01(0, std::make_pair(0, 1), 10.0); + Street s14(1, std::make_pair(1, 4), 10.0); + Street s02(2, std::make_pair(0, 2), 10.0); + Street s24(3, std::make_pair(2, 4), 10.0); + Street s03(4, std::make_pair(0, 3), 10.0); + Street s34(5, std::make_pair(3, 4), 10.0); + graph.addStreets(s01, s14, s02, s24, s03, s34); + graph.setEdgeWeight("uniform"); + graph.computeEdgeKBetweennessCentralities(3); + + // Verify all 6 edges appear in the computations + Id const edgeIds[] = {static_cast(0), + static_cast(1), + static_cast(2), + static_cast(3), + static_cast(4), + static_cast(5)}; + + std::unordered_map edgeCentralities; + for (Id const edgeId : edgeIds) { + auto const bc = graph.edge(edgeId).getAttribute("betweennessCentrality"); + REQUIRE(bc.has_value()); + edgeCentralities[edgeId] = *bc; + } + + // Edges on the paths (0,1), (1,4), (0,2), (2,4), (0,3), (3,4) contribute equally + // because all three paths are shortest and equally likely. + // Each edge should appear in 1 out of 3 paths from node 0 to node 4. + // Normalization factor: (N-1)*(N-2) = (5-1)*(5-2) = 4*3 = 12 + // Each edge appears 1 time across all 3 paths for (0,4), and it's not on paths + // from other pairs. So contribution = 1/12. + // Actually, let's verify that all edges contribute positively and K=3 >= K=2 + RoadNetwork graphK2{}; + Street k2_s01(0, std::make_pair(0, 1), 10.0); + Street k2_s14(1, std::make_pair(1, 4), 10.0); + Street k2_s02(2, std::make_pair(0, 2), 10.0); + Street k2_s24(3, std::make_pair(2, 4), 10.0); + Street k2_s03(4, std::make_pair(0, 3), 10.0); + Street k2_s34(5, std::make_pair(3, 4), 10.0); + graphK2.addStreets(k2_s01, k2_s14, k2_s02, k2_s24, k2_s03, k2_s34); + graphK2.setEdgeWeight("uniform"); + graphK2.computeEdgeKBetweennessCentralities(2); + + bool hasIncrease = false; + for (Id const edgeId : edgeIds) { + auto const bcK2 = + graphK2.edge(edgeId).getAttribute("betweennessCentrality"); + auto const bcK3 = edgeCentralities[edgeId]; + REQUIRE(bcK2.has_value()); + // K=3 should find at least one additional path for at least one edge pair + CHECK(bcK3 + 1e-12 >= *bcK2); + if (bcK3 > *bcK2 + 1e-12) { + hasIncrease = true; + } + } + // Verify that at least some edges increased from K=2 to K=3 + // (This validates the global candidate heap: without it, no 3rd path would be found) + CHECK(hasIncrease); + } +}