cds2.graph

Graph algorithms on scipy.sparse.csgraph plus PageRank.

Graph algorithms built on scipy.sparse.csgraph, plus PageRank.

DegreeResult dataclass

In/out/total degrees of each node.

Source code in src/cds2/graph.py

@dataclass(frozen=True)
class DegreeResult:
    """In/out/total degrees of each node."""

    in_degree: IntArray
    out_degree: IntArray
    total: IntArray

ComponentResult dataclass

Connected-component labeling of a graph.

Source code in src/cds2/graph.py

@dataclass(frozen=True)
class ComponentResult:
    """Connected-component labeling of a graph."""

    count: int
    labels: IntArray

CentralityResult dataclass

Per-node centrality scores.

Source code in src/cds2/graph.py

@dataclass(frozen=True)
class CentralityResult:
    """Per-node centrality scores."""

    scores: FloatArray

CommunityResult dataclass

Label-propagation community detection output.

Source code in src/cds2/graph.py

@dataclass(frozen=True)
class CommunityResult:
    """Label-propagation community detection output."""

    labels: IntArray
    n_communities: int
    modularity: float

from_edges(n, edges, weighted=False, directed=True)

Build a CSR adjacency matrix from an edge list.

Unweighted edges get weight 1. Undirected graphs mirror every edge.

Source code in src/cds2/graph.py

def from_edges(
    n: int,
    edges: Sequence[tuple[int, int] | tuple[int, int, float]],
    weighted: bool = False,
    directed: bool = True,
) -> sparse.csr_matrix:
    """Build a CSR adjacency matrix from an edge list.

    Unweighted edges get weight 1. Undirected graphs mirror every edge.
    """
    if n <= 0:
        msg = "n must be positive"
        raise ValueError(msg)
    if not edges:
        return sparse.csr_matrix((n, n), dtype=float)
    rows: list[int] = []
    cols: list[int] = []
    data: list[float] = []
    for edge in edges:
        rows.append(int(edge[0]))
        cols.append(int(edge[1]))
        data.append(float(edge[2]) if weighted and len(edge) > 2 else 1.0)
    for u, v in zip(rows, cols, strict=False):
        if not (0 <= u < n and 0 <= v < n):
            msg = f"edge ({u}, {v}) references a node outside 0..{n - 1}"
            raise ValueError(msg)
    adj = sparse.csr_matrix((data, (rows, cols)), shape=(n, n))
    if not directed:
        adj = adj.maximum(adj.T)
    return sparse.csr_matrix(adj)

degree(adj, directed=True)

Node degrees from a CSR/COO adjacency matrix.

Source code in src/cds2/graph.py

def degree(adj: object, directed: bool = True) -> DegreeResult:
    """Node degrees from a CSR/COO adjacency matrix."""
    matrix = sparse.csr_matrix(adj)
    out_degree = (
        np.asarray(matrix.sum(axis=1)).ravel().astype(np.int64)
        if _is_unweighted(matrix)
        else np.asarray((matrix != 0).sum(axis=1)).ravel().astype(np.int64)
    )
    in_degree = (
        np.asarray(matrix.sum(axis=0)).ravel().astype(np.int64)
        if _is_unweighted(matrix)
        else np.asarray((matrix != 0).sum(axis=0)).ravel().astype(np.int64)
    )
    total = in_degree + out_degree if directed else out_degree.copy()
    if not directed:
        in_degree = out_degree.copy()
        total = out_degree.copy()
    return DegreeResult(in_degree=in_degree, out_degree=out_degree, total=total)

connected_components(adj, directed=False, connection='weak')

Label the connected components of a graph.

Source code in src/cds2/graph.py

def connected_components(
    adj: object, directed: bool = False, connection: str = "weak"
) -> ComponentResult:
    """Label the connected components of a graph."""
    count, labels = cs_connected_components(
        sparse.csr_matrix(adj), directed=directed, connection=connection
    )
    return ComponentResult(count=int(count), labels=np.asarray(labels, dtype=np.int64))

