diff --git a/.github/workflows/deploy-docs.yml b/.github/workflows/deploy-docs.yml index 8921295..6e6cebd 100644 --- a/.github/workflows/deploy-docs.yml +++ b/.github/workflows/deploy-docs.yml @@ -17,16 +17,20 @@ concurrency: group: "pages" cancel-in-progress: false +env: + # Force Node 24 for JS actions still on Node 20 (deprecation 2026-06-16). + FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: "true" + jobs: build: name: Build documentation runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 with: fetch-depth: 0 - - uses: actions/setup-python@v5 + - uses: actions/setup-python@v6 with: python-version: "3.10" diff --git a/.github/workflows/generate-visualizations.yml b/.github/workflows/generate-visualizations.yml index e305f1e..b3cb63f 100644 --- a/.github/workflows/generate-visualizations.yml +++ b/.github/workflows/generate-visualizations.yml @@ -11,15 +11,19 @@ on: - "docs/**/*.rst" - ".github/workflows/generate-visualizations.yml" +env: + # Force Node 24 for JS actions still on Node 20 (deprecation 2026-06-16). + FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: "true" + jobs: generate-gallery: name: Build visual comparison gallery runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - - uses: actions/setup-python@v5 + - uses: actions/setup-python@v6 with: python-version: "3.10" diff --git a/.github/workflows/publish-to-PyPI.yml b/.github/workflows/publish-to-PyPI.yml index 16083fc..0b0e75a 100644 --- a/.github/workflows/publish-to-PyPI.yml +++ b/.github/workflows/publish-to-PyPI.yml @@ -12,14 +12,18 @@ on: permissions: contents: read +env: + # Force Node 24 for JS actions still on Node 20 (deprecation 2026-06-16). + FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: "true" + jobs: run-tests: name: Run tests runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - - uses: actions/setup-python@v5 + - uses: actions/setup-python@v6 with: python-version: "3.10" @@ -34,11 +38,11 @@ jobs: needs: run-tests runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 with: fetch-depth: 0 - - uses: actions/setup-python@v5 + - uses: actions/setup-python@v6 with: python-version: "3.10" diff --git a/.github/workflows/publish-to-conda.yml b/.github/workflows/publish-to-conda.yml index eecb10d..9e8aee8 100644 --- a/.github/workflows/publish-to-conda.yml +++ b/.github/workflows/publish-to-conda.yml @@ -11,6 +11,10 @@ on: types: [published] workflow_dispatch: +env: + # Force Node 24 for JS actions still on Node 20 (deprecation 2026-06-16). + FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: "true" + jobs: build-and-test: name: Build and test conda package @@ -19,7 +23,7 @@ jobs: max-parallel: 5 steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 with: fetch-depth: 0 diff --git a/README.md b/README.md index d61ef9a..2e395da 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,8 @@ Backbone extraction algorithms for complex networks, built on [NetworkX](https://networkx.org/). -This library provides 65 functions across 9 modules for extracting backbone -structures from weighted and unweighted networks. +This library provides 87 functions across 10 modules for extracting backbone +structures from weighted, unweighted, and higher-order (hypergraph) networks. Full documentation: https://www.brianckeegan.com/networkx_backbone/ @@ -36,6 +36,7 @@ pip install -e ".[full]" | **proximity** | Neighborhood-similarity scoring | `jaccard_backbone`, `dice_backbone`, `cosine_backbone`, `hub_promoted_index`, `hub_depressed_index`, `adamic_adar_index`, `resource_allocation_index`, `local_path_index`, and more | | **hybrid** | Combined approaches | `glab_filter` | | **bipartite** | Bipartite projection backbones | `simple_projection`, `hyper_projection`, `probs_projection`, `ycn_projection`, `sdsm`, `fdsm`, `fixedfill`, `fixedrow`, `fixedcol`, `backbone` | +| **hypergraph** | Higher-order (hypergraph) backbones | `mdl_hypergraph_backbone`, `hypergraph_compression_ratio`, `intersection_graph`, `maximal_hyperedges`, `order_filter`, `s_components`, `statistically_validated_hypergraph`, `statistically_validated_cores` | | **unweighted** | Sparsification for unweighted graphs | `sparsify`, `lspar`, `local_degree` | | **filters** | Post-hoc filtering utilities | `multigraph_to_weighted`, `threshold_filter`, `fraction_filter`, `boolean_filter`, `consensus_backbone` | | **measures** | Evaluation and comparison | `node_fraction`, `edge_fraction`, `weight_fraction`, `reachability`, `ks_degree`, `ks_weight`, `compare_backbones` | @@ -50,6 +51,17 @@ Core method families used in [netbone](https://gitlab.liris.cnrs.fr/coregraphie/ - Structural: `global_threshold_filter`, `global_sparsification`, `primary_linkage_analysis`, `edge_betweenness_filter`, `high_salience_skeleton`, `doubly_stochastic_filter`, `maximum_spanning_tree_backbone` - Hybrid: `glab_filter` +Every backbone *model* in Neal's +[Backbone 3.0](https://doi.org/10.1371/journal.pone.0349258) R package is also +covered (`disparity`, `mlf`, `lans`, `sdsm`, `fdsm`, `fixedfill`/`fixedrow`/`fixedcol`, +`bicm`, `fastball`, and the `backbone_from_*` wrappers), including its +hypergraph-projection input via `hypergraph_to_bipartite`, plus its +cross-cutting **features**: multiple-testing correction (`adjust_pvalues`, +`threshold_filter(mtc=...)`), **signed** backbones (`signed=True` adds a `sign` +edge attribute), and **SDSM-EC** edge constraints (`sdsm(prohibited=, required=)`). +See [docs/design/backbone-3.0-coverage.md](docs/design/backbone-3.0-coverage.md) +for the full coverage analysis. + ## Quick Start ```python @@ -70,6 +82,19 @@ print(f"Edges kept: {nb.edge_fraction(G, backbone):.1%}") print(f"Nodes kept: {nb.node_fraction(G, backbone):.1%}") ``` +### Significance options: multiple-testing correction and signed backbones + +```python +# Correct p-values for the number of edges tested (Bonferroni, Holm, BH/FDR, BY) +backbone = nb.threshold_filter(scored, "disparity_pvalue", 0.05, mtc="bh") + +# Signed backbone: keep significantly strong (+1) and significantly weak (-1) +# edges under a two-tailed test; read direction from the "sign" attribute +signed = nb.disparity_filter(G, signed=True) +strong = nb.threshold_filter(signed, "disparity_pvalue", 0.05, mtc="holm") +positives = [(u, v) for u, v, d in strong.edges(data=True) if d["sign"] == 1] +``` + ### Disparity filter visualization ![Disparity filter on Les Miserables](docs/_static/graph_gallery/les_miserables/disparity_filter.png) @@ -96,6 +121,35 @@ backbone = nb.threshold_filter(scored, "sdsm_pvalue", 0.05, mode="below") Projection weights follow the simple/hyper/ProbS/YCN formulations described in [Coscia & Neffke (2017)](https://arxiv.org/abs/1906.09081). +### Hypergraph backbones + +Backbone a hypergraph (a collection of arbitrary-size hyperedges) directly: + +```python +# Parameter-free MDL backbone -- prunes nested/redundant hyperedges +# (Kirkley, Felippe, Malizia & Battiston, 2026) +H = [(1, 2, 3, 4), (1, 2, 3), (2, 3, 4), (8, 9)] +result = nb.mdl_hypergraph_backbone(H) # method="auto" runs both greedy +print(result.backbone) # [frozenset({1, 2, 3, 4}), frozenset({8, 9})] +print(result.compression_ratio) # inverse compression ratio eta +# method="edge" (fastest single pass) or "node" are also available + +# Statistically validated hypergraph (Musciotto, Battiston & Mantegna, 2021) +events = [(1, 2)] * 5 + [(3, 4)] * 100 # repeats = interaction counts +svh = nb.statistically_validated_hypergraph(events, alpha=0.05) +``` + +Interoperate with the higher-order ecosystem (all optional, lazily imported), or +reuse the bipartite projection backbones via the incidence graph: + +```python +B, nodes = nb.hypergraph_to_bipartite(H) # -> NetworkX bipartite graph +scored = nb.sdsm(B, agent_nodes=nodes) # projection backbone of a hypergraph + +nb.write_hif(H, "graph.hif") # HIF interchange (xgi/HNX/HGX/HAT) +edges = nb.from_xgi(xgi_hypergraph) # xgi / hypernetx / hypergraphx / hat +``` + ### Comparing multiple methods ```python @@ -137,13 +191,21 @@ Key papers behind the implemented methods: - Coscia, M. & Neffke, F. M. (2017). [Network backboning with noisy data](https://doi.ieeecomputersociety.org/10.1109/ICDE.2017.100). *Proc. IEEE ICDE*, 425-436. - Coscia, M. & Neffke, F. M. (2017). [Network backboning with noisy data (arXiv:1906.09081)](https://arxiv.org/abs/1906.09081). +- Dianati, N. (2016). [Unwinding the hairball graph: Pruning algorithms for weighted complex networks](https://doi.org/10.1103/PhysRevE.93.012304). *Physical Review E*, 93, 012304. +- Foti, N. J., Hughes, J. M., & Rockmore, D. N. (2011). [Nonparametric sparsification of complex multiscale networks](https://doi.org/10.1371/journal.pone.0016431). *PLoS One*, 6(2), e16431. - Girvan, M., & Newman, M. E. J. (2002). [Community structure in social and biological networks](https://doi.org/10.1073/pnas.122653799). *PNAS*, 99(12), 7821-7826. +- Godard, K., & Neal, Z. P. (2022). [fastball: A fast algorithm to sample bipartite graphs with fixed degree sequences](https://doi.org/10.1093/comnet/cnac049). *J. Complex Networks*, 10(6), cnac049. - Grady, D., Thiemann, C., & Brockmann, D. (2012). [Robust classification of salient links in complex networks](https://doi.org/10.1038/ncomms1847). *Nature Communications*, 3, 864. - Hamann, M., Lindner, G., Meyerhenke, H., Staudt, C. L., and Wagner, D. (2016). [Structure-Preserving Sparsification Methods for Social Networks](https://doi.org/10.1007/s13278-016-0332-2). Social Network Analysis and Mining, 6, 22. - Simas, T., Correia, R. B., & Rocha, L. M. (2021). [The distance backbone of complex networks](https://doi.org/10.1093/comnet/cnab021). *J. Complex Networks*, 9(6), cnab021. - Neal, Z. P. (2014). [The backbone of bipartite projections](https://doi.org/10.1016/j.socnet.2014.06.001). *Social Networks*, 39, 84-97. +- Neal, Z. P., Domagalski, R., & Sagan, B. (2021). [Comparing alternatives to the fixed degree sequence model for extracting the backbone of bipartite projections](https://doi.org/10.1038/s41598-021-03238-3). *Scientific Reports*, 11, 23929. - Neal, Z. P. (2022). [backbone: An R package to extract network backbones](https://doi.org/10.1371/journal.pone.0269137). *PLoS One*, 17(5), e0269137. -- Satuluri, V., Parthasarathy, S., & Ruan, Y. (2011). [Local graph sparsification for scalable clustering](https://doi.org/10.1145/1989323.1989399). *SIGMOD*, 721-732.Serrano, M. A., Boguna, M., & Vespignani, A. (2009). [Extracting the multiscale backbone of complex weighted networks](https://doi.org/10.1073/pnas.0808904106). *PNAS*, 106(16), 6483-6488. +- Neal, Z. P. (2026). [Backbone 3.0: An R package for extracting network backbones](https://doi.org/10.1371/journal.pone.0349258). *PLoS One*, 21, e0349258. +- Neal, Z. P., & Neal, J. W. (2023). [Stochastic Degree Sequence Model with Edge Constraints (SDSM-EC) for Backbone Extraction](https://doi.org/10.1007/978-3-031-53468-3_11). *Complex Networks 12*, 127-136. +- Saracco, F., Di Clemente, R., Gabrielli, A., & Squartini, T. (2015). [Randomizing bipartite networks: the case of the World Trade Web](https://doi.org/10.1038/srep10595). *Scientific Reports*, 5, 10595. +- Satuluri, V., Parthasarathy, S., & Ruan, Y. (2011). [Local graph sparsification for scalable clustering](https://doi.org/10.1145/1989323.1989399). *SIGMOD*, 721-732. +- Serrano, M. A., Boguna, M., & Vespignani, A. (2009). [Extracting the multiscale backbone of complex weighted networks](https://doi.org/10.1073/pnas.0808904106). *PNAS*, 106(16), 6483-6488. - Van Nuffel, N., Heyndrickx, C., & Wets, G. (2010). Measuring hierarchy and reciprocity in networks. - Yassin, A., Haidar, A., Cherifi, H., Seba, H., & Togni, O. (2023). [An evaluation tool for backbone extraction techniques in weighted complex networks](https://doi.org/10.1038/s41598-023-42076-3). *Scientific Reports*, 13, 17000. - Yassin A., Cherifi, H., Seba, H., & Togni, O. (2025). [Backbone extraction through statistical edge filtering: A comparative study](https://doi.org/10.1371/journal.pone.0316141). *PLoS One*, 20(1): e0316141. diff --git a/docs/api/filters.rst b/docs/api/filters.rst index 7d00fed..fefffb4 100644 --- a/docs/api/filters.rst +++ b/docs/api/filters.rst @@ -19,3 +19,5 @@ Complexity classes are provided in each function docstring. .. autofunction:: boolean_filter .. autofunction:: consensus_backbone + +.. autofunction:: adjust_pvalues diff --git a/docs/api/hypergraph.rst b/docs/api/hypergraph.rst new file mode 100644 index 0000000..d71bdea --- /dev/null +++ b/docs/api/hypergraph.rst @@ -0,0 +1,87 @@ +Hypergraph Methods +================== + +Methods that operate directly on hypergraphs (collections of arbitrary-size +hyperedges) rather than on dyadic graphs. Inputs are plain iterables of +hyperedges (each hyperedge an iterable of node labels); the backbone is returned +as a list of :class:`frozenset` hyperedges inside a +:class:`~networkx_backbone.HypergraphBackbone` result. + +The :func:`~networkx_backbone.mdl_hypergraph_backbone` method implements the +parameter-free, information-theoretic (minimum description length) backbone of +Kirkley, Felippe, Malizia & Battiston (2026), which prunes nested and redundant +hyperedges and naturally extends to weighted hypergraphs. + +.. automodule:: networkx_backbone.hypergraph + :no-members: + +.. currentmodule:: networkx_backbone + +.. autofunction:: mdl_hypergraph_backbone + +.. autofunction:: hypergraph_compression_ratio + +.. autofunction:: intersection_graph + +.. autoclass:: HypergraphBackbone + :no-members: + +.. rubric:: Structural methods + +Purely structural hypergraph backbones and utilities with no dyadic analog. +These return plain hyperedge collections rather than a +:class:`~networkx_backbone.HypergraphBackbone`. + +.. autofunction:: maximal_hyperedges + +.. autofunction:: order_filter + +.. autofunction:: s_components + +.. rubric:: Statistical methods + +Hypothesis-testing hypergraph backbones that validate hyperedges/groups +recurring more than expected under a null model (Musciotto, Battiston & +Mantegna, 2021). These require ``scipy`` and return a +:class:`~networkx_backbone.ValidatedHypergraph`. + +.. autofunction:: statistically_validated_hypergraph + +.. autofunction:: statistically_validated_cores + +.. autoclass:: ValidatedHypergraph + :no-members: + +.. rubric:: Interoperability and ingestion + +Convert between the hyperedge-list representation and a NetworkX incidence +bipartite graph (enabling the bipartite projection backbones), the HIF +interchange format, and the ``xgi`` / ``HyperNetX`` / ``HypergraphX`` / HAT +hypergraph classes. The third-party libraries are optional and imported lazily. + +.. automodule:: networkx_backbone.hypergraph_io + :no-members: + +.. currentmodule:: networkx_backbone + +.. autofunction:: hypergraph_to_bipartite + +.. autofunction:: read_hif + +.. autofunction:: write_hif + +.. autofunction:: from_xgi + +.. autofunction:: to_xgi + +.. autofunction:: from_hypernetx + +.. autofunction:: to_hypernetx + +.. autofunction:: from_hypergraphx + +.. autofunction:: to_hypergraphx + +.. autofunction:: from_hat + +.. autofunction:: to_hat diff --git a/docs/api/index.rst b/docs/api/index.rst index 005ca9e..14c86cb 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -1,7 +1,7 @@ API Reference ============= -Complete API documentation for backbone functions across 9 modules. +Complete API documentation for backbone functions across 10 modules. API examples are standardized on ``nx.les_miserables_graph()`` for non-bipartite methods and ``nx.davis_southern_women_graph()`` for bipartite methods. @@ -32,12 +32,15 @@ an aggregate summary in :doc:`../user_guide/complexity`. * - :doc:`bipartite` - 11 - Projection backbones, fixed null models, and high-level wrappers + * - :doc:`hypergraph` + - 12 + - MDL backbone, compression ratio, intersection / s-line graph, inclusion (toplex) reduction, order filter, s-components, and statistically validated hypergraphs / cores * - :doc:`unweighted` - 3 - Sparsification for unweighted graphs (LSpar, local degree) * - :doc:`filters` - - 5 - - Post-hoc filtering utilities and graph-conversion support + - 6 + - Post-hoc filtering utilities, multiple-testing correction, and graph-conversion support * - :doc:`measures` - 7 - Evaluation measures for comparing backbones @@ -54,6 +57,7 @@ an aggregate summary in :doc:`../user_guide/complexity`. proximity hybrid bipartite + hypergraph unweighted filters measures diff --git a/docs/concepts.rst b/docs/concepts.rst index 9002cee..0c65e89 100644 --- a/docs/concepts.rst +++ b/docs/concepts.rst @@ -15,9 +15,10 @@ is a sparser graph that preserves the essential structure of the original. Taxonomy of methods ------------------- -The 65 functions in ``networkx-backbone`` are organized into nine modules based +The 87 functions in ``networkx-backbone`` are organized into ten modules based on the approach they take. The method taxonomy aligns with the categories used -in ``netbone`` (Yassin et al., 2023; https://gitlab.liris.cnrs.fr/coregraphie/netbone). +in ``netbone`` (Yassin et al., 2023; https://gitlab.liris.cnrs.fr/coregraphie/netbone), +extended with a hypergraph module for higher-order networks. Statistical methods ^^^^^^^^^^^^^^^^^^^ @@ -33,6 +34,14 @@ These methods produce a p-value or z-score for each edge. - :func:`~networkx_backbone.lans_filter` -- nonparametric empirical CDF (Foti et al., 2011) - :func:`~networkx_backbone.multiple_linkage_analysis` -- local linkage significance (Van Nuffel et al., 2010; Yassin et al., 2023) +P-values can be corrected for multiple comparisons with +:func:`~networkx_backbone.adjust_pvalues` or the ``mtc`` argument of +:func:`~networkx_backbone.threshold_filter` (Bonferroni, Holm, Hochberg, +Benjamini-Hochberg, Benjamini-Yekutieli). The ``disparity_filter``, +``marginal_likelihood_filter``, and ``lans_filter`` methods also accept +``signed=True`` for a two-tailed test that keeps significantly strong (``+1``) +and significantly weak (``-1``) edges, each tagged with a ``"sign"`` attribute. + Structural methods ^^^^^^^^^^^^^^^^^^ @@ -93,6 +102,44 @@ significant edges from bipartite graph projections: - :func:`~networkx_backbone.backbone_from_projection` / :func:`~networkx_backbone.backbone` -- high-level wrappers +Hypergraph methods +^^^^^^^^^^^^^^^^^^ + +The :mod:`~networkx_backbone.hypergraph` module backbones higher-order networks +(hypergraphs) directly, rather than dyadic graphs. Unlike the projection-based +:mod:`~networkx_backbone.bipartite` methods, which reduce a hypergraph to a +weighted pairwise graph, these methods return a sub-hypergraph (a subset of +hyperedges). + +- :func:`~networkx_backbone.mdl_hypergraph_backbone` -- parameter-free, + information-theoretic (minimum description length) backbone that prunes nested + and redundant hyperedges, with an optional weighted extension (Kirkley, + Felippe, Malizia & Battiston, 2026) +- :func:`~networkx_backbone.hypergraph_compression_ratio` -- inverse compression + ratio achieved by the MDL backbone +- :func:`~networkx_backbone.intersection_graph` -- graph (or s-line graph) + linking hyperedges that share at least *s* nodes +- :func:`~networkx_backbone.maximal_hyperedges` -- inclusion (toplex) reduction, + keeping only hyperedges not contained in another +- :func:`~networkx_backbone.order_filter` -- select hyperedges by order (size) +- :func:`~networkx_backbone.s_components` -- s-connected components of a + hypergraph +- :func:`~networkx_backbone.statistically_validated_hypergraph` and + :func:`~networkx_backbone.statistically_validated_cores` -- statistical + validation of recurring hyperedges/groups under a null model (Musciotto, + Battiston & Mantegna, 2021) + +Because a hypergraph backbone is a subset of hyperedges rather than a graph, this +module returns a :class:`~networkx_backbone.HypergraphBackbone` (or +:class:`~networkx_backbone.ValidatedHypergraph`) result instead of using the +:mod:`~networkx_backbone.filters` utilities. + +The :mod:`~networkx_backbone.hypergraph_io` helpers convert hyperedge lists to and +from a NetworkX incidence bipartite graph (so the +:mod:`~networkx_backbone.bipartite` projection backbones apply to hypergraphs), +the HIF interchange format, and the ``xgi`` / ``HyperNetX`` / ``HypergraphX`` / +HAT hypergraph classes (optional, imported lazily). + Unweighted methods ^^^^^^^^^^^^^^^^^^ diff --git a/docs/design/backbone-3.0-coverage.md b/docs/design/backbone-3.0-coverage.md new file mode 100644 index 0000000..77f71b0 --- /dev/null +++ b/docs/design/backbone-3.0-coverage.md @@ -0,0 +1,173 @@ +# Coverage evaluation: Neal's *Backbone 3.0* (PLOS ONE, 2026) + +- **Status:** Evaluation + gap proposals (no library code changed by this document) +- **Trigger:** Evaluate `networkx-backbone`'s coverage of the methods in + *Neal, Z. P. (2026), "Backbone 3.0: An R package for extracting network + backbones," PLOS ONE* (DOI 10.1371/journal.pone.0349258), propose + implementations for any gaps, and add the paper to the references. + +> Note on sources: the PLOS article and PDF are not reachable from this +> environment's network allowlist. This evaluation is based on the authoritative +> source of the R package itself — `zpneal/backbone` **v3.0.4** (DESCRIPTION, +> NAMESPACE, `NEWS.md`, and the `R/backbone_from_*` sources, fetched from +> `raw.githubusercontent.com`) — cross-checked with the paper's metadata. + +--- + +## 1. Summary + +**Every backbone *model* in Backbone 3.0 is already implemented in +`networkx-backbone`, including its hypergraph-projection capability.** The +package's exported models map one-to-one onto functions here, and this library +goes well beyond Backbone 3.0 (noise-corrected, ECM, multiscale linkage, +structural, proximity, and a full hypergraph module: MDL backbone, SVH/SVC, +interop). + +The differences are **features layered on the statistical models**, not missing +models. Backbone 3.0 adds four cross-cutting options; three are now implemented +in `networkx-backbone` and the fourth is parked: + +| Feature | What it is | Status | +|---------|-----------|--------| +| **`mtc`** multiple-testing correction | Bonferroni / Holm / Hochberg / BH (fdr) / BY adjustment of edge p-values before thresholding | ✅ `adjust_pvalues`, `threshold_filter(mtc=...)` | +| **`signed`** backbones | Two-tailed test retaining significantly *strong* (+) and significantly *weak* (−) edges, with a `sign` attribute | ✅ `signed=` on disparity/mlf/lans + sdsm/fdsm/fixed* | +| **SDSM-EC** edge constraints | `sdsm` with prohibited/required edges (structural 0s/1s; Neal & Neal 2023) | ✅ `sdsm(prohibited=, required=)` | +| **`narrative`** | Auto-generated methods text + citations for a chosen backbone | Parked (low value; utility only) | + +## 2. What Backbone 3.0 provides + +Backbone 3.0 is organized by **input type**, with a model chosen per input +(mirroring this library's own `backbone_from_weighted` / `_from_projection` / +`_from_unweighted` wrappers): + +- **Weighted networks** (`backbone_from_weighted`): models `disparity`, `lans`, + `mlf`, and `global` (with a length-2 parameter giving a *signed* global + threshold). Controlled by `alpha`, `signed`, `mtc`, `missing_as_zero`. +- **Bipartite projections *and hypergraphs*** (`backbone_from_projection`): the + input `B` is a bipartite network **or a hypergraph, as an incidence matrix**; + models `sdsm` (incl. **SDSM-EC** when structural 0s/1s are present), `fdsm`, + `fixedfill`, `fixedrow`, `fixedcol`. Controlled by `alpha`, `signed`, `mtc`, + `missing_as_zero`, `trials`. +- **Unweighted networks** (`backbone_from_unweighted`): `sparsify`-family. +- A unified `backbone()` wrapper dispatching on input class, plus `bicm`, + `fastball`, and `print`/`summary`/`plot` for backbone objects. + +New since 2.x (from `NEWS.md`): modular rewrite organized by input type; the +`backbone()` wrapper; **hypergraph** projection input; SDSM-EC structural +constraints; backbone objects with `narrative`; removal of edgelist input and the +ordinal SDSM. + +## 3. Coverage matrix (models) + +| Backbone 3.0 export / model | `networkx-backbone` equivalent | Status | +|-----------------------------|-------------------------------|--------| +| `disparity` | `disparity_filter` / `disparity` | ✅ | +| `lans` | `lans_filter` / `lans` | ✅ | +| `mlf` | `marginal_likelihood_filter` / `mlf` | ✅ | +| `global` (threshold) | `global_threshold_filter` | ✅ (one-sided; signed variant → §4.2) | +| `sdsm` | `sdsm` | ✅ (edge constraints → §4.3) | +| `fdsm` | `fdsm` | ✅ | +| `fixedfill` / `fixedrow` / `fixedcol` | `fixedfill` / `fixedrow` / `fixedcol` | ✅ | +| `bicm` | `bicm` | ✅ | +| `fastball` | `fastball` | ✅ | +| `sparsify` (unweighted) | `sparsify` / `lspar` / `local_degree` | ✅ | +| `backbone_from_weighted/_projection/_unweighted`, `backbone` | same names | ✅ | +| **hypergraph** projection input | `hypergraph_to_bipartite` → `sdsm`/`fdsm`/`fixed*` | ✅ | + +**Beyond Backbone 3.0:** `noise_corrected_filter`, `ecm_filter`, +`multiple_linkage_analysis`; the entire `structural` and `proximity` families; +and the `hypergraph` module (`mdl_hypergraph_backbone`, +`statistically_validated_hypergraph`/`_cores`, `maximal_hyperedges`, +`order_filter`, `s_components`) with HIF/xgi/HyperNetX/HypergraphX/HAT interop. + +## 4. Gap analysis and proposed implementations + +### 4.1 Multiple-testing correction (`mtc`) — ✅ implemented + +Backbone 3.0 adjusts the matrix of edge p-values with R's `p.adjust()` before +thresholding at `alpha`. `networkx-backbone` now provides `adjust_pvalues()` +(reproducing R's `p.adjust` for `bonferroni`, `holm`, `hochberg`, `bh`/`fdr`, +`by`; verified against `scipy.stats.false_discovery_control`) and an `mtc` +parameter on `threshold_filter`, so every p-value method is FDR/Bonferroni-aware +through one shared path. Original API sketch: + +```python +# networkx_backbone/filters.py +def adjust_pvalues(pvalues, method="bh"): + """Return multiplicity-adjusted p-values. + + method in {"bonferroni","holm","hochberg","hommel","bh"/"fdr","by"}. + """ + +def threshold_filter(G, score, threshold, mode="below", *, mtc="none", ...): + """When mtc != "none" and mode == "below", adjust the score attribute across + all edges with adjust_pvalues(..., mtc) before applying the threshold.""" +``` + +`adjust_pvalues` is ~30 lines of pure Python (Bonferroni/Holm step-down, +Benjamini–Hochberg/Yekutieli step-up); the existing `_bh_threshold` in +`hypergraph.py` is the BH building block. This makes every p-value method +(`disparity`, `mlf`, `lans`, `sdsm`, `fdsm`, `fixed*`) FDR/Bonferroni-aware +through one shared path, matching Backbone 3.0's `mtc` semantics. + +### 4.2 Signed backbones (`signed`) — ✅ implemented + +With `signed=TRUE`, Backbone 3.0 runs a **two-tailed** test and keeps edges that +are significantly *strong* (sign `+1`) **and** significantly *weak* (sign `−1`), +annotating each retained edge with a `sign`. + +`networkx-backbone` now exposes `signed=False` on the weighted scorers +(`disparity_filter`/`mlf`/`lans_filter`) and the projection null models +(`sdsm`/`fdsm`/`fixedfill`/`fixedrow`/`fixedcol`). Each adds the lower-tail +p-value, stores a two-sided p-value (`2·min(p_hi, p_lo)`, clipped to 1) and a +`"sign"` edge attribute (`+1` strong, `−1` weak). The `sign` flows through +`threshold_filter`/`boolean_filter` automatically (edge data is copied). +`backbone_from_weighted` and `backbone_from_projection` forward `signed`. + +### 4.3 SDSM with edge constraints (SDSM-EC) — ✅ implemented + +Backbone 3.0's `sdsm` switches to **SDSM-EC** (Neal & Neal 2023) when the +incidence matrix carries structural values: `10` = prohibited edge, `11` = +required edge. + +`networkx-backbone`'s `sdsm` now accepts `prohibited`/`required` as iterables of +`(agent, artifact)` pairs and fixes those cells' null probability to 0/1 before +the Poisson-binomial test. Defaults (`None`) leave behavior unchanged; unknown +nodes are ignored; constraints compose with `signed` and are forwarded by +`backbone_from_projection`. + +### 4.4 Narrative (`narrative`) — parked + +Backbone 3.0 can emit suggested methods text and citations for a chosen backbone. + +**Proposal (low priority).** A `describe_backbone(method, alpha=..., mtc=...)` +helper returning a citation/methods string per method. Pure formatting; no +algorithmic content. Could also populate a `narrative` field on the hypergraph +result objects. + +### 4.5 Non-gaps + +- **Hypergraph input** — covered (§3); this library additionally offers native + hypergraph backbones beyond projection. +- **`missing_as_zero`** — the projection scorers already test *all* agent pairs + (including zero co-occurrence), so absent edges are evaluated; a weighted-input + `missing_as_zero` toggle is a minor convenience if desired. +- **`print`/`summary`/`plot`** — this library returns NetworkX graphs and provides + its own `visualization` module and `measures` (`compare_backbones`). + +## 5. Status + +**`mtc`, `signed`, and SDSM-EC are implemented** (additive parameters; no change +to default behavior), closing the methodological gaps versus Backbone 3.0. Only +`narrative` (auto methods text — a utility, not an algorithm) is parked for a +future follow-up. + +## 6. References + +- Neal, Z. P. (2026). *Backbone 3.0: An R package for extracting network + backbones.* PLOS ONE. https://doi.org/10.1371/journal.pone.0349258 +- Neal, Z. P. (2022). *backbone: An R package to extract network backbones.* PLOS + ONE, 17(5), e0269137. https://doi.org/10.1371/journal.pone.0269137 +- Neal, Z. P., & Neal, J. W. (2023). *Stochastic Degree Sequence Model with Edge + Constraints (SDSM-EC) for Backbone Extraction.* Complex Networks 12, 127-136. + https://doi.org/10.1007/978-3-031-53468-3_11 diff --git a/docs/design/hypergraph-backboning.md b/docs/design/hypergraph-backboning.md new file mode 100644 index 0000000..24e8feb --- /dev/null +++ b/docs/design/hypergraph-backboning.md @@ -0,0 +1,341 @@ +# Design proposal: hypergraph network backboning + +- **Status:** Draft / for discussion +- **Scope:** Evaluation + proposed integration plan (no library code changed by this document) +- **Trigger:** Request to evaluate incorporating hypergraph network backboning + methods into `networkx-backbone`, specifically the method of *Kirkley, Felippe, + Malizia & Battiston, "Hypergraph backboning," arXiv:2606.00893 (2026)*, with + optional interoperability for `xgi`, `hypergraphx`, `HyperNetX`, and the + `Hypergraph-Analysis-Toolbox`, and an assessment of which hypergraph backbone + methods are genuinely distinct from the existing ensemble. + +> **Revision note (after reading the paper).** An earlier draft of this document +> was written without access to arXiv:2606.00893 and hypothesized that the paper +> belonged to the *statistical* hyperedge-filtering family (SVH/SVC). The paper +> has now been read in full: it is a **parameter-free, information-theoretic +> (MDL) compression** method — a *different paradigm* from statistical +> null-model testing. Sections 2, 6, and 9 are updated accordingly; the headline +> method to port is now the MDL backbone, with statistical/structural methods +> repositioned as complementary. + +--- + +## 1. Summary and recommendation + +**Recommendation: yes — port the paper's MDL hypergraph backbone as the headline +method, add optional interoperability with the major hypergraph libraries +(centered on the HIF interchange format), and offer a small set of complementary +hypergraph-native methods.** The case is now stronger than in the first draft: the +method is principled, **parameter-free** (unweighted), **handles weights** (the +paper's main novelty — no prior hypergraph filter did), needs only +`networkx` + `numpy`/`scipy`, and runs in minutes on real data. + +Hypergraph backboning organizes into three paradigms: + +| Paradigm | What it does | Output | Library status | +|----------|--------------|--------|----------------| +| **A — Projection backboning** | Hypergraph → weighted pairwise graph → null-model edge test | a graph | **already works** (verified, §5) | +| **B1 — Statistical hyperedge filtering** | Keep hyperedges over-expressed vs a null model (needs significance level α) | sub-hypergraph | not present (optional, §6) | +| **B2 — Information-theoretic (MDL) backboning ← this paper** | Compress nested/redundant hyperedges; keep the minimal "parent" set | sub-hypergraph | **port this** (§2, §9) | + +Key conclusions: + +1. **Family A already works today** with the existing `bipartite` module (§5). +2. **Most weighted-hypergraph backbones are generalizations** of methods we ship + and become reachable for free once a hypergraph can enter the pipeline (§6.1). +3. **The paper's MDL method is genuinely sui generis** — it exploits overlap and + nestedness, structural signatures unique to hypergraphs, via a global + compression objective with no dyadic analog and no significance parameter + (§2, §6.2). It is the right headline method to port. +4. **Optional interop is cheap** because all four target libraries converge on an + incidence/bipartite substrate and on the **HIF** JSON format; none needs to + become a *required* dependency (§7). + +Phasing (details in §9): **Phase 0** ingestion (HIF + per-library adapters) and +surfacing of Family A; **Phase 1** the MDL backbone (unweighted + weighted); +**Phase 2** complementary statistical (SVH/SVC) and structural (toplex, s-line) +methods. Out of scope: hypergraph neural networks; any required third-party +hypergraph dependency. + +## 2. The source paper (arXiv:2606.00893) + +*A. Kirkley, H. Felippe, F. Malizia, F. Battiston, "Hypergraph backboning" (2026).* + +**Idea.** Given a hypergraph `G` on `N` nodes (undirected hyperedges, no repeated +nodes within an edge, no multi-edges) with `L` distinct hyperedge sizes (orders), +find a backbone `B ⊆ G`: a subset of hyperedges ("parents") such that every +non-backbone hyperedge ("child") can be cheaply reconstructed from a parent it +overlaps. The best backbone is the one that **minimizes a two-part description +length** (MDL) — equivalently, that maximizes the structural redundancy +(overlap/nestedness) explained. It is **fully nonparametric** for unweighted +hypergraphs. + +**Unweighted objective.** With `log ≡ log2` and `C(n,k)` the binomial: + +- Transmit each parent `p ∈ B`: `H(p) = log L + log C(N, |p|)`; `L(B) = Σ_p H(p)`. +- Transmit each child `c` from its parent `p` (with `|p ∩ c| ≥ 1`): + `H(c|p) = log L + log min(|p|,|c|) + log C(|p|, |p∩c|) + log C(N−|p|, |c|−|p∩c|)`. +- Total: `L(G,B) = L(B) + Σ_{p∈B} Σ_{c∈∂p} H(c|p)`, where `∂p` are the children of + `p`. The optimum is `B* = argmin_B L(G,B)`. +- Equivalent reduced-mutual-information form: `L(G,B) = L(G,G) − Σ_c R(c, p(c))`, + so minimizing description length = **maximizing parent–child overlap/nestedness**. +- **Inverse compression ratio** `η = L(G,B*) / L(G,G) ∈ [0,1]` measures how + compressible (redundant) the hypergraph is (`η→0` very compressible, `η=1` none). + +**Weighted extension.** `Lw(G,B) = L(G,B) + Σ_e L(w(e), b_e)`, where `b_e∈{0,1}` is +backbone membership and weights follow a role-dependent **Poisson or Geometric** +prior under an empirical-Bayes mean constraint. A single hyperparameter +`γ ∈ (0,1]` trades off weight vs. topology: `γ=1` recovers the unweighted +objective (weights ignored); `γ→0` makes it infinitely costly to leave a +high-weight edge out of the backbone. The backbone-inclusion reward is **linear in +the weight** `w(e)`, with a closed-form weight threshold `w*` below which weight no +longer favors inclusion. (Integer weights `≥1`; continuous weights need a +resolution parameter.) **No prior hypergraph filtering method handled weights** — +this is the paper's central novelty and aligns with this library's weighted focus. + +**Optimization (Appendix D).** Exact minimization is combinatorial; the paper uses +greedy approximations on the **intersection graph** `Int(G)` (one node per +hyperedge; link two hyperedges that share ≥1 node). Parent–child assignments form a +**partition of `Int(G)` into disjoint stars** (each child has exactly one parent; +parents/children don't nest further). Two greedy schemes — "node" addition and +"edge" addition (the latter usually better) — are run and the lower description +length is kept. Greedy compression is **indistinguishable from exact** on small +samples. + +**Complexity / cost.** Bottleneck is building `Int(G)`: `O(Σ_i |G_i|²)` over node +neighborhoods `G_i = {e : i ∈ e}` (best case `O(N)`, worst `O(N|G|²)`); optional +random pair sampling gives `O(N s²)`. Empirically ~`N^1.17`, **≤6 minutes** on the +empirical corpus with a plain Python implementation. A **local variant** +(Appendix E) backbones each node neighborhood separately. + +**Inputs / outputs.** Input: a hyperedge list over `N` nodes (+ optional integer +weights). Output: the backbone sub-hypergraph `B`, the parent→children assignment +(star forest), and `η`. + +**Dependencies implied.** Only a hyperedge list and arithmetic; `Int(G)` is an +ordinary graph (build with NetworkX), and the log-binomials use +`scipy.special.gammaln`. **No third-party hypergraph library is required to +implement it.** No public reference code is cited (the authors describe a "simple +Python implementation"; datasets come from the Hypergraphx-data repository). + +## 3. Background: three paradigms + +- **A — projection backboning.** Represent `G` as an incidence (node × hyperedge) + structure, project to a node–node weighted graph, apply a degree-preserving null + model. Output is a graph. Lineage: Neal's `backbone` (Backbone 3.0, 2026), + Coscia & Neffke (2017). +- **B1 — statistical hyperedge filtering.** Keep hyperedges over-expressed vs a + configuration null, *given a significance level α*. Output is a sub-hypergraph. + Reference: Musciotto, Battiston & Mantegna (2021); impl in HGX (`get_svh`/`get_svc`). +- **B2 — information-theoretic (MDL) backboning — this paper.** Compress nested and + redundant hyperedges; keep the minimal parent set. Output is a sub-hypergraph. + **Parameter-free** (unweighted), weighted via one knob. Distinct from B1: it is a + global compression optimum, not a per-hyperedge hypothesis test, and needs no α. + +Families B1/B2 share the "sub-hypergraph output" problem that does not fit the +current graph-in/graph-out filters. + +## 4. Current library capabilities relevant to this + +- **Two idioms already in use:** *score-then-filter* (e.g. `disparity_filter` → + `threshold_filter`) and *direct boolean flag* (e.g. + `maximum_spanning_tree_backbone` → `boolean_filter` on `mst_keep`). The MDL + method is a global optimizer, so it maps onto the **boolean-flag idiom** (annotate + each hyperedge with an `mdl_keep` role), not onto per-edge p-value thresholding. +- **Dependency policy:** core `networkx`-only; `numpy`/`scipy` are the optional + `[full]` extra, imported lazily. The MDL method fits this exactly. +- **The `bipartite` module is already an incidence engine** (`_bipartite_projection_matrix`, + `sdsm`/`fdsm`/`fixed*`, `fastball`, `bicm`), so Family A is essentially built. + +## 5. Key finding: Family A already works today + +```python +import networkx as nx +import networkx_backbone as nb + +hyperedges = {"H1": [1, 2, 3], "H2": [1, 2, 3], "H3": [3, 4, 5, 6], "H4": [5, 6]} +nodes = sorted({v for members in hyperedges.values() for v in members}) + +B = nx.Graph() # incidence bipartite graph +B.add_nodes_from(nodes, bipartite=0) # nodes +B.add_nodes_from(hyperedges, bipartite=1) # hyperedges as "artifacts" +for h, members in hyperedges.items(): + for v in members: + B.add_edge(v, h) + +scored = nb.sdsm(B, agent_nodes=nodes, projection="hyper") +backbone = nb.threshold_filter(scored, "sdsm_pvalue", 0.30, mode="below") +# backbone.edges() -> [(1, 2), (4, 5), (4, 6), (5, 6)] (verified) +``` + +A large slice of "hypergraph backboning" is therefore a **latent, undocumented +capability** today (projection family). The paper's method is a *different* output +type (a sub-hypergraph) and is the new work. + +## 6. Method inventory: generalizations vs. sui generis + +### 6.1 Generalizations (reachable once a hypergraph enters the pipeline) + +| Hypergraph method | Reduces to | Notes | +|-------------------|-----------|-------| +| Hyperedge global-weight threshold | `global_threshold_filter` | trivial | +| Degree-preserving projection null (SDSM/FDSM/fixed*) | `sdsm`/`fdsm`/`fixed*` on incidence | **already works** (§5) | +| Disparity / MLF / LANS / NC / ECM on the projection | matching statistical filter | applied to the projected graph | +| Clique- or line-graph expansion + any graph backbone | existing graph methods | transform, then any current method | + +### 6.2 Sui generis methods (no clean dyadic analog) + +| Method | Paradigm | Reference / impl | Distinctness | +|--------|----------|------------------|--------------| +| **MDL hypergraph backbone (this paper)** | B2 (compression) | Kirkley+ 2026 | **Headline.** Global MDL optimum over parent/child overlap & nestedness; parameter-free; weighted via γ; sub-hypergraph output. Naive toplex/maximal-face reduction is a degenerate special case. | +| Statistically Validated Hypergraph / Cores (SVH/SVC) | B1 (statistical) | Musciotto+ 2021; HGX `get_svh`/`get_svc` | Complementary α-based alternative; group-level null-model test. | +| Toplex / inclusion (encapsulation) reduction | structural | HNX `toplexes()`, XGI `encapsulation_dag` | Cheap heuristic; subsumed by the MDL objective. | +| s-connectivity / s-line-graph backbone | structural | HNX `s_components` | Parameterized by shared-node threshold `s`; no dyadic analog. | +| Order-resolved hyperedge filtering | utility | building block | Meaningful only with variable arity. | + +**Bottom line for the user's question:** hypergraph backboning is *mostly* +generalizations of existing methods (§6.1), **but the paper's MDL method is +genuinely new** and cannot be produced by the current ensemble — both because the +objective exploits higher-order overlap/nestedness and because the output is a +sub-hypergraph. It is worth porting; SVH/SVC and the structural primitives are +worthwhile but secondary. + +## 7. Optional interoperability with hypergraph libraries + +All four target libraries converge on an incidence/bipartite substrate and **all +support the HIF JSON interchange format**, so interop is cheap and adds **no +required dependency**. + +| Library | Hypergraph type(s) | → incidence / bipartite (in) | ← construct (from our output) | HIF | +|---------|--------------------|------------------------------|-------------------------------|-----| +| **XGI** (`xgi`) | `Hypergraph`, `DiHypergraph`, `SimplicialComplex` | `xgi.to_bipartite_graph`, `xgi.to_incidence_matrix` | `xgi.from_bipartite_graph`, `xgi.from_incidence_matrix` | `xgi.read_hif`/`write_hif` | +| **HyperNetX** (`hypernetx`) | `Hypergraph` | `.bipartite()`, `.incidence_matrix()`, `.incidence_dict` | `Hypergraph.from_bipartite`, `.restrict_to_edges(keep)` | supported | +| **HypergraphX** (`hypergraphx`) | `Hypergraph`, Temporal/Directed/Multiplex | `.binary_incidence_matrix(return_mapping=True)` | `Hypergraph(edge_list=...)` | supported | +| **HAT** (`HAT`) | `Hypergraph` (tensor/incidence) | `.incidence_matrix` | `Hypergraph(incidence_matrix=...)` | import/export | + +- **`restrict_to_edges`** (HNX) is the natural way to return the MDL backbone as a + native object in each library; HGX/HNX also give SVH/toplex references. +- Design: (1) **HIF** as a stdlib-only hub (`read_hif`/`write_hif` into our + hyperedge-list form); (2) thin lazy `from_*`/`to_*` adapters (~10–20 lines each, + since the converters above already exist); (3) optional extras + `networkx-backbone[xgi|hypernetx|hypergraphx|hat]`. Core install unchanged; each + adapter imports its library lazily and errors with an install hint if absent. + +## 8. Gap analysis and challenges + +1. **No native hypergraph type in NetworkX.** Use a hyperedge-list / incidence + representation internally (§9.1); the MDL optimizer's `Int(G)` is itself a + NetworkX graph. +2. **Sub-hypergraph output.** Doesn't flow through `threshold_filter`. Mirror the + boolean-flag idiom: annotate hyperedges with an `mdl_keep` role + a + `hyperedge_filter`, and optionally return a native library object. +3. **Global optimizer, not a per-edge score.** The MDL backbone is a combinatorial + optimum (greedy), unlike independent per-edge p-values — document it as a + "direct" method like the spanning-tree/metric backbones. +4. **One knob for weights (γ).** Parameter-free unweighted; `γ` only for weighted, + default `γ=1`. Far less parameter burden than α-based methods. +5. **Dependency policy preserved.** `networkx` + `numpy`/`scipy` + (`scipy.special.gammaln`); the four libraries stay optional; HIF is stdlib-only. + +## 9. Proposed design and re-scoped plan + +### 9.1 Representation + +Internal: a **hyperedge list** (tuples/`frozenset`s) + node ordering, with +incidence matrix / incidence bipartite graph available on demand. Zero new deps; +matches the substrate every target library and HIF share. + +### 9.2 Phase 0 — ingestion + surface Family A (**implemented**) + +`networkx_backbone.hypergraph_io` provides `hypergraph_to_bipartite` (enabling the +existing SDSM/FDSM/fixed projection backbones on hypergraphs), `read_hif`/`write_hif` +(stdlib-only HIF interchange, cross-checked against xgi in tests), and lazy +`from_*`/`to_*` adapters for xgi, HyperNetX, HypergraphX, and HAT. No required +dependency is added. (A worked tutorial is still a nice-to-have follow-up.) + +### 9.3 Phase 1 — MDL hypergraph backbone (the paper) + +A `hypergraph` module implementing Kirkley+ 2026: + +```python +# networkx_backbone/hypergraph.py +def mdl_hypergraph_backbone( + hyperedges, weights=None, gamma=1.0, prior="poisson", + method="auto", # "node" | "edge" | "auto" (run both, keep lower L) + sample_pairs=None, seed=None, +): + """Return the MDL-optimal backbone (Kirkley et al. 2026). + + Returns a result with: backbone hyperedges, parent->children assignment, + description_length, and compression_ratio eta. Also annotates each hyperedge + with an `mdl_keep` boolean role. + """ + +def hyperedge_filter(scored, role="mdl_keep"): + """Boolean-flag filter -> the backbone sub-hypergraph.""" + +def hypergraph_compression_ratio(hyperedges, backbone=None, weights=None, gamma=1.0): + """Inverse compression ratio eta (Eq. 8) as an evaluation measure.""" +``` + +Building blocks: `intersection_graph(hyperedges)` (NetworkX), the parent/child MDL +terms (`gammaln`-based), greedy "node"/"edge" optimizers over the star partition, +and optional pair-sampling for very large/dense inputs. The local variant +(Appendix E) is a natural follow-up. + +### 9.4 Phase 2 — complementary methods + +Structural sui generis methods (**implemented**): `maximal_hyperedges` +(inclusion/toplex reduction), `order_filter` (order-resolved filtering), and +`intersection_graph(..., s=...)` + `s_components` (s-line graph / s-connectivity). +These are stdlib + networkx only and return plain hyperedge collections. + +Statistical sui generis methods (**implemented**): `statistically_validated_hypergraph` +(SVH) and `statistically_validated_cores` (SVC) — the Musciotto, Battiston & +Mantegna (2021) null-model validation, faithfully re-implemented from HGX +`get_svh`/`get_svc` with a lazy `scipy` dependency (binomial tail + Benjamini-Hochberg +FDR). P-values verified exactly against `scipy.stats.binom`. + +### 9.5 Testing strategy + +- Phase 0: round-trip converter tests per library + HIF; equivalence of + `hypergraph_to_bipartite + sdsm` with the manual recipe. +- Phase 1: **reproduce the paper's controlled synthetic experiments** — planted + fully-nested simplices and parent/child-with-noise hypergraphs (Figs. 1–3): + backbone recovers planted top faces; `η` behaves as reported; greedy ≈ exact on + tiny inputs; `γ=1` ≡ unweighted; `γ→0` forces high-weight edges in; determinism + via `seed`. +- Cross-validate against HGX/HNX where installed (skipped otherwise). + +## 10. Resolved questions and remaining implementation decisions + +The first-draft "open questions" are now answered by the paper: +family = **B2 (MDL compression)**; null model = **none** (information-theoretic, +parameter-free unweighted); correction = **n/a**; output = **sub-hypergraph** + star +forest + `η`; relationship to SVH/SVC = **distinct paradigm**. + +Resolved during implementation: +1. Optimizer — **both** greedy schemes implemented (`method="edge"`, `"node"`), + with `method="auto"` (the default) running both and keeping the lower + description length, per the paper. +2. Weight prior default `poisson`, `γ` default `1.0`. +3. Result object — a `HypergraphBackbone` dataclass (backbone, assignment, `η`, + description lengths). + +Remaining nice-to-have: the Appendix E local (per-node-neighbourhood) variant. + +## 11. References + +- **Kirkley, Felippe, Malizia & Battiston — Hypergraph backboning — arXiv:2606.00893 (2026).** +- Musciotto, Battiston & Mantegna — Detecting informative higher-order interactions + in statistically validated hypergraphs — Communications Physics (2021), + arXiv:2103.16484. https://www.nature.com/articles/s42005-021-00710-4 +- Backbone 3.0: An R package for extracting network backbones — PLOS One (2026). + https://journals.plos.org/plosone/article?id=10.1371%2Fjournal.pone.0349258 +- Coscia & Neffke — Network backboning with noisy data — arXiv:1906.09081. +- HIF: The hypergraph interchange format — arXiv:2507.11520; + standard: https://github.com/HIF-org/HIF-standard +- Hypergraphx (HGX) — arXiv:2303.15356. https://github.com/HGX-Team/hypergraphx +- HyperNetX (HNX) — arXiv:2310.11626. https://github.com/pnnl/HyperNetX +- XGI — https://github.com/xgi-org/xgi +- Hypergraph Analysis Toolbox (HAT) — https://github.com/Jpickard1/Hypergraph-Analysis-Toolbox diff --git a/docs/tutorials/bipartite_backbone.rst b/docs/tutorials/bipartite_backbone.rst index 4755643..4e35eac 100644 --- a/docs/tutorials/bipartite_backbone.rst +++ b/docs/tutorials/bipartite_backbone.rst @@ -134,8 +134,49 @@ Partition selection note which makes partition selection straightforward. In general, pass whichever partition you want to project as ``agent_nodes``. +Signed backbones +---------------- + +By default the null-model tests are one-tailed and retain only +*significantly strong* co-occurrences. Pass ``signed=True`` to ``sdsm``, +``fdsm``, ``fixedfill``, ``fixedrow``, or ``fixedcol`` for a **two-tailed** test: +the stored p-value becomes two-sided and each edge gains a ``"sign"`` attribute +(``+1`` significantly strong, ``-1`` significantly weak):: + + H = nb.sdsm(B, agent_nodes=women_nodes, signed=True) + backbone = nb.threshold_filter(H, "sdsm_pvalue", 0.05, mode="below") + weak = [(u, v) for u, v, d in backbone.edges(data=True) if d["sign"] == -1] + +Edge constraints (SDSM-EC) +-------------------------- + +The Stochastic Degree Sequence Model with Edge Constraints (Neal & Neal, 2023) +lets you fix some cells of the null model: ``prohibited`` edges cannot occur +(null probability 0) and ``required`` edges must occur (null probability 1). +Pass them to ``sdsm`` as iterables of ``(agent, artifact)`` pairs:: + + H = nb.sdsm( + B, + agent_nodes=women_nodes, + prohibited=[(women_nodes[0], event_nodes[0])], + required=[(women_nodes[1], event_nodes[1])], + ) + References ---------- - Coscia, M., & Neffke, F. M. (2017). *Network backboning with noisy data*. https://arxiv.org/abs/1906.09081 +- Neal, Z. P. (2014). *The backbone of bipartite projections*. Social Networks, + 39, 84-97. +- Neal, Z. P., Domagalski, R., & Sagan, B. (2021). *Comparing alternatives to the + fixed degree sequence model for extracting the backbone of bipartite + projections*. Scientific Reports, 11, 23929. +- Saracco, F., Di Clemente, R., Gabrielli, A., & Squartini, T. (2015). + *Randomizing bipartite networks: the case of the World Trade Web*. Scientific + Reports, 5, 10595. +- Godard, K., & Neal, Z. P. (2022). *fastball: A fast algorithm to sample + bipartite graphs with fixed degree sequences*. J. Complex Networks, 10(6), + cnac049. +- Neal, Z. P., & Neal, J. W. (2023). *Stochastic Degree Sequence Model with Edge + Constraints (SDSM-EC) for Backbone Extraction*. Complex Networks 12, 127-136. diff --git a/docs/tutorials/hypergraph_backbone.rst b/docs/tutorials/hypergraph_backbone.rst new file mode 100644 index 0000000..e66930f --- /dev/null +++ b/docs/tutorials/hypergraph_backbone.rst @@ -0,0 +1,150 @@ +Hypergraph Backbones +==================== + +This tutorial demonstrates backbone extraction for **hypergraphs** -- networks +whose edges (hyperedges) join an arbitrary number of nodes, capturing +higher-order interactions that pairwise graphs cannot. It covers the +information-theoretic (MDL) backbone, structural reductions, statistical +validation, and interoperability with the higher-order ecosystem. + +What is a hypergraph backbone? +------------------------------ + +Real higher-order datasets -- co-authorship teams, group chats, protein +complexes -- are often dense and redundant, with smaller interactions nested +inside larger ones. A hypergraph backbone keeps only the essential hyperedges. +Unlike the :doc:`bipartite_backbone` methods (which project a hypergraph to a +pairwise graph), these methods return a **sub-hypergraph**: a subset of the +original hyperedges. + +A hypergraph is represented as a plain iterable of hyperedges, each an iterable +of node labels:: + + import networkx_backbone as nb + + H = [ + (1, 2, 3, 4), # a 4-node interaction + (1, 2, 3), # nested inside the first hyperedge (redundant) + (2, 3, 4), # also nested + (8, 9), # a separate pairwise interaction + ] + +MDL backbone (parameter-free) +----------------------------- + +:func:`~networkx_backbone.mdl_hypergraph_backbone` implements the +information-theoretic method of Kirkley, Felippe, Malizia & Battiston (2026). It +keeps a minimal set of "parent" hyperedges from which the remaining "child" +hyperedges can be reconstructed via overlap and nestedness, by minimizing a +two-part description length. It is **fully nonparametric** for unweighted +hypergraphs:: + + result = nb.mdl_hypergraph_backbone(H) + + print(result.backbone) # [frozenset({1, 2, 3, 4}), frozenset({8, 9})] + print(result.compression_ratio) # inverse compression ratio eta in [0, 1] + print(result.assignment) # parent -> [pruned child hyperedges] + +The result is a :class:`~networkx_backbone.HypergraphBackbone`. The inverse +compression ratio ``eta`` (also available via +:func:`~networkx_backbone.hypergraph_compression_ratio`) is near 0 for highly +redundant hypergraphs and equals 1 when no compression is possible. + +Choosing the optimizer +~~~~~~~~~~~~~~~~~~~~~~~ + +Exact minimization is combinatorial, so the backbone is found with a greedy +heuristic over the intersection graph. Two schemes are available (Appendix D of +the paper): + +- ``method="edge"`` adds parent--child links by decreasing gain; +- ``method="node"`` adds backbone parents by decreasing total gain; +- ``method="auto"`` (the default) runs **both** and keeps the lower description + length, exactly as the paper does. + +:: + + edge = nb.mdl_hypergraph_backbone(H, method="edge") + node = nb.mdl_hypergraph_backbone(H, method="node") + auto = nb.mdl_hypergraph_backbone(H, method="auto") + + print(edge.description_length, node.description_length, auto.description_length) + +``"auto"`` never does worse than either scheme. Use ``"edge"`` for the fastest +single pass on large hypergraphs. + +Weighted hypergraphs +~~~~~~~~~~~~~~~~~~~~~ + +When hyperedges carry integer weights (interaction strengths), the backbone +balances topology against weight through a single knob ``gamma`` in ``(0, 1]``. +``gamma=1`` ignores weights; ``gamma`` near 0 makes it costly to leave a +high-weight hyperedge out of the backbone:: + + weights = [1, 1, 5, 1] # the {2, 3, 4} interaction is strong + result = nb.mdl_hypergraph_backbone(H, weights=weights, gamma=0.1) + +Structural reductions +--------------------- + +Several purely structural backbones have no analog in pairwise graphs: + +- :func:`~networkx_backbone.maximal_hyperedges` -- inclusion (toplex) reduction: + drop every hyperedge contained in another:: + + nb.maximal_hyperedges(H) # [frozenset({1, 2, 3, 4}), frozenset({8, 9})] + +- :func:`~networkx_backbone.order_filter` -- keep hyperedges by order (size):: + + nb.order_filter(H, min_order=3) # interactions of three or more nodes + +- :func:`~networkx_backbone.s_components` -- group hyperedges into s-connected + components (those sharing at least ``s`` nodes):: + + nb.s_components(H, s=2) + +Statistical validation (SVH / SVC) +---------------------------------- + +:func:`~networkx_backbone.statistically_validated_hypergraph` keeps hyperedges +that **recur** more often than expected under a null model preserving node +activity (Musciotto, Battiston & Mantegna, 2021). Repeated hyperedges (or +integer ``weights``) are treated as interaction counts:: + + events = [(1, 2)] * 5 + [(3, 4)] * 100 # {1, 2} co-occurs 5x; {3, 4} is background + result = nb.statistically_validated_hypergraph(events, alpha=0.05) + print(result.validated) # [frozenset({1, 2})] + print(result.pvalues[frozenset({1, 2})]) + +:func:`~networkx_backbone.statistically_validated_cores` additionally validates +significant sub-groups (cores), attributing significance to the largest +validated group. Both require ``scipy``. + +Interoperability +---------------- + +Convert a hypergraph to its incidence bipartite graph to reuse the projection +backbones (:doc:`bipartite_backbone`) such as SDSM and FDSM:: + + B, nodes = nb.hypergraph_to_bipartite(H) + scored = nb.sdsm(B, agent_nodes=nodes) + projection_backbone = nb.threshold_filter(scored, "sdsm_pvalue", 0.05, mode="below") + +Exchange hypergraphs with the wider ecosystem through the HIF interchange format +(standard library only) or the optional ``xgi``, ``HyperNetX``, ``HypergraphX``, +and HAT adapters (imported lazily):: + + nb.write_hif(H, "hypergraph.hif") + edges = nb.read_hif("hypergraph.hif") + + # With the relevant library installed: + # hg = nb.to_xgi(H); edges = nb.from_hypernetx(hnx_hypergraph) + +References +---------- + +- Kirkley, A., Felippe, H., Malizia, F., & Battiston, F. (2026). *Hypergraph + backboning*. arXiv:2606.00893. +- Musciotto, F., Battiston, F., & Mantegna, R. N. (2021). *Detecting informative + higher-order interactions in statistically validated hypergraphs*. + Communications Physics, 4, 218. diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index a070a24..77626bc 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -11,6 +11,7 @@ extraction methods. proximity_scoring structural_methods bipartite_backbone + hypergraph_backbone comparing_methods unweighted_sparsification les_miserables_benchmark diff --git a/docs/tutorials/statistical_backbone.rst b/docs/tutorials/statistical_backbone.rst index a58f037..6239939 100644 --- a/docs/tutorials/statistical_backbone.rst +++ b/docs/tutorials/statistical_backbone.rst @@ -125,3 +125,59 @@ Use :func:`~networkx_backbone.compare_backbones` to compare the results:: ef = metrics["edge_fraction"] nf = metrics["node_fraction"] print(f"{name:25s}: edges={ef:.1%}, nodes={nf:.1%}") + +Multiple-testing correction +--------------------------- + +Because a statistical backbone tests every edge, correcting for multiple +comparisons avoids retaining edges by chance. Pass ``mtc`` to +:func:`~networkx_backbone.threshold_filter` (or use +:func:`~networkx_backbone.adjust_pvalues` directly). Supported corrections match +R's ``p.adjust``: ``"bonferroni"``, ``"holm"``, ``"hochberg"``, ``"bh"``/``"fdr"`` +(Benjamini-Hochberg), and ``"by"`` (Benjamini-Yekutieli):: + + scored = nb.disparity_filter(G) + + raw = nb.threshold_filter(scored, "disparity_pvalue", 0.05) + fdr = nb.threshold_filter(scored, "disparity_pvalue", 0.05, mtc="bh") + bonf = nb.threshold_filter(scored, "disparity_pvalue", 0.05, mtc="bonferroni") + + print(f"raw={raw.number_of_edges()}, BH={fdr.number_of_edges()}, " + f"Bonferroni={bonf.number_of_edges()}") + +Signed backbones +---------------- + +By default the statistical filters run a one-tailed test that keeps only +*significantly strong* edges. With ``signed=True``, the +``disparity_filter``, ``marginal_likelihood_filter``, and ``lans_filter`` run a +**two-tailed** test: the stored p-value becomes two-sided, and each edge gains a +``"sign"`` attribute that is ``+1`` for a significantly strong edge and ``-1`` +for a significantly weak one:: + + signed = nb.disparity_filter(G, signed=True) + backbone = nb.threshold_filter(signed, "disparity_pvalue", 0.05, mtc="bh") + + positive = [(u, v) for u, v, d in backbone.edges(data=True) if d["sign"] == 1] + negative = [(u, v) for u, v, d in backbone.edges(data=True) if d["sign"] == -1] + print(f"strong (+): {len(positive)}, weak (-): {len(negative)}") + +References +---------- + +- Serrano, M. A., Boguna, M., & Vespignani, A. (2009). *Extracting the multiscale + backbone of complex weighted networks*. PNAS, 106(16), 6483-6488. +- Coscia, M., & Neffke, F. M. (2017). *Network backboning with noisy data*. + Proc. IEEE ICDE, 425-436. +- Dianati, N. (2016). *Unwinding the hairball graph: Pruning algorithms for + weighted complex networks*. Physical Review E, 93, 012304. +- Gemmetto, V., Cardillo, A., & Garlaschelli, D. (2017). *Irreducible network + backbones: unbiased graph filtering via maximum entropy*. arXiv:1706.00230. +- Foti, N. J., Hughes, J. M., & Rockmore, D. N. (2011). *Nonparametric + sparsification of complex multiscale networks*. PLOS ONE, 6(2), e16431. +- Van Nuffel, N., Heyndrickx, C., & Wets, G. (2010). *Measuring hierarchy and + reciprocity in networks*. +- Benjamini, Y., & Hochberg, Y. (1995). *Controlling the false discovery rate*. + J. Royal Statistical Society B, 57(1), 289-300. +- Benjamini, Y., & Yekutieli, D. (2001). *The control of the false discovery rate + in multiple testing under dependency*. Annals of Statistics, 29(4), 1165-1188. diff --git a/networkx_backbone/__init__.py b/networkx_backbone/__init__.py index b9f2163..768dedf 100644 --- a/networkx_backbone/__init__.py +++ b/networkx_backbone/__init__.py @@ -2,13 +2,16 @@ Backbone extraction algorithms for complex networks. This package provides algorithms for extracting backbone structures from -networks, organized into nine submodules: +networks, organized into ten submodules: - **statistical**: Hypothesis-testing methods (disparity, noise-corrected, etc.) - **structural**: Topology-based methods (threshold, spanning tree, salience, etc.) - **proximity**: Neighborhood-similarity edge scoring (Jaccard, Dice, cosine, etc.) - **hybrid**: Combined statistical/structural methods (GLAB) - **bipartite**: Bipartite projection backbones (SDSM, FDSM, fixed models, wrappers) +- **hypergraph**: Hypergraph backbones (MDL pruning, statistical validation, + structural reduction); ingestion/interoperability helpers live in + ``hypergraph_io`` (HIF, xgi, HyperNetX, HypergraphX, HAT) - **unweighted**: Sparsification for unweighted graphs (LSpar, local degree) - **filters**: Post-hoc filtering utilities (threshold, fraction, boolean, consensus) - **measures**: Evaluation measures for comparing backbones @@ -20,6 +23,8 @@ from networkx_backbone.proximity import * # noqa: F401,F403 from networkx_backbone.hybrid import * # noqa: F401,F403 from networkx_backbone.bipartite import * # noqa: F401,F403 +from networkx_backbone.hypergraph import * # noqa: F401,F403 +from networkx_backbone.hypergraph_io import * # noqa: F401,F403 from networkx_backbone.unweighted import * # noqa: F401,F403 from networkx_backbone.filters import * # noqa: F401,F403 from networkx_backbone.measures import * # noqa: F401,F403 @@ -85,6 +90,31 @@ "backbone_from_weighted", "backbone_from_unweighted", "backbone", + # Hypergraph + "intersection_graph", + "mdl_hypergraph_backbone", + "hypergraph_compression_ratio", + "HypergraphBackbone", + "maximal_hyperedges", + "order_filter", + "s_components", + "statistically_validated_hypergraph", + "statistically_validated_cores", + "ValidatedHypergraph", + "svh", + "svc", + # Hypergraph interoperability + "hypergraph_to_bipartite", + "read_hif", + "write_hif", + "from_xgi", + "to_xgi", + "from_hypernetx", + "to_hypernetx", + "from_hypergraphx", + "to_hypergraphx", + "from_hat", + "to_hat", # Unweighted "sparsify", "lspar", @@ -95,6 +125,7 @@ "fraction_filter", "boolean_filter", "consensus_backbone", + "adjust_pvalues", # Measures "node_fraction", "edge_fraction", diff --git a/networkx_backbone/bipartite.py b/networkx_backbone/bipartite.py index 4545faf..4e09641 100644 --- a/networkx_backbone/bipartite.py +++ b/networkx_backbone/bipartite.py @@ -58,6 +58,17 @@ def _validate_bipartite(B, agent_nodes): raise nx.NetworkXError("agent_nodes contains nodes not in B.") +def _combine_tails(p_hi, p_lo): + """Two-sided p-value and sign from upper/lower tails (see statistical module). + + ``sign`` is ``+1`` for a significantly strong (over-represented) co-occurrence + and ``-1`` for a significantly weak (under-represented) one. + """ + p_hi = min(max(p_hi, 0.0), 1.0) + p_lo = min(max(p_lo, 0.0), 1.0) + return min(1.0, 2.0 * min(p_hi, p_lo)), (1 if p_hi <= p_lo else -1) + + def _bipartite_projection_matrix(B, agent_nodes): """Build the co-occurrence matrix for agent nodes. @@ -374,6 +385,12 @@ def bicm(B, agent_nodes, return_labels=False): agent ``i`` and artifact ``k``. (P, agents, artifacts) : tuple Returned when ``return_labels=True``. + + References + ---------- + .. [1] Saracco, F., Di Clemente, R., Gabrielli, A., & Squartini, T. (2015). + Randomizing bipartite networks: the case of the World Trade Web. + *Scientific Reports*, 5, 10595. """ import numpy as np @@ -413,6 +430,12 @@ def fastball(matrix, n_swaps=None, seed=None): ------- randomized : np.ndarray Randomized binary matrix with preserved row/column sums. + + References + ---------- + .. [1] Godard, K., & Neal, Z. P. (2022). fastball: A fast algorithm to + sample bipartite graphs with fixed degree sequences. *Journal of Complex + Networks*, 10(6), cnac049. """ import numpy as np @@ -463,6 +486,7 @@ def fixedfill( B, agent_nodes, alpha=0.05, + signed=False, projection="simple", projection_weight="weight", projection_directed=False, @@ -503,9 +527,16 @@ def fixedfill( for i in range(na): for j in range(i + 1, na): obs = int(observed[i, j]) - pval = float(sp_stats.binom.sf(obs - 1, nf, q)) - if pval < alpha: - backbone.add_edge(agents[i], agents[j], fixedfill_pvalue=pval) + p_hi = float(sp_stats.binom.sf(obs - 1, nf, q)) + if signed: + p_lo = float(sp_stats.binom.cdf(obs, nf, q)) + pval, sign = _combine_tails(p_hi, p_lo) + if pval < alpha: + backbone.add_edge( + agents[i], agents[j], fixedfill_pvalue=pval, sign=sign + ) + elif p_hi < alpha: + backbone.add_edge(agents[i], agents[j], fixedfill_pvalue=p_hi) return _apply_projection_weights( backbone, @@ -523,6 +554,7 @@ def fixedrow( B, agent_nodes, alpha=0.05, + signed=False, projection="simple", projection_weight="weight", projection_directed=False, @@ -558,9 +590,16 @@ def fixedrow( obs = int(observed[i, j]) di = int(row_sums[i]) dj = int(row_sums[j]) - pval = float(sp_stats.hypergeom.sf(obs - 1, nf, di, dj)) - if pval < alpha: - backbone.add_edge(agents[i], agents[j], fixedrow_pvalue=pval) + p_hi = float(sp_stats.hypergeom.sf(obs - 1, nf, di, dj)) + if signed: + p_lo = float(sp_stats.hypergeom.cdf(obs, nf, di, dj)) + pval, sign = _combine_tails(p_hi, p_lo) + if pval < alpha: + backbone.add_edge( + agents[i], agents[j], fixedrow_pvalue=pval, sign=sign + ) + elif p_hi < alpha: + backbone.add_edge(agents[i], agents[j], fixedrow_pvalue=p_hi) return _apply_projection_weights( backbone, @@ -578,6 +617,7 @@ def fixedcol( B, agent_nodes, alpha=0.05, + signed=False, projection="simple", projection_weight="weight", projection_directed=False, @@ -617,12 +657,24 @@ def fixedcol( obs = float(observed[i, j]) if sigma2 > 0: z = (obs - mu) / np.sqrt(sigma2) - pval = float(1.0 - sp_stats.norm.cdf(z)) + p_hi = float(1.0 - sp_stats.norm.cdf(z)) + p_lo = float(sp_stats.norm.cdf(z)) + elif obs > mu: + p_hi, p_lo = 0.0, 1.0 + elif obs < mu: + p_hi, p_lo = 1.0, 0.0 + else: + p_hi, p_lo = 1.0, 1.0 + if signed: + pval, sign = _combine_tails(p_hi, p_lo) + if pval < alpha: + backbone.add_edge( + agents[i], agents[j], fixedcol_pvalue=pval, sign=sign + ) else: - pval = 0.0 if obs > mu else 1.0 - pval = float(max(min(pval, 1.0), 0.0)) - if pval < alpha: - backbone.add_edge(agents[i], agents[j], fixedcol_pvalue=pval) + p_hi = float(max(min(p_hi, 1.0), 0.0)) + if p_hi < alpha: + backbone.add_edge(agents[i], agents[j], fixedcol_pvalue=p_hi) return _apply_projection_weights( backbone, @@ -703,6 +755,7 @@ def backbone_from_projection( projection_directed=projection_directed, projection_max_iter=projection_max_iter, projection_tol=projection_tol, + **kwargs, ) if method_l == "fixedrow": return fixedrow( @@ -714,6 +767,7 @@ def backbone_from_projection( projection_directed=projection_directed, projection_max_iter=projection_max_iter, projection_tol=projection_tol, + **kwargs, ) if method_l == "fixedcol": return fixedcol( @@ -725,6 +779,7 @@ def backbone_from_projection( projection_directed=projection_directed, projection_max_iter=projection_max_iter, projection_tol=projection_tol, + **kwargs, ) raise ValueError( "Unknown projection method. Choose one of: " @@ -737,6 +792,8 @@ def backbone_from_weighted( method="disparity", weight="weight", alpha=0.05, + signed=False, + mtc="none", collapse_multiedges=True, edge_type_attr=None, **kwargs, @@ -747,6 +804,13 @@ def backbone_from_weighted( Parameters ---------- + signed : bool, optional (default=False) + If ``True``, use a two-tailed test (for ``disparity``/``mlf``/``lans``) + so the backbone retains significantly strong and significantly weak + edges, each annotated with a ``"sign"`` attribute. + mtc : string, optional (default="none") + Multiple-testing correction applied before thresholding (see + :func:`~networkx_backbone.adjust_pvalues`). collapse_multiedges : bool, optional (default=True) If ``True`` and ``G`` is a ``MultiGraph`` or ``MultiDiGraph``, collapse parallel edges using @@ -770,16 +834,16 @@ def backbone_from_weighted( method_l = method.lower() if method_l in ("disparity", "disparity_filter"): - scored = disparity_filter(G, weight=weight) - return threshold_filter(scored, "disparity_pvalue", alpha, mode="below") + scored = disparity_filter(G, weight=weight, signed=signed) + return threshold_filter(scored, "disparity_pvalue", alpha, mode="below", mtc=mtc) if method_l in ("mlf", "marginal_likelihood", "marginal_likelihood_filter"): - scored = marginal_likelihood_filter(G, weight=weight) - return threshold_filter(scored, "ml_pvalue", alpha, mode="below") + scored = marginal_likelihood_filter(G, weight=weight, signed=signed) + return threshold_filter(scored, "ml_pvalue", alpha, mode="below", mtc=mtc) if method_l in ("lans", "lans_filter"): - scored = lans_filter(G, weight=weight) - return threshold_filter(scored, "lans_pvalue", alpha, mode="below") + scored = lans_filter(G, weight=weight, signed=signed) + return threshold_filter(scored, "lans_pvalue", alpha, mode="below", mtc=mtc) if method_l in ("global", "global_threshold", "global_threshold_filter"): threshold = kwargs.get("threshold") @@ -864,6 +928,9 @@ def sdsm( agent_nodes, alpha=0.05, weight=None, + signed=False, + prohibited=None, + required=None, projection="simple", projection_weight="weight", projection_directed=False, @@ -891,6 +958,16 @@ def sdsm( weight : None or string, optional (default=None) Not used for SDSM (bipartite is unweighted); reserved for API consistency. + signed : bool, optional (default=False) + If ``True``, run a two-tailed test: the stored p-value becomes two-sided + and a ``"sign"`` edge attribute marks significantly strong (``+1``) + versus significantly weak (``-1``) co-occurrences. + prohibited : iterable of (agent, artifact) pairs or None, optional + SDSM-EC edge constraints (Neal & Neal, 2023): cells fixed to probability + 0 in the null model (edges that cannot occur). + required : iterable of (agent, artifact) pairs or None, optional + SDSM-EC edge constraints: cells fixed to probability 1 in the null model + (edges that must occur). projection : {"simple", "hyper", "probs", "ycn"}, optional Projection weighting assigned to each returned edge. projection_weight : str, optional @@ -919,6 +996,9 @@ def sdsm( .. [1] Neal, Z. P. (2014). The backbone of bipartite projections: Inferring relationships from co-authorship, co-sponsorship, co-attendance and other co-behaviors. *Social Networks*, 39, 84-97. + .. [2] Neal, Z. P., & Neal, J. W. (2023). Stochastic Degree Sequence Model + with Edge Constraints (SDSM-EC) for Backbone Extraction. *Complex + Networks 12*, 127-136. Examples -------- @@ -963,6 +1043,17 @@ def sdsm( P = np.outer(row_sums, col_sums) / total P = np.clip(P, 0, 1) + # SDSM-EC (Neal & Neal 2023): fix the null probability of constrained cells. + if prohibited or required: + a_idx = {v: idx for idx, v in enumerate(agents)} + f_idx = {v: idx for idx, v in enumerate(artifacts)} + for a, f in required or (): + if a in a_idx and f in f_idx: + P[a_idx[a], f_idx[f]] = 1.0 + for a, f in prohibited or (): + if a in a_idx and f in f_idx: + P[a_idx[a], f_idx[f]] = 0.0 + for i in range(na): for j in range(i + 1, na): obs = observed[i, j] @@ -976,12 +1067,23 @@ def sdsm( if sigma2 > 0: z = (obs - mu) / np.sqrt(sigma2) - pval = 1.0 - sp_stats.norm.cdf(z) + p_hi = 1.0 - sp_stats.norm.cdf(z) + p_lo = sp_stats.norm.cdf(z) + elif obs > mu: + p_hi, p_lo = 0.0, 1.0 + elif obs < mu: + p_hi, p_lo = 1.0, 0.0 else: - pval = 0.0 if obs > mu else 1.0 + p_hi, p_lo = 1.0, 1.0 - pval = float(max(min(pval, 1.0), 0.0)) - backbone.add_edge(agents[i], agents[j], sdsm_pvalue=pval) + if signed: + pval, sign = _combine_tails(p_hi, p_lo) + backbone.add_edge(agents[i], agents[j], sdsm_pvalue=pval, sign=sign) + else: + backbone.add_edge( + agents[i], agents[j], + sdsm_pvalue=float(max(min(p_hi, 1.0), 0.0)), + ) return _apply_projection_weights( backbone, @@ -1006,6 +1108,7 @@ def fdsm( alpha=0.05, trials=1000, seed=None, + signed=False, projection="simple", projection_weight="weight", projection_directed=False, @@ -1034,6 +1137,10 @@ def fdsm( Number of Monte Carlo randomisations. seed : integer, random_state, or None (default) Random seed for reproducibility. + signed : bool, optional (default=False) + If ``True``, run a two-tailed test: the stored p-value becomes two-sided + and a ``"sign"`` edge attribute marks significantly strong (``+1``) + versus significantly weak (``-1``) co-occurrences. projection : {"simple", "hyper", "probs", "ycn"}, optional Projection weighting assigned to each returned edge. projection_weight : str, optional @@ -1089,22 +1196,31 @@ def fdsm( row_sums = R.sum(axis=1) col_sums = R.sum(axis=0) - # Count how many times the random co-occurrence >= observed + # Count how often the random co-occurrence is >= observed (upper tail) and + # <= observed (lower tail, for signed backbones). exceed_count = np.zeros((na, na), dtype=int) + below_count = np.zeros((na, na), dtype=int) for _ in range(trials): R_rand = _random_bipartite_matrix(row_sums, col_sums, rng) co_rand = R_rand @ R_rand.T np.fill_diagonal(co_rand, 0) exceed_count += (co_rand >= observed).astype(int) + if signed: + below_count += (co_rand <= observed).astype(int) backbone = nx.Graph() backbone.add_nodes_from(agents) for i in range(na): for j in range(i + 1, na): - pval = exceed_count[i, j] / trials - backbone.add_edge(agents[i], agents[j], fdsm_pvalue=float(pval)) + p_hi = exceed_count[i, j] / trials + if signed: + p_lo = below_count[i, j] / trials + pval, sign = _combine_tails(p_hi, p_lo) + backbone.add_edge(agents[i], agents[j], fdsm_pvalue=pval, sign=sign) + else: + backbone.add_edge(agents[i], agents[j], fdsm_pvalue=float(p_hi)) return _apply_projection_weights( backbone, diff --git a/networkx_backbone/filters.py b/networkx_backbone/filters.py index 7c32051..c9c2598 100644 --- a/networkx_backbone/filters.py +++ b/networkx_backbone/filters.py @@ -16,9 +16,104 @@ "fraction_filter", "boolean_filter", "consensus_backbone", + "adjust_pvalues", ] +_MTC_METHODS = ("none", "bonferroni", "holm", "hochberg", "bh", "fdr", "by") + + +def adjust_pvalues(pvalues, method="bh"): + """Apply a multiple-testing correction to a collection of p-values. + + Reproduces R's ``p.adjust`` for the methods used by Neal's *Backbone 3.0* + ``mtc`` option, so statistical backbone p-values can be corrected for the + number of edges tested before thresholding. + + Parameters + ---------- + pvalues : iterable of float + The raw p-values. + method : string, optional (default="bh") + One of ``"none"``, ``"bonferroni"``, ``"holm"`` (step-down), + ``"hochberg"`` (step-up), ``"bh"``/``"fdr"`` (Benjamini-Hochberg), or + ``"by"`` (Benjamini-Yekutieli). + + Returns + ------- + adjusted : list of float + Adjusted p-values in the same order as the input, each in ``[0, 1]``. + + Raises + ------ + ValueError + If *method* is not recognised. + + References + ---------- + .. [1] Holm, S. (1979). A simple sequentially rejective multiple test + procedure. *Scandinavian Journal of Statistics*, 6(2), 65-70. + .. [2] Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple + tests of significance. *Biometrika*, 75(4), 800-802. + .. [3] Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery + rate. *J. Royal Statistical Society B*, 57(1), 289-300. + .. [4] Benjamini, Y., & Yekutieli, D. (2001). The control of the false + discovery rate in multiple testing under dependency. *Annals of + Statistics*, 29(4), 1165-1188. + + Examples + -------- + >>> from networkx_backbone import adjust_pvalues + >>> adjust_pvalues([0.01, 0.02, 0.5], method="bonferroni") + [0.03, 0.06, 1.0] + """ + method_l = method.lower() + if method_l not in _MTC_METHODS: + raise ValueError( + f"method must be one of {_MTC_METHODS}, got {method!r}" + ) + + p = [float(x) for x in pvalues] + n = len(p) + if n == 0: + return [] + if method_l == "none": + return [min(1.0, max(0.0, x)) for x in p] + if method_l == "bonferroni": + return [min(1.0, n * x) for x in p] + + if method_l == "holm": + # Step-down over ascending p-values with running maximum. + order = sorted(range(n), key=lambda k: p[k]) + out = [0.0] * n + running = 0.0 + for rank, idx in enumerate(order): + running = max(running, (n - rank) * p[idx]) + out[idx] = min(1.0, running) + return out + + # Step-up methods over descending p-values with running minimum. + order = sorted(range(n), key=lambda k: p[k], reverse=True) + out = [0.0] * n + running = float("inf") + if method_l == "hochberg": + for j, idx in enumerate(order): + running = min(running, (j + 1) * p[idx]) + out[idx] = min(1.0, running) + return out + if method_l in ("bh", "fdr"): + for j, idx in enumerate(order): + running = min(running, (n / (n - j)) * p[idx]) + out[idx] = min(1.0, running) + return out + # Benjamini-Yekutieli + c = sum(1.0 / k for k in range(1, n + 1)) + for j, idx in enumerate(order): + running = min(running, c * (n / (n - j)) * p[idx]) + out[idx] = min(1.0, running) + return out + + def multigraph_to_weighted(G, weight="weight", edge_type_attr=None): """Convert a MultiGraph/MultiDiGraph into a weighted simple graph. @@ -98,7 +193,8 @@ def multigraph_to_weighted(G, weight="weight", edge_type_attr=None): def threshold_filter( - G, score, threshold, mode="below", filter_on="edges", include_all_nodes=True + G, score, threshold, mode="below", filter_on="edges", include_all_nodes=True, + mtc="none", ): """Retain edges or nodes whose score passes a threshold test. @@ -124,6 +220,10 @@ def threshold_filter( - If ``filter_on="nodes"`` and ``True``, all retained nodes are kept even if isolated in the induced subgraph. If ``False``, retained nodes with degree 0 are removed. + mtc : string, optional (default="none") + Multiple-testing correction applied to *score* values before + thresholding (see :func:`adjust_pvalues`). Only valid with + ``mode="below"`` (p-values). Default ``"none"`` leaves scores unchanged. Returns ------- @@ -135,8 +235,9 @@ def threshold_filter( Raises ------ ValueError - If *mode* is not ``"below"`` or ``"above"``, or if *filter_on* is - not ``"edges"`` or ``"nodes"``. + If *mode* is not ``"below"`` or ``"above"``, if *filter_on* is + not ``"edges"`` or ``"nodes"``, or if *mtc* is used with + ``mode="above"``. Examples -------- @@ -151,31 +252,39 @@ def threshold_filter( """ if mode not in ("below", "above"): raise ValueError(f"mode must be 'below' or 'above', got {mode!r}") + if mtc != "none" and mode != "below": + raise ValueError("mtc correction requires mode='below' (p-values)") + + def _passes(val): + return (mode == "below" and val < threshold) or ( + mode == "above" and val >= threshold + ) if filter_on == "edges": H = G.__class__() if include_all_nodes: H.add_nodes_from(G.nodes(data=True)) - for u, v, data in G.edges(data=True): - val = data.get(score) - if val is None: - continue - if (mode == "below" and val < threshold) or ( - mode == "above" and val >= threshold - ): + scored = [ + (u, v, data) for u, v, data in G.edges(data=True) if data.get(score) is not None + ] + if mtc != "none": + values = adjust_pvalues([data[score] for _, _, data in scored], mtc) + else: + values = [data[score] for _, _, data in scored] + for (u, v, data), val in zip(scored, values): + if _passes(val): H.add_edge(u, v, **data) return H elif filter_on == "nodes": - keep = set() - for node, data in G.nodes(data=True): - val = data.get(score) - if val is None: - continue - if (mode == "below" and val < threshold) or ( - mode == "above" and val >= threshold - ): - keep.add(node) + scored = [ + (node, data) for node, data in G.nodes(data=True) if data.get(score) is not None + ] + if mtc != "none": + values = adjust_pvalues([data[score] for _, data in scored], mtc) + else: + values = [data[score] for _, data in scored] + keep = {node for (node, _), val in zip(scored, values) if _passes(val)} H = G.subgraph(keep).copy() if not include_all_nodes: isolates = list(nx.isolates(H)) diff --git a/networkx_backbone/hypergraph.py b/networkx_backbone/hypergraph.py new file mode 100644 index 0000000..bc5881d --- /dev/null +++ b/networkx_backbone/hypergraph.py @@ -0,0 +1,1130 @@ +""" +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). +hypergraph_compression_ratio + Inverse compression ratio ``eta`` achieved by the MDL backbone. +intersection_graph + Graph linking hyperedges that share at least one node. +""" + +import math +from dataclasses import dataclass, field + +import networkx as nx + +from networkx_backbone._docstrings import append_complexity_docstrings + +__all__ = [ + "intersection_graph", + "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 favours 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 + + +# --------------------------------------------------------------------------- +# Result container +# --------------------------------------------------------------------------- + + +@dataclass +class HypergraphBackbone: + """Result of :func:`mdl_hypergraph_backbone`. + + 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. + 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. + """ + + backbone: list = field(default_factory=list) + assignment: dict = 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 + + def __len__(self): + return len(self.backbone) + + @property + def n_input_hyperedges(self): + """Total hyperedges = backbone parents + all children.""" + return len(self.backbone) + sum(len(v) for v in self.assignment.values()) + + @property + def fraction_kept(self): + """Fraction of (unique) input hyperedges retained in the backbone.""" + total = self.n_input_hyperedges + return len(self.backbone) / total if total else 1.0 + + +@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 = field(default_factory=list) + pvalues: dict = field(default_factory=dict) + counts: dict = field(default_factory=dict) + alpha: float = 0.01 + method: str = "svh" + n_nodes: int = 0 + n_instances: int = 0 + + def __len__(self): + return len(self.validated) + + def __iter__(self): + return iter(self.validated) + + +# --------------------------------------------------------------------------- +# Public API +# --------------------------------------------------------------------------- + + +def intersection_graph(hyperedges, s=1, weight="overlap"): + """Build the intersection graph (or s-line graph) of a hypergraph. + + Each hyperedge becomes a node (labelled 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. + + Examples + -------- + >>> from networkx_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() + for i, e in enumerate(edges): + graph.add_node(i, members=e) + + 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 + + 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" optimiser: 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" optimiser: 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``. + """ + 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 + + while True: + chosen, best_gain = None, 1e-12 # require a strictly positive improvement + for e in range(n_edges): + if in_backbone[e]: + continue + 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 + if gain > best_gain: + best_gain, chosen = gain, e + if chosen is None: + break + in_backbone[chosen] = True + best_saving[chosen] = 0.0 + best_parent[chosen] = None + _attach_children(chosen) + + # 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]} + + +def _description_length( + edges, parent_of, weight_list, n_nodes, n_orders, weighted, gamma, mean_weight, prior +): + """Return ``(L(G, B), L(G, G))`` in bits for a given parent/child assignment.""" + dl = 0.0 + dl0 = 0.0 + for i, e in enumerate(edges): + size = len(e) + dl0 += _parent_codelength(size, n_orders, n_nodes) + if weighted: + dl0 += _weight_codelength(weight_list[i], 1, gamma, mean_weight, prior) + if i in parent_of: + p = parent_of[i] + overlap = len(e & edges[p]) + dl += _child_codelength(size, len(edges[p]), overlap, n_orders, n_nodes) + if weighted: + dl += _weight_codelength(weight_list[i], 0, gamma, mean_weight, prior) + else: + dl += _parent_codelength(size, n_orders, n_nodes) + if weighted: + dl += _weight_codelength(weight_list[i], 1, gamma, mean_weight, prior) + return dl, dl0 + + +def mdl_hypergraph_backbone( + hyperedges, + weights=None, + gamma=1.0, + prior="poisson", + method="auto", +): + """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 minimises the two-part description length + ``L(G, B) = L(B) + L(G|B)`` (Eqs (1)-(6)); equivalently it maximises 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 optimiser (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. + + Notes + ----- + Exact minimisation 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 networkx_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 + + all_nodes = set() + sizes = set() + for e in edges: + all_nodes.update(e) + sizes.add(len(e)) + n_nodes = len(all_nodes) + n_orders = len(sizes) + result.n_nodes = n_nodes + result.n_orders = n_orders + + mean_weight = sum(weight_list) / n_edges + weighted = weights is not None and any(w != 1.0 for w in weight_list) + result.weighted = weighted + + # Pairwise overlaps via node neighbourhoods (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 + + # Precompute per-edge weight terms (0 when unweighted). + if weighted: + wterm = [_weight_term(w, gamma, mean_weight, prior) for w in weight_list] + else: + 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, n_nodes) + neighbors[i].append((j, rmi)) + neighbors[j].append((i, rmi)) + + if method == "edge": + parent_of = _greedy_edge(n_edges, edges, neighbors, wterm) + elif method == "node": + parent_of = _greedy_node(n_edges, neighbors, wterm) + else: # "auto": run both and keep the lower description length. + parent_of = min( + ( + _greedy_edge(n_edges, edges, neighbors, wterm), + _greedy_node(n_edges, neighbors, wterm), + ), + key=lambda po: _description_length( + edges, po, weight_list, n_nodes, n_orders, + weighted, gamma, mean_weight, prior, + )[0], + ) + + dl, dl0 = _description_length( + edges, parent_of, weight_list, n_nodes, n_orders, + weighted, gamma, mean_weight, prior, + ) + + result.backbone = [edges[i] for i in range(n_edges) if i not in parent_of] + assignment = {} + 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 + + +def hypergraph_compression_ratio( + hyperedges, weights=None, gamma=1.0, prior="poisson" +): + """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]``. + + Examples + -------- + >>> from networkx_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 + + +def maximal_hyperedges(hyperedges): + """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. + + Examples + -------- + >>> from networkx_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(range(len(edges)), key=lambda i: (-len(edges[i]), sorted(edges[i]))) + kept = [] + for i in order: + e = edges[i] + if not any(e < bigger for bigger in kept): + kept.append(e) + return kept + + +def order_filter(hyperedges, min_order=None, max_order=None, orders=None): + """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``. + + Examples + -------- + >>> from networkx_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 + + +def s_components(hyperedges, s=1): + """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. + + Examples + -------- + >>> from networkx_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) + ] + components.sort( + key=lambda comp: (-len(comp), sorted(sorted(e) for e in comp)) + ) + return components + + +# --------------------------------------------------------------------------- +# 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 statistically_validated_hypergraph( + hyperedges, weights=None, max_order=None, alpha=0.01 +): + """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. + + Returns + ------- + result : ValidatedHypergraph + The validated hyperedges plus per-hyperedge p-values and counts. + + 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 networkx_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 + + 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 = {} + 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 = {} + order_nodes = set() + for e, w in members: + order_nodes.update(e) + for node in e: + degree[node] = degree.get(node, 0) + w + + groups = [e for e, _ in members] + pvals = [ + _svh_pvalue(w, n_instances, [degree[n] for n in e], binom) + for e, w in members + ] + 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 + + +def statistically_validated_cores( + hyperedges, weights=None, min_order=2, max_order=None, alpha=0.01 +): + """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. + max_order : int or None, optional (default=None) + Largest group size to test. ``None`` uses the largest hyperedge size. + alpha : float, optional (default=0.01) + FDR significance level. + + Returns + ------- + result : ValidatedHypergraph + The validated cores plus per-group p-values and co-occurrence counts. + + 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 networkx_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' + """ + from itertools import combinations + + from scipy.stats import binom + + edges, mult = _normalize_multiplicities(hyperedges, weights) + result = ValidatedHypergraph(alpha=alpha, method="svc") + if not edges: + return result + + n_instances = sum(mult) + degree = {} + 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) + + validated_groups = [] + for order in range(top, min_order - 1, -1): + drop = set() + for g in validated_groups: + if len(g) > order: + drop.update(frozenset(c) for c in combinations(tuple(g), order)) + + counts = {} + 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_pvalue(counts[g], n_instances, [degree[n] for n in g], binom) + for g in groups + ] + 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))." + ), + }, + "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) diff --git a/networkx_backbone/hypergraph_io.py b/networkx_backbone/hypergraph_io.py new file mode 100644 index 0000000..b139777 --- /dev/null +++ b/networkx_backbone/hypergraph_io.py @@ -0,0 +1,318 @@ +""" +Hypergraph ingestion and interoperability adapters. + +``networkx_backbone`` represents a hypergraph as a plain iterable of hyperedges +(each an iterable of node labels), which every hypergraph backbone method in +:mod:`networkx_backbone.hypergraph` consumes. This module converts that +representation to and from: + +- a **NetworkX incidence bipartite graph**, so bipartite projection backbones + (:func:`~networkx_backbone.sdsm`, :func:`~networkx_backbone.fdsm`, ...) can be + applied to hypergraphs; +- the **HIF** (Hypergraph Interchange Format) JSON standard; +- the **xgi**, **HyperNetX**, **HypergraphX**, and **Hypergraph Analysis + Toolbox (HAT)** hypergraph classes. + +The third-party libraries are imported lazily; none is a required dependency. +HIF is the recommended interchange path because all four libraries read and +write it, so it needs only the standard library. +""" + +import importlib +import json +import os + +import networkx as nx + +__all__ = [ + "hypergraph_to_bipartite", + "read_hif", + "write_hif", + "from_xgi", + "to_xgi", + "from_hypernetx", + "to_hypernetx", + "from_hypergraphx", + "to_hypergraphx", + "from_hat", + "to_hat", +] + + +def _as_hyperedges(hyperedges): + """Normalise to a list of non-empty :class:`frozenset` hyperedges (order kept).""" + out = [] + for e in hyperedges: + fs = frozenset(e) + if fs: + out.append(fs) + return out + + +def _require(module, pip_name=None): + """Import an optional dependency or raise a helpful ImportError.""" + try: + return importlib.import_module(module) + except ImportError as exc: # pragma: no cover - exercised when lib absent + raise ImportError( + f"'{module}' is required for this conversion; install it with " + f"`pip install {pip_name or module}`" + ) from exc + + +# --------------------------------------------------------------------------- +# NetworkX incidence bipartite graph +# --------------------------------------------------------------------------- + + +def hypergraph_to_bipartite(hyperedges, edge_prefix="he"): + """Convert a hypergraph to its incidence bipartite graph. + + Nodes of the hypergraph form one partition (``bipartite=0``); hyperedges form + the other (``bipartite=1``, labelled ``f"{edge_prefix}{i}"`` with a + ``"members"`` attribute). The result is exactly the input expected by the + bipartite projection backbones, so a hypergraph can be backboned with, e.g., + :func:`~networkx_backbone.sdsm` or :func:`~networkx_backbone.fdsm`. + + Parameters + ---------- + hyperedges : iterable of iterables + The hypergraph. + edge_prefix : str, optional (default="he") + Prefix for the generated hyperedge-node labels. + + Returns + ------- + B : networkx.Graph + The incidence bipartite graph. + nodes : list + The hypergraph's nodes (the agent partition), suitable to pass as + ``agent_nodes`` to the bipartite methods. + + Examples + -------- + >>> import networkx_backbone as nb + >>> B, nodes = nb.hypergraph_to_bipartite([(1, 2, 3), (2, 3, 4), (5, 6)]) + >>> import networkx as nx + >>> nx.is_bipartite(B) + True + >>> scored = nb.sdsm(B, agent_nodes=nodes) + """ + edges = _as_hyperedges(hyperedges) + nodes = sorted({v for e in edges for v in e}, key=repr) + + B = nx.Graph() + B.add_nodes_from(nodes, bipartite=0) + for i, e in enumerate(edges): + edge_node = f"{edge_prefix}{i}" + B.add_node(edge_node, bipartite=1, members=e) + for v in e: + B.add_edge(v, edge_node) + return B, nodes + + +# --------------------------------------------------------------------------- +# HIF (Hypergraph Interchange Format) +# --------------------------------------------------------------------------- + + +def write_hif(hyperedges, target=None, weights=None, network_type="undirected"): + """Serialise a hypergraph to the HIF (Hypergraph Interchange Format) standard. + + HIF is a JSON format shared by xgi, HyperNetX, HypergraphX, and HAT, making it + the most portable interchange path. + + Parameters + ---------- + hyperedges : iterable of iterables + The hypergraph. + target : str, path, file-like, or None, optional (default=None) + Where to write. If ``None``, the HIF dictionary is returned instead of + being written. + weights : iterable of numbers or None, optional (default=None) + Optional per-hyperedge weights, written to the HIF ``edges`` list. + network_type : str, optional (default="undirected") + Value for the HIF ``network-type`` field. + + Returns + ------- + result : dict or target + The HIF dictionary if *target* is ``None``; otherwise *target*. + """ + raw = list(hyperedges) + if weights is not None: + weights = list(weights) + if len(weights) != len(raw): + raise ValueError("weights must have the same length as hyperedges") + + incidences = [] + edge_records = [] + edge_id = 0 + for i, e in enumerate(raw): + members = frozenset(e) + if not members: + continue + for node in members: + incidences.append({"edge": edge_id, "node": node}) + if weights is not None: + edge_records.append({"edge": edge_id, "weight": weights[i]}) + edge_id += 1 + + hif = {"network-type": network_type, "metadata": {}, "incidences": incidences} + if edge_records: + hif["edges"] = edge_records + + if target is None: + return hif + if hasattr(target, "write"): + json.dump(hif, target) + return target + with open(target, "w") as handle: + json.dump(hif, handle) + return target + + +def read_hif(source, return_weights=False): + """Read a hypergraph from the HIF (Hypergraph Interchange Format) standard. + + Parameters + ---------- + source : dict, str, path, or file-like + A parsed HIF dictionary, a path to a HIF JSON file, or an open file. + return_weights : bool, optional (default=False) + If ``True``, also return per-hyperedge weights (defaulting to 1 for edges + without an explicit weight). + + Returns + ------- + hyperedges : list of frozenset + The hyperedges, in order of first appearance among the incidences. + weights : list of numbers + Returned only when *return_weights* is ``True``. + + Raises + ------ + TypeError + If *source* is not a dict, path, or file-like object. + """ + if isinstance(source, dict): + hif = source + elif isinstance(source, (str, os.PathLike)): + with open(source) as handle: + hif = json.load(handle) + elif hasattr(source, "read"): + hif = json.load(source) + else: + raise TypeError("source must be a dict, path, or file-like object") + + groups = {} + for incidence in hif.get("incidences", []): + groups.setdefault(incidence["edge"], set()).add(incidence["node"]) + + edge_ids = list(groups) + hyperedges = [frozenset(groups[eid]) for eid in edge_ids] + if not return_weights: + return hyperedges + + weight_map = {rec["edge"]: rec.get("weight", 1) for rec in hif.get("edges", [])} + weights = [weight_map.get(eid, 1) for eid in edge_ids] + return hyperedges, weights + + +# --------------------------------------------------------------------------- +# xgi +# --------------------------------------------------------------------------- + + +def from_xgi(H): + """Convert an :class:`xgi.Hypergraph` to a list of :class:`frozenset` hyperedges.""" + return [frozenset(members) for members in H.edges.members()] + + +def to_xgi(hyperedges): + """Convert hyperedges to an :class:`xgi.Hypergraph` (requires ``xgi``).""" + xgi = _require("xgi") + return xgi.Hypergraph([set(e) for e in _as_hyperedges(hyperedges)]) + + +# --------------------------------------------------------------------------- +# HyperNetX +# --------------------------------------------------------------------------- + + +def from_hypernetx(H): + """Convert a :class:`hypernetx.Hypergraph` to a list of :class:`frozenset` hyperedges.""" + return [frozenset(members) for members in H.incidence_dict.values()] + + +def to_hypernetx(hyperedges): + """Convert hyperedges to a :class:`hypernetx.Hypergraph` (requires ``hypernetx``).""" + hnx = _require("hypernetx") + return hnx.Hypergraph( + {i: set(e) for i, e in enumerate(_as_hyperedges(hyperedges))} + ) + + +# --------------------------------------------------------------------------- +# HypergraphX +# --------------------------------------------------------------------------- + + +def from_hypergraphx(H): + """Convert a ``hypergraphx.Hypergraph`` to a list of :class:`frozenset` hyperedges.""" + return [frozenset(e) for e in H.get_edges()] + + +def to_hypergraphx(hyperedges): + """Convert hyperedges to a ``hypergraphx.Hypergraph`` (requires ``hypergraphx``).""" + hgx = _require("hypergraphx") + return hgx.Hypergraph(edge_list=[tuple(e) for e in _as_hyperedges(hyperedges)]) + + +# --------------------------------------------------------------------------- +# Hypergraph Analysis Toolbox (HAT) +# --------------------------------------------------------------------------- + + +def _hat_incidence_matrix(H): + for attr in ("IM", "incidence_matrix", "incidenceMatrix"): + value = getattr(H, attr, None) + if value is None: + continue + return value() if callable(value) else value + raise AttributeError("could not find an incidence matrix on the HAT hypergraph") + + +def from_hat(H): + """Convert a HAT ``Hypergraph`` to a list of :class:`frozenset` hyperedges. + + Node labels are positional (row indices of the incidence matrix), as HAT is + incidence/tensor-based. + """ + np = _require("numpy") + matrix = np.asarray(_hat_incidence_matrix(H)) + hyperedges = [] + for col in range(matrix.shape[1]): + members = frozenset(int(r) for r in np.flatnonzero(matrix[:, col])) + if members: + hyperedges.append(members) + return hyperedges + + +def to_hat(hyperedges): + """Convert hyperedges to a HAT ``Hypergraph`` via an incidence matrix. + + Node labels are replaced by contiguous integer indices (HAT is positional). + Requires ``numpy`` and ``HAT``. + """ + np = _require("numpy") + hat = _require("HAT", pip_name="HypergraphAnalysisToolbox") + + edges = _as_hyperedges(hyperedges) + nodes = sorted({v for e in edges for v in e}, key=repr) + index = {v: i for i, v in enumerate(nodes)} + matrix = np.zeros((len(nodes), len(edges)), dtype=int) + for j, e in enumerate(edges): + for v in e: + matrix[index[v], j] = 1 + return hat.Hypergraph(incidence_matrix=matrix) diff --git a/networkx_backbone/statistical.py b/networkx_backbone/statistical.py index 183247c..212ba3e 100644 --- a/networkx_backbone/statistical.py +++ b/networkx_backbone/statistical.py @@ -56,12 +56,26 @@ def _validate_weights(G, weight): ) +def _combine_tails(p_hi, p_lo): + """Combine upper/lower-tail p-values into a two-sided p-value and a sign. + + Returns ``(two_sided_pvalue, sign)`` where ``sign`` is ``+1`` for a + significantly strong (over-represented) edge and ``-1`` for a significantly + weak (under-represented) edge. Used by the signed backbone variants. + """ + p_hi = min(max(p_hi, 0.0), 1.0) + p_lo = min(max(p_lo, 0.0), 1.0) + pval = min(1.0, 2.0 * min(p_hi, p_lo)) + sign = 1 if p_hi <= p_lo else -1 + return pval, sign + + # ===================================================================== # 1. Disparity filter -- Serrano et al. (2009) # ===================================================================== -def disparity_filter(G, weight="weight"): +def disparity_filter(G, weight="weight", signed=False): r"""Compute disparity filter p-values for each edge. The disparity filter [1]_ tests whether an edge's weight is @@ -84,12 +98,16 @@ def disparity_filter(G, weight="weight"): A NetworkX graph. weight : string, optional (default="weight") Edge attribute key for weights. All weights must be positive. + signed : bool, optional (default=False) + If ``True``, run a two-tailed test: the stored p-value becomes two-sided + and a ``"sign"`` edge attribute marks significantly strong (``+1``) + versus significantly weak (``-1``) edges. Returns ------- H : graph A copy of *G* (same type) with ``"disparity_pvalue"`` added as an - edge attribute. + edge attribute (and ``"sign"`` when ``signed=True``). Raises ------ @@ -125,26 +143,32 @@ def disparity_filter(G, weight="weight"): for u, v, data in H.edges(data=True): w = data[weight] if G.is_directed(): - pval = _disparity_node_pvalue(w, strength[u], degree[u]) + p_hi, p_lo = _disparity_node_tails(w, strength[u], degree[u]) else: - pval_u = _disparity_node_pvalue(w, strength[u], degree[u]) - pval_v = _disparity_node_pvalue(w, strength[v], degree[v]) - pval = min(pval_u, pval_v) - data["disparity_pvalue"] = pval + hi_u, lo_u = _disparity_node_tails(w, strength[u], degree[u]) + hi_v, lo_v = _disparity_node_tails(w, strength[v], degree[v]) + # OR rule: strong/weak if significant from either endpoint. + p_hi = min(hi_u, hi_v) + p_lo = min(lo_u, lo_v) + if signed: + data["disparity_pvalue"], data["sign"] = _combine_tails(p_hi, p_lo) + else: + data["disparity_pvalue"] = p_hi return H -def _disparity_node_pvalue(w, s, k): - """Disparity p-value from one node's perspective.""" +def _disparity_node_tails(w, s, k): + """Upper- and lower-tail disparity p-values from one node's perspective.""" if k <= 1: - return 1.0 + return 1.0, 1.0 p = min(w / s, 1.0) try: - alpha = (1.0 - p) ** (k - 1) + p_hi = (1.0 - p) ** (k - 1) except (OverflowError, ValueError): - alpha = 0.0 - return max(alpha, 0.0) + p_hi = 0.0 + p_hi = max(p_hi, 0.0) + return p_hi, max(1.0 - p_hi, 0.0) # ===================================================================== @@ -240,7 +264,7 @@ def noise_corrected_filter(G, weight="weight"): # ===================================================================== -def marginal_likelihood_filter(G, weight="weight"): +def marginal_likelihood_filter(G, weight="weight", signed=False): r"""Compute marginal likelihood p-values for each edge. The marginal likelihood filter [1]_ considers edge weights as @@ -253,6 +277,9 @@ def marginal_likelihood_filter(G, weight="weight"): A NetworkX graph. Integer weights are recommended. weight : string, optional (default="weight") Edge attribute key for weights. All weights must be positive. + signed : bool, optional (default=False) + If ``True``, run a two-tailed test and add a ``"sign"`` edge attribute + (``+1`` significantly strong, ``-1`` significantly weak). Returns ------- @@ -305,11 +332,16 @@ def marginal_likelihood_filter(G, weight="weight"): denom = W - su if denom > 0 and n_param > 0: p_param = min(sv / denom, 1.0) - pval = sp_stats.binom.sf(int(round(w)) - 1, n_param, p_param) + k = int(round(w)) + p_hi = sp_stats.binom.sf(k - 1, n_param, p_param) + p_lo = sp_stats.binom.cdf(k, n_param, p_param) else: - pval = 1.0 + p_hi, p_lo = 1.0, 1.0 - data["ml_pvalue"] = float(pval) + if signed: + data["ml_pvalue"], data["sign"] = _combine_tails(p_hi, p_lo) + else: + data["ml_pvalue"] = float(p_hi) return H @@ -451,7 +483,7 @@ def ecm_filter(G, weight="weight", max_iter=1000, tol=1e-6): # ===================================================================== -def lans_filter(G, weight="weight"): +def lans_filter(G, weight="weight", signed=False): r"""Compute LANS (Locally Adaptive Network Sparsification) p-values. LANS [1]_ is a nonparametric method that makes no distributional @@ -468,6 +500,9 @@ def lans_filter(G, weight="weight"): A NetworkX graph. weight : string, optional (default="weight") Edge attribute key for weights. All weights must be positive. + signed : bool, optional (default=False) + If ``True``, run a two-tailed test and add a ``"sign"`` edge attribute + (``+1`` significantly strong, ``-1`` significantly weak). Returns ------- @@ -515,13 +550,19 @@ def lans_filter(G, weight="weight"): if G.is_directed(): ecdf_u = _empirical_cdf(w, node_weights[u]) - pval = 1.0 - ecdf_u + p_hi = 1.0 - ecdf_u + p_lo = ecdf_u else: ecdf_u = _empirical_cdf(w, node_weights[u]) ecdf_v = _empirical_cdf(w, node_weights[v]) - pval = 1.0 - max(ecdf_u, ecdf_v) + # OR rule: strong/weak if significant from either endpoint. + p_hi = 1.0 - max(ecdf_u, ecdf_v) + p_lo = min(ecdf_u, ecdf_v) - data["lans_pvalue"] = max(pval, 0.0) + if signed: + data["lans_pvalue"], data["sign"] = _combine_tails(p_hi, p_lo) + else: + data["lans_pvalue"] = max(p_hi, 0.0) return H @@ -600,19 +641,19 @@ def _empirical_cdf(w, sorted_weights): return lo / n -def disparity(G, weight="weight"): +def disparity(G, weight="weight", signed=False): """Alias for :func:`disparity_filter`.""" - return disparity_filter(G, weight=weight) + return disparity_filter(G, weight=weight, signed=signed) -def mlf(G, weight="weight"): +def mlf(G, weight="weight", signed=False): """Alias for :func:`marginal_likelihood_filter`.""" - return marginal_likelihood_filter(G, weight=weight) + return marginal_likelihood_filter(G, weight=weight, signed=signed) -def lans(G, weight="weight"): +def lans(G, weight="weight", signed=False): """Alias for :func:`lans_filter`.""" - return lans_filter(G, weight=weight) + return lans_filter(G, weight=weight, signed=signed) _COMPLEXITY = { diff --git a/tests/test_bipartite.py b/tests/test_bipartite.py index 80cd3ee..8cfe7ee 100644 --- a/tests/test_bipartite.py +++ b/tests/test_bipartite.py @@ -307,3 +307,91 @@ def test_backbone_dispatch(davis_southern_women_graph, davis_women_nodes, weight with pytest.raises(ValueError): backbone(weighted_triangle, method="not_a_method") + + +class TestSignedProjections: + def test_sdsm_signed_adds_sign(self, davis_southern_women_graph, davis_women_nodes): + H = sdsm(davis_southern_women_graph, agent_nodes=davis_women_nodes, signed=True) + assert all("sign" in d and d["sign"] in (-1, 1) for _, _, d in H.edges(data=True)) + assert {d["sign"] for _, _, d in H.edges(data=True)} == {-1, 1} + assert all(0.0 <= d["sdsm_pvalue"] <= 1.0 for _, _, d in H.edges(data=True)) + + def test_sdsm_unsigned_has_no_sign(self, davis_southern_women_graph, davis_women_nodes): + H = sdsm(davis_southern_women_graph, agent_nodes=davis_women_nodes) + assert all("sign" not in d for _, _, d in H.edges(data=True)) + + def test_fdsm_signed_deterministic(self, davis_southern_women_graph, davis_women_nodes): + kw = dict(agent_nodes=davis_women_nodes, trials=200, seed=7, signed=True) + a = fdsm(davis_southern_women_graph, **kw) + b = fdsm(davis_southern_women_graph, **kw) + assert all("sign" in d for _, _, d in a.edges(data=True)) + assert {(u, v): d["sign"] for u, v, d in a.edges(data=True)} == { + (u, v): d["sign"] for u, v, d in b.edges(data=True) + } + + @pytest.mark.parametrize("model", [fixedrow, fixedcol, fixedfill]) + def test_fixed_models_signed(self, model, davis_southern_women_graph, davis_women_nodes): + H = model(davis_southern_women_graph, davis_women_nodes, alpha=0.3, signed=True) + assert H.number_of_edges() > 0 + assert all("sign" in d and d["sign"] in (-1, 1) for _, _, d in H.edges(data=True)) + + def test_backbone_from_projection_signed_passthrough( + self, davis_southern_women_graph, davis_women_nodes + ): + H = backbone_from_projection( + davis_southern_women_graph, + davis_women_nodes, + method="fixedrow", + alpha=0.3, + signed=True, + ) + assert all("sign" in d for _, _, d in H.edges(data=True)) + + +class TestSDSMEdgeConstraints: + def test_constraints_change_null(self, davis_southern_women_graph, davis_women_nodes): + B = davis_southern_women_graph + women = davis_women_nodes + artifact = next(n for n, d in B.nodes(data=True) if d["bipartite"] == 1) + cell = (women[0], artifact) + + base = sdsm(B, agent_nodes=women) + prohibited = sdsm(B, agent_nodes=women, prohibited=[cell]) + required = sdsm(B, agent_nodes=women, required=[cell]) + + diff_pro = sum( + 1 + for u, v, d in prohibited.edges(data=True) + if abs(d["sdsm_pvalue"] - base[u][v]["sdsm_pvalue"]) > 1e-9 + ) + diff_req = sum( + 1 + for u, v, d in required.edges(data=True) + if abs(d["sdsm_pvalue"] - base[u][v]["sdsm_pvalue"]) > 1e-9 + ) + assert diff_pro > 0 and diff_req > 0 + + def test_unknown_constraint_nodes_ignored( + self, davis_southern_women_graph, davis_women_nodes + ): + B = davis_southern_women_graph + base = sdsm(B, agent_nodes=davis_women_nodes) + constrained = sdsm( + B, agent_nodes=davis_women_nodes, prohibited=[("missing", "absent")] + ) + assert base.number_of_edges() == constrained.number_of_edges() + + def test_constraints_compose_with_signed( + self, davis_southern_women_graph, davis_women_nodes + ): + artifact = next( + n for n, d in davis_southern_women_graph.nodes(data=True) + if d["bipartite"] == 1 + ) + H = sdsm( + davis_southern_women_graph, + agent_nodes=davis_women_nodes, + signed=True, + required=[(davis_women_nodes[0], artifact)], + ) + assert all("sign" in d for _, _, d in H.edges(data=True)) diff --git a/tests/test_filters.py b/tests/test_filters.py index e7af27b..38889bd 100644 --- a/tests/test_filters.py +++ b/tests/test_filters.py @@ -4,6 +4,7 @@ import pytest from networkx_backbone import ( + adjust_pvalues, boolean_filter, consensus_backbone, disparity_filter, @@ -13,6 +14,74 @@ ) +class TestAdjustPvalues: + PS = [0.001, 0.008, 0.012, 0.03, 0.04, 0.2, 0.5, 0.9] + + def test_bonferroni(self): + assert adjust_pvalues([0.01, 0.02, 0.5], "bonferroni") == [0.03, 0.06, 1.0] + + def test_holm_step_down(self): + assert adjust_pvalues([0.01, 0.02, 0.5], "holm") == pytest.approx( + [0.03, 0.04, 0.5] + ) + + def test_none_clips_to_unit_interval(self): + assert adjust_pvalues([0.2, 1.5, -0.1], "none") == [0.2, 1.0, 0.0] + + def test_order_preserved_and_bounded(self): + adj = adjust_pvalues(self.PS, "bh") + assert len(adj) == len(self.PS) + assert all(0.0 <= a <= 1.0 for a in adj) + + @pytest.mark.parametrize("method", ["bh", "fdr", "by"]) + def test_matches_scipy(self, method): + sp = pytest.importorskip("scipy.stats") + import numpy as np + + expected = sp.false_discovery_control( + self.PS, method="by" if method == "by" else "bh" + ) + assert np.allclose(adjust_pvalues(self.PS, method), expected) + + def test_correction_is_conservative(self): + raw = self.PS + for method in ("bonferroni", "holm", "hochberg", "bh", "by"): + adj = adjust_pvalues(raw, method) + assert all(a >= r - 1e-12 for a, r in zip(adj, raw)) + + def test_empty(self): + assert adjust_pvalues([], "bh") == [] + + def test_invalid_method_raises(self): + with pytest.raises(ValueError): + adjust_pvalues([0.1], "hommel") + + +class TestThresholdFilterMTC: + def test_mtc_reduces_or_equal_edges(self): + G = nx.les_miserables_graph() + H = disparity_filter(G) + raw = threshold_filter(H, "disparity_pvalue", 0.05).number_of_edges() + for method in ("bonferroni", "bh"): + corrected = threshold_filter( + H, "disparity_pvalue", 0.05, mtc=method + ).number_of_edges() + assert corrected <= raw + + def test_mtc_requires_below_mode(self): + G = nx.les_miserables_graph() + H = disparity_filter(G) + with pytest.raises(ValueError): + threshold_filter(H, "disparity_pvalue", 0.05, mode="above", mtc="bh") + + def test_mtc_none_matches_plain(self): + G = nx.les_miserables_graph() + H = disparity_filter(G) + a = threshold_filter(H, "disparity_pvalue", 0.3).number_of_edges() + b = threshold_filter(H, "disparity_pvalue", 0.3, mtc="none").number_of_edges() + assert a == b + + class TestMultigraphToWeighted: def test_multigraph_parallel_edge_count(self): G = nx.MultiGraph() diff --git a/tests/test_hypergraph.py b/tests/test_hypergraph.py new file mode 100644 index 0000000..fdd6f81 --- /dev/null +++ b/tests/test_hypergraph.py @@ -0,0 +1,474 @@ +"""Tests for hypergraph backbone extraction (MDL method, Kirkley et al. 2026).""" + +import itertools + +import pytest + +from networkx_backbone import ( + HypergraphBackbone, + ValidatedHypergraph, + hypergraph_compression_ratio, + intersection_graph, + maximal_hyperedges, + mdl_hypergraph_backbone, + order_filter, + s_components, + statistically_validated_cores, + statistically_validated_hypergraph, + svc, + svh, +) +from networkx_backbone.hypergraph import ( + _child_codelength, + _parent_codelength, + _reduced_mutual_information, +) + + +# --------------------------------------------------------------------------- +# Fixtures / helpers +# --------------------------------------------------------------------------- + + +def simplex_with_subfaces(nodes, min_order=4): + """A top face plus all of its nested sub-hyperedges down to ``min_order``.""" + nodes = list(nodes) + edges = [tuple(nodes)] + for k in range(min_order, len(nodes)): + edges.extend(itertools.combinations(nodes, k)) + return edges + + +# --------------------------------------------------------------------------- +# Information-theoretic primitives +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "size_p,size_c,overlap", + [(5, 4, 4), (5, 3, 3), (5, 2, 2), (4, 4, 3), (6, 3, 2)], +) +def test_rmi_symmetric_and_identity(size_p, size_c, overlap): + n_nodes = 12 + r_pc = _reduced_mutual_information(size_p, size_c, overlap, n_nodes) + r_cp = _reduced_mutual_information(size_c, size_p, overlap, n_nodes) + assert r_pc == pytest.approx(r_cp) # R is symmetric (Eqs. 9-10) + + # R(c, p) == H(c) - H(c|p) (Eq. 9), using L=1 so the log L term cancels. + h_c = _parent_codelength(size_c, 1, n_nodes) + h_c_given_p = _child_codelength(size_c, size_p, overlap, 1, n_nodes) + assert r_pc == pytest.approx(h_c - h_c_given_p) + + +def test_rmi_positive_for_nesting_when_universe_large(): + # Nested c ⊂ p gives positive RMI once N is sufficiently larger than |p|. + assert _reduced_mutual_information(5, 4, 4, 20) > 0 + # ...and can be non-positive when the universe is tiny. + assert _reduced_mutual_information(5, 4, 4, 5) <= 0 + + +# --------------------------------------------------------------------------- +# intersection_graph +# --------------------------------------------------------------------------- + + +def test_intersection_graph_basic(): + I = intersection_graph([(1, 2, 3), (2, 3, 4), (5, 6)]) + assert I.number_of_nodes() == 3 + assert I.number_of_edges() == 1 # only the first two overlap + assert I[0][1]["overlap"] == 2 + assert I.nodes[0]["members"] == frozenset({1, 2, 3}) + + +def test_intersection_graph_disjoint_has_no_edges(): + I = intersection_graph([(0, 1), (2, 3), (4, 5)]) + assert I.number_of_nodes() == 3 + assert I.number_of_edges() == 0 + + +def test_intersection_graph_s_threshold(): + G = [(1, 2, 3), (2, 3, 4), (5, 6)] + assert intersection_graph(G, s=1).number_of_edges() == 1 + assert intersection_graph(G, s=2).number_of_edges() == 1 # overlap is exactly 2 + assert intersection_graph(G, s=3).number_of_edges() == 0 + + +def test_intersection_graph_invalid_s_raises(): + with pytest.raises(ValueError): + intersection_graph([(1, 2, 3)], s=0) + + +# --------------------------------------------------------------------------- +# Structural methods: inclusion reduction, order filter, s-components +# --------------------------------------------------------------------------- + + +def test_maximal_hyperedges_removes_subsets(): + result = maximal_hyperedges([(1, 2, 3), (1, 2), (2, 3), (4, 5), (4, 5, 6)]) + assert set(result) == {frozenset({1, 2, 3}), frozenset({4, 5, 6})} + # Ordered by decreasing size. + assert [len(e) for e in result] == sorted((len(e) for e in result), reverse=True) + + +def test_maximal_hyperedges_nested_chain(): + assert maximal_hyperedges([(1,), (1, 2), (1, 2, 3)]) == [frozenset({1, 2, 3})] + + +def test_maximal_hyperedges_all_maximal(): + result = maximal_hyperedges([(1, 2), (3, 4), (5, 6)]) + assert set(result) == {frozenset({1, 2}), frozenset({3, 4}), frozenset({5, 6})} + + +def test_maximal_hyperedges_merges_duplicates(): + assert maximal_hyperedges([(1, 2, 3), (3, 2, 1)]) == [frozenset({1, 2, 3})] + + +def test_order_filter_min_max_and_orders(): + G = [(1, 2), (1, 2, 3), (1, 2, 3, 4)] + assert set(order_filter(G, min_order=3)) == { + frozenset({1, 2, 3}), + frozenset({1, 2, 3, 4}), + } + assert order_filter(G, max_order=2) == [frozenset({1, 2})] + assert order_filter(G, min_order=3, max_order=3) == [frozenset({1, 2, 3})] + assert set(order_filter(G, orders=[2, 4])) == { + frozenset({1, 2}), + frozenset({1, 2, 3, 4}), + } + + +def test_order_filter_invalid_range_raises(): + with pytest.raises(ValueError): + order_filter([(1, 2, 3)], min_order=4, max_order=2) + + +def test_s_components_threshold(): + G = [(1, 2, 3), (2, 3, 4), (5, 6, 7)] + # s=1 and s=2: first two hyperedges connect; the third is isolated. + for s in (1, 2): + comps = s_components(G, s=s) + assert [len(c) for c in comps] == [2, 1] + # s=3: no pair shares 3 nodes, so every hyperedge is its own component. + assert [len(c) for c in s_components(G, s=3)] == [1, 1, 1] + + +def test_s_components_chain(): + comps = s_components([(1, 2), (2, 3), (3, 4), (10, 11)], s=1) + assert [len(c) for c in comps] == [3, 1] + # The big component holds the chain; the disjoint pair stands alone. + assert frozenset({10, 11}) in comps[1] + + +# --------------------------------------------------------------------------- +# Core MDL backbone behaviour +# --------------------------------------------------------------------------- + + +def test_recovers_top_faces_of_nested_simplices(): + # Two disjoint 5-simplices, each with all size-4 sub-faces. The MDL backbone + # should keep exactly the two top faces and prune every redundant sub-face. + a, b = range(0, 5), range(5, 10) + G = simplex_with_subfaces(a) + simplex_with_subfaces(b) + + result = mdl_hypergraph_backbone(G) + + assert isinstance(result, HypergraphBackbone) + assert len(result.backbone) == 2 + assert frozenset(a) in result.backbone + assert frozenset(b) in result.backbone + assert result.compression_ratio < 1.0 + assert result.n_nodes == 10 + # Every pruned sub-face is recorded as a child of a retained parent. + assert result.n_input_hyperedges == len(G) + + +def test_downward_closure_compresses_and_keeps_tops(): + a, b = range(0, 5), range(5, 10) + G = simplex_with_subfaces(a, min_order=2) + simplex_with_subfaces(b, min_order=2) + + result = mdl_hypergraph_backbone(G) + + assert frozenset(a) in result.backbone + assert frozenset(b) in result.backbone + assert result.fraction_kept < 0.5 # substantial sparsification + assert 0.0 <= result.compression_ratio < 1.0 + + +def test_disjoint_hyperedges_are_incompressible(): + G = [(0, 1), (2, 3), (4, 5)] + result = mdl_hypergraph_backbone(G) + assert len(result.backbone) == 3 + assert result.assignment == {} + assert result.compression_ratio == pytest.approx(1.0) + + +def test_description_length_never_exceeds_baseline(): + G = simplex_with_subfaces(range(0, 6), min_order=2) + result = mdl_hypergraph_backbone(G) + assert result.description_length <= result.baseline_description_length + assert 0.0 <= result.compression_ratio <= 1.0 + + +def test_compression_ratio_helper_matches_backbone(): + G = simplex_with_subfaces(range(0, 5)) + simplex_with_subfaces(range(5, 10)) + eta = hypergraph_compression_ratio(G) + assert eta == pytest.approx(mdl_hypergraph_backbone(G).compression_ratio) + assert 0.0 <= eta <= 1.0 + + +def test_deterministic(): + G = simplex_with_subfaces(range(0, 5), min_order=2) + simplex_with_subfaces( + range(5, 10), min_order=2 + ) + r1 = mdl_hypergraph_backbone(G) + r2 = mdl_hypergraph_backbone(G) + assert set(r1.backbone) == set(r2.backbone) + + +# --------------------------------------------------------------------------- +# Weighted model +# --------------------------------------------------------------------------- + + +def _family_a_with_filler(): + """Family A (top + one size-4 child) plus a disjoint filler simplex. + + The filler raises the node count so the size-4 child of A is topologically + redundant (pruned) in the unweighted backbone. + """ + A = [(0, 1, 2, 3, 4), (0, 1, 2, 3)] + filler = simplex_with_subfaces(range(5, 10)) + return A, filler + + +def test_gamma_one_equivalent_to_unweighted(): + A, filler = _family_a_with_filler() + G = A + filler + weights = [5.0, 5.0] + [5.0] * len(filler) # uniform weights + + unweighted = mdl_hypergraph_backbone(G) + g1 = mdl_hypergraph_backbone(G, weights=weights, gamma=1.0) + + assert g1.weighted is True + assert set(g1.backbone) == set(unweighted.backbone) + + +def test_low_gamma_forces_high_weight_hyperedge_into_backbone(): + A, filler = _family_a_with_filler() + G = A + filler + weights = [1.0, 100.0] + [1.0] * len(filler) # the size-4 child is heavy + + unweighted = mdl_hypergraph_backbone(G) + weighted = mdl_hypergraph_backbone(G, weights=weights, gamma=0.01) + + child = frozenset({0, 1, 2, 3}) + assert child not in unweighted.backbone # topologically redundant + assert child in weighted.backbone # weight keeps it in the backbone + + +def test_geometric_prior_runs(): + A, filler = _family_a_with_filler() + G = A + filler + weights = [1.0, 100.0] + [1.0] * len(filler) + result = mdl_hypergraph_backbone(G, weights=weights, gamma=0.05, prior="geometric") + assert frozenset({0, 1, 2, 3}) in result.backbone + + +# --------------------------------------------------------------------------- +# Edge cases and input handling +# --------------------------------------------------------------------------- + + +def test_empty_hypergraph(): + result = mdl_hypergraph_backbone([]) + assert len(result.backbone) == 0 + assert result.compression_ratio == pytest.approx(1.0) + + +def test_single_hyperedge(): + result = mdl_hypergraph_backbone([(1, 2, 3)]) + assert result.backbone == [frozenset({1, 2, 3})] + assert result.assignment == {} + assert result.compression_ratio == pytest.approx(1.0) + + +def test_duplicate_hyperedges_are_merged(): + result = mdl_hypergraph_backbone([(1, 2, 3), (1, 2, 3), (3, 2, 1)]) + assert result.backbone == [frozenset({1, 2, 3})] + assert result.n_input_hyperedges == 1 + + +def test_repeated_nodes_within_hyperedge_ignored(): + result = mdl_hypergraph_backbone([(1, 1, 2, 2, 3)]) + assert result.backbone == [frozenset({1, 2, 3})] + + +def test_duplicate_weighted_hyperedges_sum_weights(): + # Two copies of the same hyperedge merge; their weights add (3 + 4 = 7 > 1), + # so the merged hyperedge is treated as weighted. + result = mdl_hypergraph_backbone( + [(1, 2, 3), (1, 2, 3)], weights=[3.0, 4.0], gamma=0.5 + ) + assert result.backbone == [frozenset({1, 2, 3})] + assert result.weighted is True + + +# --------------------------------------------------------------------------- +# Validation +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("gamma", [0.0, -0.1, 1.5]) +def test_invalid_gamma_raises(gamma): + with pytest.raises(ValueError): + mdl_hypergraph_backbone([(1, 2, 3)], weights=[2.0], gamma=gamma) + + +def test_invalid_prior_raises(): + with pytest.raises(ValueError): + mdl_hypergraph_backbone([(1, 2, 3)], weights=[2.0], prior="bogus") + + +def test_invalid_method_raises(): + with pytest.raises(ValueError): + mdl_hypergraph_backbone([(1, 2, 3)], method="bogus") + + +@pytest.mark.parametrize("method", ["edge", "node", "auto"]) +def test_mdl_methods_recover_nested_top_faces(method): + a, b = range(0, 5), range(5, 10) + G = simplex_with_subfaces(a) + simplex_with_subfaces(b) + result = mdl_hypergraph_backbone(G, method=method) + assert frozenset(a) in result.backbone + assert frozenset(b) in result.backbone + assert result.description_length <= result.baseline_description_length + + +def test_mdl_auto_is_best_of_edge_and_node(): + # On a downward-closed simplex the two greedy schemes differ; "auto" must + # match the lower-description-length one. + G = simplex_with_subfaces(range(0, 6), min_order=2) + edge = mdl_hypergraph_backbone(G, method="edge").description_length + node = mdl_hypergraph_backbone(G, method="node").description_length + auto = mdl_hypergraph_backbone(G, method="auto").description_length + assert auto == pytest.approx(min(edge, node)) + + +def test_mdl_node_weighting_keeps_high_weight_child(): + # The gamma weight knob must work under the node optimiser too. + A = [(0, 1, 2, 3, 4), (0, 1, 2, 3)] + filler = simplex_with_subfaces(range(5, 10)) + G = A + filler + weights = [1.0, 100.0] + [1.0] * len(filler) + result = mdl_hypergraph_backbone(G, weights=weights, gamma=0.01, method="node") + assert frozenset({0, 1, 2, 3}) in result.backbone + + +def test_weights_length_mismatch_raises(): + with pytest.raises(ValueError): + mdl_hypergraph_backbone([(1, 2, 3), (2, 3, 4)], weights=[1.0]) + + +def test_weight_below_one_raises(): + with pytest.raises(ValueError): + mdl_hypergraph_backbone([(1, 2, 3)], weights=[0.5]) + + +# --------------------------------------------------------------------------- +# Statistically validated hypergraphs (SVH) and cores (SVC) +# --------------------------------------------------------------------------- + + +def test_svh_pvalue_matches_binomial(): + from scipy.stats import binom + + edges = [(1, 2)] * 5 + [(1, 3), (2, 4), (3, 4), (5, 6)] + result = statistically_validated_hypergraph(edges) + # order-2 instances N=9; node 1 and 2 each have degree 6; {1,2} co-occurs 5x. + expected = binom.sf(5 - 1, 9, (6 / 9) * (6 / 9)) + assert result.pvalues[frozenset({1, 2})] == pytest.approx(expected) + + +def test_svh_validates_over_represented_pair(): + # {1,2} recurs 5x; a high-multiplicity background pair inflates the instance + # count so the planted pair is significant while the background is not. + edges = [(1, 2)] * 5 + [(3, 4)] * 100 + result = statistically_validated_hypergraph(edges, alpha=0.05) + assert isinstance(result, ValidatedHypergraph) + assert frozenset({1, 2}) in result.validated + assert frozenset({3, 4}) not in result.validated # expected background + assert result.method == "svh" + + +def test_svh_weights_equal_repeated_edges(): + repeated = statistically_validated_hypergraph([(1, 2)] * 5 + [(3, 4)] * 100) + weighted = statistically_validated_hypergraph( + [(1, 2), (3, 4)], weights=[5, 100] + ) + assert weighted.pvalues[frozenset({1, 2})] == pytest.approx( + repeated.pvalues[frozenset({1, 2})] + ) + assert set(weighted.validated) == set(repeated.validated) + + +def test_svh_alpha_monotonic(): + edges = [(1, 2)] * 5 + [(3, 4)] * 100 + strict = set(svh(edges, alpha=0.001).validated) + loose = set(svh(edges, alpha=0.5).validated) + assert strict <= loose + + +def test_svh_single_occurrences_not_validated(): + result = statistically_validated_hypergraph([(1, 2), (3, 4), (5, 6)], alpha=0.05) + assert result.validated == [] + + +def test_svh_max_order_restricts_tests(): + edges = [(1, 2)] * 5 + [(3, 4)] * 100 + [(5, 6, 7)] * 5 + result = statistically_validated_hypergraph(edges, max_order=2) + assert all(len(e) == 2 for e in result.pvalues) + + +def test_svh_empty(): + result = statistically_validated_hypergraph([]) + assert len(result) == 0 + assert list(result) == [] + + +def test_svc_validates_core_and_drops_subgroups(): + edges = [(1, 2, 3)] * 5 + [(7, 8)] * 200 + result = statistically_validated_cores(edges, alpha=0.05) + assert result.method == "svc" + assert frozenset({1, 2, 3}) in result.validated + # Sub-pairs of a validated core are not separately validated. + assert frozenset({1, 2}) not in result.validated + assert frozenset({2, 3}) not in result.validated + + +def test_svc_counts_subgroup_cooccurrence(): + # {1,2} appears inside both a triangle and a separate pair: co-occurrence 2. + edges = [(1, 2, 3), (1, 2)] + result = statistically_validated_cores(edges, min_order=2) + assert result.counts[frozenset({1, 2})] == 2 + + +@pytest.mark.parametrize("func", [statistically_validated_hypergraph, statistically_validated_cores]) +def test_svh_svc_integer_multiplicity_required(func): + with pytest.raises(ValueError): + func([(1, 2), (3, 4)], weights=[1.5, 2]) + + +@pytest.mark.parametrize("func", [statistically_validated_hypergraph, statistically_validated_cores]) +def test_svh_svc_weights_length_checked(func): + with pytest.raises(ValueError): + func([(1, 2), (3, 4)], weights=[1]) + + +def test_svh_svc_aliases(): + edges = [(1, 2)] * 5 + [(3, 4)] * 100 + assert set(svh(edges).validated) == set( + statistically_validated_hypergraph(edges).validated + ) + assert set(svc(edges).validated) == set( + statistically_validated_cores(edges).validated + ) diff --git a/tests/test_hypergraph_io.py b/tests/test_hypergraph_io.py new file mode 100644 index 0000000..e1ecec5 --- /dev/null +++ b/tests/test_hypergraph_io.py @@ -0,0 +1,114 @@ +"""Tests for hypergraph ingestion / interoperability adapters.""" + +import importlib + +import networkx as nx +import pytest + +import networkx_backbone as nb + +HYPEREDGES = [(1, 2, 3), (2, 3, 4), (5, 6)] +EXPECTED = {frozenset(e) for e in HYPEREDGES} + + +# --------------------------------------------------------------------------- +# Incidence bipartite graph +# --------------------------------------------------------------------------- + + +def test_hypergraph_to_bipartite_structure(): + B, nodes = nb.hypergraph_to_bipartite(HYPEREDGES) + assert nx.is_bipartite(B) + assert nodes == [1, 2, 3, 4, 5, 6] + node_part = {n for n, d in B.nodes(data=True) if d["bipartite"] == 0} + edge_part = {n for n, d in B.nodes(data=True) if d["bipartite"] == 1} + assert node_part == set(nodes) + assert len(edge_part) == 3 + # Each hyperedge node records its members. + assert any(B.nodes[e]["members"] == frozenset({1, 2, 3}) for e in edge_part) + + +def test_hypergraph_to_bipartite_feeds_sdsm(): + B, nodes = nb.hypergraph_to_bipartite(HYPEREDGES) + scored = nb.sdsm(B, agent_nodes=nodes) + assert all("sdsm_pvalue" in d for _, _, d in scored.edges(data=True)) + + +# --------------------------------------------------------------------------- +# HIF +# --------------------------------------------------------------------------- + + +def test_hif_roundtrip_dict(): + hif = nb.write_hif(HYPEREDGES) + assert hif["network-type"] == "undirected" + assert set(nb.read_hif(hif)) == EXPECTED + + +def test_hif_roundtrip_with_weights(): + hif = nb.write_hif([(1, 2), (3, 4)], weights=[5, 2]) + edges, weights = nb.read_hif(hif, return_weights=True) + assert dict(zip((tuple(sorted(e)) for e in edges), weights)) == { + (1, 2): 5, + (3, 4): 2, + } + + +def test_hif_roundtrip_file(tmp_path): + path = tmp_path / "h.json" + nb.write_hif(HYPEREDGES, path) + assert set(nb.read_hif(path)) == EXPECTED + + +def test_write_hif_weight_length_checked(): + with pytest.raises(ValueError): + nb.write_hif([(1, 2), (3, 4)], weights=[1]) + + +def test_read_hif_bad_source_raises(): + with pytest.raises(TypeError): + nb.read_hif(12345) + + +# --------------------------------------------------------------------------- +# Cross-library interoperability (run when a library is available) +# --------------------------------------------------------------------------- + + +def test_hif_crosscheck_with_xgi(tmp_path): + xgi = pytest.importorskip("xgi") + path = tmp_path / "h.json" + # xgi reads HIF we wrote. + nb.write_hif(HYPEREDGES, path) + Hx = xgi.read_hif(path) + assert {frozenset(m) for m in Hx.edges.members()} == EXPECTED + # We read HIF xgi wrote. + xgi.write_hif(xgi.Hypergraph([list(e) for e in HYPEREDGES]), path) + assert set(nb.read_hif(path)) == EXPECTED + + +ADAPTERS = [ + ("xgi", nb.to_xgi, nb.from_xgi), + ("hypernetx", nb.to_hypernetx, nb.from_hypernetx), + ("hypergraphx", nb.to_hypergraphx, nb.from_hypergraphx), + ("HAT", nb.to_hat, nb.from_hat), +] + + +@pytest.mark.parametrize("module,to_fn,from_fn", ADAPTERS, ids=[a[0] for a in ADAPTERS]) +def test_adapter_roundtrip_or_importerror(module, to_fn, from_fn): + try: + importlib.import_module(module) + except ImportError: + # Absent (or broken) optional dependency raises a helpful ImportError. + with pytest.raises(ImportError, match=module): + to_fn(HYPEREDGES) + return + + obj = to_fn(HYPEREDGES) + recovered = from_fn(obj) + if module == "HAT": + # HAT is positional: labels become indices, so compare the order multiset. + assert sorted(len(e) for e in recovered) == sorted(len(set(e)) for e in HYPEREDGES) + else: + assert set(recovered) == EXPECTED diff --git a/tests/test_les_miserables.py b/tests/test_les_miserables.py index d63c1f6..cf1d6da 100644 --- a/tests/test_les_miserables.py +++ b/tests/test_les_miserables.py @@ -27,7 +27,7 @@ def _warn_if_no_edge_reduction(method_name, original_edges, filtered_edges): "disparity_filter", lambda G: nb.disparity_filter(G), lambda H: nb.threshold_filter(H, "disparity_pvalue", DEFAULT_PVALUE, mode="below"), - 247, + 9, ), ( "noise_corrected_filter", diff --git a/tests/test_statistical.py b/tests/test_statistical.py index c685329..e518028 100644 --- a/tests/test_statistical.py +++ b/tests/test_statistical.py @@ -137,3 +137,46 @@ def test_aliases_match_filter_outputs(self, weighted_triangle): assert "disparity_pvalue" in h_disparity[u][v] assert "lans_pvalue" in h_lans[u][v] assert "ml_pvalue" in h_mlf[u][v] + + +class TestSignedBackbones: + @staticmethod + def _hub_star(): + # Hub "c" with one very strong edge and several weak edges. + G = nx.Graph() + G.add_weighted_edges_from( + [("c", 1, 100.0), ("c", 2, 1.0), ("c", 3, 1.0), ("c", 4, 1.0)] + ) + return G + + def test_disparity_signed_marks_strong_and_weak(self): + G = self._hub_star() + H = disparity_filter(G, signed=True) + assert H["c"][1]["sign"] == 1 # strong edge + assert H["c"][2]["sign"] == -1 # weak edge + assert all(0.0 <= d["disparity_pvalue"] <= 1.0 for _, _, d in H.edges(data=True)) + + def test_unsigned_unchanged(self): + G = self._hub_star() + H = disparity_filter(G) + assert all("sign" not in d for _, _, d in H.edges(data=True)) + # one-tailed p-value preserved: (1 - w/s)^(k-1) for the strong edge + assert H["c"][1]["disparity_pvalue"] == pytest.approx((1 - 100 / 103) ** 3) + + @pytest.mark.parametrize( + "scorer,attr", + [(disparity_filter, "disparity_pvalue"), (mlf, "ml_pvalue"), (lans, "lans_pvalue")], + ) + def test_signed_adds_sign_and_two_sided_pvalue(self, scorer, attr): + G = nx.les_miserables_graph() + H = scorer(G, signed=True) + assert all("sign" in d and d["sign"] in (-1, 1) for _, _, d in H.edges(data=True)) + assert all(0.0 <= d[attr] <= 1.0 for _, _, d in H.edges(data=True)) + + def test_signed_pvalue_is_two_sided(self): + # Two-sided p-value is twice the smaller one-tailed tail (clipped to 1). + G = self._hub_star() + unsigned = disparity_filter(G) + signed = disparity_filter(G, signed=True) + p_hi = unsigned["c"][1]["disparity_pvalue"] + assert signed["c"][1]["disparity_pvalue"] == pytest.approx(min(1.0, 2 * p_hi))