Repository: xgi-org/xgi

Issue:
uniform_erdos_renyi_hypergraph produces way too many edges
Unless I'm seriously mistaken as to how this is supposed to work, the number of edges produced by `uniform_erdos_renyi_hypergraph` is off by a factor of $m!$.

```python
import xgi
import matplotlib.pyplot as plt

p = 1e-6
n = 1000
m = 3
plt.hist([len(xgi.uniform_erdos_renyi_hypergraph(n, m, p, p_type="prob").edges) for i in range(1000)],
    label="Num. of edges")
plt.axvline(math.comb(n, m)*p, label="Exp. number of edges")
plt.legend()
```
![download](https://user-images.githubusercontent.com/699404/233212915-e45b9b1b-026d-4f6b-9389-36a8b62aa351.png)

Possible solutions: lower `q`  by $m!$ or even better, use [nth combination](https://docs.python.org/3/library/itertools.html#itertools-recipes:~:text=def%20nth_combination) instead of n-th product when defining `_index_to_edge`.

Hints from issue comments:
@leotrs: Hi @arashbm thanks for your report. We'll look into this shortly!@maximelucas: Thanks @arashbm ! We'll have a look.

We also have the plan to move than function to `random` and rename it `fast_random_hypergraph` and have it deal with both uniform and non-uniform hypergraphs (see https://github.com/xgi-org/xgi/pull/321#issuecomment-1497974162).

In the meantime, you can instead use `xgi.random_hypergraph()` with the argument `order`, which does the same thing in a different way (see [docs](https://xgi.readthedocs.io/en/stable/api/generators/xgi.generators.random.html#xgi.generators.random.random_hypergraph)).@arashbm: @maximelucas The combinatorial explosion of the mask array makes it impossible to use in some scenarios: e.g.
```pycon
>>> import xgi
>>> xgi.random_hypergraph(1000, order=5, ps=[1e-12])  # => should create a network with ~1350 edges
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 9.72 PiB for an array with shape (1368173298991500,) and data type float64
```

But that also has an easy solution of not pre-allocating mask:

```python
import xgi
import numpy as np
import itertools

def random_hypergraph(N, ps, order=None, seed=None):
    rng = np.random.default_rng(seed)

    if order is not None:
        if len(ps) != 1:
            raise ValueError("ps must contain a single element if order is an int")

    if (np.any(np.array(ps) < 0)) or (np.any(np.array(ps) > 1)):
        raise ValueError("All elements of ps must be between 0 and 1 included.")

    nodes = range(N)
    hyperedges = []

    for i, p in enumerate(ps):

        if order is not None:
            d = order
        else:
            d = i + 1  # order, ps[0] is prob of edges (d=1)

        potential_edges = itertools.combinations(nodes, d + 1)
        edges_to_add = [e for e in potential_edges if rng.random() <= p]

        hyperedges += edges_to_add

    H = xgi.empty_hypergraph()
    H.add_nodes_from(nodes)
    H.add_edges_from(hyperedges)

    return H
```

This doesn't fail but this still takes $O({{N}\choose{m}})$ time to run. The original algorithm (with geometric "jumps") is still superior.@nwlandry: Hi @arashbm ! Thanks for raising this issue! I added this function to XGI with my most recent ArXiv paper (https://arxiv.org/abs/2302.13967, see appendix). The problem is because we are iterating over the Cartesian product, not the unique combinations. I did this because the index _to_edge function is much slower when getting a combination from an index instead on a product from an index. If the hypergraph is sparse enough, than it should be fine though. I will try to dig up the combinations version of index_to_edge and we can chat more!@arashbm: @nwlandry That makes total sense. My current WIP implementation is also [based on Cartesian product](https://github.com/reticula-network/reticula/blob/42ed11ebd3cccb16af2ad66dd803ba07d9940cb2/src/random_networks/uniform_hypergraph.tpp#L53-L58) as this help alleviate the issue with fixed sized integers, but I have to [reject](https://github.com/reticula-network/reticula/blob/42ed11ebd3cccb16af2ad66dd803ba07d9940cb2/src/random_networks/uniform_hypergraph.tpp#L63-L64) any candidate edge (node combinations) that are not strictly monotonically increasing order, which makes for a lot of wasted CPU cycles.

Nice paper BTW. I quite liked the presentation in SIAM NS.@nwlandry: Thank you!! :)

Okay, I ran the following:

```python3
from scipy.special import comb
import random


def _index_to_edge_comb(i, n, k):
    """
    returns the i-th combination of k numbers chosen from 0,1,2,...,n-1
    """
    c = []
    r = i
    j = -1
    for s in range(1,k+1):
        cs = j+1
        while r - comb(n-1-cs,k-s, exact=True) > 0:
            r -= comb(n-1-cs,k-s, exact=True)
            cs += 1
        c.append(cs)
        j = cs
    return c

def _index_to_edge(index, n, m):
    """Generate a hyperedge given an index in the list of possible edges.

    Parameters
    ----------
    index : int > 0
        The index of the hyperedge in the list of all possible hyperedges.
    n : int > 0
        The number of nodes
    m : int > 0
        The hyperedge size.

    Returns
    -------
    list
        The reconstructed hyperedge

    See Also
    --------
    _index_to_edge_partition

    References
    ----------
    https://stackoverflow.com/questions/53834707/element-at-index-in-itertools-product
    """
    return [(index // (n**r) % n) for r in range(m - 1, -1, -1)]

# Run timeit on the results:
n = 1000
m = 3

max_prod = n**m
max_comb = comb(n, m, exact=True)

%timeit _index_to_edge(random.randrange(max_prod), n, m)
%timeit _index_to_edge_comb(random.randrange(max_comb), n, m)
```

The results were 1.02 µs ± 5.64 ns per loop when iterating over the cartesian product and 315 µs ± 6.63 µs per loop when iterating over unique combinations. Of course, there is a factor of $\approx m!$ fewer indices to loop over as you previously mentioned, but even when accounting for this difference, iterating over the combinations is ~50x slower.@nwlandry: I think that at the very least, the documentation could be improved to explain exactly what the probability is of (i.e., cartesian product instead of combinations). Any other things that you would suggest?@arashbm: I think it's a problem of naming and expectations. In dyadic undirected Erdős-Rényi $G(n, p)$ networks, edges $(i, j)$ and $(j, i)$ don't get a separate Bernoulli trial each with probability $p$. But if I use `uniform_erdos_renyi_hypergraph(n, 2, p, p_type="prob")` I get twice as many edges as what i expect, $p \frac{n (n-1)}{2}$.

Of course one alternative is to just call it something else and add to the documentation and hope that that's enough distinction to avoid confusion.

Another alternative is to continue with the Cartesian as it is but divide $p$ by $m!$ before proceeding. This works but has the downside of creating duplicate edges, a behaviour that is again inconsistent with what people would expect from something called `erdos_renyi`.

You can continue with the Cartesian and keep `p` as it is, but replace the `len(set(...)) == m` check to something more restrictive e.g., reject all edges except the ones that have a strictly ascending order of nodes. This has the $m!$ overhead of extra loops but no duplicates.

I would say no amount of documentation is going to fix API usability problems. I'm not even sure renaming would create enough distinction. But in any case I'm not in a position to suggest if you should break compatibility or not.@maximelucas: > The combinatorial explosion of the mask array makes it impossible to use in some scenarios
> ```
> numpy.core._exceptions._ArrayMemoryError: Unable to allocate 9.72 PiB for an array with shape
> ```

Ah that's a fair point @arashbm, thanks!

> But that also has an easy solution of not pre-allocating mask:  

Yea that's actually how we used to do it. We changed to using a numpy array to make it slightly faster, but I didn't think of the problem you mentioned. I'll change it back. 
@maximelucas: I opened #361  related to this, is there anything else for which this issue should remain open? @nwlandry @nwlandry: Yeah, I don't think I've addressed this yet, but I will prioritize getting this done.

Solution patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 6b7c83639..610dddc64 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -13,12 +13,12 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
+        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v4
 
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 
@@ -112,7 +112,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu, macos, windows]
-        python-version: ['3.8', '3.9', '3.10', '3.11', '3.12']
+        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
     steps:
     - uses: actions/checkout@v4
 
diff --git a/README.md b/README.md
index 672aeb0c5..4a2575e37 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 <img src='https://github.com/xgi-org/xgi/blob/main/logo/logo.svg' alt="XGI" width='50%'></img>
 
 [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
-[![Supports Python versions 3.8 and above.](https://img.shields.io/badge/Python%20versions%20supported-3.8%2B-forest
+[![Supports Python versions 3.9 and above.](https://img.shields.io/badge/Python%20versions%20supported-3.9%2B-forest
 )](https://www.repostatus.org/#active)
 [![Test Status](https://github.com/xgi-org/xgi/workflows/test/badge.svg?branch=main)](https://github.com/xgi-org/xgi/actions?query=workflow%3A%22test%22)
 [![codecov](https://codecov.io/gh/xgi-org/xgi/branch/main/graph/badge.svg?token=BI6TX2WDSG)](https://codecov.io/gh/xgi-org/xgi)
@@ -42,7 +42,7 @@ Comple**X** **G**roup **I**nteractions (**XGI**) is a Python package for higher-
 * Draw higher-order networks in a variety of **visually striking ways** (See our [gallery](https://xgi.readthedocs.io/en/stable/gallery.html) for several examples.)
 
 ## Installation<a id="installation"></a>
-XGI runs on Python 3.8 or higher.
+XGI runs on Python 3.9 or higher.
 
 To install the latest version of XGI, run the following command:
 ```sh
diff --git a/docs/source/api/algorithms.rst b/docs/source/api/algorithms.rst
index 41a5a91ce..b6ebf1218 100644
--- a/docs/source/api/algorithms.rst
+++ b/docs/source/api/algorithms.rst
@@ -11,5 +11,6 @@ algorithms package
    ~xgi.algorithms.centrality
    ~xgi.algorithms.clustering
    ~xgi.algorithms.connected
+   ~xgi.algorithms.properties
    ~xgi.algorithms.shortest_path
-   ~xgi.algorithms.properties
\ No newline at end of file
+   ~xgi.algorithms.simpliciality
\ No newline at end of file
diff --git a/docs/source/api/algorithms/xgi.algorithms.simpliciality.rst b/docs/source/api/algorithms/xgi.algorithms.simpliciality.rst
new file mode 100644
index 000000000..e4d2e4a5c
--- /dev/null
+++ b/docs/source/api/algorithms/xgi.algorithms.simpliciality.rst
@@ -0,0 +1,14 @@
+xgi.algorithms.simpliciality
+============================
+
+.. currentmodule:: xgi.algorithms.simpliciality
+
+.. automodule:: xgi.algorithms.simpliciality
+   
+   .. rubric:: Functions
+   
+   .. autofunction:: edit_simpliciality
+   .. autofunction:: simplicial_edit_distance
+   .. autofunction:: face_edit_simpliciality
+   .. autofunction:: mean_face_edit_distance
+   .. autofunction:: simplicial_fraction
\ No newline at end of file
diff --git a/docs/source/api/stats/xgi.stats.nodestats.rst b/docs/source/api/stats/xgi.stats.nodestats.rst
index 23f0a66e4..7d9195039 100644
--- a/docs/source/api/stats/xgi.stats.nodestats.rst
+++ b/docs/source/api/stats/xgi.stats.nodestats.rst
@@ -15,4 +15,7 @@
    .. autofunction:: h_eigenvector_centrality
    .. autofunction:: local_clustering_coefficient
    .. autofunction:: node_edge_centrality
-   .. autofunction:: two_node_clustering_coefficient
\ No newline at end of file
+   .. autofunction:: two_node_clustering_coefficient
+   .. autofunction:: local_simplicial_fraction
+   .. autofunction:: local_edit_simpliciality
+   .. autofunction:: local_face_edit_simpliciality
\ No newline at end of file
diff --git a/docs/source/api_reference.rst b/docs/source/api_reference.rst
index 51777d2c4..c103f902b 100644
--- a/docs/source/api_reference.rst
+++ b/docs/source/api_reference.rst
@@ -20,8 +20,14 @@ XGI is organized into the following subpackages:
 
 .. note::
 
-   This page describes release |version| of XGI. Find out more about this release and previous releases in the `release notes <https://github.com/xgi-org/xgi/releases>`_.
-   
+   .. only:: stable_version
+
+      This page describes release |version| of XGI. Find out more about this release and previous releases in the `release notes <https://github.com/xgi-org/xgi/releases>`_.
+      
+   .. only:: dev_version
+      
+      This page describes the latest XGI development with features that may not be released yet. To see officially released features, visit the `API reference <https://xgi.readthedocs.io/en/stable/api_reference.html>`_ for the most recent stable version.
+
 .. toctree::
    :maxdepth: 2
 
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 3056bbc41..180ae8492 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -130,11 +130,12 @@
 else:
     raise Exception(f"Error: HTTP response {r.status_code}")
 
-
 rst_epilog = f"""
 .. role:: raw-html(raw)
    :format: html
 .. |release_announcement| replace:: :raw-html:`<a href={release_url}><button type="button" class="version-button">XGI {release_version} released! {release_date}</button></a>`
+.. |release_date| replace:: {release_date}
+.. |release_version| replace:: {release_version}
 """
 
 today = release_date
@@ -200,6 +201,12 @@
 else:
     version_match = "stable"
 
+
+if version_match == "stable":
+    tags.add("stable_version")
+elif version_match == "dev":
+    tags.add("dev_version")
+
 # documentation.
 html_theme_options = {
     "logo": {
diff --git a/docs/source/installing.rst b/docs/source/installing.rst
index 2c848b9cd..367d91f9b 100644
--- a/docs/source/installing.rst
+++ b/docs/source/installing.rst
@@ -13,7 +13,7 @@ To install XGI, execute the following in the command line:
    pip install xgi
 
 
-XGI was developed and tested for Python 3.8-3.12 on Mac OS, Windows, and Ubuntu.
+XGI was developed and tested for Python 3.9-3.13 on Mac OS, Windows, and Ubuntu.
 
 Once installed, go directly to the `User Guides <user_guides.html>`_ to get started!
 
diff --git a/docs/source/using-xgi.rst b/docs/source/using-xgi.rst
index d2c2ec82f..dbd32d0d3 100644
--- a/docs/source/using-xgi.rst
+++ b/docs/source/using-xgi.rst
@@ -27,6 +27,11 @@ Nicholas W. Landry, Ilya Amburg, Mirah Shi, and Sinan Aksoy, "Filtering higher-o
 :bdg-link-primary-line:`Paper <https://doi.org/10.1088/2632-072X/ad253a>`
 :bdg-link-primary-line:`Code <https://github.com/nwlandry/filtering-higher-order-datasets>`
 
+Nicholas W. Landry, William Thompson, Laurent Hébert-Dufresne, and Jean-Gabriel Young, "Reconstructing networks from simple and complex contagions", Physical Review E **110**, L042301 (2024).
+
+:bdg-link-primary-line:`Paper <https://doi.org/10.1103/PhysRevE.110.L042301>`
+:bdg-link-primary-line:`Code <https://github.com/nwlandry/complex-network-reconstruction>`
+
 Nicholas W. Landry, Jean-Gabriel Young, and Nicole Eikmeier, "The simpliciality of higher-order networks", *EPJ Data Science* **13**, 17 (2024).
 
 :bdg-link-primary-line:`Paper <https://doi.org/10.1140/epjds/s13688-024-00458-1>`
@@ -37,6 +42,11 @@ Leonie Neuhäuser, Michael Scholkemper, Francesco Tudisco, and Michael T. Schaub
 :bdg-link-primary-line:`Paper <https://doi.org/10.1126/sciadv.adh4053>`
 :bdg-link-primary-line:`Code <https://git.rwth-aachen.de/netsci/publication-2023-learning-the-effective-order-of-a-hypergraph-dynamical-system>`
 
+Yuanzhao Zhang, Per Sebastian Skardal, Federico Battiston, Giovanni Petri, and Maxime Lucas, "Deeper but smaller: Higher-order interactions increase linear stability but shrink basins", *Science Advances* **10**,eado8049 (2024).
+
+:bdg-link-primary-line:`Paper <https://doi.org/10.1126/sciadv.ado8049>`
+:bdg-link-primary-line:`Code <https://github.com/maximelucas/basins_and_triangles>`
+
 2023
 ----
 
@@ -81,12 +91,6 @@ Robin Delabays, Giulia De Pasquale, Florian Dörfler, and Yuanzhao Zhang, "Hyper
 :bdg-link-primary-line:`Paper <https://arxiv.org/abs/2402.00078>`
 :bdg-link-primary-line:`Code <https://github.com/TaylorBasedHypergraphInference/THIS>`
 
-Nicholas W. Landry, William Thompson, Laurent Hébert-Dufresne, and Jean-Gabriel Young, "Complex contagions can outperform simple contagions for network reconstruction with dense networks or saturated dynamics", arXiv:2405.00129
-
-:bdg-link-primary-line:`Paper <https://arxiv.org/abs/2405.00129>`
-:bdg-link-primary-line:`Code <https://github.com/nwlandry/complex-network-reconstruction>`
-
-
 Maxime Lucas, Luca Gallo, Arsham Ghavasieh, Federico Battiston, and Manlio De Domenico, "Functional reducibility of higher-order networks", arXiv:2404.08547 (2024).
 
 :bdg-link-primary-line:`Paper <https://arxiv.org/abs/2404.08547>`
@@ -104,10 +108,6 @@ Gonzalo Contreras-Aso, Cristian Pérez-Corral, and Miguel Romance, "Uplifting ed
 :bdg-link-primary-line:`Paper <https://arxiv.org/abs/2310.20335>`
 :bdg-link-primary-line:`Code <https://github.com/LaComarca-Lab/non-uniform-hypergraphs>`
 
-Yuanzhao Zhang, Per Sebastian Skardal, Federico Battiston, Giovanni Petri, and Maxime Lucas, "Deeper but smaller: Higher-order interactions increase linear stability but shrink basins", arXiv:2309.16581 (2023).
-
-:bdg-link-primary-line:`Paper <https://arxiv.org/abs/2309.16581>`
-:bdg-link-primary-line:`Code <https://github.com/maximelucas/basins_and_triangles>`
 
 Theses
 ======
diff --git a/pyproject.toml b/pyproject.toml
index bb25b33a3..66fd47954 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ name = "xgi"
 dynamic = ["version"]
 description = "XGI is a Python package for higher-order networks."
 readme = "README.md"
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 license = {file = "LICENSE.md"}
 authors = [
     {name = "XGI Developers", email = "nicholas.landry@virginia.edu"}
@@ -30,11 +30,11 @@ classifiers = [
     "License :: OSI Approved :: BSD License",
     "Operating System :: OS Independent",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Programming Language :: Python :: 3 :: Only",
     "Topic :: Software Development :: Libraries :: Python Modules",
     "Topic :: Scientific/Engineering :: Bio-Informatics",
diff --git a/xgi/algorithms/__init__.py b/xgi/algorithms/__init__.py
index 4ccac56e7..c5986ea25 100644
--- a/xgi/algorithms/__init__.py
+++ b/xgi/algorithms/__init__.py
@@ -1,7 +1,16 @@
-from . import assortativity, centrality, clustering, connected, properties
+from . import (
+    assortativity,
+    centrality,
+    clustering,
+    connected,
+    properties,
+    shortest_path,
+    simpliciality,
+)
 from .assortativity import *
 from .centrality import *
 from .clustering import *
 from .connected import *
 from .properties import *
 from .shortest_path import *
+from .simpliciality import *
diff --git a/xgi/algorithms/simpliciality.py b/xgi/algorithms/simpliciality.py
new file mode 100644
index 000000000..3110a9049
--- /dev/null
+++ b/xgi/algorithms/simpliciality.py
@@ -0,0 +1,508 @@
+from itertools import chain, combinations
+
+import numpy as np
+from scipy.special import binom
+
+from ..core import Hypergraph
+from ..utils import Trie
+
+__all__ = [
+    "edit_simpliciality",
+    "simplicial_edit_distance",
+    "face_edit_simpliciality",
+    "mean_face_edit_distance",
+    "simplicial_fraction",
+]
+
+
+def edit_simpliciality(H, min_size=2, exclude_min_size=True):
+    """Computes the edit simpliciality.
+
+    The fraction of sub-edges contained when compared to a simplicial complex.
+
+    Parameters
+    ----------
+    H : xgi.Hypergraph
+        The hypergraph of interest
+    min_size: int, optional
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces. For more details, see
+        the Notes below. By default, 2.
+    exclude_min_size : bool, optional
+        Whether to exclude minimal simplices when counting simplices.
+        For more detailed information, see the Notes below. By default, True.
+
+    Returns
+    -------
+    float
+        The edit simpliciality
+
+    See Also
+    --------
+    simplicial_edit_distance
+
+    Notes
+    -----
+    1. The formal definition of a simplicial complex can be unnecessarily
+    strict when used to represent perfect inclusion structures.
+    By definition, a simplex always contains singletons
+    (edges comprising a single node) and the empty set.
+    Several datasets will not include such interactions by construction.
+    To circumvent this issue, we use a relaxed definition of
+    downward closure that excludes edges of a certain size or smaller
+    wherever it makes sense. By default, we set the minimum size
+    to be 2 since some datasets do not contain singletons.
+
+    2. Hyperedges we call “minimal faces” may significantly skew the
+    simpliciality of a dataset. The minimal faces of a hypergraph :math:`H`
+    are the edges of the minimal size, i.e., :math:`|e| = \min(K)`, where :math:`K`
+    is the set of sizes that we consider based on note 1.
+    (In a traditional simplicial complex, the minimal faces are singletons).
+    With the size restrictions in place, the minimal faces of a hypergraph
+    are always simplices because, by definition, there are no smaller edges
+    for these edges to include. When measuring the simpliciality of a dataset,
+    it is most meaningful to focus on the faces for which inclusion is possible,
+    and so, by default, we exclude these minimal faces when counting potential
+    simplices.
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    return 1 - simplicial_edit_distance(
+        H, min_size=min_size, exclude_min_size=exclude_min_size
+    )
+
+
+def simplicial_edit_distance(H, min_size=2, exclude_min_size=True, normalize=True):
+    """Computes the simplicial edit distance.
+
+    The number (or fraction) of sub-edges needed to be added
+    to a hypergraph to make it a simplicial complex.
+
+    Parameters
+    ----------
+    H : xgi.Hypergraph
+        The hypergraph of interest
+    min_size: int, optional
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces. For more details, see
+        the Notes below. By default, 2.
+    exclude_min_size : bool, optional
+        Whether to exclude minimal simplices when counting simplices.
+        For more detailed information, see the Notes below. By default, True.
+
+    normalize : bool, optional
+        Whether to normalize by the total number of edges
+
+    Returns
+    -------
+    float
+        The edit simpliciality
+
+    See Also
+    --------
+    edit_simpliciality
+
+    Notes
+    -----
+    1. The formal definition of a simplicial complex can be unnecessarily
+    strict when used to represent perfect inclusion structures.
+    By definition, a simplex always contains singletons
+    (edges comprising a single node) and the empty set.
+    Several datasets will not include such interactions by construction.
+    To circumvent this issue, we use a relaxed definition of
+    downward closure that excludes edges of a certain size or smaller
+    wherever it makes sense. By default, we set the minimum size
+    to be 2 since some datasets do not contain singletons.
+
+    2. Hyperedges we call “minimal faces” may significantly skew the
+    simpliciality of a dataset. The minimal faces of a hypergraph :math:`H`
+    are the edges of the minimal size, i.e., :math:`|e| = \min(K)`, where :math:`K`
+    is the set of sizes that we consider based on note 1.
+    (In a traditional simplicial complex, the minimal faces are singletons).
+    With the size restrictions in place, the minimal faces of a hypergraph
+    are always simplices because, by definition, there are no smaller edges
+    for these edges to include. When measuring the simpliciality of a dataset,
+    it is most meaningful to focus on the faces for which inclusion is possible,
+    and so, by default, we exclude these minimal faces when counting potential
+    simplices.
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    edges = H.edges.filterby("size", min_size, "geq").members()
+
+    t = Trie()
+    t.build_trie(edges)
+
+    maxH = Hypergraph(
+        H.edges.maximal()
+        .filterby("size", min_size + exclude_min_size, "geq")
+        .members(dtype=dict)
+    )
+    if not maxH.edges:
+        return np.nan
+    
+    id_to_num = dict(zip(maxH.edges, range(maxH.num_edges)))
+
+    ms = 0
+    for id1, e in maxH.edges.members(dtype=dict).items():
+        redundant_missing_faces = set()
+        for id2 in maxH.edges.neighbors(id1):
+            if id_to_num[id2] < id_to_num[id1]:
+                c = maxH._edge[id2].intersection(e)
+                if len(c) >= min_size:
+                    redundant_missing_faces.update(_missing_subfaces(t, c, min_size))
+
+                    # we don't have to worry about the intersection being a max face
+                    # because a) there are no multiedges and b) these are all maximal
+                    # faces so no inclusions.
+                    if not t.search(c):
+                        redundant_missing_faces.add(frozenset(c))
+
+        mf = _count_missing_subfaces(t, e, min_size)
+        rmf = len(redundant_missing_faces)
+        ms += mf - rmf
+
+    if normalize:
+        s = len(edges)
+        mf = maxH.num_edges
+        if s - mf + ms > 0:
+            return ms / (s - mf + ms)
+        else:
+            return np.nan
+    else:
+        return ms
+
+
+def face_edit_simpliciality(H, min_size=2, exclude_min_size=True):
+    """Computes the face edit simpliciality.
+
+    The average fraction of sub-edges contained in a hyperedge
+    relative to a simplex.
+
+    Parameters
+    ----------
+    H : xgi.Hypergraph
+        The hypergraph of interest
+    min_size: int, optional
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces. For more details, see
+        the Notes below. By default, 2.
+    exclude_min_size : bool, optional
+        Whether to exclude minimal simplices when counting simplices.
+        For more detailed information, see the Notes below. By default, True.
+
+    Returns
+    -------
+    float
+        The face edit simpliciality
+
+    See Also
+    --------
+    mean_face_edit_distance
+
+    Notes
+    -----
+    1. The formal definition of a simplicial complex can be unnecessarily
+    strict when used to represent perfect inclusion structures.
+    By definition, a simplex always contains singletons
+    (edges comprising a single node) and the empty set.
+    Several datasets will not include such interactions by construction.
+    To circumvent this issue, we use a relaxed definition of
+    downward closure that excludes edges of a certain size or smaller
+    wherever it makes sense. By default, we set the minimum size
+    to be 2 since some datasets do not contain singletons.
+
+    2. Hyperedges we call “minimal faces” may significantly skew the
+    simpliciality of a dataset. The minimal faces of a hypergraph :math:`H`
+    are the edges of the minimal size, i.e., :math:`|e| = \min(K)`, where :math:`K`
+    is the set of sizes that we consider based on note 1.
+    (In a traditional simplicial complex, the minimal faces are singletons).
+    With the size restrictions in place, the minimal faces of a hypergraph
+    are always simplices because, by definition, there are no smaller edges
+    for these edges to include. When measuring the simpliciality of a dataset,
+    it is most meaningful to focus on the faces for which inclusion is possible,
+    and so, by default, we exclude these minimal faces when counting potential
+    simplices.
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    return 1 - mean_face_edit_distance(
+        H, min_size=min_size, exclude_min_size=exclude_min_size
+    )
+
+
+def mean_face_edit_distance(H, min_size=2, exclude_min_size=True, normalize=True):
+    """Computes the mean face edit distance
+
+    The average number (or fraction) of sub-edges needed to be added to make
+    a hyperedge a simplex.
+
+    Parameters
+    ----------
+    H : Hypergraph
+        The hypergraph of interest
+    min_size: int, optional
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces. For more details, see
+        the Notes below. By default, 2.
+    exclude_min_size : bool, optional
+        Whether to exclude minimal simplices when counting simplices.
+        For more detailed information, see the Notes below. By default, True.
+    normalize : bool, optional
+        Whether to normalize the face edit distance, by default True
+
+    Returns
+    -------
+    float
+        The mean face edit distance
+
+    See Also
+    --------
+    face_edit_simpliciality
+
+    Notes
+    -----
+    1. The formal definition of a simplicial complex can be unnecessarily
+    strict when used to represent perfect inclusion structures.
+    By definition, a simplex always contains singletons
+    (edges comprising a single node) and the empty set.
+    Several datasets will not include such interactions by construction.
+    To circumvent this issue, we use a relaxed definition of
+    downward closure that excludes edges of a certain size or smaller
+    wherever it makes sense. By default, we set the minimum size
+    to be 2 since some datasets do not contain singletons.
+
+    2. Hyperedges we call “minimal faces” may significantly skew the
+    simpliciality of a dataset. The minimal faces of a hypergraph :math:`H`
+    are the edges of the minimal size, i.e., :math:`|e| = \min(K)`, where :math:`K`
+    is the set of sizes that we consider based on note 1.
+    (In a traditional simplicial complex, the minimal faces are singletons).
+    With the size restrictions in place, the minimal faces of a hypergraph
+    are always simplices because, by definition, there are no smaller edges
+    for these edges to include. When measuring the simpliciality of a dataset,
+    it is most meaningful to focus on the faces for which inclusion is possible,
+    and so, by default, we exclude these minimal faces when counting potential
+    simplices.
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    t = Trie()
+    t.build_trie(H.edges.filterby("size", min_size, "geq").members())
+
+    max_faces = (
+        H.edges.maximal().filterby("size", min_size + exclude_min_size, "geq").members()
+    )
+    avg_d = 0
+    for e in max_faces:
+        if len(e) >= min_size:
+            d = _count_missing_subfaces(t, e, min_size=min_size)  # missing subfaces
+            m = _max_number_of_subfaces(min_size, len(e))
+            if normalize and m != 0:
+                d *= 1.0 / m
+            avg_d += d / len(max_faces)
+    return avg_d
+
+
+def simplicial_fraction(H, min_size=2, exclude_min_size=True):
+    """Computing the simplicial fraction for a hypergraph.
+
+    What fraction of the hyperedges are simplices?
+
+    Parameters
+    ----------
+    H : Hypergraph
+        The hypergraph of interest
+    min_size: int, optional
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces. For more details, see
+        the Notes below. By default, 2.
+    exclude_min_size : bool, optional
+        Whether to exclude minimal simplices when counting simplices.
+        For more detailed information, see the Notes below. By default, True.
+
+    Returns
+    -------
+    float
+        The simplicial fraction
+
+    Notes
+    -----
+    1. The formal definition of a simplicial complex can be unnecessarily
+    strict when used to represent perfect inclusion structures.
+    By definition, a simplex always contains singletons
+    (edges comprising a single node) and the empty set.
+    Several datasets will not include such interactions by construction.
+    To circumvent this issue, we use a relaxed definition of
+    downward closure that excludes edges of a certain size or smaller
+    wherever it makes sense. By default, we set the minimum size
+    to be 2 since some datasets do not contain singletons.
+
+    2. Hyperedges we call “minimal faces” may significantly skew the
+    simpliciality of a dataset. The minimal faces of a hypergraph :math:`H`
+    are the edges of the minimal size, i.e., :math:`|e| = \min(K)`, where :math:`K`
+    is the set of sizes that we consider based on note 1.
+    (In a traditional simplicial complex, the minimal faces are singletons).
+    With the size restrictions in place, the minimal faces of a hypergraph
+    are always simplices because, by definition, there are no smaller edges
+    for these edges to include. When measuring the simpliciality of a dataset,
+    it is most meaningful to focus on the faces for which inclusion is possible,
+    and so, by default, we exclude these minimal faces when counting potential
+    simplices.
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    try:
+        ns = _count_simplices(H, min_size, exclude_min_size)
+        ps = _potential_simplices(H, min_size, exclude_min_size)
+        return ns / ps
+    except ZeroDivisionError:
+        return np.nan
+
+
+#### Helper functions
+def _powerset(iterable, min_size=1, max_size=None):
+    """Generates a modified powerset.
+
+    User can specify the maximum and minimum size
+    of the sets in the powerset.
+
+    Parameters
+    ----------
+    iterable : iterable
+        The set for which to compute the powerset.
+    min_size: int, default: 1
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces.
+    max_size : int, default: None.
+        The maximum size to include when computing
+        the power set. When max_size=None, it generates
+        the powerset including the edge itself.
+
+    Returns
+    -------
+    itertools.chain
+        a generator of the sets in the powerset.
+    """
+    s = iterable
+    if max_size is None:
+        max_size = len(s)
+
+    return chain.from_iterable(
+        combinations(s, r) for r in range(min_size, max_size + 1)
+    )
+
+
+def _count_missing_subfaces(t, face, min_size=1):
+    """Computing the edit distance for a single face.
+
+    Parameters
+    ----------
+    t : Trie
+        The trie representing the hypergraph
+    face : iterable
+        The edge for which to find the edit distance
+    min_size: int, default: 1
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces.
+
+    Returns
+    -------
+    int
+        The edit distance
+    """
+    sub_edges = list(_powerset(face, min_size=min_size, max_size=len(face) - 1))
+    count = 0
+    for e in sub_edges:
+        if not t.search(e):
+            count += 1
+    return count
+
+
+def _missing_subfaces(t, face, min_size=1):
+    """Computing the edit distance for a single face.
+
+    Parameters
+    ----------
+    t : Trie
+        The trie representing the hypergraph
+    face : iterable
+        The edge for which to find the edit distance
+    min_size: int, default: 1
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces.
+
+    Returns
+    -------
+    int
+        The edit distance
+    """
+    sub_edges = list(_powerset(face, min_size=min_size, max_size=len(face) - 1))
+    ms = set()
+    for e in sub_edges:
+        if not t.search(e):
+            ms.add(frozenset(e))
+    return ms
+
+
+def _max_number_of_subfaces(min_size, max_size):
+    d = 2**max_size - 2  # subtract 2 for the face itself and the empty set
+    for i in range(1, min_size):
+        d -= binom(max_size, i)
+    return int(d)
+
+
+def _potential_simplices(H, min_size=2, exclude_min_size=True):
+    # record total number of hyperedges that are potential simplices
+    return len(H.edges.filterby("size", min_size + exclude_min_size, "geq"))
+
+
+def _count_simplices(H, min_size=2, exclude_min_size=True):
+    # build trie data structure
+    t = Trie()
+    all_edges = H.edges.members()
+    t.build_trie(all_edges)
+
+    edges = H.edges.filterby("size", min_size + exclude_min_size, "geq").members()
+
+    # for each hyperedge, determine if it's a simplex
+    count = 0
+    # The following loop is embarassingly parallel, so parallelize to increase speed would be good
+    for e in edges:
+        if _is_simplex(t, e, min_size):
+            count += 1
+    return count
+
+
+def _is_simplex(t, edge, min_size=2):
+    for e in _powerset(edge, min_size):
+        if not t.search(e):
+            return False
+    return True
diff --git a/xgi/core/dihypergraph.py b/xgi/core/dihypergraph.py
index fbcbc707e..6a4ce2aed 100644
--- a/xgi/core/dihypergraph.py
+++ b/xgi/core/dihypergraph.py
@@ -1039,8 +1039,7 @@ def cleanup(self, isolates=False, relabel=True, in_place=True):
 
             convert_labels_to_integers(_DH, in_place=True)
 
-        if not in_place:
-            return _DH
+        return _DH
 
     def freeze(self):
         """Method for freezing a dihypergraph which prevents it from being modified
diff --git a/xgi/core/hypergraph.py b/xgi/core/hypergraph.py
index 134e77257..329c86222 100644
--- a/xgi/core/hypergraph.py
+++ b/xgi/core/hypergraph.py
@@ -1487,8 +1487,7 @@ def cleanup(
 
             convert_labels_to_integers(_H, in_place=True)
 
-        if not in_place:
-            return _H
+        return _H
 
     def freeze(self):
         """Method for freezing a hypergraph which prevents it from being modified
diff --git a/xgi/core/simplicialcomplex.py b/xgi/core/simplicialcomplex.py
index 09885251e..e6f417e40 100644
--- a/xgi/core/simplicialcomplex.py
+++ b/xgi/core/simplicialcomplex.py
@@ -848,8 +848,7 @@ def cleanup(self, isolates=False, connected=True, relabel=True, in_place=True):
 
             convert_labels_to_integers(_S, in_place=True)
 
-        if not in_place:
-            return _S
+        return _S
 
     def freeze(self):
         """Method for freezing a simplicial complex
diff --git a/xgi/generators/uniform.py b/xgi/generators/uniform.py
index c010e8ad9..a5216fd10 100644
--- a/xgi/generators/uniform.py
+++ b/xgi/generators/uniform.py
@@ -7,9 +7,10 @@
 from functools import reduce
 
 import numpy as np
+from scipy.special import comb
 
 from ..exception import XGIError
-from .classic import empty_hypergraph
+from .classic import complete_hypergraph, empty_hypergraph
 
 __all__ = [
     "uniform_hypergraph_configuration_model",
@@ -107,6 +108,12 @@ def uniform_hypergraph_configuration_model(k, m, seed=None):
 def uniform_HSBM(n, m, p, sizes, seed=None):
     """Create a uniform hypergraph stochastic block model (HSBM).
 
+    This uses a fast method for generating hyperedges
+    so that instead of the algorithm being of complexity
+    :math:`\mathcal{O}(N^m)`, it can be as fast as
+    :math:`\mathcal{O}(m(N + |E|))`. See the references
+    for more details.
+
     Parameters
     ----------
     n : int
@@ -140,11 +147,21 @@ def uniform_HSBM(n, m, p, sizes, seed=None):
     --------
     uniform_HPPM
 
+    Notes
+    -----
+    Because XGI only stores edges as sets, when self-loops occur,
+    they become smaller edges (for example, the edge (0, 0, 0)
+    will be mapped to {0}). However, because this is explicitly
+    a *uniform* method, we discard these edges so that this is the case.
+    For sparse networks, this is a rare occurrence and this method offers
+    an order of magnitude speedup.
+
     References
     ----------
-    Nicholas W. Landry and Juan G. Restrepo.
-    "Polarization in hypergraphs with community structure."
-    Preprint, 2023. https://doi.org/10.48550/arXiv.2302.13967
+    Nicholas W. Landry and Juan G. Restrepo,
+    "Opinion disparity in hypergraphs with community structure",
+    Phys. Rev. E **108**, 034311 (2024).
+    https://doi.org/10.1103/PhysRevE.108.034311
     """
 
     # Check if dimensions match
@@ -156,7 +173,7 @@ def uniform_HSBM(n, m, p, sizes, seed=None):
     if len(set(np.shape(p))) != 1:
         raise XGIError("'p' must be a square tensor.")
     if np.max(p) > 1 or np.min(p) < 0:
-        raise XGIError("Entries of 'p' not in [0,1].")
+        raise XGIError("Entries of 'p' not in [0, 1].")
     if np.sum(sizes) != n:
         raise XGIError("Sum of sizes does not match n")
 
@@ -190,6 +207,12 @@ def uniform_HSBM(n, m, p, sizes, seed=None):
             while index < max_index:
                 indices = _index_to_edge_partition(index, partition_sizes, m)
                 e = {partition[block[i]][indices[i]] for i in range(m)}
+                # edge ids are not guaranteed to be unique
+                # and when casting to a set, they will no
+                # longer be of size m.
+                # for instance (0, 0, 0) becomes {0}
+                # if we accept these edges, the hypergraph
+                # will not longer be uniform, so we discard them.
                 if len(e) == m:
                     H.add_edge(e)
                 index += np.random.geometric(p[block])
@@ -199,6 +222,12 @@ def uniform_HSBM(n, m, p, sizes, seed=None):
 def uniform_HPPM(n, m, k, epsilon, rho=0.5, seed=None):
     """Construct the m-uniform hypergraph planted partition model (m-HPPM)
 
+    This uses a fast method for generating hyperedges
+    so that instead of the algorithm being of complexity
+    :math:`\mathcal{O}(N^m)`, it can be as fast as
+    :math:`\mathcal{O}(m(N + |E|))`. See the references
+    for more details.
+
     Parameters
     ----------
     n : int > 0
@@ -230,11 +259,21 @@ def uniform_HPPM(n, m, k, epsilon, rho=0.5, seed=None):
     --------
     uniform_HSBM
 
+    Notes
+    -----
+    Because XGI only stores edges as sets, when self-loops occur,
+    they become smaller edges (for example, the edge (0, 0, 0)
+    will be mapped to {0}). However, because this is explicitly
+    a *uniform* method, we discard these edges so that this is the case.
+    For sparse networks, this is a rare occurrence and this method offers
+    an order of magnitude speedup.
+
     References
     ----------
-    Nicholas W. Landry and Juan G. Restrepo.
-    "Polarization in hypergraphs with community structure."
-    Preprint, 2023. https://doi.org/10.48550/arXiv.2302.13967
+    Nicholas W. Landry and Juan G. Restrepo,
+    "Opinion disparity in hypergraphs with community structure",
+    Phys. Rev. E **108**, 034311 (2024).
+    https://doi.org/10.1103/PhysRevE.108.034311
     """
 
     if rho < 0 or rho > 1:
@@ -260,12 +299,18 @@ def uniform_HPPM(n, m, k, epsilon, rho=0.5, seed=None):
     return uniform_HSBM(n, m, p, sizes, seed=seed)
 
 
-def uniform_erdos_renyi_hypergraph(n, m, p, p_type="degree", seed=None):
+def uniform_erdos_renyi_hypergraph(n, m, p, p_type="prob", multiedges=False, seed=None):
     """Generate an m-uniform Erdős–Rényi hypergraph
 
     This creates a hypergraph with `n` nodes where
-    hyperedges of size `m` are created at random to
-    obtain a mean degree of `k`.
+    hyperedges of size `m` are created at random 
+    with probability (or to obtain a mean degree of) `p`.
+
+    This uses a fast method for generating hyperedges
+    so that instead of the algorithm being of complexity
+    :math:`\mathcal{O}(N^m)`, it can be as fast as
+    :math:`\mathcal{O}(m(N + |E|))`. See the references
+    for more details.
 
     Parameters
     ----------
@@ -273,11 +318,21 @@ def uniform_erdos_renyi_hypergraph(n, m, p, p_type="degree", seed=None):
         Number of nodes
     m : int > 0
         Hyperedge size
-    p : float or int > 0
-        Mean expected degree if p_type="degree" and
-        probability of an m-hyperedge if p_type="prob"
-    p_type : str
-        "degree" or "prob", by default "degree"
+    p : float or int >= 0
+        Probability of an m-hyperedge if p_type="prob" and
+        mean expected degree if p_type="degree"
+    p_type : str, optional
+        Determines the way p is interpreted (see p for detail).
+        Valid options are "prob" or "degree", by default "prob"
+    multiedges : bool, optional
+        Whether or not to allow multiedges. If True, there
+        can be significant speedups but at the cost of creating
+        (potentially unwanted) artifacts. When multiedges=True,
+        it treats each edge permutation as distinct, which can
+        lead to multiedges, especially for dense hypergraphs.
+        For sparse hypergraphs, however, this is unlikely to
+        be the case.
+        By default, False.
     seed : integer or None (default)
         The seed for the random number generator
 
@@ -286,10 +341,30 @@ def uniform_erdos_renyi_hypergraph(n, m, p, p_type="degree", seed=None):
     Hypergraph
         The Erdos Renyi hypergraph
 
-
     See Also
     --------
     ~xgi.generators.random.random_hypergraph
+
+    Notes
+    -----
+    When `multiedges=True`, there is the possibility of generating
+    (potentially unwanted) artifacts, like multiedges and loopy
+    hyperedges which are not present in the original Erdos-Renyi model.
+    Because hypergraphs in XGI only store edges as sets, when an edge
+    such as (0, 0, 0) is generated will be mapped to {0}. However,
+    because this is explicitly a *uniform* method, we discard these edges.
+    For sparse networks, this is a rare occurrence and allowing these
+    artifacts offers an order of magnitude speedup. Although allowing
+    loopy hyperedges and multiedges is not the default behavior, this
+    (as well as the associated performance boost) is enabled by setting
+    `multiedges=True`.
+
+    References
+    ----------
+    Nicholas W. Landry and Juan G. Restrepo,
+    "Opinion disparity in hypergraphs with community structure",
+    Phys. Rev. E **108**, 034311 (2024).
+    https://doi.org/10.1103/PhysRevE.108.034311
     """
     if seed is not None:
         np.random.seed(seed)
@@ -299,27 +374,70 @@ def uniform_erdos_renyi_hypergraph(n, m, p, p_type="degree", seed=None):
     H.add_nodes_from(node_labels)
 
     if p_type == "degree":
-        q = p / (m * n ** (m - 1))  # wiring probability
+        if multiedges:
+            q = p / (m * n ** (m - 1))  # wiring probability
+        else:
+            q = p * n / (m * comb(n, m))
     elif p_type == "prob":
         q = p
     else:
         raise XGIError("Invalid p_type!")
 
     if q > 1 or q < 0:
-        raise XGIError("Probability not in [0,1].")
+        raise XGIError("Probability not in [0, 1].")
+
+    if q == 1 and not multiedges:
+        return complete_hypergraph(n, order=m - 1)
+    if q == 0:
+        H = empty_hypergraph()
+        H.add_nodes_from(range(n))
+        return H
 
     index = np.random.geometric(q) - 1  # -1 b/c zero indexing
-    max_index = n**m
-    while index < max_index:
-        e = set(_index_to_edge(index, n, m))
+    if multiedges:
+        max_index = n**m
+        f = _index_to_edge_prod
+    else:
+        max_index = comb(n, m, exact=True)
+        f = _index_to_edge_comb
+
+    while index <= max_index:
+        e = set(f(index, n, m))
+        # if f corresponds to _index_to_edge_prod,
+        # edge ids are not guaranteed to be unique
+        # and when casting to a set, they will no
+        # longer be of size m.
+        # for instance (0, 0, 0) becomes {0}
+        # if we accept these edges, the hypergraph
+        # will not longer be uniform, so we discard them.
         if len(e) == m:
             H.add_edge(e)
+        # We no longer subtract 1 because if we did, the minimum
+        # value of the right-hand side would be zero, meaning that
+        # we sample the same index multiple times.
         index += np.random.geometric(q)
     return H
 
 
-def _index_to_edge(index, n, m):
-    """Generate a hyperedge given an index in the list of possible edges.
+def _index_to_edge_prod(index, n, m):
+    """Generate a hyperedge from an index given the
+    number of nodes and size of hyperedges.
+
+    In this method, it treats each edge permutation as distinct, which can
+    lead to multiedges and loopy edges, especially for dense hypergraphs.
+    Imagine that there is a hypergraph with 4 nodes and an edge size of 3.
+    We write out each edge (allowing duplicate entries) incrementing the last entry first,
+    followed by the second-to-last entry and so on, with each edge corresponding to an index
+    starting at zero. For example, (0, 0, 0) has index 0, (0, 0, 1) has index 1,
+    (0, 0, 2) has index 2, (0, 0, 3) has index 3, (0, 1, 0) has index 4, and so on.
+    This function will, for instance,
+    return (0, 0, 3) for index 3, network size 4, and edge size 3.
+
+    Because hypergraphs in XGI only store edges as sets, the edge (0, 0, 0),
+    for example, generated by this function will be mapped to {0}. However,
+    because this is explicitly a *uniform* method, we discard these edges.
+    For sparse networks, this is a rare occurrence and this method offers
+    an order of magnitude speedup.
 
     Parameters
     ----------
@@ -333,11 +451,12 @@ def _index_to_edge(index, n, m):
     Returns
     -------
     list
-        The reconstructed hyperedge
+        The hyperedge to which that index corresponds
 
     See Also
     --------
     _index_to_edge_partition
+    _index_to_edge_comb
 
     References
     ----------
@@ -346,13 +465,23 @@ def _index_to_edge(index, n, m):
     return [(index // (n**r) % n) for r in range(m - 1, -1, -1)]
 
 
-def _index_to_edge_partition(index, partition_sizes, m):
-    """Generate a hyperedge given an index in the list of possible edges
-    and a partition of community labels.
+def _index_to_edge_comb(index, n, m):
+    """Generate a hyperedge from an index given the number of nodes and size of hyperedges.
+
+    Imagine that there is a hypergraph with 4 nodes and an edge size of 3.
+    We write out each edge incrementing the last entry first, followed by the
+    second-to-last entry and so on, with each edge corresponding to an index
+    starting at zero. For example, (0, 1, 2) has index 0, (0, 1, 3) has index 0,
+    (0, 2, 3) has index 2, and (1, 2, 3) has index 3. This function will, for instance,
+    return (0, 2, 3) for index 2, network size 4, and edge size 3.
+
+
+    In this function, we prohibit multiedges, so each edge corresponds to a
+    unique index.
 
     Parameters
     ----------
-    index : int > 0
+    index : int >= 0
         The index of the hyperedge in the list of all possible hyperedges.
     n : int > 0
         The number of nodes
@@ -362,11 +491,65 @@ def _index_to_edge_partition(index, partition_sizes, m):
     Returns
     -------
     list
-        The reconstructed hyperedge
+        The hyperedge to which that index corresponds
+
+    See Also
+    --------
+    _index_to_edge_partition
+    _index_to_edge_prod
+
+    References
+    ----------
+    https://math.stackexchange.com/questions/1227409/indexing-all-combinations-without-making-list
+    """
+    c = []
+    r = index + 1  # makes it zero indexed
+    j = -1
+    for s in range(1, m + 1):
+        cs = j + 1
+        while r - comb(n - 1 - cs, m - s, exact=True) > 0:
+            r -= comb(n - 1 - cs, m - s, exact=True)
+            cs += 1
+        c.append(cs)
+        j = cs
+    return c
+
+
+def _index_to_edge_partition(index, partition_sizes, m):
+    """Generate a hyperedge from an index given the
+    number of nodes, size of hyperedges, and community sizes.
+
+    Imagine that there is a hypergraph with 10 nodes, an edge size of 3,
+    and two communities, the first of size 8 and the second of size 2.
+    We start out by specifying which community each node belongs to
+    and index into each community. For example, suppose the nodes
+    belong to communities 1, 1, and 2. Thene write out each edge
+    (allowing duplicate entries) incrementing the last entry first,
+    followed by the second-to-last entry and so on, with each edge
+    corresponding to an index starting at zero. For example, (0, 0, 0) has index 0,
+    (0, 0, 1) has index 1, (0, 1, 0) has index 2, (0, 1, 1) has index 3,
+    (0, 2, 0) has index 4, and so on. These are indices in each partition,
+    however, and we need the original labels of each node in each partition
+    to recover the nodes in each edge.
+
+    Parameters
+    ----------
+    index : int > 0
+        The index of the hyperedge in the list of all possible hyperedges.
+    partition_sizes : list or numpy array
+        The sizes of the partitions to which the nodes belong.
+    m : int > 0
+        The hyperedge size.
+
+    Returns
+    -------
+    list
+        The indices in each partition to which that index corresponds
 
     See Also
     --------
-    _index_to_edge
+    _index_to_edge_prod
+    _index_to_edge_comb
 
     """
     try:
diff --git a/xgi/readwrite/hif.py b/xgi/readwrite/hif.py
index c3de88487..9d00ee1f9 100644
--- a/xgi/readwrite/hif.py
+++ b/xgi/readwrite/hif.py
@@ -44,13 +44,13 @@ def write_hif(G, path):
     isolates = set(G.nodes.isolates())
     nodes_with_attrs = set(n for n in G.nodes if G.nodes[n])
     for n in isolates.union(nodes_with_attrs):
-        attr = {"attr": G.nodes[n]} if G.nodes[n] else {}
+        attr = {"attrs": G.nodes[n]} if G.nodes[n] else {}
         data["nodes"].append(IDDict({"node": n}) + attr)
 
     empty = set(G.edges.empty())
     edges_with_attrs = set(e for e in G.edges if G.edges[e])
     for e in empty.union(edges_with_attrs):
-        attr = {"attr": G.edges[e]} if G.edges[e] else {}
+        attr = {"attrs": G.edges[e]} if G.edges[e] else {}
         data["edges"].append(IDDict({"edge": e}) + attr)
 
     # hyperedge dict
@@ -164,8 +164,8 @@ def _convert_id(i, idtype):
     if "nodes" in data:
         for record in data["nodes"]:
             n = _convert_id(record["node"], nodetype)
-            if "attr" in record:
-                attr = record["attr"]
+            if "attrs" in record:
+                attr = record["attrs"]
             else:
                 attr = {}
 
@@ -178,8 +178,8 @@ def _convert_id(i, idtype):
     if "edges" in data:
         for record in data["edges"]:
             e = _convert_id(record["edge"], edgetype)
-            if "attr" in record:
-                attr = record["attr"]
+            if "attrs" in record:
+                attr = record["attrs"]
             else:
                 attr = {}
             if e not in G._edge:
diff --git a/xgi/stats/nodestats.py b/xgi/stats/nodestats.py
index a1fc91fc9..b09a871e9 100644
--- a/xgi/stats/nodestats.py
+++ b/xgi/stats/nodestats.py
@@ -493,3 +493,117 @@ def katz_centrality(net, bunch, cutoff=100):
     """
     c = xgi.katz_centrality(net, cutoff=cutoff)
     return {n: c[n] for n in c if n in bunch}
+
+
+def local_simplicial_fraction(net, bunch, min_size=2, exclude_min_size=True):
+    """The local simplicial fraction.
+
+    Parameters
+    ----------
+    net : xgi.Hypergraph
+        The network.
+    bunch : Iterable
+        Nodes in `net`.
+    min_size: int, default: 2
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces.
+    exclude_min_size : bool, optional
+        Whether to include minimal simplices when counting simplices, by default True
+
+    Returns
+    -------
+    dict
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    s = dict()
+    for n in bunch:
+        nbrs = net.nodes.neighbors(n)
+        if len(nbrs) == 0:
+            s[n] = np.nan
+        else:
+            nbrs.add(n)
+            sh = xgi.subhypergraph(net, nodes=nbrs)
+            s[n] = xgi.simplicial_fraction(sh, min_size, exclude_min_size)
+    return s
+
+
+def local_edit_simpliciality(net, bunch, min_size=2, exclude_min_size=True):
+    """The local edit simpliciality.
+
+    Parameters
+    ----------
+    net : xgi.Hypergraph
+        The network.
+    bunch : Iterable
+        Nodes in `net`.
+    min_size: int, default: 2
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces.
+    exclude_min_size : bool, optional
+        Whether to include minimal simplices when counting simplices, by default True
+
+    Returns
+    -------
+    dict
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    s = dict()
+    for n in bunch:
+        nbrs = net.nodes.neighbors(n)
+        if len(nbrs) == 0:
+            s[n] = np.nan
+        else:
+            nbrs.add(n)
+            sh = xgi.subhypergraph(net, nodes=nbrs)
+            s[n] = xgi.edit_simpliciality(sh, min_size, exclude_min_size)
+    return s
+
+
+def local_face_edit_simpliciality(net, bunch, min_size=2, exclude_min_size=True):
+    """The local face edit simpliciality.
+
+    Parameters
+    ----------
+    net : xgi.Hypergraph
+        The network.
+    bunch : Iterable
+        Nodes in `net`.
+    min_size: int, default: 2
+        The minimum hyperedge size to include when
+        calculating whether a hyperedge is a simplex
+        by counting subfaces.
+    exclude_min_size : bool, optional
+        Whether to include minimal simplices when counting simplices, by default True
+
+    Returns
+    -------
+    dict
+
+    References
+    ----------
+    "The simpliciality of higher-order order networks"
+    by Nicholas Landry, Jean-Gabriel Young, and Nicole Eikmeier,
+    *EPJ Data Science* **13**, 17 (2024).
+    """
+    s = dict()
+    for n in bunch:
+        nbrs = net.nodes.neighbors(n)
+        if len(nbrs) == 0:
+            s[n] = np.nan
+        else:
+            nbrs.add(n)
+            sh = xgi.subhypergraph(net, nodes=nbrs)
+            s[n] = xgi.face_edit_simpliciality(sh, min_size, exclude_min_size)
+    return s
diff --git a/xgi/utils/__init__.py b/xgi/utils/__init__.py
index 8c3aad65a..570aa121e 100644
--- a/xgi/utils/__init__.py
+++ b/xgi/utils/__init__.py
@@ -1,2 +1,3 @@
-from . import utilities
+from . import trie, utilities
+from .trie import *
 from .utilities import *
diff --git a/xgi/utils/trie.py b/xgi/utils/trie.py
new file mode 100644
index 000000000..b95bf9baf
--- /dev/null
+++ b/xgi/utils/trie.py
@@ -0,0 +1,44 @@
+"""Trie implementation for simpliciality.
+
+This is not meant to be user-facing, but rather
+contains auxiliary functions for the simpliciality
+methods.
+"""
+
+# This Trie implementation comes from user Ajay Rawat, https://stackoverflow.com/questions/11015320/how-to-create-a-trie-in-python
+
+__all__ = ["Trie"]
+
+
+class TrieNode:
+    def __init__(self):
+        # Dict: Key = letter, Item = TrieNode
+        self.children = {}
+        self.end = False
+
+
+class Trie:
+    def __init__(self):
+        self.root = TrieNode()
+
+    def build_trie(self, words):
+        for word in words:
+            self.insert(word)
+
+    def insert(self, word):
+        node = self.root
+        for char in sorted(word):
+            if char not in node.children:
+                node.children[char] = TrieNode()
+            node = node.children[char]
+        node.end = True
+
+    def search(self, word):
+        node = self.root
+        for char in sorted(word):
+            if char in node.children:
+                node = node.children[char]
+            else:
+                return False
+
+        return node.end


Test patch:
diff --git a/tests/algorithms/test_simpliciality.py b/tests/algorithms/test_simpliciality.py
new file mode 100644
index 000000000..f93d21deb
--- /dev/null
+++ b/tests/algorithms/test_simpliciality.py
@@ -0,0 +1,382 @@
+import numpy as np
+
+import xgi
+
+
+def test_edit_simpliciality(
+    sc1_with_singletons,
+    h_missing_one_singleton,
+    h_missing_one_link,
+    h_links_and_triangles2,
+    h1,
+):
+    # simplicial complex
+    es = xgi.edit_simpliciality(sc1_with_singletons)
+    assert es == 1.0
+
+    es = xgi.edit_simpliciality(sc1_with_singletons, min_size=1)
+    assert es == 1.0
+
+    es = xgi.edit_simpliciality(sc1_with_singletons, min_size=1, exclude_min_size=False)
+    assert es == 1.0
+
+    # h1
+    es = xgi.edit_simpliciality(h_missing_one_singleton)
+    assert es == 1.0
+
+    es = xgi.edit_simpliciality(h_missing_one_singleton, min_size=1)
+    assert np.allclose(es, 5 / 6)
+
+    es = xgi.edit_simpliciality(
+        h_missing_one_singleton, min_size=1, exclude_min_size=False
+    )
+    assert np.allclose(es, 5 / 6)
+
+    # h2
+    es = xgi.edit_simpliciality(h_missing_one_link)
+    assert np.allclose(es, 2 / 3)
+
+    es = xgi.edit_simpliciality(h_missing_one_link, min_size=1)
+    assert np.allclose(es, 5 / 6)
+
+    # links and triangles 2
+    es = xgi.edit_simpliciality(h_links_and_triangles2)
+    assert np.allclose(es, 2 / 3)
+
+    es = xgi.edit_simpliciality(h_links_and_triangles2, min_size=1)
+    assert np.allclose(es, 1 / 3)
+
+    es = xgi.edit_simpliciality(h_links_and_triangles2, exclude_min_size=False)
+    assert np.allclose(es, 3 / 5)
+
+    # test h1
+    es = xgi.edit_simpliciality(h1)
+    s = 4
+    m = 4 + 10
+    mf = 3
+    assert np.allclose(es, (s - mf) / (m + s - mf))
+
+    es = xgi.edit_simpliciality(h1, min_size=1)
+    s = 4
+    m = 4 + 10 + 7
+    mf = 3
+    assert np.allclose(es, (s - mf) / (m + s - mf))
+
+    es = xgi.edit_simpliciality(h1, exclude_min_size=False)
+    s = 4
+    m = 4 + 10
+    mf = 3
+    assert np.allclose(es, (s - mf) / (m + s - mf))
+
+    H = xgi.Hypergraph([[1, 2, 3, 4], [3, 4, 5, 6], [2, 3, 6, 7], [1, 2], [3, 4], [2, 3], [3, 6], [3, 4]])
+    es = xgi.edit_simpliciality(H, exclude_min_size=False)
+    assert np.allclose(es, 0.1785714285714286)
+
+
+def test_simplicial_edit_distance(
+    sc1_with_singletons,
+    h_missing_one_singleton,
+    h_missing_one_link,
+    h_links_and_triangles2,
+    h1,
+):
+    # simplicial complex
+    sed = xgi.simplicial_edit_distance(sc1_with_singletons)
+    assert sed == 0.0
+
+    sed = xgi.simplicial_edit_distance(sc1_with_singletons, min_size=1)
+    assert sed == 0.0
+
+    sed = xgi.simplicial_edit_distance(
+        sc1_with_singletons, min_size=1, exclude_min_size=False
+    )
+    assert sed == 0.0
+
+    # h1
+    sed = xgi.simplicial_edit_distance(h_missing_one_singleton)
+    assert sed == 0.0
+
+    sed = xgi.simplicial_edit_distance(h_missing_one_singleton, min_size=1)
+    assert np.allclose(sed, 1 / 6)
+
+    sed = xgi.simplicial_edit_distance(
+        h_missing_one_singleton, min_size=1, exclude_min_size=False
+    )
+    assert np.allclose(sed, 1 / 6)
+
+    # h2
+    sed = xgi.simplicial_edit_distance(h_missing_one_link)
+    assert np.allclose(sed, 1 / 3)
+
+    sed = xgi.simplicial_edit_distance(h_missing_one_link, min_size=1)
+    assert np.allclose(sed, 1 / 6)
+
+    # links and triangles 2
+    sed = xgi.simplicial_edit_distance(h_links_and_triangles2)
+    assert np.allclose(sed, 1 / 3)
+
+    sed = xgi.simplicial_edit_distance(h_links_and_triangles2, min_size=1)
+    assert np.allclose(sed, 2 / 3)
+
+    sed = xgi.simplicial_edit_distance(h_links_and_triangles2, exclude_min_size=False)
+    assert np.allclose(sed, 2 / 5)
+
+    # test h1
+    sed = xgi.simplicial_edit_distance(h1)
+    s = 4
+    m = 4 + 10
+    mf = 3
+    assert np.allclose(sed, m / (m + s - mf))
+
+    sed = xgi.simplicial_edit_distance(h1, min_size=1)
+    s = 4
+    m = 4 + 10 + 7
+    mf = 3
+    assert np.allclose(sed, m / (m + s - mf))
+
+    sed = xgi.simplicial_edit_distance(h1, exclude_min_size=False)
+    s = 4
+    m = 4 + 10
+    mf = 3
+    assert np.allclose(sed, m / (m + s - mf))
+
+    sed = xgi.simplicial_edit_distance(h1, exclude_min_size=False, normalize=False)
+    assert np.allclose(sed, m)
+
+
+def test_face_edit_simpliciality(
+    sc1_with_singletons,
+    h_missing_one_singleton,
+    h_missing_one_link,
+    h_links_and_triangles2,
+):
+    # simplicial complex
+    fes = xgi.face_edit_simpliciality(sc1_with_singletons)
+    assert fes == 1.0
+
+    fes = xgi.face_edit_simpliciality(sc1_with_singletons, min_size=1)
+    assert fes == 1.0
+
+    fes = xgi.face_edit_simpliciality(
+        sc1_with_singletons, min_size=1, exclude_min_size=False
+    )
+    assert fes == 1.0
+
+    # h1
+    fes = xgi.face_edit_simpliciality(h_missing_one_singleton)
+    assert fes == 1.0
+
+    fes = xgi.face_edit_simpliciality(h_missing_one_singleton, min_size=1)
+    assert np.allclose(fes, 5 / 6)
+
+    fes = xgi.face_edit_simpliciality(
+        h_missing_one_singleton, min_size=1, exclude_min_size=False
+    )
+    assert np.allclose(fes, 5 / 6)
+
+    # h2
+    fes = xgi.face_edit_simpliciality(h_missing_one_link)
+    assert np.allclose(fes, 2 / 3)
+
+    fes = xgi.face_edit_simpliciality(h_missing_one_link, min_size=1)
+    assert np.allclose(fes, 5 / 6)
+
+    # links and triangles 2
+    fes = xgi.face_edit_simpliciality(h_links_and_triangles2)
+    assert np.allclose(fes, 2 / 3)
+
+    fes = xgi.face_edit_simpliciality(h_links_and_triangles2, min_size=1)
+    assert np.allclose(fes, 2 / 9)
+
+    fes = xgi.face_edit_simpliciality(h_links_and_triangles2, exclude_min_size=False)
+    assert np.allclose(fes, 7 / 9)
+
+
+def test_mean_face_edit_distance(
+    sc1_with_singletons,
+    h_missing_one_singleton,
+    h_missing_one_link,
+    h_links_and_triangles2,
+):
+    # simplicial complex
+    mfed = xgi.mean_face_edit_distance(sc1_with_singletons)
+    assert mfed == 0.0
+
+    mfed = xgi.mean_face_edit_distance(sc1_with_singletons, min_size=1)
+    assert mfed == 0.0
+
+    mfed = xgi.mean_face_edit_distance(
+        sc1_with_singletons, min_size=1, exclude_min_size=False
+    )
+    assert mfed == 0.0
+
+    # h1
+    mfed = xgi.mean_face_edit_distance(h_missing_one_singleton)
+    assert mfed == 0.0
+
+    mfed = xgi.mean_face_edit_distance(h_missing_one_singleton, min_size=1)
+    assert np.allclose(mfed, 1 / 6)
+
+    mfed = xgi.mean_face_edit_distance(
+        h_missing_one_singleton, min_size=1, exclude_min_size=False
+    )
+    assert np.allclose(mfed, 1 / 6)
+
+    # h2
+    mfed = xgi.mean_face_edit_distance(h_missing_one_link)
+    assert np.allclose(mfed, 1 / 3)
+
+    mfed = xgi.mean_face_edit_distance(h_missing_one_link, min_size=1)
+    assert np.allclose(mfed, 1 / 6)
+
+    # links and triangles 2
+    mfed = xgi.mean_face_edit_distance(h_links_and_triangles2)
+    assert np.allclose(mfed, 1 / 3)
+
+    mfed = xgi.mean_face_edit_distance(h_links_and_triangles2, min_size=1)
+    assert np.allclose(mfed, 7 / 9)
+
+    mfed = xgi.mean_face_edit_distance(h_links_and_triangles2, exclude_min_size=False)
+    assert np.allclose(mfed, 2 / 9)
+
+
+def test_simplicial_fraction(
+    sc1_with_singletons, h_missing_one_singleton, h_missing_one_link
+):
+    # simplicial complex
+    sf = xgi.simplicial_fraction(sc1_with_singletons)
+    assert sf == 1.0
+
+    sf = xgi.simplicial_fraction(sc1_with_singletons, min_size=1)
+    assert sf == 1.0
+
+    sf = xgi.simplicial_fraction(
+        sc1_with_singletons, min_size=1, exclude_min_size=False
+    )
+    assert sf == 1.0
+
+    # h1
+    sf = xgi.simplicial_fraction(h_missing_one_singleton)
+    assert sf == 1.0
+
+    sf = xgi.simplicial_fraction(h_missing_one_singleton, min_size=1)
+    assert sf == 1 / 4
+
+    sf = xgi.simplicial_fraction(
+        h_missing_one_singleton, min_size=1, exclude_min_size=False
+    )
+    assert sf == 0.5
+
+    # h2
+    sf = xgi.simplicial_fraction(h_missing_one_link)
+    assert sf == 0
+
+    sf = xgi.simplicial_fraction(h_missing_one_link, min_size=1)
+    assert sf == 2 / 3
+
+
+def test_is_simplex(sc1_with_singletons, h_missing_one_singleton):
+    t = xgi.Trie()
+    edges = sc1_with_singletons.edges.members()
+    t.build_trie(edges)
+
+    is_simplex = xgi.algorithms.simpliciality._is_simplex
+
+    assert is_simplex(t, {1, 2, 3}, min_size=2)
+    assert is_simplex(t, {1, 2, 3}, min_size=1)
+    assert is_simplex(t, {1, 2}, min_size=1)
+
+    t = xgi.Trie()
+    edges = h_missing_one_singleton.edges.members()
+    t.build_trie(edges)
+
+    assert is_simplex(t, {1, 2, 3})
+    assert not is_simplex(t, {1, 2, 3}, min_size=1)
+    assert not is_simplex(t, {2, 3}, min_size=1)
+    assert is_simplex(t, {1, 2}, min_size=1)
+
+
+def test_count_simplices(sc1_with_singletons, h_missing_one_singleton):
+    count_simplices = xgi.algorithms.simpliciality._count_simplices
+
+    ns = count_simplices(sc1_with_singletons)
+    assert ns == 1
+
+    ns = count_simplices(sc1_with_singletons, min_size=1)
+    assert ns == 4
+
+    ns = count_simplices(sc1_with_singletons, min_size=1, exclude_min_size=False)
+    assert ns == 7
+
+    ns = count_simplices(h_missing_one_singleton, min_size=1)
+    assert ns == 1
+
+
+def test_potential_simplices(sc1_with_singletons, h_missing_one_link):
+    potential_simplices = xgi.algorithms.simpliciality._potential_simplices
+
+    ps = potential_simplices(sc1_with_singletons)
+    assert ps == 1
+
+    ps = potential_simplices(sc1_with_singletons, min_size=1)
+    assert ps == 4
+
+    ps = potential_simplices(sc1_with_singletons, min_size=1, exclude_min_size=False)
+    assert ps == 7
+
+    ps = potential_simplices(h_missing_one_link, min_size=1)
+    assert ps == 3
+
+
+def test_powerset():
+    powerset = xgi.algorithms.simpliciality._powerset
+    a = {1, 2, 3}
+
+    # test default behavior
+    subsets = {frozenset(s) for s in powerset(a, min_size=1)}
+    assert subsets == {
+        frozenset({1}),
+        frozenset({2}),
+        frozenset({3}),
+        frozenset({1, 2}),
+        frozenset({1, 3}),
+        frozenset({2, 3}),
+        frozenset({1, 2, 3}),
+    }
+
+    subsets = {frozenset(s) for s in powerset(a, min_size=2)}
+    assert subsets == {
+        frozenset({1, 2}),
+        frozenset({1, 3}),
+        frozenset({2, 3}),
+        frozenset({1, 2, 3}),
+    }
+
+    subsets = {frozenset(s) for s in powerset(a, max_size=2)}
+    assert subsets == {
+        frozenset({1}),
+        frozenset({2}),
+        frozenset({3}),
+        frozenset({1, 2}),
+        frozenset({1, 3}),
+        frozenset({2, 3}),
+    }
+
+
+def test_count_missing_subfaces(h_missing_one_link):
+    count_missing_subfaces = xgi.algorithms.simpliciality._count_missing_subfaces
+    t = xgi.Trie()
+    t.build_trie(h_missing_one_link.edges.members())
+    assert count_missing_subfaces(t, {1}, min_size=2) == 0
+    assert count_missing_subfaces(t, {2, 3}, min_size=2) == 0
+    assert count_missing_subfaces(t, {2, 3}) == 0
+    assert count_missing_subfaces(t, {1, 2, 3}) == 1
+    assert count_missing_subfaces(t, {1, 2, 3}, min_size=2) == 1
+
+
+def test_max_number_of_subfaces():
+    max_number_of_subfaces = xgi.algorithms.simpliciality._max_number_of_subfaces
+    assert max_number_of_subfaces(1, 3) == 6
+    assert max_number_of_subfaces(2, 3) == 3
+    assert max_number_of_subfaces(1, 4) == 14
+    assert max_number_of_subfaces(2, 4) == 10
diff --git a/tests/conftest.py b/tests/conftest.py
index 453aea26e..e3d20b255 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -270,3 +270,36 @@ def dihyperwithattrs(diedgelist2, attr0, attr1, attr2, attr3, attr4, attr5):
     H.add_edges_from(diedgelist2)
     H.set_edge_attributes({0: attr3, 1: attr4, 2: attr5})
     return H
+
+
+### Fixtures for simpliciality
+
+
+@pytest.fixture
+def sc1_with_singletons():
+    return xgi.Hypergraph([{1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}])
+
+
+@pytest.fixture
+def h_missing_one_singleton():
+    return xgi.Hypergraph([{1}, {2}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}])
+
+
+@pytest.fixture
+def h_missing_one_link():
+    return xgi.Hypergraph([{1}, {2}, {3}, {1, 3}, {2, 3}, {1, 2, 3}])
+
+
+@pytest.fixture
+def h_links_and_triangles():
+    return xgi.Hypergraph([{1, 3}, {2, 3}, {1, 2, 3}])
+
+
+@pytest.fixture
+def h_links_and_triangles2():
+    return xgi.Hypergraph([{1, 3}, {2, 3}, {1, 2, 3}, {1, 4}, {2, 3, 4}, {2, 4}])
+
+
+@pytest.fixture
+def h1():
+    return xgi.Hypergraph([{1, 2, 3}, {2, 3, 4, 5}, {5, 6, 7}, {5, 6}])
diff --git a/tests/generators/test_uniform.py b/tests/generators/test_uniform.py
index ec0b4118b..ce0e7f0b8 100644
--- a/tests/generators/test_uniform.py
+++ b/tests/generators/test_uniform.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+from scipy.special import comb
 
 import xgi
 from xgi.exception import XGIError
@@ -150,14 +151,29 @@ def test_uniform_HPPM():
 def test_uniform_erdos_renyi_hypergraph():
     m = 2
     n = 10
-    k = 2
-    H1 = xgi.uniform_erdos_renyi_hypergraph(n, m, k, seed=0)
+    p = 1
+    H1 = xgi.uniform_erdos_renyi_hypergraph(n, m, p, seed=0)
+    ne1 = H1.num_edges
+    H1.merge_duplicate_edges(rename="tuple")
+    ne2 = H1.num_edges
+    assert ne1 == ne2
+    assert ne1 == comb(n, m)
 
     assert H1.num_nodes == 10
     assert xgi.unique_edge_sizes(H1) == [2]
 
+    H2 = xgi.uniform_erdos_renyi_hypergraph(n, m, p, seed=0, multiedges=True)
+    print(H2.edges)
+    ne1 = H2.num_edges
+    H2.merge_duplicate_edges()
+    ne2 = H2.num_edges
+    assert ne1 != ne2
+    assert ne1 == n**m - n  # remove loopy edges
+
     # test that the seed works
-    H2 = xgi.uniform_erdos_renyi_hypergraph(n, m, k, seed=0)
+    p = 0.1
+    H1 = xgi.uniform_erdos_renyi_hypergraph(n, m, p, seed=0)
+    H2 = xgi.uniform_erdos_renyi_hypergraph(n, m, p, seed=0)
 
     assert H1.edges.members(dtype=dict) == H2.edges.members(dtype=dict)
 
diff --git a/tests/readwrite/test_hif.py b/tests/readwrite/test_hif.py
index 3b33ea490..62ae6173b 100644
--- a/tests/readwrite/test_hif.py
+++ b/tests/readwrite/test_hif.py
@@ -58,19 +58,19 @@ def test_to_hif(
     assert jsondata["metadata"] == {"name": "test"}
 
     nodes = [
-        {"node": 1, "attr": {"color": "red", "name": "horse"}},
-        {"node": 2, "attr": {"color": "blue", "name": "pony"}},
-        {"node": 3, "attr": {"color": "yellow", "name": "zebra"}},
-        {"node": 4, "attr": {"color": "red", "name": "orangutan", "age": 20}},
-        {"node": 5, "attr": {"color": "blue", "name": "fish", "age": 2}},
+        {"node": 1, "attrs": {"color": "red", "name": "horse"}},
+        {"node": 2, "attrs": {"color": "blue", "name": "pony"}},
+        {"node": 3, "attrs": {"color": "yellow", "name": "zebra"}},
+        {"node": 4, "attrs": {"color": "red", "name": "orangutan", "age": 20}},
+        {"node": 5, "attrs": {"color": "blue", "name": "fish", "age": 2}},
     ]
 
     edges = [
-        {"edge": 0, "attr": {"color": "blue"}},
-        {"edge": 1, "attr": {"color": "red", "weight": 2}},
-        {"edge": 2, "attr": {"color": "yellow"}},
-        {"edge": 3, "attr": {"color": "purple"}},
-        {"edge": 4, "attr": {"color": "purple", "name": "test"}},
+        {"edge": 0, "attrs": {"color": "blue"}},
+        {"edge": 1, "attrs": {"color": "red", "weight": 2}},
+        {"edge": 2, "attrs": {"color": "yellow"}},
+        {"edge": 3, "attrs": {"color": "purple"}},
+        {"edge": 4, "attrs": {"color": "purple", "name": "test"}},
     ]
 
     incidences = [
diff --git a/tests/stats/test_nodestats.py b/tests/stats/test_nodestats.py
index 8770b9e53..ddb2401f0 100644
--- a/tests/stats/test_nodestats.py
+++ b/tests/stats/test_nodestats.py
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 
 import xgi
@@ -117,3 +118,173 @@ def test_attrs(hyperwithattrs, attr1, attr2, attr3, attr4, attr5):
 
     with pytest.raises(ValueError):
         H.nodes.attrs(-1).asdict()
+
+
+def test_local_edit_simpliciality(
+    sc1_with_singletons,
+    h_links_and_triangles2,
+):
+    # no nodes
+    H = xgi.Hypergraph()
+
+    assert H.local_edit_simpliciality() == dict()
+    assert H.nodes.local_edit_simpliciality().asdict() == dict()
+
+    val = sc1_with_singletons.nodes.local_edit_simpliciality.asdict()
+    true_val = {
+        1: 1.0,
+        2: 1.0,
+        3: 1.0,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_edit_simpliciality.asdict()
+    true_val = {
+        1: 2 / 3,
+        2: 2 / 3,
+        3: 2 / 3,
+        4: 2 / 3,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_edit_simpliciality(min_size=1).asdict()
+    true_val = {
+        1: 1 / 3,
+        2: 1 / 3,
+        3: 1 / 3,
+        4: 1 / 3,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_edit_simpliciality(
+        min_size=1, exclude_min_size=False
+    ).asdict()
+    true_val = {
+        1: 1 / 3,
+        2: 1 / 3,
+        3: 1 / 3,
+        4: 1 / 3,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+
+def test_local_face_edit_simpliciality(
+    sc1_with_singletons,
+    h_links_and_triangles2,
+):
+    # no nodes
+    H = xgi.Hypergraph()
+
+    assert H.local_face_edit_simpliciality() == dict()
+    assert H.nodes.local_face_edit_simpliciality().asdict() == dict()
+
+    val = sc1_with_singletons.nodes.local_face_edit_simpliciality.asdict()
+    true_val = {
+        1: 1.0,
+        2: 1.0,
+        3: 1.0,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_face_edit_simpliciality.asdict()
+    true_val = {
+        1: 2 / 3,
+        2: 2 / 3,
+        3: 2 / 3,
+        4: 2 / 3,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_face_edit_simpliciality(
+        min_size=1
+    ).asdict()
+    true_val = {
+        1: 0.2222222222222222,
+        2: 0.2222222222222222,
+        3: 0.2222222222222222,
+        4: 0.2222222222222222,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_face_edit_simpliciality(
+        min_size=1, exclude_min_size=False
+    ).asdict()
+    true_val = {
+        1: 0.2222222222222222,
+        2: 0.2222222222222222,
+        3: 0.2222222222222222,
+        4: 0.2222222222222222,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+
+def test_local_simplicial_fraction(
+    sc1_with_singletons,
+    h_links_and_triangles2,
+):
+    # no nodes
+    H = xgi.Hypergraph()
+
+    assert H.local_simplicial_fraction() == dict()
+    assert H.nodes.local_simplicial_fraction().asdict() == dict()
+
+    val = sc1_with_singletons.nodes.local_simplicial_fraction.asdict()
+    true_val = {
+        1: 1.0,
+        2: 1.0,
+        3: 1.0,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_simplicial_fraction.asdict()
+    true_val = {
+        1: 0.0,
+        2: 0.0,
+        3: 0.0,
+        4: 0.0,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_simplicial_fraction(min_size=1).asdict()
+    true_val = {
+        1: 0.0,
+        2: 0.0,
+        3: 0.0,
+        4: 0.0,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)
+
+    val = h_links_and_triangles2.nodes.local_simplicial_fraction(
+        min_size=1, exclude_min_size=False
+    ).asdict()
+    true_val = {
+        1: 0.0,
+        2: 0.0,
+        3: 0.0,
+        4: 0.0,
+    }
+    for n in true_val:
+        assert np.allclose(val[n], true_val[n])
+    assert sorted(val) == sorted(true_val)