shortest_paths(adj, indices=None, method='dijkstra', unweighted=False)

Pairwise shortest paths from indices to all nodes (Dijkstra or Bellman-Ford).

Source code in src/cds2/graph.py

def shortest_paths(
    adj: object, indices: IntArray | None = None, method: str = "dijkstra", unweighted: bool = False
) -> ShortestPaths:
    """Pairwise shortest paths from ``indices`` to all nodes (Dijkstra or Bellman-Ford)."""
    matrix = sparse.csr_matrix(adj)
    source = None if indices is None else np.asarray(indices, dtype=np.int32)
    if method == "dijkstra":
        dist: Any = cs_dijkstra(matrix, indices=source, unweighted=unweighted)
    elif method == "bellman_ford":
        dist = bellman_ford(matrix, indices=source)
    else:
        msg = f"unsupported shortest-path method: {method!r}"
        raise ValueError(msg)
    return ShortestPaths(method=method, distances=np.asarray(dist, dtype=float))

single_source_shortest_paths(adj, source, method='dijkstra')

Distance vector from one node to all others.

Source code in src/cds2/graph.py

def single_source_shortest_paths(adj: object, source: int, method: str = "dijkstra") -> FloatArray:
    """Distance vector from one node to all others."""
    result = shortest_paths(adj, indices=np.array([source]), method=method)
    return np.asarray(result.distances[0], dtype=float)

floyd_warshall_paths(adj)

All-pairs shortest paths via the Floyd-Warshall algorithm.

Source code in src/cds2/graph.py

def floyd_warshall_paths(adj: object) -> ShortestPaths:
    """All-pairs shortest paths via the Floyd-Warshall algorithm."""
    dist = cs_floyd_warshall(sparse.csr_matrix(adj))
    return ShortestPaths(method="floyd_warshall", distances=np.asarray(dist))

bellman_ford_paths(adj, indices=None)

Shortest paths tolerating negative weights (no negative cycles).

Source code in src/cds2/graph.py

def bellman_ford_paths(adj: object, indices: IntArray | None = None) -> ShortestPaths:
    """Shortest paths tolerating negative weights (no negative cycles)."""
    return shortest_paths(adj, indices=indices, method="bellman_ford")

minimum_spanning_forest(adj)

Minimum spanning tree/forest of an undirected weighted graph.

Source code in src/cds2/graph.py

def minimum_spanning_forest(adj: object) -> sparse.csr_matrix:
    """Minimum spanning tree/forest of an undirected weighted graph."""
    symmetric = sparse.csr_matrix(adj).maximum(sparse.csr_matrix(adj).T)
    upper = sparse.triu(symmetric, format="csr")
    return sparse.csr_matrix(minimum_spanning_tree(upper))

pagerank(adj, damping=0.85, max_iter=100, tol=1e-10)

PageRank scores for a finite square graph with non-negative weights.

Source code in src/cds2/graph.py

