Draft PR: Add modularity and modularity_adata functions to scanpy.metrics

src/scanpy/metrics/__init__.py

@@ -3,7 +3,7 @@
 from __future__ import annotations
 
 from ._gearys_c import gearys_c
-from ._metrics import confusion_matrix
+from ._metrics import confusion_matrix, modularity, modularity_adata
 from ._morans_i import morans_i
 
-__all__ = ["confusion_matrix", "gearys_c", "morans_i"]
+__all__ = ["confusion_matrix", "gearys_c", "modularity", "modularity_adata", "morans_i"]

src/scanpy/metrics/_metrics.py

@@ -2,15 +2,22 @@
 
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, cast
 
 import numpy as np
 import pandas as pd
 from natsort import natsorted
 from pandas.api.types import CategoricalDtype
+from scipy.sparse import coo_matrix
+
+from .._compat import CSRBase, SpBase
 
 if TYPE_CHECKING:
     from collections.abc import Sequence
+    from typing import Literal
+
+    from anndata import AnnData
+    from numpy.typing import ArrayLike
 
 
 def confusion_matrix(
@@ -60,7 +67,9 @@ def confusion_matrix(
     orig, new = pd.Series(orig), pd.Series(new)
     assert len(orig) == len(new)
 
-    unique_labels = pd.unique(np.concatenate((orig.values, new.values)))
+    unique_labels = pd.unique(
+        np.concatenate((np.asarray(orig.values), np.asarray(new.values)))
+    )
 
     # Compute
     mtx = _confusion_matrix(orig, new, labels=unique_labels)
@@ -89,3 +98,95 @@ def confusion_matrix(
     df = df.loc[np.array(orig_idx), np.array(new_idx)]
 
     return df
+
+
+def modularity(
+    connectivities: ArrayLike | SpBase,
+    labels: pd.Series | ArrayLike,
+    mode: Literal["UNDIRECTED", "DIRECTED"] = "UNDIRECTED",
+) -> float:
+    # accepting both dense or spare matrices as the connectivity graph
+    # setting mode between directed and undirected
+    """Compute the modularity of a graph given its connectivities and labels.
+
+    Parameters
+    ----------
+    connectivities: array-like or sparse matrix
+        Weighted adjacency matrix representing the graph. Can be a dense NumPy array or a sparse CSR matrix.
+    labels: array-like or pandas.Series
+        Cluster labels for each node in the graph.
+    mode: str
+        The mode of the graph. Can be "UNDIRECTED" or "DIRECTED". Default is "UNDIRECTED".
+
+    Returns
+    -------
+    float
+        The modularity of the graph based on the provided clustering.
+    """
+    try:
+        # try to import igraph in case the user wants to calculate modularity
+        # not in the main module to avoid import errors
+        import igraph as ig
+    except ImportError as e:
+        msg = "igraph is require for computing modularity"
+        raise ImportError(msg) from e
+    if isinstance(connectivities, CSRBase | SpBase):
+        # check if the connectivities is a sparse matrix
+        coo = coo_matrix(connectivities)
+        edges = list(zip(coo.row, coo.col, strict=False))
+        # converting to the coo format to extract the edges and weights
+        # storing only non-zero elements and their indices
+        weights = coo.data.tolist()
+        graph = ig.Graph(edges=edges, directed=mode == "DIRECTED")
+        graph.es["weight"] = weights
+    else:
+        # if the graph is dense, creates it directly using igraph's adjacency matrix
+        dense_array = np.asarray(connectivities)
+        igraph_mode = ig.ADJ_UNDIRECTED if mode == "UNDIRECTED" else ig.ADJ_DIRECTED
+        graph = ig.Graph.Weighted_Adjacency(dense_array.tolist(), mode=igraph_mode)
+    # cluster labels to integer codes required by igraph
+    labels = pd.Categorical(np.asarray(labels)).codes
+
+    return graph.modularity(labels)
+
+
+def modularity_adata(
+    adata: AnnData,
+    *,
+    labels: str | ArrayLike = "leiden",
+    obsp: str = "connectivities",
+    mode: Literal["UNDIRECTED", "DIRECTED"] = "UNDIRECTED",
+) -> float:
+    # default to leiden labels and connectivities as it is more common
+    """Compute modularity from an AnnData object using stored graph and clustering labels.
+
+    Parameters
+    ----------
+    adata: AnnData
+        The AnnData object containing the data.
+    labels: str or array-like
+        The key in adata.obs that contains the cluster labels.
+    obsp: str
+        The key in adata.obsp that contains the connectivities.
+
+    Returns
+    -------
+    float
+        The modularity of the graph based on the provided clustering.
+    """
+    # if labels is a key in adata.obs, get the values from adata.obs
+    # otherwise, assume it is an array-like object
+    label_array = adata.obs[labels] if isinstance(labels, str) else labels
+    connectivities = adata.obsp[obsp]
+
+    if isinstance(connectivities, CSRBase):
+        # converting to dense if it is a CSR matrix
+        dense = connectivities.toarray()
+    else:
+        toarray = getattr(connectivities, "toarray", None)
+        dense = toarray() if callable(toarray) else connectivities
+
+    if isinstance(dense, pd.DataFrame):
+        dense = dense.values
+    dense = cast("ArrayLike", dense)
+    return modularity(dense, label_array, mode=mode)

src/scanpy/neighbors/__init__.py

@@ -21,9 +21,11 @@
 from .._utils import NeighborsView, _doc_params, get_literal_vals
 from . import _connectivity
 from ._common import (
+    _get_indices_distances_from_dense_matrix,
     _get_indices_distances_from_sparse_matrix,
     _get_sparse_matrix_from_indices_distances,
 )
+from ._connectivity import umap
 from ._doc import doc_n_pcs, doc_use_rep
 from ._types import _KnownTransformer, _Method
 
@@ -74,6 +76,7 @@ def neighbors(  # noqa: PLR0913
     n_neighbors: int = 15,
     n_pcs: int | None = None,
     *,
+    distance_matrix: np.ndarray | None = None,
     use_rep: str | None = None,
     knn: bool = True,
     method: _Method = "umap",
@@ -186,6 +189,13 @@ def neighbors(  # noqa: PLR0913
     :doc:`/how-to/knn-transformers`
 
     """
+    # if distance_matrix is not None:
+    #     return neighbors_from_distance(
+    #         adata,
+    #         distance_matrix,
+    #         n_neighbors=n_neighbors,
+    #         method=method,
+    #     )
     start = logg.info("computing neighbors")
     adata = adata.copy() if copy else adata
     if adata.is_view:  # we shouldn't need this here...
@@ -248,6 +258,65 @@ def neighbors(  # noqa: PLR0913
     return adata if copy else None
 
 
+# def neighbors_from_distance(
+#     adata: AnnData,
+#     distance_matrix: np.ndarray,
+#     n_neighbors: int = 15,
+#     method: Literal["umap", "gauss"] = "umap",  # default to umap
+# ) -> None:
+#     # computes the neighborhood graph from a precomputed distance matrix
+#     # only umap is supported
+#     # skipping PCA and distance computation and goes straight to the graph
+#     """Compute neighbors from a precomputer distance matrix.
+
+#     Parameters
+#     ----------
+#     adata : AnnData
+#         Annotated data matrix.
+#     distance_matrix : np.ndarray
+#         Precomputed dense or sparse distance matrix.
+#     n_neighbors : int
+#         Number of nearest neighbors to use in the graph.
+#     method : str
+#         Method to use for computing the graph. Currently only 'umap' is supported.
+#     """
+#     if isinstance(distance_matrix, SpBase):
+#         distance_matrix = np.asarray(sparse.csr_matrix(distance_matrix).toarray())
+
+#     knn_indices, knn_distances = _get_indices_distances_from_dense_matrix(
+#         distance_matrix, n_neighbors
+#     )
+#     if method == "umap":
+#         connectivities = umap(
+#             knn_indices,
+#             knn_distances,
+#             n_obs=adata.n_obs,
+#             n_neighbors=n_neighbors,
+#         )
+#     elif method == "gauss":
+#         connectivities = _connectivity.gauss(
+#             sparse.csr_matrix(distance_matrix),
+#             n_neighbors,
+#             knn=True,
+#         )
+#     else:
+#         msg = f"Method {method} not implemented."
+#         raise NotImplementedError(msg)
+
+#     adata.obsp["connectivities"] = connectivities
+#     adata.obsp["distances"] = sparse.csr_matrix(distance_matrix)
+#     adata.uns["neighbors"] = {
+#         "connectivities_key": "connectivities",
+#         "distances_key": "distances",
+#         "params": {
+#             "n_neighbors": n_neighbors,
+#             "method": method,
+#             "random_state": 0,
+#             "metric": "euclidean",
+#         },
+#     }
+
+
 class FlatTree(NamedTuple):  # noqa: D101
     hyperplanes: None
     offsets: None

tests/test_metrics.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import importlib.util
 import warnings
 from functools import partial
 from string import ascii_letters
@@ -10,8 +11,10 @@
 import pytest
 import threadpoolctl
 from scipy import sparse
+from scipy.sparse import csr_matrix  # noqa: TID251
 
 import scanpy as sc
+from scanpy.metrics import modularity, modularity_adata
 from testing.scanpy._helpers.data import pbmc68k_reduced
 from testing.scanpy._pytest.params import ARRAY_TYPES
 
@@ -196,3 +199,123 @@ def test_confusion_matrix_api():
     pd.testing.assert_frame_equal(
         expected, sc.metrics.confusion_matrix(data["a"], "b", data)
     )
+
+
+# importing igraph, louvain, leiden if available
+HAS_IGRAPH = importlib.util.find_spec("igraph") is not None
+HAS_LOUVAIN = importlib.util.find_spec("louvain") is not None
+HAS_LEIDEN = importlib.util.find_spec("leidenalg") is not None
+
+
+# Test 1: Sample graph with clear community structure (dense & sparse, directed & undirected)
+@pytest.mark.parametrize(
+    "mode", ["UNDIRECTED", "DIRECTED"], ids=["undirected", "directed"]
+)
+@pytest.mark.parametrize("use_sparse", [False, True], ids=["sparse", "dense"])
+def test_modularity_sample_structure(mode, use_sparse):
+    if HAS_IGRAPH is False:
+        pytest.skip("igraph is not installed")
+    # 4 node adjacency matrix with two separate 2-node communities
+    mat = np.array(
+        [
+            [1, 1, 0, 0],
+            [1, 1, 0, 0],
+            [0, 0, 1, 1],
+            [0, 0, 1, 1],
+        ]
+    )
+    labels = ["A", "A", "B", "B"]
+    adj = csr_matrix(mat) if use_sparse else mat
+    score = modularity(adj, labels, mode=mode)
+
+    # Modularity should be between 0 and 1
+    assert 0 <= score <= 1
+
+
+# Test 2: Edge case when all nodes belong to the same community/cluster
+def test_modularity_single_community():
+    if HAS_IGRAPH is False:
+        pytest.skip("igraph is not installed")
+    # fully connected graph sample
+    adj = np.ones((4, 4)) - np.eye(4)
+    labels = ["A", "A", "A", "A"]
+    score = modularity(adj, labels)
+
+    # modularity should be 0
+    assert score == pytest.approx(0.0, rel=1e-6)
+
+
+# Test 3: Invalad input, labels length does not match adjacency matrix size
+def test_modularity_invalid_labels():
+    if HAS_IGRAPH is False:
+        pytest.skip("igraph is not installed")
+    from igraph._igraph import InternalError
+
+    adj = np.eye(4)
+    labels = ["A", "A", "B"]
+    with pytest.raises(
+        InternalError,
+        match="Membership vector size differs",
+    ):
+        modularity(adj, labels)
+
+
+# Test 4: Pass both Louvain and Leiden clustering algorithms
+@pytest.mark.parametrize("cluster_method", ["louvain", "leiden"])
+def test_modularity_adata_multiple_clusterings(cluster_method):
+    if HAS_IGRAPH is False:
+        pytest.skip("igraph is not installed")
+    if cluster_method == "louvain" and HAS_LOUVAIN is False:
+        pytest.skip("louvain is not installed")
+    if cluster_method == "leiden" and HAS_LEIDEN is False:
+        pytest.skip("leiden is not installed")
+    # Loading PBMC Data and compute PCA and neighbors graph
+    adata = sc.datasets.pbmc3k()
+    sc.pp.pca(adata)
+    sc.pp.neighbors(adata)
+    # Compute modularity using both Louvain and Leiden clustering
+    if cluster_method == "louvain":
+        sc.tl.louvain(adata)
+        score_louvain = modularity_adata(
+            adata, labels="louvain", obsp="connectivities", mode="UNDIRECTED"
+        )
+        # Score should be between 0 and 1
+        assert 0 <= score_louvain <= 1
+    if cluster_method == "leiden":
+        sc.tl.leiden(adata)
+        score_leiden = modularity_adata(
+            adata, labels="leiden", obsp="connectivities", mode="UNDIRECTED"
+        )
+        # Score should be between 0 and 1
+        assert 0 <= score_leiden <= 1
+
+
+# Test 5: Modularity should be the same no matter the order of the labels
+def test_modularity_order():
+    if HAS_IGRAPH is False:
+        pytest.skip("igraph is not installed")
+    adj = np.array(
+        [
+            [1, 1, 0, 0],
+            [1, 1, 0, 0],
+            [0, 0, 1, 1],
+            [0, 0, 1, 1],
+        ]
+    )
+    labels1 = ["A", "A", "B", "B"]
+    labels2 = ["B", "B", "A", "A"]
+    score_1 = modularity(adj, labels1)
+    score_2 = modularity(adj, labels2)
+    assert score_1 == score_2
+
+
+# Test 6: Modularity on disconnected graph lke edge-case behavior in some algorithms
+def test_modularity_disconnected_graph():
+    if HAS_IGRAPH is False:
+        pytest.skip("igraph is not installed")
+    adj = np.zeros((4, 4))
+    labels = ["A", "B", "C", "D"]
+    score = modularity(adj, labels)
+
+    # Modularity should be undefined for disconnected graphs
+    assert np.isnan(score)