"""
Hypergraph backbone methods.
These methods operate directly on hypergraphs (collections of hyperedges, where
each hyperedge is a set of nodes of arbitrary size) rather than on dyadic graphs.
A hyperedge is given as any iterable of node labels; a hypergraph is an iterable
of hyperedges, for example ``[(1, 2, 3), (2, 4), (1, 2, 3, 5)]``. NetworkX has
no native hypergraph type, so inputs and outputs use plain Python collections
(the backbone is returned as a list of :class:`frozenset` hyperedges).
Methods
-------
mdl_hypergraph_backbone
Parameter-free information-theoretic (MDL) backbone that prunes nested and
redundant hyperedges (Kirkley, Felippe, Malizia & Battiston, 2026).
local_mdl_hypergraph_backbone
Local variant that backbones each node neighborhood separately and
aggregates the results (Kirkley et al., 2026, Appendix E).
hypergraph_compression_ratio
Inverse compression ratio ``eta`` achieved by the MDL backbone.
intersection_graph
Graph linking hyperedges that share at least one node.
"""
from __future__ import annotations
import math
from dataclasses import dataclass, field
from typing import TYPE_CHECKING
import networkx as nx
from nx_backbone._docstrings import append_complexity_docstrings
if TYPE_CHECKING:
from collections.abc import Hashable, Iterable, Iterator
from nx_backbone._typing import Hyperedge, Hyperedges
__all__ = [
"intersection_graph",
"mdl_hypergraph_backbone",
"local_mdl_hypergraph_backbone",
"hypergraph_compression_ratio",
"HypergraphBackbone",
"maximal_hyperedges",
"order_filter",
"s_components",
"statistically_validated_hypergraph",
"statistically_validated_cores",
"ValidatedHypergraph",
"svh",
"svc",
]
_LN2 = math.log(2.0)
# ---------------------------------------------------------------------------
# Information-theoretic primitives (all codelengths in bits, log base 2)
# ---------------------------------------------------------------------------
def _log2(x):
return math.log2(x)
def _log2_factorial(n):
"""log2(n!) via the log-gamma function (stdlib, no numpy/scipy needed)."""
return math.lgamma(n + 1.0) / _LN2
def _log2_binom(n, k):
"""log2 of the binomial coefficient C(n, k)."""
if k < 0 or k > n or n < 0:
return float("-inf")
return (
math.lgamma(n + 1.0) - math.lgamma(k + 1.0) - math.lgamma(n - k + 1.0)
) / _LN2
def _parent_codelength(size, n_orders, n_nodes):
"""H(p): bits to transmit a parent (backbone) hyperedge. Eq. (1)."""
return _log2(n_orders) + _log2_binom(n_nodes, size)
def _child_codelength(size_c, size_p, overlap, n_orders, n_nodes):
"""H(c|p): bits to transmit a child hyperedge from its parent. Eq. (3)."""
return (
_log2(n_orders)
+ _log2(min(size_p, size_c))
+ _log2_binom(size_p, overlap)
+ _log2_binom(n_nodes - size_p, size_c - overlap)
)
def _reduced_mutual_information(size_p, size_c, overlap, n_nodes):
"""R(c, p): reduced mutual information between hyperedges. Eq. (11).
Symmetric in the two hyperedge sizes. Larger overlap / nestedness gives a
larger value, so the description length favors assigning highly overlapping
hyperedges as children of a shared parent.
"""
p_not_c = size_p - overlap
c_not_p = size_c - overlap
rest = n_nodes - size_p - c_not_p # = n_nodes - |p ∪ c|
log2_multinomial = (
math.lgamma(n_nodes + 1.0)
- math.lgamma(overlap + 1.0)
- math.lgamma(p_not_c + 1.0)
- math.lgamma(c_not_p + 1.0)
- math.lgamma(rest + 1.0)
) / _LN2
return (
_log2_binom(n_nodes, size_p)
+ _log2_binom(n_nodes, size_c)
- log2_multinomial
- _log2(min(size_p, size_c))
)
# ---------------------------------------------------------------------------
# Weight model (empirical-Bayes Poisson / Geometric prior). Sec. III.
# ---------------------------------------------------------------------------
def _expected_weight(role, gamma, mean_weight):
"""Expected weight mu_b for a parent (role=1) or child (role=0). Eq. (15)."""
return 1.0 + 2.0 * (gamma ** (1 - role)) * (mean_weight - 1.0) / (1.0 + gamma)
def _weight_codelength(weight, role, gamma, mean_weight, prior):
"""L(w, b): bits to transmit a hyperedge weight given its role. Eqs (16)-(17)."""
mu = _expected_weight(role, gamma, mean_weight)
if prior == "poisson":
if mu - 1.0 <= 0.0:
return 0.0
return (
(mu - 1.0) / _LN2
- (weight - 1.0) * _log2(mu - 1.0)
+ _log2_factorial(weight - 1.0)
)
if prior == "geometric":
if mu <= 1.0:
return 0.0
return _log2(mu) + (weight - 1.0) * _log2(mu / (mu - 1.0))
raise ValueError(f"prior must be 'poisson' or 'geometric', got {prior!r}")
def _weight_term(weight, gamma, mean_weight, prior):
"""L(w, 1) - L(w, 0): weight-dependent reward for keeping an edge as a parent.
Negative for high-weight hyperedges when ``gamma < 1``, which discourages
demoting them to children (i.e. encourages keeping them in the backbone).
Zero when ``gamma == 1`` or the hypergraph is effectively unweighted.
"""
return _weight_codelength(
weight, 1, gamma, mean_weight, prior
) - _weight_codelength(weight, 0, gamma, mean_weight, prior)
# ---------------------------------------------------------------------------
# Input handling
# ---------------------------------------------------------------------------
def _normalize_hyperedges(hyperedges, weights):
"""Validate input and collapse duplicate / empty hyperedges.
Returns ``(edges, weights)`` where ``edges`` is a list of unique
:class:`frozenset` hyperedges and ``weights`` a parallel list of floats.
Duplicate hyperedges are merged; their weights are summed.
"""
raw = [frozenset(e) for e in hyperedges]
if weights is not None:
weights = [float(w) for w in weights]
if len(weights) != len(raw):
raise ValueError("weights must have the same length as hyperedges")
if any(w < 1.0 for w in weights):
raise ValueError("weights must be >= 1")
edges = []
out_weights = []
index = {}
for i, fs in enumerate(raw):
if len(fs) == 0:
continue
w = weights[i] if weights is not None else 1.0
if fs in index:
if weights is not None:
out_weights[index[fs]] += w
else:
index[fs] = len(edges)
edges.append(fs)
out_weights.append(w)
return edges, out_weights
def _sorted_labels(labels):
"""Sort node labels deterministically, tolerating mixed / unorderable types.
Uses the labels' natural order when they are mutually comparable (so
integer labels sort numerically) and falls back to ``repr`` otherwise.
"""
labels = list(labels)
try:
return sorted(labels)
except TypeError:
return sorted(labels, key=repr)
def _sorted_by(items, key):
"""``sorted(items, key=key)``, retrying on ``repr``-ed labels on TypeError.
*key* is called as ``key(item, sort_labels)`` where ``sort_labels`` is a
function returning a canonical sorted list of an iterable of node labels.
"""
items = list(items)
try:
return sorted(items, key=lambda it: key(it, sorted))
except TypeError:
return sorted(items, key=lambda it: key(it, lambda xs: sorted(map(repr, xs))))
# ---------------------------------------------------------------------------
# Result container
# ---------------------------------------------------------------------------
[docs]
@dataclass
class HypergraphBackbone:
"""Result of :func:`mdl_hypergraph_backbone` or its local variant.
Attributes
----------
backbone : list of frozenset
The retained ("parent") hyperedges forming the structural backbone.
assignment : dict
Mapping from each parent hyperedge to the list of child hyperedges it
encodes (its star in the intersection graph). Parents with no children
are absent from this mapping. Empty for a local backbone, which has no
single global star partition (see *local_backbones*).
n_nodes : int
Number of distinct nodes in the input hypergraph.
n_orders : int
Number of distinct hyperedge sizes (orders) in the input.
weighted : bool
Whether edge weights influenced the backbone.
gamma : float
Weight/topology trade-off used (only meaningful when ``weighted``).
prior : str
Weight prior used (``"poisson"`` or ``"geometric"``).
description_length : float
Description length ``L(G, B*)`` of the input given the backbone (bits).
baseline_description_length : float
Description length ``L(G, G)`` with no backbone (bits).
compression_ratio : float
Inverse compression ratio ``eta = L(G, B*) / L(G, G)`` in ``[0, 1]``.
Smaller means more redundant structure was removed. For a local
backbone this is ``eta^(local)`` (Eq. (E4)).
n_hyperedges : int or None
Number of distinct (deduplicated, non-empty) input hyperedges, or
``None`` if unknown (then derived from *backbone* and *assignment*).
aggregation : str or None
For a local backbone, the rule used to combine the neighborhood
backbones (``"union"``, ``"majority"`` or ``"all"``); ``None`` for the
global backbone.
local_backbones : dict
For a local backbone, a mapping from each node ``i`` to its
neighborhood backbone ``B_i*`` (a list of frozenset), the MDL optimum
of the hyperedges containing ``i`` (Eq. (E1)). Empty otherwise.
local_description_lengths : dict
For a local backbone, a mapping from each node ``i`` to the
neighborhood optimum ``L(G_i, B_i*)`` in bits. Empty otherwise.
"""
backbone: list[Hyperedge] = field(default_factory=list)
assignment: dict[Hyperedge, list[Hyperedge]] = field(default_factory=dict)
n_nodes: int = 0
n_orders: int = 0
weighted: bool = False
gamma: float = 1.0
prior: str = "poisson"
description_length: float = 0.0
baseline_description_length: float = 0.0
compression_ratio: float = 1.0
n_hyperedges: int | None = None
aggregation: str | None = None
local_backbones: dict[Hashable, list[Hyperedge]] = field(default_factory=dict)
local_description_lengths: dict[Hashable, float] = field(default_factory=dict)
def __len__(self) -> int:
return len(self.backbone)
@property
def n_input_hyperedges(self) -> int:
"""Total distinct input hyperedges (backbone parents + all children)."""
if self.n_hyperedges is not None:
return self.n_hyperedges
return len(self.backbone) + sum(len(v) for v in self.assignment.values())
@property
def fraction_kept(self) -> float:
"""Fraction of (unique) input hyperedges retained in the backbone."""
total = self.n_input_hyperedges
return len(self.backbone) / total if total else 1.0
[docs]
@dataclass
class ValidatedHypergraph:
"""Result of :func:`statistically_validated_hypergraph` / cores.
Attributes
----------
validated : list of frozenset
The statistically validated hyperedges (SVH) or cores (SVC) -- the
sub-hypergraph that survived FDR validation.
pvalues : dict
Mapping from each tested hyperedge/group to its p-value.
counts : dict
Mapping from each tested hyperedge/group to its observed co-occurrence
count (multiplicity for SVH; number of containing instances for SVC).
alpha : float
Significance level used for FDR validation.
method : str
``"svh"`` or ``"svc"``.
n_nodes : int
Number of distinct nodes in the input hypergraph.
n_instances : int
Total number of hyperedge instances (sum of multiplicities).
"""
validated: list[Hyperedge] = field(default_factory=list)
pvalues: dict[Hyperedge, float] = field(default_factory=dict)
counts: dict[Hyperedge, int] = field(default_factory=dict)
alpha: float = 0.01
method: str = "svh"
n_nodes: int = 0
n_instances: int = 0
def __len__(self) -> int:
return len(self.validated)
def __iter__(self) -> Iterator[Hyperedge]:
return iter(self.validated)
# ---------------------------------------------------------------------------
# Public API
# ---------------------------------------------------------------------------
[docs]
def intersection_graph(
hyperedges: Hyperedges, s: int = 1, weight: str = "overlap"
) -> nx.Graph:
"""Build the intersection graph (or s-line graph) of a hypergraph.
Each hyperedge becomes a node (labeled by its integer index after
deduplication); two hyperedges are linked when they share at least *s*
nodes. With ``s=1`` this is the intersection graph over which MDL
parent--child relationships are formed (Kirkley et al. 2026, Appendix D);
with ``s > 1`` it is the *s-line graph*, the basis of s-connectivity in
higher-order networks.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Each hyperedge is an iterable of node labels.
s : int, optional (default=1)
Minimum shared-node count ``|e_i ∩ e_j|`` for two hyperedges to be
linked. Must be ``>= 1``.
weight : string, optional (default="overlap")
Edge attribute name used to store the overlap size ``|e_i ∩ e_j|``.
Returns
-------
I : networkx.Graph
Graph on hyperedge indices. Each node has a ``"members"`` attribute
(the hyperedge as a :class:`frozenset`); each edge stores the overlap
size under *weight*.
Raises
------
ValueError
If *s* is less than 1.
See Also
--------
s_components : Connected components of the s-line graph.
hypergraph_to_bipartite : Node--hyperedge incidence graph.
Examples
--------
>>> from nx_backbone import intersection_graph
>>> I = intersection_graph([(1, 2, 3), (2, 3, 4), (5, 6)])
>>> I.number_of_nodes()
3
>>> I[0][1]["overlap"]
2
>>> intersection_graph([(1, 2, 3), (2, 3, 4), (5, 6)], s=3).number_of_edges()
0
"""
if s < 1:
raise ValueError(f"s must be >= 1, got {s}")
edges, _ = _normalize_hyperedges(hyperedges, None)
graph: nx.Graph = nx.Graph()
for i, e in enumerate(edges):
graph.add_node(i, members=e)
node_to_edges: dict[Hashable, list[int]] = {}
for i, e in enumerate(edges):
for v in e:
node_to_edges.setdefault(v, []).append(i)
overlaps: dict[tuple[int, int], int] = {}
for incident in node_to_edges.values():
m = len(incident)
if m < 2:
continue
for a in range(m):
for b in range(a + 1, m):
i, j = incident[a], incident[b]
key = (i, j) if i < j else (j, i)
overlaps[key] = overlaps.get(key, 0) + 1
for (i, j), o in overlaps.items():
if o >= s:
graph.add_edge(i, j, **{weight: o})
return graph
def _greedy_edge(n_edges, edges, neighbors, wterm):
"""Greedy "edge" optimizer: accept parent->child links by decreasing gain.
Builds a star partition of the intersection graph by greedily accepting the
highest-gain parent--child assignment that preserves the star structure
(Kirkley et al. 2026, Appendix D). Returns ``parent_of`` (child index ->
parent index).
"""
candidates = []
for i in range(n_edges):
si = len(edges[i])
for j, rmi in neighbors[i]:
gain = rmi + wterm[i] # i becomes a child of parent j
if gain > 0.0:
candidates.append((gain, len(edges[j]), si, i, j))
# Highest gain first; prefer the larger hyperedge as parent on ties.
candidates.sort(key=lambda t: (-t[0], -t[1], t[2], t[3], t[4]))
UNDECIDED, PARENT, CHILD = 0, 1, 2
role = [UNDECIDED] * n_edges
parent_of = {}
for _gain, _psize, _csize, c, p in candidates:
if role[c] != UNDECIDED or role[p] == CHILD:
continue
role[c] = CHILD
parent_of[c] = p
if role[p] == UNDECIDED:
role[p] = PARENT
return parent_of
def _greedy_node(n_edges, neighbors, wterm):
"""Greedy "node" optimizer: add backbone parents by decreasing total gain.
The alternative scheme of Kirkley et al. 2026 (Appendix D): repeatedly add to
the backbone the hyperedge whose promotion to a parent most increases the
total parent--child savings (a facility-location-style greedy), then force any
still-uncovered hyperedge to be a parent. Returns ``parent_of``.
The marginal gain of promoting a hyperedge never increases as the backbone
grows (each child's best saving only rises, and backbone members drop out of
the sum), so the argmax is found with a *lazy* max-heap of stale gains
(CELF): only the top entry is re-evaluated each step. Ties are broken by
the smallest hyperedge index, exactly as an exhaustive scan would, so the
result is identical to recomputing every gain each round.
"""
import heapq
in_backbone = [False] * n_edges
best_saving = [0.0] * n_edges
best_parent = [None] * n_edges
def _attach_children(parent):
for c, rmi in neighbors[parent]:
if in_backbone[c]:
continue
gain = rmi + wterm[c]
if gain > best_saving[c]:
best_saving[c] = gain
best_parent[c] = parent
def _gain(e):
gain = -best_saving[e]
for c, rmi in neighbors[e]:
if in_backbone[c]:
continue
delta = (rmi + wterm[c]) - best_saving[c]
if delta > 0.0:
gain += delta
return gain
min_gain = 1e-12 # require a strictly positive improvement
heap = [(-_gain(e), e) for e in range(n_edges)]
heapq.heapify(heap)
while heap:
neg_stale, e = heap[0]
if -neg_stale <= min_gain:
break # every true gain is bounded by its stale value
fresh = _gain(e)
if fresh != -neg_stale:
heapq.heapreplace(heap, (-fresh, e))
continue
heapq.heappop(heap)
in_backbone[e] = True
best_saving[e] = 0.0
best_parent[e] = None
_attach_children(e)
# Force-cover any hyperedge still without a parent (no overlap with backbone).
for e in range(n_edges):
if not in_backbone[e] and best_parent[e] is None:
in_backbone[e] = True
_attach_children(e)
return {e: best_parent[e] for e in range(n_edges) if not in_backbone[e]}
class _MDLFit:
"""The MDL objective of one hypergraph, ready to optimize and evaluate.
Holds everything the description length depends on -- the node count
``N``, the number of distinct orders ``L``, the empirical mean weight and
the reduced mutual information of every overlapping pair -- for the
deduplicated hyperedge list *edges*. Used both for the whole hypergraph
(:func:`mdl_hypergraph_backbone`) and for each node neighborhood ``G_i``
considered in isolation (:func:`local_mdl_hypergraph_backbone`).
"""
def __init__(self, edges, weight_list, weights_given, gamma, prior):
self.edges = edges
self.weight_list = weight_list
self.gamma = gamma
self.prior = prior
n_edges = len(edges)
self.n_edges = n_edges
all_nodes = set()
sizes = set()
for e in edges:
all_nodes.update(e)
sizes.add(len(e))
self.n_nodes = len(all_nodes)
self.n_orders = len(sizes)
self.mean_weight = sum(weight_list) / n_edges
self.weighted = weights_given and any(w != 1.0 for w in weight_list)
# Pairwise overlaps via node neighborhoods (Appendix D).
node_to_edges = {}
for i, e in enumerate(edges):
for v in e:
node_to_edges.setdefault(v, []).append(i)
overlaps = {}
for incident in node_to_edges.values():
m = len(incident)
if m < 2:
continue
for a in range(m):
for b in range(a + 1, m):
i, j = incident[a], incident[b]
key = (i, j) if i < j else (j, i)
overlaps[key] = overlaps.get(key, 0) + 1
# Per-edge weight terms (0 when unweighted).
if self.weighted:
self.wterm = [
_weight_term(w, gamma, self.mean_weight, prior) for w in weight_list
]
else:
self.wterm = [0.0] * n_edges
# Intersection-graph adjacency with reduced mutual information per pair.
neighbors = [[] for _ in range(n_edges)]
for (i, j), o in overlaps.items():
rmi = _reduced_mutual_information(
len(edges[i]), len(edges[j]), o, self.n_nodes
)
neighbors[i].append((j, rmi))
neighbors[j].append((i, rmi))
self.neighbors = neighbors
def optimize(self, method):
"""Return ``parent_of`` (child index -> parent index) from a greedy."""
if method == "edge":
return _greedy_edge(self.n_edges, self.edges, self.neighbors, self.wterm)
if method == "node":
return _greedy_node(self.n_edges, self.neighbors, self.wterm)
# "auto": run both and keep the lower description length.
return min(
(
_greedy_edge(self.n_edges, self.edges, self.neighbors, self.wterm),
_greedy_node(self.n_edges, self.neighbors, self.wterm),
),
key=self.description_length,
)
def description_length(self, parent_of):
"""``L(G, B)`` in bits for a parent/child assignment (Eqs (5), (13))."""
edges = self.edges
dl = 0.0
for i, e in enumerate(edges):
role = 0 if i in parent_of else 1
if role:
dl += _parent_codelength(len(e), self.n_orders, self.n_nodes)
else:
p = edges[parent_of[i]]
dl += _child_codelength(
len(e), len(p), len(e & p), self.n_orders, self.n_nodes
)
if self.weighted:
dl += _weight_codelength(
self.weight_list[i], role, self.gamma, self.mean_weight, self.prior
)
return dl
def baseline_description_length(self):
"""``L(G, G)`` in bits: every hyperedge transmitted as a parent."""
return self.description_length({})
def assign_children(self, parents):
"""Best ``parent_of`` for a fixed backbone (set of parent indices).
Each non-backbone hyperedge becomes the child of the overlapping
backbone hyperedge with the largest reduced mutual information, which
minimizes its child codelength ``H(c|p)`` (Eq. (12)); ties go to the
smallest index. Returns ``None`` if some hyperedge overlaps no parent.
"""
parent_of = {}
for c in range(self.n_edges):
if c in parents:
continue
best = None
for p, rmi in self.neighbors[c]:
if p in parents and (
best is None or rmi > best[0] or (rmi == best[0] and p < best[1])
):
best = (rmi, p)
if best is None:
return None
parent_of[c] = best[1]
return parent_of
[docs]
def mdl_hypergraph_backbone(
hyperedges: Hyperedges,
weights: Iterable[float] | None = None,
gamma: float = 1.0,
prior: str = "poisson",
method: str = "auto",
) -> HypergraphBackbone:
"""Extract a hypergraph backbone via minimum description length (MDL).
Implements the parameter-free information-theoretic backboning method of
Kirkley, Felippe, Malizia & Battiston [1]_. The backbone ``B`` is a subset
of hyperedges ("parents") chosen so that the remaining "child" hyperedges can
be cheaply reconstructed from a parent they overlap, exploiting the nested
and redundant structure unique to higher-order networks. The optimal
backbone minimizes the two-part description length
``L(G, B) = L(B) + L(G|B)`` (Eqs (1)-(6)); equivalently it maximizes the
total parent--child reduced mutual information (Eq. (12)).
The method is **fully nonparametric** for unweighted hypergraphs. Edge
weights are incorporated through an empirical-Bayes prior with a single knob
*gamma* trading off weight against topology (Sec. III).
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Each hyperedge is an iterable of node labels (treated as
a set; repeated nodes within a hyperedge are ignored). Duplicate
hyperedges are merged.
weights : iterable of numbers or None, optional (default=None)
Optional per-hyperedge weights (each ``>= 1``), parallel to
*hyperedges*. ``None`` means an unweighted hypergraph.
gamma : float, optional (default=1.0)
Weight/topology trade-off in ``(0, 1]`` (Eq. (14)). ``gamma=1`` ignores
weights (recovers the unweighted backbone); ``gamma -> 0`` makes it
increasingly costly to leave a high-weight hyperedge out of the backbone.
Ignored when *weights* is ``None``.
prior : {"poisson", "geometric"}, optional (default="poisson")
Weight prior family (Eqs (16)-(17)). Ignored when *weights* is ``None``.
method : {"auto", "edge", "node"}, optional (default="auto")
Greedy optimizer (Kirkley et al. 2026, Appendix D). ``"edge"`` adds
parent--child links by decreasing gain; ``"node"`` adds backbone parents
by decreasing total gain; ``"auto"`` runs both and keeps the lower
description length (the procedure used in [1]_).
Returns
-------
result : HypergraphBackbone
The backbone hyperedges, the parent--child assignment, and compression
diagnostics (see :class:`HypergraphBackbone`).
Raises
------
ValueError
If *gamma* is not in ``(0, 1]``, *prior* is unknown, *weights* has the
wrong length or contains a value ``< 1``, or *method* is unsupported.
See Also
--------
hypergraph_compression_ratio : Only the compression ratio ``eta``.
maximal_hyperedges : Purely set-inclusion (toplex) pruning.
statistically_validated_hypergraph : Statistical hyperedge filtering.
Notes
-----
Exact minimization is combinatorial; this uses the greedy heuristics of
[1]_, which forms a maximum-reward partition of the intersection graph into
disjoint stars (each child attached to a single parent). On small inputs the
greedy compression matches exhaustive search.
References
----------
.. [1] Kirkley, A., Felippe, H., Malizia, F., & Battiston, F. (2026).
Hypergraph backboning. arXiv:2606.00893.
Examples
--------
>>> from nx_backbone import mdl_hypergraph_backbone
>>> # A 4-node hyperedge with two nested (redundant) sub-hyperedges.
>>> G = [(1, 2, 3, 4), (1, 2, 3), (2, 3, 4), (8, 9)]
>>> result = mdl_hypergraph_backbone(G)
>>> frozenset({1, 2, 3, 4}) in result.backbone
True
>>> result.compression_ratio <= 1.0
True
"""
if not 0.0 < gamma <= 1.0:
raise ValueError(f"gamma must be in (0, 1], got {gamma}")
if prior not in ("poisson", "geometric"):
raise ValueError(f"prior must be 'poisson' or 'geometric', got {prior!r}")
if method not in ("edge", "node", "auto"):
raise ValueError(f"method must be 'edge', 'node', or 'auto', got {method!r}")
edges, weight_list = _normalize_hyperedges(hyperedges, weights)
n_edges = len(edges)
result = HypergraphBackbone(gamma=gamma, prior=prior)
if n_edges == 0:
return result
fit = _MDLFit(edges, weight_list, weights is not None, gamma, prior)
parent_of = fit.optimize(method)
dl = fit.description_length(parent_of)
dl0 = fit.baseline_description_length()
result.n_nodes = fit.n_nodes
result.n_orders = fit.n_orders
result.weighted = fit.weighted
result.n_hyperedges = n_edges
result.backbone = [edges[i] for i in range(n_edges) if i not in parent_of]
assignment: dict[Hyperedge, list[Hyperedge]] = {}
for c, p in parent_of.items():
assignment.setdefault(edges[p], []).append(edges[c])
result.assignment = assignment
result.description_length = dl
result.baseline_description_length = dl0
result.compression_ratio = dl / dl0 if dl0 > 0 else 1.0
return result
[docs]
def hypergraph_compression_ratio(
hyperedges: Hyperedges,
weights: Iterable[float] | None = None,
gamma: float = 1.0,
prior: str = "poisson",
) -> float:
"""Inverse compression ratio ``eta`` of the MDL backbone (Eq. (8)).
A convenience wrapper returning only
:attr:`HypergraphBackbone.compression_ratio` from
:func:`mdl_hypergraph_backbone`. ``eta`` lies in ``[0, 1]``; values near 0
indicate a highly redundant (compressible) hypergraph, while ``eta = 1``
indicates no compressible structure.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph.
weights : iterable of numbers or None, optional (default=None)
Optional per-hyperedge weights.
gamma : float, optional (default=1.0)
Weight/topology trade-off (see :func:`mdl_hypergraph_backbone`).
prior : {"poisson", "geometric"}, optional (default="poisson")
Weight prior family.
Returns
-------
eta : float
Inverse compression ratio in ``[0, 1]``.
Raises
------
ValueError
For invalid *weights*, *gamma*, or *prior* (see
:func:`mdl_hypergraph_backbone`).
See Also
--------
mdl_hypergraph_backbone : The full backbone and diagnostics.
Examples
--------
>>> from nx_backbone import hypergraph_compression_ratio
>>> eta = hypergraph_compression_ratio([(1, 2, 3, 4), (1, 2, 3), (2, 3, 4)])
>>> 0.0 <= eta <= 1.0
True
"""
return mdl_hypergraph_backbone(
hyperedges, weights=weights, gamma=gamma, prior=prior
).compression_ratio
_AGGREGATIONS = ("union", "majority", "all")
[docs]
def local_mdl_hypergraph_backbone(
hyperedges: Hyperedges,
weights: Iterable[float] | None = None,
gamma: float = 1.0,
prior: str = "poisson",
aggregation: str = "majority",
method: str = "auto",
) -> HypergraphBackbone:
"""Extract the local (per-node neighborhood) MDL hypergraph backbone.
The local variant of :func:`mdl_hypergraph_backbone` from Appendix E of
Kirkley, Felippe, Malizia & Battiston [1]_, the higher-order analogue of
the disparity filter. Each node ``i`` has a *neighborhood*
``G_i = {e in G : i in e}``, the hyperedges containing it. The MDL
objective is minimized on every neighborhood separately, "considered in
isolation" (Eq. (E1)),
.. math:: B_i^* = \\arg\\min_{B_i \\subseteq G_i} \\mathcal{L}(G_i, B_i),
and the neighborhood backbones ``B_i*`` are combined into one backbone
``B^(local)`` by an *aggregation* rule. Because each hyperedge is judged
only against the other hyperedges of its own nodes, the backbone keeps
hyperedges that are informative locally even when the global backbone
encodes them as children of a hyperedge elsewhere; local backbones are
therefore denser and compress less than the global one [1]_.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Each hyperedge is an iterable of node labels (treated
as a set). Duplicate hyperedges are merged.
weights : iterable of numbers or None, optional (default=None)
Optional per-hyperedge weights (each ``>= 1``), parallel to
*hyperedges*. As in [1]_, the weighted objective ``L_w`` (Eq. (13))
then replaces ``L`` in every neighborhood.
gamma : float, optional (default=1.0)
Weight/topology trade-off in ``(0, 1]``; see
:func:`mdl_hypergraph_backbone`. Ignored when *weights* is ``None``.
prior : {"poisson", "geometric"}, optional (default="poisson")
Weight prior family. Ignored when *weights* is ``None``.
aggregation : {"majority", "union", "all"}, optional (default="majority")
How the neighborhood backbones are combined. A hyperedge ``e`` is
kept if it belongs to ``B_i*`` for
- ``"union"``: at least one of its nodes (Eq. (E2)); every node then
stays in the backbone, but the backbone tends to be dense;
- ``"majority"``: at least half of its nodes (the rule used for
Table II of [1]_; see Notes);
- ``"all"``: every one of its nodes (the most conservative rule,
often very sparse).
method : {"auto", "edge", "node"}, optional (default="auto")
Greedy optimizer applied to each neighborhood; see
:func:`mdl_hypergraph_backbone`.
Returns
-------
result : HypergraphBackbone
The local backbone ``B^(local)`` (in input order) with
``description_length`` = Eq. (E3), ``baseline_description_length`` =
``sum_i L(G_i, G_i)`` and ``compression_ratio`` = ``eta^(local)``
(Eq. (E4)). Per-node results are in ``local_backbones`` (each
``B_i*``) and ``local_description_lengths`` (each
``L(G_i, B_i*)``); ``aggregation`` records the rule. ``assignment``
is empty, since there is no single global parent--child partition.
Raises
------
ValueError
If *aggregation* is unknown, or for any invalid argument of
:func:`mdl_hypergraph_backbone`.
See Also
--------
mdl_hypergraph_backbone : The global MDL backbone.
disparity_filter : The dyadic local backbone this generalises.
Notes
-----
Every hyperedge of ``G_i`` contains ``i``, so all pairs of hyperedges in a
neighborhood overlap and any non-empty subset of ``G_i`` is a valid
backbone of it. The neighborhood objectives involve ``O(sum_i |G_i|^2)``
hyperedge pairs, as many as the global intersection graph is built from.
Appendix E of [1]_ leaves a few details implicit; this implementation
makes the following choices.
- "Considered in isolation" is read literally: the node count ``N``, the
number of distinct orders ``L`` and, for weighted input, the mean weight
of the empirical-Bayes prior are those of ``G_i``, not of ``G``.
- *Majority* means at least half of a hyperedge's nodes. The paper notes
that, with the union rule, the majority rule is "consistent with the
disparity filter convention for pairwise edges in which the edge is kept
[...] if it belongs to either node's neighborhood backbone"; for a
two-node hyperedge only the at-least-half reading does that.
- In Eq. (E3), ``L(G_i, B_i* ∩ B^(local))`` assigns each remaining
hyperedge of ``G_i`` to the parent in ``B_i* ∩ B^(local)`` that encodes
it most cheaply. If that intersection is empty (possible only under
``"majority"`` or ``"all"``), the neighborhood is charged its
uncompressed cost ``L(G_i, G_i)``.
- The denominator ``L^(local)(G, G)`` of Eq. (E4) is taken as
``sum_i L(G_i, G_i)``, the cost with no compression in any
neighborhood (read literally, ``B_i* ∩ G = B_i*`` would make
``eta^(local)`` identically 1 under the union rule). With the union
rule ``eta^(local) <= 1``; with ``"majority"`` or ``"all"`` it can
exceed 1 when dropping a neighborhood's parents forces costlier
encodings of its remaining hyperedges.
References
----------
.. [1] Kirkley, A., Felippe, H., Malizia, F., & Battiston, F. (2026).
Hypergraph backboning. arXiv:2606.00893.
Examples
--------
A core hyperedge ``{0, ..., 5}`` with two nested sub-hyperedges and two
"satellites", each joining four core nodes to three peripheral nodes of
its own. The global backbone encodes the satellite ``{0, 1, 2, 3, 6, 7,
8}`` as a child, but it is the only hyperedge in the neighborhoods of
nodes 6, 7 and 8, so the local backbone keeps it (while still pruning the
nested sub-hyperedges):
>>> from nx_backbone import (
... local_mdl_hypergraph_backbone,
... mdl_hypergraph_backbone,
... )
>>> G = [
... (0, 1, 2, 3, 4, 5),
... (0, 1, 2, 3, 4),
... (0, 1, 2, 3, 5),
... (0, 1, 2, 3, 6, 7, 8),
... (2, 3, 4, 5, 9, 10, 11),
... ]
>>> satellite = frozenset({0, 1, 2, 3, 6, 7, 8})
>>> satellite in mdl_hypergraph_backbone(G).backbone
False
>>> local = local_mdl_hypergraph_backbone(G) # majority rule
>>> satellite in local.backbone
True
>>> len(local.backbone)
3
>>> local.local_backbones[6] == [satellite]
True
"""
if aggregation not in _AGGREGATIONS:
raise ValueError(
f"aggregation must be 'union', 'majority', or 'all', got {aggregation!r}"
)
if not 0.0 < gamma <= 1.0:
raise ValueError(f"gamma must be in (0, 1], got {gamma}")
if prior not in ("poisson", "geometric"):
raise ValueError(f"prior must be 'poisson' or 'geometric', got {prior!r}")
if method not in ("edge", "node", "auto"):
raise ValueError(f"method must be 'edge', 'node', or 'auto', got {method!r}")
edges, weight_list = _normalize_hyperedges(hyperedges, weights)
n_edges = len(edges)
result = HypergraphBackbone(gamma=gamma, prior=prior, aggregation=aggregation)
if n_edges == 0:
return result
node_to_edges: dict[Hashable, list[int]] = {}
sizes = set()
for k, e in enumerate(edges):
sizes.add(len(e))
for v in e:
node_to_edges.setdefault(v, []).append(k)
result.n_nodes = len(node_to_edges)
result.n_orders = len(sizes)
result.n_hyperedges = n_edges
# Eq. (E1): the MDL backbone of every neighborhood, in isolation.
votes = [0] * n_edges
neighborhoods = []
for v in _sorted_labels(node_to_edges):
members = node_to_edges[v] # global indices, increasing
fit = _MDLFit(
[edges[k] for k in members],
[weight_list[k] for k in members],
weights is not None,
gamma,
prior,
)
result.weighted = result.weighted or fit.weighted
parents = {j for j in range(fit.n_edges) if j not in fit.optimize(method)}
for j in parents:
votes[members[j]] += 1
result.local_backbones[v] = [edges[members[j]] for j in sorted(parents)]
result.local_description_lengths[v] = fit.description_length(
fit.assign_children(parents)
)
neighborhoods.append((members, fit, parents))
# Aggregate the neighborhood backbones (Eq. (E2) and the paper's
# majority and all-nodes rules).
if aggregation == "union":
keep = [votes[k] >= 1 for k in range(n_edges)]
elif aggregation == "majority":
keep = [2 * votes[k] >= len(edges[k]) for k in range(n_edges)]
else:
keep = [votes[k] == len(edges[k]) for k in range(n_edges)]
result.backbone = [edges[k] for k in range(n_edges) if keep[k]]
# Eqs (E3)-(E4): neighborhood codelengths restricted to B^(local).
dl = 0.0
dl0 = 0.0
for members, fit, parents in neighborhoods:
baseline = fit.baseline_description_length()
dl0 += baseline
kept = {j for j in parents if keep[members[j]]}
dl += fit.description_length(fit.assign_children(kept)) if kept else baseline
result.description_length = dl
result.baseline_description_length = dl0
result.compression_ratio = dl / dl0 if dl0 > 0 else 1.0
return result
[docs]
def maximal_hyperedges(hyperedges: Hyperedges) -> list[Hyperedge]:
"""Inclusion (toplex) reduction: keep only the maximal hyperedges.
Removes every hyperedge that is a strict subset of another hyperedge,
retaining the *toplexes* -- hyperedges not contained in any other. This is
the simplest structural hypergraph backbone, pruning nested redundancy
purely by set inclusion (cf. ``HyperNetX``'s ``toplexes``). For a richer,
information-theoretic treatment of nested *and* overlapping redundancy, see
:func:`mdl_hypergraph_backbone`.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Duplicate hyperedges are merged.
Returns
-------
maximal : list of frozenset
The maximal hyperedges, ordered by decreasing size (ties broken by
sorted node labels; labels of mixed, unorderable types are compared by
``repr``).
See Also
--------
mdl_hypergraph_backbone : Information-theoretic pruning of nested *and*
overlapping redundancy.
order_filter : Keep hyperedges by size.
Examples
--------
>>> from nx_backbone import maximal_hyperedges
>>> sorted(map(sorted, maximal_hyperedges([(1, 2, 3), (1, 2), (2, 3), (4, 5)])))
[[1, 2, 3], [4, 5]]
"""
edges, _ = _normalize_hyperedges(hyperedges, None)
order = _sorted_by(
range(len(edges)), lambda i, srt: (-len(edges[i]), srt(edges[i]))
)
kept: list[Hyperedge] = []
for i in order:
e = edges[i]
if not any(e < bigger for bigger in kept):
kept.append(e)
return kept
[docs]
def order_filter(
hyperedges: Hyperedges,
min_order: int | None = None,
max_order: int | None = None,
orders: Iterable[int] | None = None,
) -> list[Hyperedge]:
"""Keep hyperedges whose order (size) falls in a range or set.
Order-resolved filtering has no analog in dyadic graphs, where every edge
has order 2; in a hypergraph it selects interactions at chosen scales.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Duplicate hyperedges are merged.
min_order : int or None, optional (default=None)
Keep hyperedges with size ``>= min_order``.
max_order : int or None, optional (default=None)
Keep hyperedges with size ``<= max_order``.
orders : iterable of int or None, optional (default=None)
If given, keep only hyperedges whose size is in this set (applied in
addition to *min_order*/*max_order*).
Returns
-------
selected : list of frozenset
The retained hyperedges, in first-occurrence order.
Raises
------
ValueError
If *min_order* and *max_order* are both given and ``min_order > max_order``.
See Also
--------
maximal_hyperedges : Structural (set-inclusion) pruning.
s_components : Group hyperedges by s-connectivity.
Examples
--------
>>> from nx_backbone import order_filter
>>> sorted(map(sorted, order_filter([(1, 2), (1, 2, 3), (1, 2, 3, 4)], min_order=3)))
[[1, 2, 3], [1, 2, 3, 4]]
"""
if min_order is not None and max_order is not None and min_order > max_order:
raise ValueError("min_order must not exceed max_order")
order_set = set(orders) if orders is not None else None
edges, _ = _normalize_hyperedges(hyperedges, None)
selected = []
for e in edges:
k = len(e)
if min_order is not None and k < min_order:
continue
if max_order is not None and k > max_order:
continue
if order_set is not None and k not in order_set:
continue
selected.append(e)
return selected
[docs]
def s_components(hyperedges: Hyperedges, s: int = 1) -> list[list[Hyperedge]]:
"""Group hyperedges into s-connected components.
Two hyperedges are *s-adjacent* when they share at least *s* nodes; an
s-component is a connected component of the resulting s-line graph
(:func:`intersection_graph` with the same *s*). s-connectivity is a
higher-order notion with no dyadic counterpart and underlies s-centrality
and s-distance analyses of hypergraphs.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Duplicate hyperedges are merged.
s : int, optional (default=1)
Minimum shared-node count for s-adjacency. Must be ``>= 1``.
Returns
-------
components : list of list of frozenset
Each inner list is the hyperedges of one s-connected component,
ordered by decreasing component size (ties broken by sorted node
labels; labels of mixed, unorderable types are compared by ``repr``).
Raises
------
ValueError
If *s* is less than 1.
See Also
--------
intersection_graph : The s-line graph whose components are returned.
maximal_hyperedges : Structural (set-inclusion) pruning.
Examples
--------
>>> from nx_backbone import s_components
>>> comps = s_components([(1, 2, 3), (2, 3, 4), (5, 6, 7)], s=2)
>>> [len(c) for c in comps]
[2, 1]
"""
edges, _ = _normalize_hyperedges(hyperedges, None)
graph = intersection_graph(edges, s=s)
components = [
[edges[i] for i in sorted(component)]
for component in nx.connected_components(graph)
]
return _sorted_by(
components,
lambda comp, srt: (-len(comp), sorted(srt(e) for e in comp)),
)
# ---------------------------------------------------------------------------
# Statistically validated hypergraphs (Musciotto, Battiston & Mantegna 2021)
# ---------------------------------------------------------------------------
def _normalize_multiplicities(hyperedges, weights):
"""Collapse duplicate hyperedges into integer multiplicities.
Unlike :func:`_normalize_hyperedges` (which treats the hypergraph as a set),
this counts repeated occurrences: with ``weights=None`` a hyperedge appearing
``r`` times has multiplicity ``r``; with explicit *weights* the integer
weights of duplicates are summed. Multiplicities are interaction counts for
the statistical filters.
"""
raw = [frozenset(e) for e in hyperedges]
if weights is not None:
weights = list(weights)
if len(weights) != len(raw):
raise ValueError("weights must have the same length as hyperedges")
edges = []
mult = []
index = {}
for i, fs in enumerate(raw):
if len(fs) == 0:
continue
if weights is not None:
w = weights[i]
if not (w >= 1 and float(w).is_integer()):
raise ValueError(
"statistical hypergraph filters require integer multiplicities >= 1"
)
w = int(w)
else:
w = 1
if fs in index:
mult[index[fs]] += w
else:
index[fs] = len(edges)
edges.append(fs)
mult.append(w)
return edges, mult
def _bh_threshold(pvalues, alpha, n_possible):
"""Benjamini-Hochberg FDR threshold with per-rank increment alpha/n_possible.
Returns the largest ``i * alpha / n_possible`` such that the i-th smallest
p-value is below it (0.0 if none), matching Tumminello et al. / HGX.
"""
n = len(pvalues)
if n == 0:
return 0.0
bonf = alpha / n_possible if n_possible > 0 else alpha
threshold = 0.0
for rank, p in enumerate(sorted(pvalues), start=1):
kv = rank * bonf
if p < kv:
threshold = kv
return threshold
def _svh_pvalue(observed, n_instances, degrees, binom):
"""Upper-tail p-value P(X >= observed) with X ~ Binomial(N, prod d_i / N)."""
p = 1.0
for d in degrees:
p *= d / n_instances
return float(binom.sf(observed - 1, n_instances, p))
def _svh_pvalues(observed, n_instances, degree_lists, binom):
"""Vectorised :func:`_svh_pvalue` over parallel lists of counts / degrees."""
import numpy as np
probs = []
for degrees in degree_lists:
p = 1.0
for d in degrees:
p *= d / n_instances
probs.append(p)
sf = binom.sf(np.asarray(observed) - 1, n_instances, np.asarray(probs))
return [float(v) for v in np.atleast_1d(sf)]
def _check_alpha(alpha):
if not 0.0 < alpha <= 1.0:
raise ValueError(f"alpha must be in (0, 1], got {alpha}")
def _count_subgroups(sizes, min_order, top, limit):
"""Number of sub-groups of order ``min_order..top`` over hyperedge *sizes*.
Stops counting (returning a value ``> limit``) as soon as *limit* is passed.
"""
total = 0
for size in sizes:
for k in range(min_order, min(top, size) + 1):
total += math.comb(size, k)
if total > limit:
return total
return total
[docs]
def statistically_validated_hypergraph(
hyperedges: Hyperedges,
weights: Iterable[int] | None = None,
max_order: int | None = None,
alpha: float = 0.01,
) -> ValidatedHypergraph:
"""Extract the Statistically Validated Hypergraph (SVH).
Keeps the observed hyperedges that recur (co-occur) significantly more often
than expected under a null model preserving node activity, following
Musciotto, Battiston & Mantegna [1]_ (the method implemented as ``get_svh``
in Hypergraphx). Each hyperedge of order ``k`` is tested independently per
order: with ``N`` order-``k`` instances and node activities ``d_i`` (number
of order-``k`` instances containing node ``i``), the probability of seeing a
group co-occur at least ``n`` times is ``P(X >= n)`` for
``X ~ Binomial(N, prod_i d_i / N)``. P-values are validated with a
Benjamini-Hochberg FDR at level *alpha* (corrected for the number of
possible order-``k`` hyperedges).
Unlike :func:`mdl_hypergraph_backbone` (an information-theoretic, parameter-
free method), this is a statistical hypothesis test requiring a significance
level, and is most informative for **weighted** hypergraphs whose weights are
integer interaction multiplicities.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Duplicate hyperedges are merged (multiplicities summed).
weights : iterable of int or None, optional (default=None)
Per-hyperedge integer multiplicities (interaction counts). ``None``
treats every hyperedge as occurring once.
max_order : int or None, optional (default=None)
Only test hyperedges up to this order (size). ``None`` tests all orders.
alpha : float, optional (default=0.01)
FDR significance level, in ``(0, 1]``.
Returns
-------
result : ValidatedHypergraph
The validated hyperedges plus per-hyperedge p-values and counts.
Raises
------
ValueError
If *alpha* is not in ``(0, 1]``, or *weights* has the wrong length or
contains a value that is not an integer ``>= 1``.
See Also
--------
statistically_validated_cores : Validates sub-groups (cores) as well.
mdl_hypergraph_backbone : Parameter-free information-theoretic backbone.
References
----------
.. [1] Musciotto, F., Battiston, F., & Mantegna, R. N. (2021). Detecting
informative higher-order interactions in statistically validated
hypergraphs. *Communications Physics*, 4, 218.
Examples
--------
>>> from nx_backbone import statistically_validated_hypergraph
>>> edges = [(1, 2), (1, 2), (1, 2), (1, 3), (2, 4), (5, 6)]
>>> result = statistically_validated_hypergraph(edges)
>>> isinstance(result.validated, list)
True
"""
from scipy.stats import binom
_check_alpha(alpha)
edges, mult = _normalize_multiplicities(hyperedges, weights)
result = ValidatedHypergraph(alpha=alpha, method="svh")
if not edges:
return result
all_nodes = set()
for e in edges:
all_nodes.update(e)
result.n_nodes = len(all_nodes)
result.n_instances = sum(mult)
by_order: dict[int, list[tuple[Hyperedge, int]]] = {}
for e, w in zip(edges, mult):
by_order.setdefault(len(e), []).append((e, w))
for order, members in by_order.items():
if order < 2 or (max_order is not None and order > max_order):
continue
n_instances = sum(w for _, w in members)
degree: dict[Hashable, int] = {}
order_nodes: set[Hashable] = set()
for e, w in members:
order_nodes.update(e)
for node in e:
degree[node] = degree.get(node, 0) + w
pvals = _svh_pvalues(
[w for _, w in members],
n_instances,
[[degree[n] for n in e] for e, _ in members],
binom,
)
n_possible = math.comb(len(order_nodes), order)
threshold = _bh_threshold(pvals, alpha, n_possible)
for (e, w), p in zip(members, pvals):
result.pvalues[e] = p
result.counts[e] = w
if p < threshold:
result.validated.append(e)
return result
_SVC_MAX_GROUPS = 1_000_000
[docs]
def statistically_validated_cores(
hyperedges: Hyperedges,
weights: Iterable[int] | None = None,
min_order: int = 2,
max_order: int | None = None,
alpha: float = 0.01,
max_groups: int | None = _SVC_MAX_GROUPS,
) -> ValidatedHypergraph:
"""Extract the Statistically Validated Cores (SVC).
A complement to :func:`statistically_validated_hypergraph` that validates
significant *groups* (cores) of nodes, including sub-groups that are not
themselves present as a single hyperedge (the ``get_svc`` method of
Hypergraphx, built on [1]_). Orders are processed from high to low; once a
core is validated, its sub-combinations are not re-tested at lower orders, so
significance is attributed to the largest validated group. The co-occurrence
of a group is the number of hyperedge instances (of any order) containing it,
tested against ``Binomial(N, prod_i d_i / N)`` with global node activities,
and validated with the same Benjamini-Hochberg FDR as SVH.
Parameters
----------
hyperedges : iterable of iterables
The hypergraph. Duplicate hyperedges are merged (multiplicities summed).
weights : iterable of int or None, optional (default=None)
Per-hyperedge integer multiplicities. ``None`` treats each as occurring once.
min_order : int, optional (default=2)
Smallest group size to test (``>= 1``).
max_order : int or None, optional (default=None)
Largest group size to test. ``None`` uses the largest hyperedge size.
Every sub-group of every hyperedge up to this size is enumerated -- up to
``2**k`` per hyperedge of size ``k`` -- so set it when the hypergraph
contains large hyperedges.
alpha : float, optional (default=0.01)
FDR significance level, in ``(0, 1]``.
max_groups : int or None, optional (default=1_000_000)
Safety guard: if the number of candidate sub-groups to enumerate (summed
over the distinct hyperedges, for orders ``min_order..max_order``)
exceeds this, a :class:`ValueError` is raised *before* any work is done
instead of running for an impractically long time. ``None`` disables
the guard.
Returns
-------
result : ValidatedHypergraph
The validated cores plus per-group p-values and co-occurrence counts.
Raises
------
ValueError
If *alpha* is not in ``(0, 1]``, *min_order* is less than 1, *weights*
is invalid, or the enumeration would exceed *max_groups* (lower
*max_order*, or raise / disable *max_groups*).
See Also
--------
statistically_validated_hypergraph : Validates only observed hyperedges.
maximal_hyperedges : Purely structural pruning of nested hyperedges.
References
----------
.. [1] Musciotto, F., Battiston, F., & Mantegna, R. N. (2021). Detecting
informative higher-order interactions in statistically validated
hypergraphs. *Communications Physics*, 4, 218.
Examples
--------
>>> from nx_backbone import statistically_validated_cores
>>> edges = [(1, 2, 3), (1, 2, 3), (1, 2, 3), (1, 4), (2, 5)]
>>> result = statistically_validated_cores(edges)
>>> result.method
'svc'
>>> # A single 24-node hyperedge has ~16.7 million sub-groups: bound the order.
>>> big = [tuple(range(24)), (0, 1), (0, 1), (0, 1)]
>>> len(statistically_validated_cores(big, max_order=3).pvalues)
2300
"""
from itertools import combinations
from scipy.stats import binom
_check_alpha(alpha)
if min_order < 1:
raise ValueError(f"min_order must be >= 1, got {min_order}")
edges, mult = _normalize_multiplicities(hyperedges, weights)
result = ValidatedHypergraph(alpha=alpha, method="svc")
if not edges:
return result
n_instances = sum(mult)
degree: dict[Hashable, int] = {}
all_nodes = set()
for e, w in zip(edges, mult):
all_nodes.update(e)
for node in e:
degree[node] = degree.get(node, 0) + w
result.n_nodes = len(all_nodes)
result.n_instances = n_instances
largest = max(len(e) for e in edges)
top = largest if max_order is None else min(max_order, largest)
if max_groups is not None:
n_groups = _count_subgroups((len(e) for e in edges), min_order, top, max_groups)
if n_groups > max_groups:
raise ValueError(
f"statistically_validated_cores would enumerate more than "
f"{max_groups} candidate sub-groups (largest hyperedge has "
f"{largest} nodes, testing orders {min_order}..{top}); pass a "
"smaller max_order, or raise max_groups (None disables this "
"guard)"
)
validated_groups: list[Hyperedge] = []
for order in range(top, min_order - 1, -1):
drop: set[Hyperedge] = set()
for g in validated_groups:
if len(g) > order:
drop.update(frozenset(c) for c in combinations(tuple(g), order))
counts: dict[Hyperedge, int] = {}
for e, w in zip(edges, mult):
if len(e) >= order:
for c in combinations(tuple(e), order):
fs = frozenset(c)
if fs not in drop:
counts[fs] = counts.get(fs, 0) + w
if not counts:
continue
groups = list(counts)
pvals = _svh_pvalues(
[counts[g] for g in groups],
n_instances,
[[degree[n] for n in g] for g in groups],
binom,
)
n_possible = math.comb(len(all_nodes), order)
threshold = _bh_threshold(pvals, alpha, n_possible)
for g, p in zip(groups, pvals):
result.pvalues[g] = p
result.counts[g] = counts[g]
if p < threshold:
result.validated.append(g)
validated_groups.append(g)
return result
# Short aliases
svh = statistically_validated_hypergraph
svc = statistically_validated_cores
_COMPLEXITY = {
"intersection_graph": {
"time": "O(sum_i |G_i|^2)",
"space": "O(m + P)",
"notes": "G_i=hyperedges incident to node i, m=hyperedges, P=overlapping pairs.",
},
"mdl_hypergraph_backbone": {
"time": "O(sum_i |G_i|^2 + P log P)",
"space": "O(m + P)",
"notes": (
"P=overlapping pairs. method='edge' sorts candidates; method='node'/"
"'auto' run a facility-location greedy up to O(m*(m+P))."
),
},
"local_mdl_hypergraph_backbone": {
"time": "O(sum_i k |G_i|^2 log |G_i|)",
"space": "O(sum_i |G_i|^2)",
"notes": (
"G_i=hyperedges containing node i (all pairs overlap), k=mean "
"hyperedge size. method='node'/'auto' run a facility-location "
"greedy per neighborhood, up to O(sum_i |G_i|^3)."
),
},
"hypergraph_compression_ratio": {
"time": "O(sum_i |G_i|^2 + P log P)",
"space": "O(m + P)",
},
"maximal_hyperedges": {
"time": "O(m^2 * k)",
"space": "O(m)",
"notes": "m=hyperedges, k=mean hyperedge size.",
},
"order_filter": {
"time": "O(m)",
"space": "O(m)",
},
"s_components": {
"time": "O(sum_i |G_i|^2)",
"space": "O(m + P)",
"notes": "P=overlapping pairs in the s-line graph.",
},
"statistically_validated_hypergraph": {
"time": "O(sum_e |e| + m log m)",
"space": "O(m + n)",
"notes": "m=hyperedges, n=nodes; per-order binomial tests with FDR.",
},
"statistically_validated_cores": {
"time": "O(sum_e 2^|e|)",
"space": "O(G)",
"notes": "Enumerates sub-groups per order; G=number of distinct sub-groups.",
},
}
append_complexity_docstrings(globals(), _COMPLEXITY)