Implementation of $S^1$ model (networkx#6858)

Add a geometric_soft_configuration_graph graph generator function implementing the S1 model.

---------
Co-authored-by: Ross Barnowski <rossbar@berkeley.edu>
Co-authored-by: Mridul Seth <seth.mridul@gmail.com>
Co-authored-by: Dan Schult <dschult@colgate.edu>

doc/reference/generators.rst

@@ -189,6 +189,7 @@ Geometric
    soft_random_geometric_graph
    thresholded_random_geometric_graph
    waxman_graph
+   geometric_soft_configuration_graph
 
 Line Graph
 ----------

networkx/generators/geometric.py

@@ -16,6 +16,7 @@
     "soft_random_geometric_graph",
     "thresholded_random_geometric_graph",
     "waxman_graph",
+    "geometric_soft_configuration_graph",
 ]
 
 
@@ -844,3 +845,201 @@ def thresholded_random_geometric_graph(
     )
     G.add_edges_from(edges)
     return G
+
+
+@py_random_state(5)
+@nx._dispatch(graphs=None)
+def geometric_soft_configuration_graph(
+    *, beta, n=None, gamma=None, mean_degree=None, kappas=None, seed=None
+):
+    r"""Returns a random graph from the geometric soft configuration model.
+
+    The $\mathbb{S}^1$ model [1]_ is the geometric soft configuration model
+    which is able to explain many fundamental features of real networks such as
+    small-world property, heteregenous degree distributions, high level of
+    clustering, and self-similarity.
+
+    In the geometric soft configuration model, a node $i$ is assigned two hidden
+    variables: a hidden degree $\kappa_i$, quantifying its popularity, influence,
+    or importance, and an angular position $\theta_i$ in a circle abstracting the
+    similarity space, where angular distances between nodes are a proxy for their
+    similarity. Focusing on the angular position, this model is often called
+    the $\mathbb{S}^1$ model (a one-dimensional sphere). The circle's radius is
+    adjusted to $R = N/2\pi$, where $N$ is the number of nodes, so that the density
+    is set to 1 without loss of generality.
+
+    The connection probability between any pair of nodes increases with
+    the product of their hidden degrees (i.e., their combined popularities),
+    and decreases with the angular distance between the two nodes.
+    Specifically, nodes $i$ and $j$ are connected with the probability
+
+    $p_{ij} = \frac{1}{1 + \frac{d_{ij}^\beta}{\left(\mu \kappa_i \kappa_j\right)^{\max(1, \beta)}}}$
+
+    where $d_{ij} = R\Delta\theta_{ij}$ is the arc length of the circle between
+    nodes $i$ and $j$ separated by an angular distance $\Delta\theta_{ij}$.
+    Parameters $\mu$ and $\beta$ (also called inverse temperature) control the
+    average degree and the clustering coefficient, respectively.
+
+    It can be shown [2]_ that the model undergoes a structural phase transition
+    at $\beta=1$ so that for $\beta<1$ networks are unclustered in the thermodynamic
+    limit (when $N\to \infty$) whereas for $\beta>1$ the ensemble generates
+    networks with finite clustering coefficient.
+
+    The $\mathbb{S}^1$ model can be expressed as a purely geometric model
+    $\mathbb{H}^2$ in the hyperbolic plane [3]_ by mapping the hidden degree of
+    each node into a radial coordinate as
+
+    $r_i = \hat{R} - \frac{2 \max(1, \beta)}{\beta \zeta} \ln \left(\frac{\kappa_i}{\kappa_0}\right)$
+
+    where $\hat{R}$ is the radius of the hyperbolic disk and $\zeta$ is the curvature,
+
+    $\hat{R} = \frac{2}{\zeta} \ln \left(\frac{N}{\pi}\right)
+    - \frac{2\max(1, \beta)}{\beta \zeta} \ln (\mu \kappa_0^2)$
+
+    The connection probability then reads
+
+    $p_{ij} = \frac{1}{1 + \exp\left({\frac{\beta\zeta}{2} (x_{ij} - \hat{R})}\right)}$
+
+    where
+
+    $x_{ij} = r_i + r_j + \frac{2}{\zeta} \ln \frac{\Delta\theta_{ij}}{2}$
+
+    is a good approximation of the hyperbolic distance between two nodes separated
+    by an angular distance $\Delta\theta_{ij}$ with radial coordinates $r_i$ and $r_j$.
+    For $\beta > 1$, the curvature $\zeta = 1$, for $\beta < 1$, $\zeta = \beta^{-1}$.
+
+
+    Parameters
+    ----------
+    Either `n`, `gamma`, `mean_degree` are provided or `kappas`. The values of
+    `n`, `gamma`, `mean_degree` (if provided) are used to construct a random
+    kappa-dict keyed by node with values sampled from a power-law distribution.
+
+    beta : positive number
+        Inverse temperature, controlling the clustering coefficient.
+    n : int (default: None)
+        Size of the network (number of nodes).
+        If not provided, `kappas` must be provided and holds the nodes.
+    gamma : float (default: None)
+        Exponent of the power-law distribution for hidden degrees `kappas`.
+        If not provided, `kappas` must be provided directly.
+    mean_degree : float (default: None)
+        The mean degree in the network.
+        If not provided, `kappas` must be provided directly.
+    kappas : dict (default: None)
+        A dict keyed by node to its hidden degree value.
+        If not provided, random values are computed based on a power-law
+        distribution using `n`, `gamma` and `mean_degree`.
+    seed : int, random_state, or None (default)
+        Indicator of random number generation state.
+        See :ref:`Randomness<randomness>`.
+
+    Returns
+    -------
+    Graph
+        A random geometric soft configuration graph (undirected with no self-loops).
+        Each node has three node-attributes:
+
+        - ``kappa`` that represents the hidden degree.
+
+        - ``theta`` the position in the similarity space ($\mathbb{S}^1$) which is
+          also the angular position in the hyperbolic plane.
+
+        - ``radius`` the radial position in the hyperbolic plane
+          (based on the hidden degree).
+
+
+    Examples
+    --------
+    Generate a network with specified parameters:
+
+    >>> G = nx.geometric_soft_configuration_graph(beta=1.5, n=100, gamma=2.7, mean_degree=5)
+
+    Create a geometric soft configuration graph with 100 nodes. The $\beta$ parameter
+    is set to 1.5 and the exponent of the powerlaw distribution of the hidden
+    degrees is 2.7 with mean value of 5.
+
+    Generate a network with predefined hidden degrees:
+
+    >>> kappas = {i: 10 for i in range(100)}
+    >>> G = nx.geometric_soft_configuration_graph(beta=2.5, kappas=kappas)
+
+    Create a geometric soft configuration graph with 100 nodes. The $\beta$ parameter
+    is set to 2.5 and all nodes with hidden degree $\kappa=10$.
+
+
+    References
+    ----------
+    .. [1] Serrano, M. Á., Krioukov, D., & Boguñá, M. (2008). Self-similarity
+       of complex networks and hidden metric spaces. Physical review letters, 100(7), 078701.
+
+    .. [2] van der Kolk, J., Serrano, M. Á., & Boguñá, M. (2022). An anomalous
+       topological phase transition in spatial random graphs. Communications Physics, 5(1), 245.
+
+    .. [3] Krioukov, D., Papadopoulos, F., Kitsak, M., Vahdat, A., & Boguná, M. (2010).
+       Hyperbolic geometry of complex networks. Physical Review E, 82(3), 036106.
+
+    """
+    if beta <= 0:
+        raise nx.NetworkXError("The parameter beta cannot be smaller or equal to 0.")
+
+    if kappas is not None:
+        if not all((n is None, gamma is None, mean_degree is None)):
+            raise nx.NetworkXError(
+                "When kappas is input, n, gamma and mean_degree must not be."
+            )
+
+        n = len(kappas)
+        mean_degree = sum(kappas) / len(kappas)
+    else:
+        if any((n is None, gamma is None, mean_degree is None)):
+            raise nx.NetworkXError(
+                "Please provide either kappas, or all 3 of: n, gamma and mean_degree."
+            )
+
+        # Generate `n` hidden degrees from a powerlaw distribution
+        # with given exponent `gamma` and mean value `mean_degree`
+        gam_ratio = (gamma - 2) / (gamma - 1)
+        kappa_0 = mean_degree * gam_ratio * (1 - 1 / n) / (1 - 1 / n**gam_ratio)
+        base = 1 - 1 / n
+        power = 1 / (1 - gamma)
+        kappas = {i: kappa_0 * (1 - seed.random() * base) ** power for i in range(n)}
+
+    G = nx.Graph()
+    R = n / (2 * math.pi)
+
+    # Approximate values for mu in the thermodynamic limit (when n -> infinity)
+    if beta > 1:
+        mu = beta * math.sin(math.pi / beta) / (2 * math.pi * mean_degree)
+    elif beta == 1:
+        mu = 1 / (2 * mean_degree * math.log(n))
+    else:
+        mu = (1 - beta) / (2**beta * mean_degree * n ** (1 - beta))
+
+    # Generate random positions on a circle
+    thetas = {k: seed.uniform(0, 2 * math.pi) for k in kappas}
+
+    for u in kappas:
+        for v in list(G):
+            angle = math.pi - math.fabs(math.pi - math.fabs(thetas[u] - thetas[v]))
+            dij = math.pow(R * angle, beta)
+            mu_kappas = math.pow(mu * kappas[u] * kappas[v], max(1, beta))
+            p_ij = 1 / (1 + dij / mu_kappas)
+
+            # Create an edge with a certain connection probability
+            if seed.random() < p_ij:
+                G.add_edge(u, v)
+        G.add_node(u)
+
+    nx.set_node_attributes(G, thetas, "theta")
+    nx.set_node_attributes(G, kappas, "kappa")
+
+    # Map hidden degrees into the radial coordiantes
+    zeta = 1 if beta > 1 else 1 / beta
+    kappa_min = min(kappas.values())
+    R_c = 2 * max(1, beta) / (beta * zeta)
+    R_hat = (2 / zeta) * math.log(n / math.pi) - R_c * math.log(mu * kappa_min)
+    radii = {node: R_hat - R_c * math.log(kappa) for node, kappa in kappas.items()}
+    nx.set_node_attributes(G, radii, "radius")
+
+    return G

