dijkstra keeping ties

README.md

@@ -8,12 +8,12 @@
   * [x] shortest path
   * [x] A*
   * [x] Yen's k-shortest
-  * [ ] all-pairs shortest path
+  * [x] all-pairs shortest path
   * [x] single-source shortest path
   * [x] minimum spanning tree
   * [ ] random walking
   * [x] degree centrality
-  * [ ] closeness centrality
+  * [x] closeness centrality
   * [ ] betweenness centrality
   * [ ] pagerank
   * [ ] triangle counting

src/algo/all_pairs_shortest_path.rs

@@ -1,8 +1,68 @@
 use std::cmp::Reverse;
+use std::collections::BTreeMap;
 
+use anyhow::{anyhow, bail};
 use ordered_float::OrderedFloat;
 use priority_queue::PriorityQueue;
 use rayon::prelude::*;
+use smartstring::{LazyCompact, SmartString};
+
+use crate::algo::AlgoImpl;
+use crate::data::expr::Expr;
+use crate::data::program::{MagicAlgoRuleArg, MagicSymbol};
+use crate::data::tuple::Tuple;
+use crate::data::value::DataValue;
+use crate::runtime::derived::DerivedRelStore;
+use crate::runtime::transact::SessionTx;
+
+pub(crate) struct ClosenessCentrality;
+
+impl AlgoImpl for ClosenessCentrality {
+    fn run(
+        &mut self,
+        tx: &SessionTx,
+        rels: &[MagicAlgoRuleArg],
+        opts: &BTreeMap<SmartString<LazyCompact>, Expr>,
+        stores: &BTreeMap<MagicSymbol, DerivedRelStore>,
+        out: &DerivedRelStore,
+    ) -> anyhow::Result<()> {
+        let edges = rels
+            .get(0)
+            .ok_or_else(|| anyhow!("'closeness_centrality' requires edges relation"))?;
+        let undirected = match opts.get("undirected") {
+            None => false,
+            Some(Expr::Const(DataValue::Bool(b))) => *b,
+            Some(v) => bail!(
+                "option 'undirected' for 'closeness_centrality' requires a boolean, got {:?}",
+                v
+            ),
+        };
+
+        let (graph, indices, _inv_indices, _) =
+            edges.convert_edge_to_weighted_graph(undirected, false, tx, stores)?;
+
+        let n = graph.len();
+        if n == 0 {
+            return Ok(());
+        }
+        let res: Vec<_> = (0..n)
+            .into_par_iter()
+            .map(|start| -> f64 {
+                let distances = dijkstra_cost_only(&graph, start);
+                let total_dist: f64 = distances.iter().filter(|d| d.is_finite()).cloned().sum();
+                let nc: f64 = distances.iter().filter(|d| d.is_finite()).count() as f64;
+                nc * nc / total_dist / (n - 1) as f64
+            })
+            .collect();
+        for (idx, centrality) in res.into_iter().enumerate() {
+            out.put(
+                Tuple(vec![indices[idx].clone(), DataValue::from(centrality)]),
+                0,
+            );
+        }
+        Ok(())
+    }
+}
 
 pub(crate) fn dijkstra_cost_only(edges: &[Vec<(usize, f64)>], start: usize) -> Vec<f64> {
     let mut distance = vec![f64::INFINITY; edges.len()];
@@ -28,10 +88,3 @@ pub(crate) fn dijkstra_cost_only(edges: &[Vec<(usize, f64)>], start: usize) -> V
 
     distance
 }
-
-pub(crate) fn apsp(edges: &[Vec<(usize, f64)>]) -> Vec<Vec<f64>> {
-    (0..edges.len())
-        .into_par_iter()
-        .map(|start| dijkstra_cost_only(edges, start))
-        .collect()
-}

src/algo/mod.rs

@@ -3,6 +3,7 @@ use std::collections::BTreeMap;
 use anyhow::{anyhow, bail, ensure, Result};
 use itertools::Itertools;
 use smartstring::{LazyCompact, SmartString};
+use crate::algo::all_pairs_shortest_path::ClosenessCentrality;
 
 use crate::algo::astar::ShortestPathAStar;
 use crate::algo::bfs::Bfs;
