Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions docs/references.bib
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,17 @@ @article{omnipath
doi = {10.1038/nmeth.4077}
}

@article{rustworkx,
author = {Treinish, Matthew and Carvalho, Ivan and Tsilimigkounakis, Georgios and Sá, Nahum},
title = {rustworkx: A High-Performance Graph Library for Python},
journal = {Journal of Open Source Software},
volume = {7},
number = {79},
pages = {3968},
year = {2022},
doi = {10.21105/joss.03968},
}

@InProceedings{networkx,
author = {Aric A. Hagberg and Daniel A. Schult and Pieter J. Swart},
title = {Exploring Network Structure, Dynamics, and Function using NetworkX},
Expand Down
1 change: 1 addition & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ dependencies = [
"pillow>=8",
"pooch>=1.6",
"pyyaml>=6",
"rustworkx>=0.18",
"scanpy>=1.9.3",
"scikit-image>=0.25",
# due to https://github.com/scikit-image/scikit-image/issues/6850 breaks rescale ufunc
Expand Down
2 changes: 1 addition & 1 deletion src/squidpy/_docs.py
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,7 @@ def decorator2(obj: Any) -> Any:
scalebar_kwargs
Keyword arguments for :class:`matplotlib_scalebar.ScaleBar`.
edges_kwargs
Keyword arguments for :func:`networkx.drawing.nx_pylab.draw_networkx_edges`.
Keyword arguments for :func:`networkx.drawing.nx_pylab.draw_networkx_edges` :cite:`networkx`.
kwargs
Keyword arguments for :func:`matplotlib.pyplot.scatter` or :func:`matplotlib.pyplot.imshow`.
"""
Expand Down
81 changes: 74 additions & 7 deletions src/squidpy/gr/_nhood.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,16 @@
from functools import partial
from typing import Any, NamedTuple

import networkx as nx
import numba.types as nt
import numpy as np
import pandas as pd
import rustworkx as rx
from anndata import AnnData
from numba import njit
from numba import njit, prange
from numpy.typing import NDArray
from pandas import CategoricalDtype
from scanpy import logging as logg
from scipy.sparse import csr_matrix
from spatialdata import SpatialData

from squidpy._constants._constants import Centrality
Expand Down Expand Up @@ -266,7 +267,7 @@ def centrality_scores(
%(table_key)s
%(cluster_key)s
score
Centrality measures as described in :mod:`networkx.algorithms.centrality` :cite:`networkx`.
Group centrality measures as implemented in ``rustworkx`` :cite:`rustworkx`.
If `None`, use all the options below. Valid options are:

- `{c.CLOSENESS.s!r}` - measure of how close the group is to other nodes.
Expand Down Expand Up @@ -296,19 +297,23 @@ def centrality_scores(

centralities = [Centrality(c) for c in centrality]

graph = nx.Graph(adata.obsp[connectivity_key])
# a rustworkx graph mirrors the undirected connectivity graph for the group closeness/degree
# measures; a symmetric, self-loop-free CSR feeds the clustering-coefficient kernel.
graph, adj = _build_graph(adata.obsp[connectivity_key])

cat = adata.obs[cluster_key].cat.categories.values
clusters = adata.obs[cluster_key].values

fun_dict = {}
for c in centralities:
if c == Centrality.CLOSENESS:
fun_dict[c.s] = partial(nx.algorithms.centrality.group_closeness_centrality, graph)
fun_dict[c.s] = partial(rx.group_closeness_centrality, graph)
elif c == Centrality.DEGREE:
fun_dict[c.s] = partial(nx.algorithms.centrality.group_degree_centrality, graph)
fun_dict[c.s] = partial(rx.group_degree_centrality, graph)
elif c == Centrality.CLUSTERING:
fun_dict[c.s] = partial(nx.algorithms.cluster.average_clustering, graph)
# average the per-node clustering coefficients over the group (0 if the group is empty).
node_clustering = _local_clustering(adj.indptr, adj.indices, adj.shape[0])
fun_dict[c.s] = lambda idx, cc=node_clustering: float(cc[idx].mean()) if len(idx) else 0.0
else:
raise NotImplementedError(f"Centrality `{c}` is not yet implemented.")

Expand Down Expand Up @@ -424,6 +429,68 @@ def _interaction_matrix(
return output


def _build_graph(conn: Any) -> tuple[rx.PyGraph, csr_matrix]:
"""Build the graph representations used by :func:`centrality_scores`.

Returns a :class:`rustworkx.PyGraph` mirroring the undirected connectivity graph
(used by the group closeness/degree measures) and a symmetric, self-loop-free,
index-sorted CSR matrix feeding the clustering-coefficient kernel.
"""
from scipy.sparse import triu

adj = csr_matrix(conn)
# undirected, unweighted, no self-loops: matches ``networkx.Graph(conn)`` topology.
adj = (adj + adj.T).tocsr()
adj.setdiag(0)
adj.eliminate_zeros()
adj.sort_indices() # the clustering kernel merges neighbor lists, which must be sorted.

n = adj.shape[0]
graph = rx.PyGraph(multigraph=False)
graph.add_nodes_from(range(n))
# the strict upper triangle lists each undirected edge exactly once.
rows, cols = triu(adj, k=1).nonzero()
graph.add_edges_from_no_data([(int(i), int(j)) for i, j in zip(rows, cols, strict=True)])
return graph, adj


@njit(parallel=True, cache=True)
def _local_clustering(indptr: NDArrayA, indices: NDArrayA, n: int) -> NDArrayA:
"""Local clustering coefficient per node over a symmetric, index-sorted CSR graph.

``c_v = 2 * T(v) / (k_v * (k_v - 1))`` where ``T(v)`` is the number of edges among the
neighbors of ``v`` and ``k_v`` its degree; ``c_v = 0`` when ``k_v < 2``. Triangles are
counted by intersecting sorted neighbor lists, so only existing edges are ever visited
(no length-2-path matrix is materialized). Matches :func:`networkx.clustering`.
"""
out = np.zeros(n, dtype=np.float64)
for v in prange(n):
start = indptr[v]
end = indptr[v + 1]
k = end - start
if k < 2:
continue
# summing |N(u) ∩ N(v)| over u in N(v) counts each neighbor-neighbor edge twice,
# so it already equals 2 * T(v) -> c_v = that sum / (k * (k - 1)).
two_triangles = 0
for a in range(start, end):
u = indices[a]
i = start
j = indptr[u]
u_end = indptr[u + 1]
while i < end and j < u_end:
if indices[i] == indices[j]:
two_triangles += 1
i += 1
j += 1
elif indices[i] < indices[j]:
i += 1
else:
j += 1
out[v] = two_triangles / (k * (k - 1))
return out


def _centrality_scores_helper(
cat: Iterable[Any],
clusters: Sequence[str],
Expand Down
17 changes: 17 additions & 0 deletions tests/graph/test_nhood.py
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,23 @@ def test_centrality_scores(nhood_data: AnnData):
assert adata.uns[key]["closeness_centrality"].dtype == np.dtype("float64")


def test_centrality_scores_networkx_parity(nhood_data: AnnData):
# centrality_scores swapped networkx for rustworkx (+ a numba clustering kernel); pin the
# numeric parity of all three group measures against networkx (still a dependency).
import networkx as nx

adata = nhood_data
df = centrality_scores(adata, cluster_key=_CK, connectivity_key="spatial", copy=True)

graph = nx.Graph(adata.obsp["spatial_connectivities"])
clusters = adata.obs[_CK].values
for cat in df.index:
idx = list(np.where(clusters == cat)[0])
np.testing.assert_allclose(df.loc[cat, "closeness_centrality"], nx.group_closeness_centrality(graph, idx))
np.testing.assert_allclose(df.loc[cat, "degree_centrality"], nx.group_degree_centrality(graph, idx))
np.testing.assert_allclose(df.loc[cat, "average_clustering"], nx.average_clustering(graph, idx))


@pytest.mark.parametrize("copy", [True, False])
def test_interaction_matrix_copy(nhood_data: AnnData, copy: bool):
adata = nhood_data
Expand Down
Loading