Add bc attribute for nodes and edges

src/dsf/base/Edge.hpp

@@ -15,7 +15,8 @@ namespace dsf {
     geometry::PolyLine m_geometry;
     Id m_id;
     std::pair<Id, Id> m_nodePair;
-    std::optional<double> m_weight;
+    std::optional<double> m_betweennessCentrality{std::nullopt};
+    std::optional<double> m_weight{std::nullopt};
     double m_angle;
 
     void m_setAngle(geometry::Point srcNodeCoordinates,
@@ -35,6 +36,11 @@ namespace dsf {
     /// @brief Set the edge's geometry
     /// @param geometry dsf::geometry::PolyLine The edge's geometry, a vector of pairs of doubles representing the coordinates of the edge's geometry
     void setGeometry(geometry::PolyLine geometry);
+    /// @brief Set the edge's betweenness centrality
+    /// @param betweennessCentrality The edge's betweenness centrality
+    inline void setBetweennessCentrality(double const betweennessCentrality) {
+      m_betweennessCentrality = betweennessCentrality;
+    }
     /// @brief Set the edge's weight
     /// @param weight The edge's weight
     /// @throws std::invalid_argument if the weight is less or equal to 0
@@ -56,6 +62,9 @@ namespace dsf {
     /// @brief Get the edge's geometry
     /// @return dsf::geometry::PolyLine The edge's geometry, a vector of pairs of doubles representing the coordinates of the edge's geometry
     inline auto const& geometry() const { return m_geometry; }
+    /// @brief Get the edge's betweenness centrality
+    /// @return std::optional<double> The edge's betweenness centrality
+    inline auto const& betweennessCentrality() const { return m_betweennessCentrality; }
 
     /// @brief Get the edge's angle, in radians, between the source and target nodes
     /// @return double The edge's angle, in radians

src/dsf/base/Network.hpp

@@ -1,7 +1,14 @@
 #pragma once
 
 #include <cassert>
+#include <cmath>
+#include <limits>
+#include <queue>
+#include <stack>
+#include <type_traits>
 #include <unordered_map>
+#include <unordered_set>
+#include <vector>
 
 #include "Edge.hpp"
 #include "Node.hpp"
@@ -88,6 +95,32 @@ namespace dsf {
     template <typename TEdge>
       requires(std::is_base_of_v<edge_t, TEdge>)
     TEdge& edge(Id edgeId);
+
+    /// @brief Compute betweenness centralities for all nodes using Brandes' algorithm
+    /// @tparam WeightFunc A callable type that takes a const reference to a unique_ptr<edge_t> and returns a double representing the edge weight
+    /// @param getEdgeWeight A callable that takes a const reference to a unique_ptr<edge_t> and returns a double (must be positive)
+    /// @details Implements Brandes' algorithm for directed weighted graphs.
+    ///          The computed centrality for each node v is:
+    ///          C_B(v) = sum_{s != v != t} sigma_st(v) / sigma_st
+    ///          where sigma_st is the number of shortest paths from s to t,
+    ///          and sigma_st(v) is the number of those paths passing through v.
+    ///          Results are stored via Node::setBetweennessCentrality.
+    template <typename WeightFunc>
+      requires(std::is_invocable_r_v<double, WeightFunc, std::unique_ptr<edge_t> const&>)
+    void computeBetweennessCentralities(WeightFunc getEdgeWeight);
+
+    /// @brief Compute edge betweenness centralities for all edges using Brandes' algorithm
+    /// @tparam WeightFunc A callable type that takes a const reference to a unique_ptr<edge_t> and returns a double representing the edge weight
+    /// @param getEdgeWeight A callable that takes a const reference to a unique_ptr<edge_t> and returns a double (must be positive)
+    /// @details Implements Brandes' algorithm for directed weighted graphs.
+    ///          The computed centrality for each edge e is:
+    ///          C_B(e) = sum_{s != t} sigma_st(e) / sigma_st
+    ///          where sigma_st is the number of shortest paths from s to t,
+    ///          and sigma_st(e) is the number of those paths using edge e.
+    ///          Results are stored via Edge::setBetweennessCentrality.
+    template <typename WeightFunc>
+      requires(std::is_invocable_r_v<double, WeightFunc, std::unique_ptr<edge_t> const&>)
+    void computeEdgeBetweennessCentralities(WeightFunc getEdgeWeight);
   };
   template <typename node_t, typename edge_t>
     requires(std::is_base_of_v<Node, node_t> && std::is_base_of_v<Edge, edge_t>)
@@ -242,4 +275,176 @@ namespace dsf {
   TEdge& Network<node_t, edge_t>::edge(Id edgeId) {
     return dynamic_cast<TEdge&>(*edge(edgeId));
   }
+
+  template <typename node_t, typename edge_t>
+    requires(std::is_base_of_v<Node, node_t> && std::is_base_of_v<Edge, edge_t>)
+  template <typename WeightFunc>
+    requires(std::is_invocable_r_v<double, WeightFunc, std::unique_ptr<edge_t> const&>)
+  void Network<node_t, edge_t>::computeBetweennessCentralities(WeightFunc getEdgeWeight) {
+    // Initialize all node betweenness centralities to 0
+    for (auto& [nodeId, pNode] : m_nodes) {
+      pNode->setBetweennessCentrality(0.0);
+    }
+
+    // Brandes' algorithm: run single-source Dijkstra from each node
+    for (auto const& [sourceId, sourceNode] : m_nodes) {
+      std::stack<Id> S;  // nodes in order of non-increasing distance
+      std::unordered_map<Id, std::vector<Id>> P;  // predecessors on shortest paths
+      std::unordered_map<Id, double> sigma;       // number of shortest paths
+      std::unordered_map<Id, double> dist;        // distance from source
+
+      for (auto const& [nId, _] : m_nodes) {
+        P[nId] = {};
+        sigma[nId] = 0.0;
+        dist[nId] = std::numeric_limits<double>::infinity();
+      }
+      sigma[sourceId] = 1.0;
+      dist[sourceId] = 0.0;
+
+      // Min-heap priority queue: (distance, nodeId)
+      std::priority_queue<std::pair<double, Id>,
+                          std::vector<std::pair<double, Id>>,
+                          std::greater<>>
+          pq;
+      pq.push({0.0, sourceId});
+
+      std::unordered_set<Id> visited;
+
+      while (!pq.empty()) {
+        auto [d, v] = pq.top();
+        pq.pop();
+
+        if (visited.contains(v)) {
+          continue;
+        }
+        visited.insert(v);
+        S.push(v);
+
+        for (auto const& edgeId : m_nodes.at(v)->outgoingEdges()) {
+          auto const& pEdge = m_edges.at(edgeId);
+          Id w = pEdge->target();
+          if (visited.contains(w)) {
+            continue;
+          }
+          double edgeWeight = getEdgeWeight(pEdge);
+          double newDist = dist[v] + edgeWeight;
+
+          if (newDist < dist[w]) {
+            dist[w] = newDist;
+            sigma[w] = sigma[v];
+            P[w] = {v};
+            pq.push({newDist, w});
+          } else if (std::abs(newDist - dist[w]) < 1e-12 * std::max(1.0, dist[w])) {
+            sigma[w] += sigma[v];
+            P[w].push_back(v);
+          }
+        }
+      }
+
+      // Dependency accumulation (backward pass)
+      std::unordered_map<Id, double> delta;
+      for (auto const& [nId, _] : m_nodes) {
+        delta[nId] = 0.0;
+      }
+
+      while (!S.empty()) {
+        Id w = S.top();
+        S.pop();
+        for (Id v : P[w]) {
+          delta[v] += (sigma[v] / sigma[w]) * (1.0 + delta[w]);
+        }
+        if (w != sourceId) {
+          auto currentBC = m_nodes.at(w)->betweennessCentrality();
+          m_nodes.at(w)->setBetweennessCentrality(*currentBC + delta[w]);
+        }
+      }
+    }
+  }
+
+  template <typename node_t, typename edge_t>
+    requires(std::is_base_of_v<Node, node_t> && std::is_base_of_v<Edge, edge_t>)
+  template <typename WeightFunc>
+    requires(std::is_invocable_r_v<double, WeightFunc, std::unique_ptr<edge_t> const&>)
+  void Network<node_t, edge_t>::computeEdgeBetweennessCentralities(
+      WeightFunc getEdgeWeight) {
+    // Initialize all edge betweenness centralities to 0
+    for (auto& [edgeId, pEdge] : m_edges) {
+      pEdge->setBetweennessCentrality(0.0);
+    }
+
+    // Brandes' algorithm: run single-source Dijkstra from each node
+    for (auto const& [sourceId, sourceNode] : m_nodes) {
+      std::stack<Id> S;  // nodes in order of non-increasing distance
+      // predecessors: P[w] = list of (predecessor node id, edge id from pred to w)
+      std::unordered_map<Id, std::vector<std::pair<Id, Id>>> P;
+      std::unordered_map<Id, double> sigma;  // number of shortest paths
+      std::unordered_map<Id, double> dist;   // distance from source
+
+      for (auto const& [nId, _] : m_nodes) {
+        P[nId] = {};
+        sigma[nId] = 0.0;
+        dist[nId] = std::numeric_limits<double>::infinity();
+      }
+      sigma[sourceId] = 1.0;
+      dist[sourceId] = 0.0;
+
+      // Min-heap priority queue: (distance, nodeId)
+      std::priority_queue<std::pair<double, Id>,
+                          std::vector<std::pair<double, Id>>,
+                          std::greater<>>
+          pq;
+      pq.push({0.0, sourceId});
+
+      std::unordered_set<Id> visited;
+
+      while (!pq.empty()) {
+        auto [d, v] = pq.top();
+        pq.pop();
+
+        if (visited.contains(v)) {
+          continue;
+        }
+        visited.insert(v);
+        S.push(v);
+
+        for (auto const& eId : m_nodes.at(v)->outgoingEdges()) {
+          auto const& pEdge = m_edges.at(eId);
+          Id w = pEdge->target();
+          if (visited.contains(w)) {
+            continue;
+          }
+          double edgeWeight = getEdgeWeight(pEdge);
+          double newDist = dist[v] + edgeWeight;
+
+          if (newDist < dist[w]) {
+            dist[w] = newDist;
+            sigma[w] = sigma[v];
+            P[w] = {{v, eId}};
+            pq.push({newDist, w});
+          } else if (std::abs(newDist - dist[w]) < 1e-12 * std::max(1.0, dist[w])) {
+            sigma[w] += sigma[v];
+            P[w].push_back({v, eId});
+          }
+        }
+      }
+
+      // Dependency accumulation (backward pass)
+      std::unordered_map<Id, double> delta;
+      for (auto const& [nId, _] : m_nodes) {
+        delta[nId] = 0.0;
+      }
+
+      while (!S.empty()) {
+        Id w = S.top();
+        S.pop();
+        for (auto const& [v, eId] : P[w]) {
+          double c = (sigma[v] / sigma[w]) * (1.0 + delta[w]);
+          delta[v] += c;
+          // Accumulate edge betweenness
+          auto currentBC = m_edges.at(eId)->betweennessCentrality();
+          m_edges.at(eId)->setBetweennessCentrality(*currentBC + c);
+        }
+      }
+    }
+  }
 }  // namespace dsf

src/dsf/base/Node.hpp

@@ -31,6 +31,7 @@ namespace dsf {
     std::string m_name;
     std::vector<Id> m_ingoingEdges;
     std::vector<Id> m_outgoingEdges;
+    std::optional<double> m_betweennessCentrality{std::nullopt};
 
   public:
     /// @brief Construct a new Node object with capacity 1
@@ -47,7 +48,8 @@ namespace dsf {
           m_geometry{other.m_geometry},
           m_name{other.m_name},
           m_ingoingEdges{other.m_ingoingEdges},
-          m_outgoingEdges{other.m_outgoingEdges} {}
+          m_outgoingEdges{other.m_outgoingEdges},
+          m_betweennessCentrality{other.m_betweennessCentrality} {}
     virtual ~Node() = default;
 
     Node& operator=(Node const& other) {
@@ -57,6 +59,7 @@ namespace dsf {
         m_name = other.m_name;
         m_ingoingEdges = other.m_ingoingEdges;
         m_outgoingEdges = other.m_outgoingEdges;
+        m_betweennessCentrality = other.m_betweennessCentrality;
       }
       return *this;
     }
@@ -98,24 +101,31 @@ namespace dsf {
       }
       m_outgoingEdges.push_back(edgeId);
     }
+    /// @brief Set the node's betweenness centrality
+    /// @param betweennessCentrality The node's betweenness centrality
+    inline void setBetweennessCentrality(double const betweennessCentrality) noexcept {
+      m_betweennessCentrality = betweennessCentrality;
+    }
 
     /// @brief Get the node's id
     /// @return Id The node's id
-    inline Id id() const { return m_id; }
+    inline auto id() const { return m_id; }
     /// @brief Get the node's geometry
     /// @return std::optional<geometry::Point> A geometry::Point
-    inline std::optional<geometry::Point> const& geometry() const noexcept {
-      return m_geometry;
-    }
+    inline auto const& geometry() const noexcept { return m_geometry; }
     /// @brief Get the node's name
     /// @return std::string The node's name
-    inline std::string const& name() const noexcept { return m_name; }
-
-    inline std::vector<Id> const& ingoingEdges() const noexcept { return m_ingoingEdges; }
-    inline std::vector<Id> const& outgoingEdges() const noexcept {
-      return m_outgoingEdges;
+    inline auto const& name() const noexcept { return m_name; }
+    /// @brief Get the node's ingoing edges
+    /// @return std::vector<Id> A vector of the node's ingoing edge ids
+    inline auto const& ingoingEdges() const noexcept { return m_ingoingEdges; }
+    /// @brief Get the node's outgoing edges
+    /// @return std::vector<Id> A vector of the node's outgoing edge ids
+    inline auto const& outgoingEdges() const noexcept { return m_outgoingEdges; }
+    /// @brief Get the node's betweenness centrality
+    /// @return std::optional<double> The node's betweenness centrality, or std::nullopt if not set
+    inline auto const& betweennessCentrality() const noexcept {
+      return m_betweennessCentrality;
     }
-
-    virtual bool isStation() const noexcept { return false; }
   };
 };  // namespace dsf