@@ -62,6 +63,7 @@ impl AlgoHandle {
     pub(crate) fn arity(&self) -> Result<usize> {
         Ok(match &self.name.0 as &str {
             "degree_centrality" => 4,
+            "closeness_centrality" => 2,
             "depth_first_search" | "dfs" => 1,
             "breadth_first_search" | "bfs" => 1,
             "shortest_path_dijkstra" => 4,
@@ -80,6 +82,7 @@ impl AlgoHandle {
     pub(crate) fn get_impl(&self) -> Result<Box<dyn AlgoImpl>> {
         Ok(match &self.name.0 as &str {
             "degree_centrality" => Box::new(DegreeCentrality),
+            "closeness_centrality" => Box::new(ClosenessCentrality),
             "depth_first_search" | "dfs" => Box::new(Dfs),
             "breadth_first_search" | "bfs" => Box::new(Bfs),
             "shortest_path_dijkstra" => Box::new(ShortestPathDijkstra),

src/algo/shortest_path_dijkstra.rs

@@ -4,10 +4,10 @@ use std::iter;
 
 use anyhow::{anyhow, bail, Result};
 use itertools::Itertools;
-use smartstring::{LazyCompact, SmartString};
-
 use ordered_float::OrderedFloat;
 use priority_queue::PriorityQueue;
+use smallvec::{smallvec, SmallVec};
+use smartstring::{LazyCompact, SmartString};
 
 use crate::algo::AlgoImpl;
 use crate::data::expr::Expr;
@@ -43,6 +43,14 @@ impl AlgoImpl for ShortestPathDijkstra {
                 v
             ),
         };
+        let keep_ties = match opts.get("keep_ties") {
+            None => false,
+            Some(Expr::Const(DataValue::Bool(b))) => *b,
+            Some(v) => bail!(
+                "option 'keep_ties' for 'shortest_path_dijkstra' requires a boolean, got {:?}",
+                v
+            ),
+        };
 
         let (graph, indices, inv_indices, _) =
             edges.convert_edge_to_weighted_graph(undirected, false, tx, stores)?;
@@ -80,9 +88,17 @@ impl AlgoImpl for ShortestPathDijkstra {
             let res = if let Some(tn) = &termination_nodes {
                 if tn.len() == 1 {
                     let single = Some(*tn.iter().next().unwrap());
-                    dijkstra(&graph, start, &single, &(), &())
+                    if keep_ties {
+                        dijkstra_keep_ties(&graph, start, &single, &(), &())
+                    } else {
+                        dijkstra(&graph, start, &single, &(), &())
+                    }
                 } else {
-                    dijkstra(&graph, start, tn, &(), &())
+                    if keep_ties {
+                        dijkstra_keep_ties(&graph, start, tn, &(), &())
+                    } else {
+                        dijkstra(&graph, start, tn, &(), &())
+                    }
                 }
             } else {
                 dijkstra(&graph, start, &(), &(), &())
@@ -273,3 +289,91 @@ pub(crate) fn dijkstra<FE: ForbiddenEdge, FN: ForbiddenNode, G: Goal + Clone>(
 
     ret
 }
+
+pub(crate) fn dijkstra_keep_ties<FE: ForbiddenEdge, FN: ForbiddenNode, G: Goal + Clone>(
+    edges: &[Vec<(usize, f64)>],
+    start: usize,
+    goals: &G,
+    forbidden_edges: &FE,
+    forbidden_nodes: &FN,
+) -> Vec<(usize, f64, Vec<usize>)> {
+    let mut distance = vec![f64::INFINITY; edges.len()];
+    let mut pq = PriorityQueue::new();
+    let mut back_pointers: Vec<SmallVec<[usize; 1]>> = vec![smallvec![]; edges.len()];
+    distance[start] = 0.;
+    pq.push(start, Reverse(OrderedFloat(0.)));
+    let mut goals_remaining = goals.clone();
+
+    while let Some((node, Reverse(OrderedFloat(cost)))) = pq.pop() {
+        if cost > distance[node] {
+            continue;
+        }
+
+        for (nxt_node, path_weight) in &edges[node] {
+            if forbidden_nodes.is_forbidden(*nxt_node) {
+                continue;
+            }
+            if forbidden_edges.is_forbidden(node, *nxt_node) {
+                continue;
+            }
+            let nxt_cost = cost + *path_weight;
+            if nxt_cost < distance[*nxt_node] {
+                pq.push_increase(*nxt_node, Reverse(OrderedFloat(nxt_cost)));
+                distance[*nxt_node] = nxt_cost;
+                back_pointers[*nxt_node].clear();
+                back_pointers[*nxt_node].push(node);
+            } else if nxt_cost == distance[*nxt_node] {
+                pq.push_increase(*nxt_node, Reverse(OrderedFloat(nxt_cost)));
+                back_pointers[*nxt_node].push(node);
+            }
+        }
+
+        goals_remaining.visit(node);
+        if goals_remaining.is_exhausted() {
+            break;
+        }
+    }
+
+    let ret = goals
+        .iter(edges.len())
+        .flat_map(|target| {
+            let cost = distance[target];
+            if !cost.is_finite() {
+                vec![(target, cost, vec![])]
+            } else {
+                struct CollectPath {
+                    collected: Vec<(usize, f64, Vec<usize>)>,
+                }
+
+                impl CollectPath {
+                    fn collect(
+                        &mut self,
+                        chain: &[usize],
+                        start: usize,
+                        target: usize,
+                        cost: f64,
+                        back_pointers: &[SmallVec<[usize; 1]>],
+                    ) {
+                        let last = chain.last().unwrap();
+                        let prevs = &back_pointers[*last];
+                        for nxt in prevs {
+                            let mut ret = chain.to_vec();
+                            ret.push(*nxt);
+                            if *nxt == start {
+                                ret.reverse();
+                                self.collected.push((target, cost, ret));
+                            } else {
+                                self.collect(&ret, start, target, cost, back_pointers)
+                            }
+                        }
+                    }
+                }
+                let mut cp = CollectPath { collected: vec![] };
+                cp.collect(&[target], start, target, cost, &back_pointers);
+                cp.collected
+            }
+        })
+        .collect_vec();
+
+    ret
+}