def pagerank(
    adj: object, damping: float = 0.85, max_iter: int = 100, tol: float = 1e-10
) -> FloatArray:
    """PageRank scores for a finite square graph with non-negative weights."""
    coo = sparse.coo_matrix(adj, dtype=np.float64)
    if coo.shape[0] != coo.shape[1]:
        raise ValueError("adjacency matrix must be square")
    if not np.isfinite(damping) or not (0.0 < damping < 1.0):
        raise ValueError("damping must be strictly between 0 and 1")
    if max_iter < 1:
        raise ValueError("max_iter must be at least 1")
    if not np.isfinite(tol) or tol <= 0.0:
        raise ValueError("tol must be a positive finite number")
    coo.sum_duplicates()
    weights = np.asarray(coo.data, dtype=np.float64)
    if not bool(np.all(np.isfinite(weights))):
        raise ValueError("edge weights must be finite")
    if bool(np.any(weights < 0.0)):
        raise ValueError("PageRank requires non-negative edge weights")
    n = coo.shape[0]
    if n == 0:
        return np.zeros(0, dtype=float)
    source_nodes = coo.row.astype(np.int64, copy=False)
    out_degree = np.bincount(source_nodes, weights=weights, minlength=n).astype(float)
    normalized = weights / np.where(out_degree == 0.0, 1.0, out_degree)[source_nodes]
    order = np.argsort(coo.col, kind="stable")
    sorted_targets = coo.col[order]
    follow_indices = np.ascontiguousarray(source_nodes[order], dtype=np.int64)
    follow_data = np.ascontiguousarray(normalized[order], dtype=np.float64)
    follow_indptr = np.ascontiguousarray(
        np.searchsorted(sorted_targets, np.arange(n + 1)).astype(np.int64)
    )
    dangling_indices = np.ascontiguousarray(np.flatnonzero(out_degree == 0), dtype=np.int64)

    if _HAS_PR_KERNEL and _pr_kernel is not None:
        rank_buffer, _iterations = _pr_kernel.iterate(
            follow_indptr,
            follow_indices,
            follow_data,
            n,
            damping,
            dangling_indices,
            max_iter,
            tol,
        )
        rank_vec = np.frombuffer(rank_buffer, dtype=np.float64).copy()
    else:
        rank_vec = np.full(n, 1.0 / n)
        follow_matrix = sparse.csr_matrix(
            (follow_data, follow_indices, follow_indptr), shape=(n, n)
        )
        teleport = (1.0 - damping) / n
        for _ in range(max_iter):
            dangling_mass = float(rank_vec.take(dangling_indices).sum())
            new_rank = damping * (follow_matrix @ rank_vec)
            new_rank += damping * dangling_mass / n
            new_rank += teleport
            delta = float(np.abs(new_rank - rank_vec).max())
            rank_vec = np.asarray(new_rank, dtype=np.float64)
            if delta < tol:
                break
    total = float(rank_vec.sum())
    if not np.isfinite(total) or total <= 0.0:
        raise FloatingPointError("PageRank iteration produced an invalid probability vector")
    return np.asarray(rank_vec / total, dtype=float)

topological_order(n, edges)

Kahn topological order of a DAG; raises ValueError when a cycle exists.

Source code in src/cds2/graph.py

def topological_order(n: int, edges: list[tuple[int, int]]) -> list[int]:
    """Kahn topological order of a DAG; raises ValueError when a cycle exists."""
    if n <= 0:
        msg = "n must be positive"
        raise ValueError(msg)
    adjacency: dict[int, list[int]] = {i: [] for i in range(n)}
    indegree = np.zeros(n, dtype=np.int64)
    for u, v in edges:
        adjacency[u].append(v)
        indegree[v] += 1
    queue: deque[int] = deque(int(i) for i in np.flatnonzero(indegree == 0))
    order: list[int] = []
    while queue:
        node = queue.popleft()
        order.append(node)
        for neighbor in adjacency[node]:
            indegree[neighbor] -= 1
            if indegree[neighbor] == 0:
                queue.append(neighbor)
    if len(order) != n:
        msg = "cycle detected: graph is not a DAG"
        raise ValueError(msg)
    return order

betweenness_centrality(adj, normalized=True)

Brandes betweenness centrality on the unweighted successor graph.

Directed graphs use n=1 pairs; undirected graphs divide by two.

Source code in src/cds2/graph.py

