Update BC computation to address normalization edge conditions (#5105)

networkx/networkx#7949 implements an update to normalization that was originally described in #4941.  We pushed an update in #4974 that addressed the specific examples that the user identified in the original cugraph issue.  However, while @eriknw was exploring updates to networkx to address this, he identified a few more edge conditions that needed to be satisfied.  This PR addresses those remaining edge conditions.

Note that the python tests comparing results to networkx are still disabled.  These can't be re-enabled until networkx/networkx#7949 is included in a networkx release.

Closes #5006 
Closes #5107

Authors:
  - Chuck Hastings (https://github.com/ChuckHastings)
  - Rick Ratzel (https://github.com/rlratzel)

Approvers:
  - Seunghwa Kang (https://github.com/seunghwak)
  - Joseph Nke (https://github.com/jnke2016)
  - Rick Ratzel (https://github.com/rlratzel)

URL: #5105

Author: ChuckHastings

cpp/src/centrality/betweenness_centrality_impl.cuh

@@ -540,37 +540,57 @@ rmm::device_uvector<weight_t> betweenness_centrality(
                               do_expensive_check);
   }
 
-  std::optional<weight_t> scale_factor{std::nullopt};
+  std::optional<weight_t> scale_nonsource{std::nullopt};
+  std::optional<weight_t> scale_source{std::nullopt};
+
+  weight_t num_vertices = static_cast<weight_t>(graph_view.number_of_vertices());
+  if (!include_endpoints) num_vertices = num_vertices - 1;
+
+  if ((static_cast<edge_t>(num_sources) == num_vertices) || include_endpoints) {
+    if (normalized) {
+      scale_nonsource = static_cast<weight_t>(num_sources * (num_vertices - 1));
+    } else if (graph_view.is_symmetric()) {
+      scale_nonsource =
+        static_cast<weight_t>(num_sources * 2) / static_cast<weight_t>(num_vertices);
+    } else {
+      scale_nonsource = static_cast<weight_t>(num_sources) / static_cast<weight_t>(num_vertices);
+    }
 
-  if (normalized) {
-    if (include_endpoints) {
-      if (graph_view.number_of_vertices() >= 2) {
-        scale_factor = static_cast<weight_t>(
-          std::min(static_cast<vertex_t>(num_sources), graph_view.number_of_vertices()) *
-          (graph_view.number_of_vertices() - 1));
-      }
-    } else if (graph_view.number_of_vertices() > 2) {
-      scale_factor = static_cast<weight_t>(
-        std::min(static_cast<vertex_t>(num_sources), graph_view.number_of_vertices() - 1) *
-        (graph_view.number_of_vertices() - 2));
+    scale_source = scale_nonsource;
+  } else if (normalized) {
+    scale_nonsource = static_cast<weight_t>(num_sources) * (num_vertices - 1);
+    scale_source    = static_cast<weight_t>(num_sources - 1) * (num_vertices - 1);
+  } else {
+    scale_nonsource = static_cast<weight_t>(num_sources) / num_vertices;
+    scale_source    = static_cast<weight_t>(num_sources - 1) / num_vertices;
+
+    if (graph_view.is_symmetric()) {
+      *scale_nonsource *= 2;
+      *scale_source *= 2;
     }
-  } else if (num_sources < static_cast<size_t>(graph_view.number_of_vertices())) {
-    if ((graph_view.number_of_vertices() > 1) && (num_sources > 0))
-      scale_factor =
-        (graph_view.is_symmetric() ? weight_t{2} : weight_t{1}) *
-        static_cast<weight_t>(num_sources) /
-        (include_endpoints ? static_cast<weight_t>(graph_view.number_of_vertices())
-                           : static_cast<weight_t>(graph_view.number_of_vertices() - 1));
-  } else if (graph_view.is_symmetric()) {
-    scale_factor = weight_t{2};
   }
 
-  if (scale_factor) {
-    thrust::transform(handle.get_thrust_policy(),
-                      centralities.begin(),
-                      centralities.end(),
-                      centralities.begin(),
-                      [sf = *scale_factor] __device__(auto centrality) { return centrality / sf; });
+  if (scale_nonsource) {
+    auto iter = thrust::make_zip_iterator(
+      thrust::make_counting_iterator(graph_view.local_vertex_partition_range_first()),
+      centralities.begin());
+
+    thrust::transform(
+      handle.get_thrust_policy(),
+      iter,
+      iter + centralities.size(),
+      centralities.begin(),
+      [nonsource = *scale_nonsource,
+       source    = *scale_source,
+       vertices_begin,
+       vertices_end] __device__(auto t) {
+        vertex_t v          = thrust::get<0>(t);
+        weight_t centrality = thrust::get<1>(t);
+
+        return (thrust::find(thrust::seq, vertices_begin, vertices_end, v) == vertices_end)
+                 ? centrality / nonsource
+                 : centrality / source;
+      });
   }
 
   return centralities;

cpp/tests/c_api/betweenness_centrality_test.c

@@ -113,9 +113,16 @@ int generic_betweenness_centrality_test(vertex_t* h_src,
   TEST_ASSERT(test_ret_value, ret_code == CUGRAPH_SUCCESS, "copy_to_host failed.");
 
   for (int i = 0; (i < num_vertices) && (test_ret_value == 0); ++i) {
-    TEST_ASSERT(test_ret_value,
-                nearlyEqual(h_result[h_vertices[i]], h_centralities[i], 0.0001),
-                "centralities results don't match");
+    if (isnan(h_result[h_vertices[i]])) {
+      TEST_ASSERT(test_ret_value, isnan(h_centralities[i]), "expected NaN, got a non-NaN value");
+    } else {
+      if (!nearlyEqual(h_result[h_vertices[i]], h_centralities[i], 0.0001))
+        printf("  expected: %g, got %g\n", h_result[h_vertices[i]], h_centralities[i]);
+
+      TEST_ASSERT(test_ret_value,
+                  nearlyEqual(h_result[h_vertices[i]], h_centralities[i], 0.0001),
+                  "centralities results don't match");
+    }
   }
 
   cugraph_centrality_result_free(p_result);
@@ -169,7 +176,7 @@ int test_betweenness_centrality_specific_normalized()
   weight_t h_wgt[] = {
     0.1f, 2.1f, 1.1f, 5.1f, 3.1f, 4.1f, 7.2f, 3.2f, 0.1f, 2.1f, 1.1f, 5.1f, 3.1f, 4.1f, 7.2f, 3.2f};
   vertex_t h_seeds[]  = {0, 3};
-  weight_t h_result[] = {0, 0.395833, 0.16667, 0.0833333, 0.0416667, 0.0625};
+  weight_t h_result[] = {0, 0.395833, 0.166667, 0.166667, 0.0416667, 0.0625};
 
   return generic_betweenness_centrality_test(h_src,
                                              h_dst,
@@ -197,7 +204,7 @@ int test_betweenness_centrality_specific_unnormalized()
   weight_t h_wgt[] = {
     0.1f, 2.1f, 1.1f, 5.1f, 3.1f, 4.1f, 7.2f, 3.2f, 0.1f, 2.1f, 1.1f, 5.1f, 3.1f, 4.1f, 7.2f, 3.2f};
   vertex_t h_seeds[]  = {0, 3};
-  weight_t h_result[] = {0, 7.91667, 3.33333, 1.666667, 0.833333, 1.25};
+  weight_t h_result[] = {0, 7.91667, 3.33333, 3.33333, 0.833333, 1.25};
 
   return generic_betweenness_centrality_test(h_src,
                                              h_dst,
@@ -312,17 +319,17 @@ int test_issue_4941()
     {TRUE, TRUE, FALSE, 0, {1.0, 0.4, 0.4, 0.4, 0.4}},
     {TRUE, TRUE, FALSE, 1, {1.0, 1.0, 0.25, 0.25, 0.25}},
     {TRUE, FALSE, TRUE, 0, {1.0, 0.0, 0.0, 0.0, 0.0}},
-    {TRUE, FALSE, TRUE, 1, {1.0, 0.0, 0.0, 0.0, 0.0}},
+    {TRUE, FALSE, TRUE, 1, {1.0, NAN, 0.0, 0.0, 0.0}},
     {TRUE, FALSE, FALSE, 0, {1.0, 0.0, 0.0, 0.0, 0.0}},
-    {TRUE, FALSE, FALSE, 1, {1.0, 0.0, 0.0, 0.0, 0.0}},
+    {TRUE, FALSE, FALSE, 1, {1.0, NAN, 0.0, 0.0, 0.0}},
     {FALSE, TRUE, TRUE, 0, {20.0, 8.0, 8.0, 8.0, 8.0}},
     {FALSE, TRUE, TRUE, 1, {20.0, 20.0, 5.0, 5.0, 5.0}},
     {FALSE, TRUE, FALSE, 0, {10.0, 4.0, 4.0, 4.0, 4.0}},
     {FALSE, TRUE, FALSE, 1, {10.0, 10.0, 2.5, 2.5, 2.5}},
     {FALSE, FALSE, TRUE, 0, {12.0, 0.0, 0.0, 0.0, 0.0}},
-    {FALSE, FALSE, TRUE, 1, {12.0, 0.0, 0.0, 0.0, 0.0}},
+    {FALSE, FALSE, TRUE, 1, {12, NAN, 0.0, 0.0, 0.0}},
     {FALSE, FALSE, FALSE, 0, {6.0, 0.0, 0.0, 0.0, 0.0}},
-    {FALSE, FALSE, FALSE, 1, {6.0, 0.0, 0.0, 0.0, 0.0}},
+    {FALSE, FALSE, FALSE, 1, {6.0, NAN, 0.0, 0.0, 0.0}},
   };
 
   int test_result = 0;

cpp/tests/centrality/betweenness_centrality_reference.hpp

@@ -130,41 +130,54 @@ void ref_edge_accumulation(std::vector<weight_t>& result,
   }
 }
 
-template <typename result_t>
+template <typename vertex_t, typename result_t>
 void reference_rescale(result_t* result,
+                       vertex_t const* sources,
                        bool directed,
                        bool normalize,
                        bool endpoints,
                        size_t const number_of_vertices,
                        size_t const number_of_sources)
 {
-  result_t rescale_factor            = static_cast<result_t>(1);
   result_t casted_number_of_sources  = static_cast<result_t>(number_of_sources);
   result_t casted_number_of_vertices = static_cast<result_t>(number_of_vertices);
+  if (!endpoints) casted_number_of_vertices = casted_number_of_vertices - 1;
 
-  if (normalize) {
-    if (number_of_vertices > 2) {
-      if (endpoints) {
-        rescale_factor /=
-          (number_of_sources > 0 ? casted_number_of_sources
-                                 : casted_number_of_vertices * (casted_number_of_vertices - 1));
-      } else {
-        rescale_factor /= (number_of_sources > 0
-                             ? casted_number_of_sources
-                             : (casted_number_of_vertices - 1) * (casted_number_of_vertices - 2));
-      }
+  if ((number_of_sources == number_of_vertices) || endpoints) {
+    result_t rescale_factor = static_cast<result_t>(1);
+
+    if (normalize) {
+      rescale_factor = result_t{1} / (casted_number_of_sources * (casted_number_of_vertices - 1));
+    } else if (!directed) {
+      rescale_factor = casted_number_of_vertices / (2 * casted_number_of_sources);
+    } else {
+      rescale_factor = casted_number_of_vertices / casted_number_of_sources;
     }
-  } else if (number_of_sources < number_of_vertices) {
-    rescale_factor = (endpoints ? casted_number_of_vertices : casted_number_of_vertices - 1) /
-                     (directed ? casted_number_of_sources : 2 * casted_number_of_sources);
-  } else if (!directed) {
-    rescale_factor = 2;
-  }
 
-  if (rescale_factor != result_t{1}) {
-    for (auto idx = 0; idx < number_of_vertices; ++idx) {
+    for (vertex_t idx = 0; idx < number_of_vertices; ++idx) {
       result[idx] *= rescale_factor;
     }
+  } else {
+    result_t rescale_source     = static_cast<result_t>(1);
+    result_t rescale_non_source = static_cast<result_t>(1);
+
+    if (normalize) {
+      rescale_source     = 1 / ((casted_number_of_sources - 1) * (casted_number_of_vertices - 1));
+      rescale_non_source = 1 / (casted_number_of_sources * (casted_number_of_vertices - 1));
+    } else if (directed) {
+      rescale_source     = casted_number_of_vertices / (casted_number_of_sources - 1);
+      rescale_non_source = casted_number_of_vertices / casted_number_of_sources;
+    } else {
+      rescale_source     = casted_number_of_vertices / (2 * (casted_number_of_sources - 1));
+      rescale_non_source = casted_number_of_vertices / (2 * casted_number_of_sources);
+    }
+
+    for (vertex_t idx = 0; idx < number_of_vertices; ++idx) {
+      if (std::find(sources, sources + number_of_sources, idx) == (sources + number_of_sources))
+        result[idx] *= rescale_non_source;
+      else
+        result[idx] *= rescale_source;
+    }
   }
 }
 
@@ -235,8 +248,13 @@ std::vector<weight_t> betweenness_centrality_reference(
     }
   }
 
-  reference_rescale(
-    result.data(), directed, normalize, include_endpoints, offsets.size() - 1, seeds.size());
+  reference_rescale(result.data(),
+                    seeds.data(),
+                    directed,
+                    normalize,
+                    include_endpoints,
+                    offsets.size() - 1,
+                    seeds.size());
 
   return result;
 }

python/cugraph/cugraph/tests/centrality/test_batch_betweenness_centrality_mg.py

@@ -15,6 +15,7 @@
 
 import pytest
 import numpy as np
+import networkx as nx
 
 from cugraph.dask.common.mg_utils import is_single_gpu
 from cugraph.datasets import karate
@@ -55,7 +56,10 @@ def setup_function():
 # =============================================================================
 
 
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.mg
 @pytest.mark.skipif(is_single_gpu(), reason="skipping MG testing on Single GPU system")
 @pytest.mark.parametrize("dataset", DATASETS)

python/cugraph/cugraph/tests/centrality/test_batch_edge_betweenness_centrality_mg.py

@@ -15,6 +15,7 @@
 
 import pytest
 import numpy as np
+import networkx as nx
 
 from cugraph.dask.common.mg_utils import is_single_gpu
 from cugraph.datasets import karate, netscience
@@ -53,7 +54,10 @@ def setup_function():
 
 
 # FIXME: Fails for directed = False(bc score twice as much) and normalized = True.
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.mg
 @pytest.mark.skipif(is_single_gpu(), reason="skipping MG testing on Single GPU system")
 @pytest.mark.parametrize("dataset", DATASETS)

python/cugraph/cugraph/tests/centrality/test_betweenness_centrality.py

@@ -19,7 +19,6 @@
 import networkx as nx
 
 import cudf
-import cupy
 import cugraph
 from cugraph.datasets import karate_disjoint
 from cugraph.testing import utils, SMALL_DATASETS
@@ -102,7 +101,7 @@ def calc_betweenness_centrality(
         Contains 'vertex' and  'cu_bc' 'ref_bc' columns,  where 'cu_bc'
         and 'ref_bc' are the two betweenness centrality scores to compare.
         The dataframe is expected to be sorted based on 'vertex', so that we
-        can use cupy.isclose to compare the scores.
+        can use np.isclose to compare the scores.
     """
     G = None
     Gnx = None
@@ -289,8 +288,15 @@ def _calc_bc_full(G, Gnx, normalized, weight, endpoints, k, seed, result_dtype):
 # i.e: sorted_df[idx][first_key] should be compared to
 #      sorted_df[idx][second_key]
 def compare_scores(sorted_df, first_key, second_key, epsilon=DEFAULT_EPSILON):
+    # Compare with numpy and pandas since presence of NaNs in cudf Series
+    # results in "ValueError: CuPy currently does not support masked arrays."
     errors = sorted_df[
-        ~cupy.isclose(sorted_df[first_key], sorted_df[second_key], rtol=epsilon)
+        ~np.isclose(
+            sorted_df[first_key].to_pandas(),
+            sorted_df[second_key].to_pandas(),
+            rtol=epsilon,
+            equal_nan=True,
+        )
     ]
     num_errors = len(errors)
     if num_errors > 0:
@@ -305,7 +311,10 @@ def compare_scores(sorted_df, first_key, second_key, epsilon=DEFAULT_EPSILON):
 # =============================================================================
 # Tests
 # =============================================================================
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.sg
 @pytest.mark.parametrize("graph_file", SMALL_DATASETS)
 @pytest.mark.parametrize("directed", [False, True])
@@ -542,17 +551,17 @@ def test_betweenness_centrality_nx(graph_file, directed, edgevals):
         (True, True, False, None, {0: 1.0, 1: 0.4, 2: 0.4, 3: 0.4, 4: 0.4}),
         (True, True, False, 1, {0: 1.0, 1: 1.0, 2: 0.25, 3: 0.25, 4: 0.25}),
         (True, False, True, None, {0: 1.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
-        (True, False, True, 1, {0: 1.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
+        (True, False, True, 1, {0: 1.0, 1: np.nan, 2: 0.0, 3: 0.0, 4: 0.0}),
         (True, False, False, None, {0: 1.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
-        (True, False, False, 1, {0: 1.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
+        (True, False, False, 1, {0: 1.0, 1: np.nan, 2: 0.0, 3: 0.0, 4: 0.0}),
         (False, True, True, None, {0: 20.0, 1: 8.0, 2: 8.0, 3: 8.0, 4: 8.0}),
         (False, True, True, 1, {0: 20.0, 1: 20.0, 2: 5.0, 3: 5.0, 4: 5.0}),
         (False, True, False, None, {0: 10.0, 1: 4.0, 2: 4.0, 3: 4.0, 4: 4.0}),
         (False, True, False, 1, {0: 10.0, 1: 10.0, 2: 2.5, 3: 2.5, 4: 2.5}),
         (False, False, True, None, {0: 12.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
-        (False, False, True, 1, {0: 12.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
+        (False, False, True, 1, {0: 12.0, 1: np.nan, 2: 0.0, 3: 0.0, 4: 0.0}),
         (False, False, False, None, {0: 6.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
-        (False, False, False, 1, {0: 6.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0}),
+        (False, False, False, 1, {0: 6.0, 1: np.nan, 2: 0.0, 3: 0.0, 4: 0.0}),
     ],
 )
 def test_scale_with_k_on_star_graph(normalized, endpoints, is_directed, k, expected):

python/cugraph/cugraph/tests/centrality/test_edge_betweenness_centrality.py

@@ -312,7 +312,10 @@ def generate_upper_triangle(dataframe):
     return dataframe
 
 
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.sg
 @pytest.mark.parametrize("graph_file", SMALL_DATASETS)
 @pytest.mark.parametrize("directed", DIRECTED_GRAPH_OPTIONS)
@@ -343,7 +346,10 @@ def test_edge_betweenness_centrality(
     compare_scores(sorted_df, first_key="cu_bc", second_key="ref_bc")
 
 
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.sg
 @pytest.mark.parametrize("graph_file", SMALL_DATASETS)
 @pytest.mark.parametrize("directed", DIRECTED_GRAPH_OPTIONS)
@@ -383,7 +389,10 @@ def test_edge_betweenness_centrality_k_full(
 #       the function operating the comparison inside is first proceeding
 #       to a random sampling over the number of vertices (thus direct offsets)
 #       in the graph structure instead of actual vertices identifiers
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.sg
 @pytest.mark.parametrize("graph_file", [karate_disjoint])
 @pytest.mark.parametrize("directed", DIRECTED_GRAPH_OPTIONS)
@@ -487,7 +496,10 @@ def test_edge_betweenness_invalid_dtype(
         compare_scores(sorted_df, first_key="cu_bc", second_key="ref_bc")
 
 
-@pytest.mark.skip(reason="https://github.com/networkx/networkx/pull/7908")
+@pytest.mark.skipif(
+    float(".".join(nx.__version__.split(".")[:2])) < 3.5,
+    reason="Requires networkx >= 3.5",
+)
 @pytest.mark.sg
 @pytest.mark.parametrize("graph_file", SMALL_DATASETS)
 @pytest.mark.parametrize("directed", DIRECTED_GRAPH_OPTIONS)