networkx/generators/tests/test_geometric.py

@@ -366,3 +366,127 @@ def test_geometric_edges_raises_no_pos():
     msg = "all nodes. must have a '"
     with pytest.raises(nx.NetworkXError, match=msg):
         nx.geometric_edges(G, radius=1)
+
+
+def test_number_of_nodes_S1():
+    G = nx.geometric_soft_configuration_graph(
+        beta=1.5, n=100, gamma=2.7, mean_degree=10, seed=42
+    )
+    assert len(G) == 100
+
+
+def test_set_attributes_S1():
+    G = nx.geometric_soft_configuration_graph(
+        beta=1.5, n=100, gamma=2.7, mean_degree=10, seed=42
+    )
+    kappas = nx.get_node_attributes(G, "kappa")
+    assert len(kappas) == 100
+    thetas = nx.get_node_attributes(G, "theta")
+    assert len(thetas) == 100
+    radii = nx.get_node_attributes(G, "radius")
+    assert len(radii) == 100
+
+
+def test_mean_kappas_S1():
+    G = nx.geometric_soft_configuration_graph(
+        beta=2.5, n=5000, gamma=2.7, mean_degree=10, seed=42
+    )
+    kappas = nx.get_node_attributes(G, "kappa")
+    mean_kappas = sum(kappas.values()) / len(kappas)
+    assert math.fabs(mean_kappas - 10) < 0.5
+
+
+def test_mean_degree_S1():
+    G = nx.geometric_soft_configuration_graph(
+        beta=2.5, n=5000, gamma=2.7, mean_degree=10, seed=42
+    )
+    degrees = dict(G.degree())
+    mean_degree = sum(degrees.values()) / len(degrees)
+    assert math.fabs(mean_degree - 10) < 1
+
+
+def test_dict_kappas_S1():
+    kappas = {i: 10 for i in range(1000)}
+    G = nx.geometric_soft_configuration_graph(beta=1, kappas=kappas)
+    assert len(G) == 1000
+    kappas = nx.get_node_attributes(G, "kappa")
+    assert all(kappa == 10 for kappa in kappas.values())
+
+
+def test_beta_clustering_S1():
+    G1 = nx.geometric_soft_configuration_graph(
+        beta=1.5, n=100, gamma=3.5, mean_degree=10, seed=42
+    )
+    G2 = nx.geometric_soft_configuration_graph(
+        beta=3.0, n=100, gamma=3.5, mean_degree=10, seed=42
+    )
+    assert nx.average_clustering(G1) < nx.average_clustering(G2)
+
+
+def test_wrong_parameters_S1():
+    with pytest.raises(
+        nx.NetworkXError,
+        match="Please provide either kappas, or all 3 of: n, gamma and mean_degree.",
+    ):
+        G = nx.geometric_soft_configuration_graph(
+            beta=1.5, gamma=3.5, mean_degree=10, seed=42
+        )
+
+    with pytest.raises(
+        nx.NetworkXError,
+        match="When kappas is input, n, gamma and mean_degree must not be.",
+    ):
+        kappas = {i: 10 for i in range(1000)}
+        G = nx.geometric_soft_configuration_graph(
+            beta=1.5, kappas=kappas, gamma=2.3, seed=42
+        )
+
+    with pytest.raises(
+        nx.NetworkXError,
+        match="Please provide either kappas, or all 3 of: n, gamma and mean_degree.",
+    ):
+        G = nx.geometric_soft_configuration_graph(beta=1.5, seed=42)
+
+
+def test_negative_beta_S1():
+    with pytest.raises(
+        nx.NetworkXError, match="The parameter beta cannot be smaller or equal to 0."
+    ):
+        G = nx.geometric_soft_configuration_graph(
+            beta=-1, n=100, gamma=2.3, mean_degree=10, seed=42
+        )
+
+
+def test_non_zero_clustering_beta_lower_one_S1():
+    G = nx.geometric_soft_configuration_graph(
+        beta=0.5, n=100, gamma=3.5, mean_degree=10, seed=42
+    )
+    assert nx.average_clustering(G) > 0
+
+
+def test_mean_degree_influence_on_connectivity_S1():
+    low_mean_degree = 2
+    high_mean_degree = 20
+    G_low = nx.geometric_soft_configuration_graph(
+        beta=1.2, n=100, gamma=2.7, mean_degree=low_mean_degree, seed=42
+    )
+    G_high = nx.geometric_soft_configuration_graph(
+        beta=1.2, n=100, gamma=2.7, mean_degree=high_mean_degree, seed=42
+    )
+    assert nx.number_connected_components(G_low) > nx.number_connected_components(
+        G_high
+    )
+
+
+def test_compare_mean_kappas_different_gammas_S1():
+    G1 = nx.geometric_soft_configuration_graph(
+        beta=1.5, n=2000, gamma=2.7, mean_degree=20, seed=42
+    )
+    G2 = nx.geometric_soft_configuration_graph(
+        beta=1.5, n=2000, gamma=3.5, mean_degree=20, seed=42
+    )
+    kappas1 = nx.get_node_attributes(G1, "kappa")
+    mean_kappas1 = sum(kappas1.values()) / len(kappas1)
+    kappas2 = nx.get_node_attributes(G2, "kappa")
+    mean_kappas2 = sum(kappas2.values()) / len(kappas2)
+    assert math.fabs(mean_kappas1 - mean_kappas2) < 1