def betweenness_centrality(adj: object, normalized: bool = True) -> FloatArray:
    """Brandes betweenness centrality on the unweighted successor graph.

    Directed graphs use n=1 pairs; undirected graphs divide by two.
    """
    matrix = sparse.csr_matrix(adj)
    n = matrix.shape[0]
    if n == 0:
        return np.zeros(0)
    successors_list: list[list[int]] = [[] for _ in range(n)]
    coo = matrix.tocoo()
    for u, v in zip(coo.row.tolist(), coo.col.tolist(), strict=False):
        if u != v and matrix[u, v] != 0:
            successors_list[u].append(v)

    betweenness = np.zeros(n)
    for source in range(n):
        stack: deque[int] = deque()
        predecessors: list[list[int]] = [[] for _ in range(n)]
        path_counts = np.zeros(n)
        path_counts[source] = 1.0
        distances = np.full(n, -1.0)
        distances[source] = 0.0
        visit_queue: deque[int] = deque([source])
        while visit_queue:
            node = visit_queue.popleft()
            stack.append(node)
            for neighbor in successors_list[node]:
                if distances[neighbor] < 0.0:
                    distances[neighbor] = distances[node] + 1.0
                    visit_queue.append(neighbor)
                if distances[neighbor] == distances[node] + 1.0:
                    path_counts[neighbor] += path_counts[node]
                    predecessors[neighbor].append(node)
        dependencies = np.zeros(n)
        while stack:
            target = stack.pop()
            for predecessor in predecessors[target]:
                share = path_counts[predecessor] / path_counts[target]
                dependencies[predecessor] += share * (1.0 + dependencies[target])
            if target != source:
                betweenness[target] += dependencies[target]

    # Undirected graphs are counted twice (s->t and t->s), so halve once.
    # Normalized scale is (n-1)(n-2)/2 for undirected, (n-1)(n-2) for directed.
    # The old code used (n-1)(n-2) for both and then halved again, capping an
    # undirected star centre at 0.5 where the standard definition gives 1.0.
    is_undirected = not directed_flag(matrix)
    if is_undirected:
        betweenness /= 2.0
    if normalized:
        if n <= 2:
            scale = 1.0
        elif is_undirected:
            scale = (n - 1) * (n - 2) / 2.0
        else:
            scale = (n - 1) * (n - 2)
        betweenness /= scale
    return np.asarray(betweenness, dtype=float)

directed_flag(matrix)

True when the matrix is not symmetric.

Source code in src/cds2/graph.py

def directed_flag(matrix: sparse.csr_matrix) -> bool:
    """True when the matrix is not symmetric."""
    difference = abs(matrix - matrix.T)
    return float(difference.sum()) > 1e-12

closeness_centrality(adj, wf_improved=True)

Wasserman-Faust closeness centrality from single-source distances.

Source code in src/cds2/graph.py

def closeness_centrality(adj: object, wf_improved: bool = True) -> FloatArray:
    """Wasserman-Faust closeness centrality from single-source distances."""
    matrix = sparse.csr_matrix(adj)
    n = matrix.shape[0]
    if n == 0:
        return np.zeros(0)
    distances = cs_dijkstra(matrix, unweighted=_is_unweighted(matrix))
    scores = np.zeros(n)
    for node in range(n):
        row = distances[node]
        reachable = row[np.isfinite(row) & (row > 0)]
        total = float(reachable.sum())
        if total > 0:
            raw_score = reachable.size / total
            if wf_improved and reachable.size < n - 1:
                raw_score *= reachable.size / (n - 1)
            scores[node] = raw_score
    return scores

eigenvector_centrality(adj, max_iter=200, tol=1e-10, directed=None)

Eigenvector centrality of the largest (Perron) eigenvalue.

Power iteration with an oscillation guard: bipartite graphs that fail to settle fall back to a Lanczos solve for the dominant eigenvector.

Source code in src/cds2/graph.py

def eigenvector_centrality(
    adj: object, max_iter: int = 200, tol: float = 1e-10, directed: bool | None = None
) -> FloatArray:
    """Eigenvector centrality of the largest (Perron) eigenvalue.

    Power iteration with an oscillation guard: bipartite graphs that fail to
    settle fall back to a Lanczos solve for the dominant eigenvector.
    """
    coo = sparse.coo_matrix(adj)
    n = coo.shape[0]
    if n == 0:
        return np.zeros(0)
    is_directed = directed if directed is not None else directed_flag(sparse.csr_matrix(adj))
    matrix = sparse.csr_matrix(coo)
    if is_directed:
        matrix = matrix.maximum(matrix.T).tocsr()
    vector = np.full(n, 1.0 / n)
    for _ in range(max_iter):
        new_vector = np.asarray(matrix @ vector, dtype=float).ravel()
        norm_value = float(np.linalg.norm(new_vector, 1))
        if norm_value == 0:
            return np.zeros(n)
        new_vector /= norm_value
        delta = float(np.abs(new_vector - vector).max())
        vector = new_vector
        if delta < tol:
            break
    else:
        from cds2.sparse import largest_eigenpairs

        result = largest_eigenpairs(matrix, k=1)
        vector = np.abs(np.asarray(result.eigenvectors[:, 0], dtype=float)).ravel()
        total = float(vector.sum())
        if total == 0:  # pragma: no cover - defensive: Lanczos returns a unit vector
            return np.zeros(n)
        return np.asarray(vector / total, dtype=float)
    total = vector.sum()
    return np.asarray(vector / total, dtype=float)

detect_communities(adj, max_sweeps=50, seed=None)

Label-propagation communities with modularity as the quality score.

Deterministic given seed; ties in the neighbour-label vote are broken by the rng so repeated calls with the same seed reproduce exactly.

Source code in src/cds2/graph.py

def detect_communities(
    adj: object, max_sweeps: int = 50, seed: int | None = None
) -> CommunityResult:
    """Label-propagation communities with modularity as the quality score.

    Deterministic given ``seed``; ties in the neighbour-label vote are broken
    by the rng so repeated calls with the same seed reproduce exactly.
    """
    matrix = sparse.csr_matrix(adj)
    n = matrix.shape[0]
    if n == 0:
        return CommunityResult(labels=np.zeros(0, dtype=np.int64), n_communities=0, modularity=0.0)
    neighbours = _undirected_neighbours(matrix)
    weights = np.asarray(np.abs(matrix).sum(axis=1)).ravel()
    total_weight = float(weights.sum())
    if total_weight == 0:
        labels = np.arange(n, dtype=np.int64)
        return CommunityResult(labels=labels, n_communities=n, modularity=0.0)

    rng = np.random.default_rng(seed)
    labels = np.arange(n, dtype=np.int64)
    order = rng.permutation(n)
    for _ in range(max_sweeps):
        changed = False
        for node in order:
            if not neighbours[node]:
                continue
            votes: dict[int, float] = {}
            for neighbour in neighbours[node]:
                weight = float(matrix[node, neighbour])
                key = int(labels[neighbour])
                votes[key] = votes.get(key, 0.0) + weight
            best_weight = max(votes.values())
            winners = sorted(
                label for label, weight in votes.items() if abs(weight - best_weight) < 1e-15
            )
            chosen = int(rng.choice(winners)) if len(winners) > 1 else winners[0]
            if chosen != labels[node]:
                labels[node] = chosen
                changed = True
        if not changed:
            break
    return CommunityResult(
        labels=labels,
        n_communities=int(np.unique(labels).size),
        modularity=float(modularity(adj, labels)),
    )

modularity(adj, labels)

Newman modularity Q of a partition on an undirected graph.

Source code in src/cds2/graph.py

def modularity(adj: object, labels: Sequence[int] | IntArray) -> float:
    """Newman modularity Q of a partition on an undirected graph."""
    matrix = sparse.csr_matrix(adj)
    symmetric = matrix.maximum(matrix.T).tocsr()
    n = symmetric.shape[0]
    label_array = np.asarray(labels, dtype=np.int64)
    if label_array.size != n:
        msg = "one label per node is required"
        raise ValueError(msg)
    degrees = np.asarray(symmetric.sum(axis=1)).ravel()
    two_m = float(degrees.sum())
    if two_m == 0:
        return 0.0

    same = label_array[:, None] == label_array[None, :]
    # Both sums count every undirected edge twice, matching the standard
    # form Q = sum_c [ e_c / m - (d_c / 2m)^2 ].
    edge_within = float(np.asarray(symmetric[same]).sum())
    degree_within = float((degrees[:, None] * degrees[None, :] * same).sum())
    return edge_within / two_m - degree_within / two_m**2