~tfardet/nngt-developers

NNGT: Generation - Add clustered graph v1 SUPERSEDED

~tfardet
~tfardet: 1
 Generation - Add clustered graph

 9 files changed, 757 insertions(+), 81 deletions(-)
#567906 .build.yml success
Export patchset (mbox)
How do I use this?

Copy & paste the following snippet into your terminal to import this patchset into git:

curl -s https://lists.sr.ht/~tfardet/nngt-developers/patches/24436/mbox | git am -3
Learn more about email & git
View this thread in the archives

[PATCH NNGT] Generation - Add clustered graph Export this patch

~tfardet
From: Tanguy Fardet <tanguyfardet@protonmail.com>

Add new `sparse_clustered` function and associated rewiring scheme.
Check empty edges for delete functions.
Slightly improved attribute type support.
---
 nngt/core/gt_graph.py                 |  24 +-
 nngt/core/ig_graph.py                 |  25 +-
 nngt/core/nngt_graph.py               |  68 ++--
 nngt/core/nx_graph.py                 |  21 +-
 nngt/generation/graph_connectivity.py | 484 +++++++++++++++++++++++++-
 nngt/generation/rewiring.py           |  30 +-
 nngt/lib/connect_tools.py             |  92 +++++
 testing/test_generation2.py           |  42 +++
 testing/test_rewire.py                |  52 +++
 9 files changed, 757 insertions(+), 81 deletions(-)

diff --git a/nngt/core/gt_graph.py b/nngt/core/gt_graph.py
index 00f8577..2dbebb0 100755
--- a/nngt/core/gt_graph.py
+++ b/nngt/core/gt_graph.py
@@ -90,10 +90,10 @@ class _GtNProperty(BaseProperty):
        g     = self.parent()._graph

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
                dtype = int
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
                dtype = float
            elif value_type == "string":
@@ -307,9 +307,9 @@ class _GtEProperty(BaseProperty):
            self._num_values_set[name] = num_edges

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
            elif value_type == "string":
                val = ""
@@ -829,17 +829,17 @@ class _GtGraph(GraphInterface):

    def delete_edges(self, edges):
        ''' Remove a list of edges '''
        g = self._graph
        if len(edges):
            g = self._graph

        Edge = g.edge
            Edge = g.edge

        if nonstring_container(edges[0]):
            # fast loop
            [self._graph.remove_edge(Edge(*e)) for e in edges]
        else:
            self._graph.remove_edge(Edge(*edges))
            if nonstring_container(edges[0]):
                # fast loop
                [self._graph.remove_edge(Edge(*e)) for e in edges]
            else:
                self._graph.remove_edge(Edge(*edges))

        if len(edges):
            self._eattr.edges_deleted()

            self._edges_deleted = True
diff --git a/nngt/core/ig_graph.py b/nngt/core/ig_graph.py
index 93599fe..7b5f4af 100755
--- a/nngt/core/ig_graph.py
+++ b/nngt/core/ig_graph.py
@@ -79,9 +79,9 @@ class _IgNProperty(BaseProperty):
        g = self.parent()._graph

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
            elif value_type == "string":
                val = ""
@@ -196,9 +196,9 @@ class _IgEProperty(BaseProperty):
        g = self.parent()._graph

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
            elif value_type == "string":
                val = ""
@@ -608,16 +608,17 @@ class _IGraph(GraphInterface):

    def delete_edges(self, edges):
        ''' Remove a list of edges '''
        if nonstring_container(edges[0]):
            if isinstance(edges[0], tuple):
                self._graph.delete_edges(edges)
        if len(edges):
            if nonstring_container(edges[0]):
                if isinstance(edges[0], tuple):
                    self._graph.delete_edges(edges)
                else:
                    self._graph.delete_edges([tuple(e) for e in edges])
            else:
                self._graph.delete_edges([tuple(e) for e in edges])
        else:
            self._graph.delete_edges([edges])
                self._graph.delete_edges([edges])

        for key in self._eattr:
            self._eattr._num_values_set[key] = self.edge_nb()
            for key in self._eattr:
                self._eattr._num_values_set[key] = self.edge_nb()

    def clear_all_edges(self):
        ''' Remove all edges from the graph '''
diff --git a/nngt/core/nngt_graph.py b/nngt/core/nngt_graph.py
index d3853dd..8d0238f 100644
--- a/nngt/core/nngt_graph.py
+++ b/nngt/core/nngt_graph.py
@@ -74,11 +74,12 @@ class _NProperty(BaseProperty):

    def new_attribute(self, name, value_type, values=None, val=None):
        dtype = object

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
                dtype = int
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
                dtype = float
            elif value_type == "string":
@@ -278,9 +279,9 @@ class _EProperty(BaseProperty):
            self._num_values_set[name] = num_edges

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
            elif value_type == "string":
                val = ""
@@ -869,48 +870,49 @@ class _NNGTGraph(GraphInterface):

    def delete_edges(self, edges):
        ''' Remove a list of edges '''
        g = self._graph
        if len(edges):
            g = self._graph

        old_enum = len(g._unique)
            old_enum = len(g._unique)

        if not nonstring_container(edges[0]):
            edges = [edges]
            if not nonstring_container(edges[0]):
                edges = [edges]

        if not isinstance(edges[0], tuple):
            edges = [tuple(e) for e in edges]
            if not isinstance(edges[0], tuple):
                edges = [tuple(e) for e in edges]

        # get edge ids
        e_to_eid = g._unique
            # get edge ids
            e_to_eid = g._unique

        eids = {e_to_eid[e] for e in edges}
            eids = {e_to_eid[e] for e in edges}

        # remove
        directed = g._directed
            # remove
            directed = g._directed

        for e in edges:
            if e in g._unique:
                del g._unique[e]
            for e in edges:
                if e in g._unique:
                    del g._unique[e]

            if not directed:
                if e in g._edges:
                    del g._edges[e]
                    del g._edges[e[::-1]]
                if not directed:
                    if e in g._edges:
                        del g._edges[e]
                        del g._edges[e[::-1]]

            # update in and out degrees
            s, t = e
                # update in and out degrees
                s, t = e

            g._in_deg[t] -= 1
            g._out_deg[s] -= 1
                g._in_deg[t] -= 1
                g._out_deg[s] -= 1

        # reset eids
        for i, e in enumerate(g._unique):
            g._unique[e] = i
            # reset eids
            for i, e in enumerate(g._unique):
                g._unique[e] = i

            if not directed:
                e._edges[e] = i
                e._edges[e[::-1]] = i
                if not directed:
                    e._edges[e] = i
                    e._edges[e[::-1]] = i

        self._eattr.edges_deleted(eids)
            self._eattr.edges_deleted(eids)

    def clear_all_edges(self):
        g = self._graph
diff --git a/nngt/core/nx_graph.py b/nngt/core/nx_graph.py
index 43e6d3d..542bf89 100755
--- a/nngt/core/nx_graph.py
+++ b/nngt/core/nx_graph.py
@@ -84,9 +84,9 @@ class _NxNProperty(BaseProperty):
        g = self.parent()._graph

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
            elif value_type == "string":
                val = ""
@@ -225,9 +225,9 @@ class _NxEProperty(BaseProperty):
        g = self.parent()._graph

        if val is None:
            if value_type == "int":
            if value_type in ("int", "integer"):
                val = int(0)
            elif value_type == "double":
            elif value_type in ("double", "float"):
                val = np.NaN
            elif value_type == "string":
                val = ""
@@ -736,13 +736,14 @@ class _NxGraph(GraphInterface):

    def delete_edges(self, edges):
        ''' Remove a list of edges '''
        if nonstring_container(edges[0]):
            self._graph.remove_edges_from(edges)
        else:
            self._graph.remove_edge(*edges)
        if len(edges):
            if nonstring_container(edges[0]):
                self._graph.remove_edges_from(edges)
            else:
                self._graph.remove_edge(*edges)

        for key in self._eattr:
            self._eattr._num_values_set[key] = self.edge_nb()
            for key in self._eattr:
                self._eattr._num_values_set[key] = self.edge_nb()

    def clear_all_edges(self):
        ''' Remove all edges from the graph '''
diff --git a/nngt/generation/graph_connectivity.py b/nngt/generation/graph_connectivity.py
index 4643d6f..c57f467 100755
--- a/nngt/generation/graph_connectivity.py
+++ b/nngt/generation/graph_connectivity.py
@@ -30,7 +30,8 @@ import numpy as np

import nngt
from nngt.geometry.geom_utils import conversion_magnitude
from nngt.lib.connect_tools import _set_options
from nngt.lib.connect_tools import (_set_options, _connect_ccs,
                                    _independent_edges)
from nngt.lib.logger import _log_message
from nngt.lib.test_functions import (mpi_checker, mpi_random, deprecated,
                                     on_master_process)
@@ -45,8 +46,9 @@ __all__ = [
    'from_degree_list',
    'gaussian_degree',
    'newman_watts',
    'random_scale_free',
    'price_scale_free',
    'random_scale_free',
    'sparse_clustered',
    'watts_strogatz',
]

@@ -190,17 +192,14 @@ def from_degree_list(degrees, degree_type='in', weighted=True,
        The list of degrees for each node in the graph.
    degree_type : str, optional (default: 'in')
        The type of the fixed degree, among ``'in'``, ``'out'`` or ``'total'``.
        @todo `'total'` not implemented yet.
    nodes : int, optional (default: None)
        The number of nodes in the graph.
    weighted : bool, optional (default: True)
        Whether the graph edges have weights.
    directed : bool, optional (default: True)
        @todo: only for directed graphs for now. Whether the graph is directed
        or not.
        Whether the graph is directed or not.
    multigraph : bool, optional (default: False)
        Whether the graph can contain multiple edges between two
        nodes.
        Whether the graph can contain multiple edges between two nodes.
    name : string, optional (default: "ER")
        Name of the created graph.
    shape : :class:`~nngt.geometry.Shape`, optional (default: None)
@@ -269,10 +268,6 @@ def fixed_degree(degree, degree_type='in', nodes=0, reciprocity=-1.,
        The value of the constant degree.
    degree_type : str, optional (default: 'in')
        The type of the fixed degree, among ``'in'``, ``'out'`` or ``'total'``.

        @todo
            `'total'` not implemented yet.

    nodes : int, optional (default: None)
        The number of nodes in the graph.
    reciprocity : double, optional (default: -1 to let it free)
@@ -989,6 +984,472 @@ def watts_strogatz(coord_nb, proba_shortcut=None, reciprocity_circular=1.,
    return graph_nw


def sparse_clustered(c, nodes=0, edges=None, avg_deg=None, connected=True,
                     rtol=None, exact_edge_nb=False, weighted=True,
                     directed=True, multigraph=False, name="FC", shape=None,
                     positions=None, population=None, from_graph=None,
                     **kwargs):
    r'''
    Generate a sparse random graph with given average clustering coefficient
    and degree.

    .. versionadded:: 2.5

    The original algorithm is adapted from [newman-clustered-2003]_ and leads
    to a graph with approximate clustering coefficient and number of edges.

    .. warning::

        This algorithm can only give reasonable results for sparse graphs and
        will raise an error if the requested graph density is above `c`.

    Nodes are distributed among :math:`\mu` overlapping groups of size
    :math:`\nu` and, each time two nodes belong to a common group, they are
    connected with a probability :math:`p = c`.

    For sparse graphs, the average (in/out-)degree can be approximmated as
    :math:`k = \mu p(\nu - 1)`, and the average clustering as:

    .. math::

        C^{(u)} = \frac{p \left[ p(\nu - 1) - 1 \right]}{k - 1}

    for undirected graphs and

    .. math::

        C^{(d)} = \frac{p\mu \left[ p(2\nu - 3) - 1 \right]}{2k - 1 - p}

    for all clustering modes in directed graphs.

    From these relations, we compute :math:`\mu` and :math:`\nu` as:

    .. math::

        \nu^{(u)} = 1 + \frac{1}{p} + \frac{C^{(u)} (k - 1)}{p^2}

    or

    .. math::

        \nu^{(d)} = \frac{3}{2} + \frac{1}{2p} +
                    \frac{C^{(u)} (2k - 1 - p)}{2p^2}

    and

    .. math::

        \mu = \frac{k}{p (\nu - 1)}.

    Parameters
    ----------
    c : float
        Desired value for the final average clustering in the graph.
    nodes : int, optional (default: None)
        The number of nodes in the graph.
    edges : int, optional
        The number of edges between the nodes
    avg_deg : double, optional
        Average degree of the neurons given by `edges / nodes`.
    connected : bool, optional (default: True)
        Whether the resulting graph must be connected (True) or may be
        unconnected (False).
    rtol : float, optional (default: not checked)
        Tolerance on the relative error between the target clustering `c` and
        the actual clustering of the final graph.
        If the algorithm leads to a relative error greater than `rtol`, then
        an error is raised.
    exact_edge_nb : bool, optional (default: False)
        Whether the final graph should have precisely the number of edges
        required.
    weighted : bool, optional (default: True)
        Whether the graph edges have weights.
    directed : bool, optional (default: True)
        Whether the graph should be directed or not.
    multigraph : bool, optional (default: False)
        Whether the graph can contain multiple edges between two nodes.
    name : string, optional (default: "ER")
        Name of the created graph.
    shape : :class:`~nngt.geometry.Shape`, optional (default: None)
        Shape of the neurons' environment.
    positions : :class:`numpy.ndarray`, optional (default: None)
        A 2D or 3D array containing the positions of the neurons in space.
    population : :class:`~nngt.NeuralPop`, optional (default: None)
        Population of neurons defining their biological properties (to create a
        :class:`~nngt.Network`).
    from_graph : :class:`Graph` or subclass, optional (default: None)
        Initial graph whose nodes are to be connected.
    **kwargs : keyword arguments
        If connected is True, `method` can be passed to define how the
        components should be connected among themselves. Available methods are
        "sequential", where the components are connected sequentially, forming
        a long thread and increasing the graph's diameter, "star", where all
        components are connected to the largest one, and "random", where
        components are connected randomly.
        If not provided, defaults to "random".

    Note
    ----
    `nodes` is required unless `from_graph` or `population` is provided.
    If `from_graph` is provided, all preexistent edges in the
    object will be deleted before the new connectivity is implemented.

    Returns
    -------
    graph_fc : :class:`~nngt.Graph`, or subclass
        A new generated graph or the modified `from_graph`.

    References
    ----------

    .. [newman-clustered-2003] Newman, M. E. J. Properties of Highly
        Clustered Networks, Phys. Rev. E 2003 68 (2).
        :doi:`10.1103/PhysRevE.68.026121`, :arxiv:`cond-mat/0303183`.
    '''
    # get method
    method = "star"

    if "method" in kwargs:
        method = kwargs.pop("method")

    # set node number and library graph
    graph_fc = from_graph

    if graph_fc is not None:
        nodes = graph_fc.node_nb()
        graph_fc.clear_all_edges()
    else:
        nodes = population.size if population is not None else nodes
        graph_fc = nngt.Graph(
            name=name, nodes=nodes, directed=directed, **kwargs)

    _set_options(graph_fc, population, shape, positions)

    # compute average degree and set the number of nodes
    if (avg_deg is None and edges is None) or None not in (avg_deg, edges):
        raise ValueError("Either `avg_deg` or `edges` should be provided.")

    # add edges
    ia_edges = None

    if nodes > 1:
        k = edges / nodes if avg_deg is None else avg_deg

        num_edges = int(k*nodes) if edges is None else edges

        if k/nodes > c:
            raise ValueError('Clustering cannot be lower than avg_deg/nodes, '
                             'or, equivalently, edges/nodes**2. Got '
                             'c = {} against {}.'.format(c, k/nodes))

        def func_nu(p, directed):
            if directed:
                return np.around(
                    1.5 + 1/(2*p) + (2*k - 1 - p) * c / (2*p**2)).astype(int)

            return np.around(1 + 1/p + (k - 1) * c / p**2).astype(int)

        def func_mu(nu, p):
            return k / (p * (nu - 1))

        p = c
        nu = func_nu(p, directed)

        if nu >= nodes:
            raise ValueError('Theoretical group size computed was '
                             'greater than the number of nodes.')

        mu = func_mu(nu, p)

        # assign the neurons to groups
        n_to_g = {}
        g_to_n = {}

        rng = nngt._rng

        gid = 0

        # to minimize the number of groups, we start with distributed groups
        for i in range(int(nodes / nu)):
            g_to_n[gid] = list(range(nu*i, nu*(i+1)))
            for j in range(nu*i, nu*(i+1)):
                if j in n_to_g:
                    n_to_g[j].add(gid)
                else:
                    n_to_g[j] = {gid}

            gid += 1

        if int(nodes / nu) * nu < nodes:
            s = set(range(nu*gid, nodes))

            # to minimize the number of disconnected components, we are going
            # to fill the remaining nodes of s with previous nodes (one from
            # each previous set until s is full) and assign a new random node
            # to each of these previous sets

            for g, ids in g_to_n.items():
                # get a node from nodes
                u = ids[0]
                group = set(ids) - {u}

                n_to_g[u] -= {g}

                while len(group) < nu:
                    v = rng.choice(nodes, 1)[0]

                    group.add(v)

                    if len(group) == nu:
                        if v in n_to_g:
                            n_to_g[v] = n_to_g[v].union({g})
                        else:
                            n_to_g[v] = {g}

                assert len(group) == nu

                g_to_n[g] = list(group)

                s.add(u)

                if len(s) == nu:
                    break

            g_to_n[gid] = list(s)

            for j in s:
                if j in n_to_g:
                    n_to_g[j].add(gid)
                else:
                    n_to_g[j] = {gid}

            gid += 1

        # check that all nodes belong to a group
        assert len(n_to_g) == nodes

        if nu < nodes or mu > 1:
            while gid < mu*nodes/nu:
                ids = rng.choice(nodes, nu, replace=False)

                g_to_n[gid] = ids

                for i in ids:
                    if i in n_to_g:
                        n_to_g[i].add(gid)
                    else:
                        n_to_g[i] = {gid}

                gid += 1

        # generate the edges
        edges = set()

        for group, ids in g_to_n.items():
            probas = rng.random(nu*(nu - 1))

            count = 0

            for i, u in enumerate(ids):
                targets = list(ids[i + 1:])

                if directed:
                    targets += list(ids[:i])
                for v in targets:
                    prob = probas[count]
                    count += 1

                    if directed:
                        if prob < p and (u, v) not in edges:
                            edges.add((u, v))
                    elif prob < p:
                        if (u, v) not in edges and (v, u) not in edges:
                            edges.add((u, v))

        graph_fc.new_edges(list(edges))

        # final optional checks for connectedness and edge number
        bridges = set()

        cc, cmean = None, None

        if connected:
            # make the graph fully connected
            cc = nngt.analysis.local_clustering(graph_fc)

            cmean = cc.mean()

            disconnected = True

            while disconnected:
                cids, hist = nngt.analysis.connected_components(graph_fc)

                disconnected = len(hist) > 1

                if not disconnected:
                    break

                # connect nodes from separate connected components
                if method == "star":
                    idx = np.argmax(hist)

                    for i in (set(range(len(hist))) - {idx}):
                        _connect_ccs(graph_fc, cids, idx, i, cc, c, cmean,
                                     edges, bridges)

                        if directed:
                            _connect_ccs(graph_fc, cids, i, idx, cc, c, cmean,
                                         edges, bridges)
                elif method == "sequential":
                    for i in range(1, len(hist)):
                        _connect_ccs(graph_fc, cids, i - 1, i, cc, c, cmean,
                                     edges, bridges)

                    if directed:
                        _connect_ccs(graph_fc, cids, i, 0, cc, c, cmean, edges,
                                     bridges)
                elif method == "random":
                    n = len(hist)

                    # almost surely connected random graph requires
                    # n*log(n) connections so we ask for half that to use as
                    # few edges as possible while limiting the number of loops
                    # since each iteration requires to compute the CCs
                    mult = max(1, int(0.5*np.ceil(np.log(n))))
                    sources = list(range(n))*mult

                    if directed and n > 2:
                        sources *= 2

                    targets = sources[::-1]

                    if n > 2:
                        rng.shuffle(targets)
                        
                    cset = set()

                    for i in sources:
                        j = i

                        while i == j or (i, j) in cset:
                            j = rng.integers(0, len(hist))

                        _connect_ccs(graph_fc, cids, i, j, cc, c, cmean,
                                     edges, bridges)

                        cset.add((i, j))
                else:
                    raise ValueError("Invalid `method`: '" + method + "'.")
        else:
            cc = nngt.analysis.local_clustering(graph_fc)

            cmean = cc.mean()

        if exact_edge_nb:
            proba = 1 - cc if cmean < c else cc
            proba /= proba.sum()

            while graph_fc.edge_nb() < num_edges:
                # add edges
                sources, targets = None, None

                missing = num_edges - graph_fc.edge_nb()

                if missing > 1:
                    sources = rng.choice(nodes, missing, p=proba)
                    targets = sources.copy()

                    rng.shuffle(targets)
                else:
                    sources = [rng.choice(nodes)]
                    targets = [rng.choice(nodes)]

                new_edges = []

                for u, v in zip(sources, targets):
                    cond = (u, v) not in edges

                    if not directed:
                        cond *= (v, u) not in edges

                    if u != v and cond:
                        new_edges.append((u, v))
                        edges.add((u, v))

                graph_fc.new_edges(new_edges)

            check_ind_edges = True
            count = 0

            while graph_fc.edge_nb() > num_edges and count < 1000:
                # remove edges
                remove = graph_fc.edge_nb() - num_edges

                sources, targets, chosen = [], [], None

                # remove edges that do not belong to triangles if possible
                if check_ind_edges:
                    ind_edges = \
                        _independent_edges(graph_fc).difference(bridges)

                    if ind_edges:
                        ind_edges = np.asarray(list(ind_edges))

                        if len(ind_edges) > remove:
                            chosen = rng.choice(len(ind_edges), remove,
                                                replace=False)
                        else:
                            chosen = slice(None)

                        sources = list(ind_edges[chosen, 0])
                        targets = list(ind_edges[chosen, 1])
                    else:
                        check_ind_edges = False

                if len(sources) < remove:
                    remaining = remove - len(sources)
                    proba = np.square(
                        graph_fc.get_degrees("out")).astype(float)
                    proba /= proba.sum()
                    sources.extend(rng.choice(nodes, remaining, p=proba))

                    proba = np.square(graph_fc.get_degrees("in")).astype(float)
                    proba /= proba.sum()
                    targets.extend(rng.choice(nodes, remaining, p=proba))

                delete = []

                for u, v in zip(sources, targets):
                    if (u, v) in edges and (u, v) not in bridges:
                        delete.append((u, v))
                        edges.discard((u, v))
                    elif not directed:
                        if (v, u) in edges and (v, u) not in bridges:
                            delete.append((v, u))
                            edges.discard((v, u))

                graph_fc.delete_edges(delete)

                count += 1

            if count == 1000:
                raise RuntimeError("Algorithm did not converge, try with "
                                   "exact_edge_nb=False to avoid that issue.")

    if rtol is not None:
        cmean = nngt.analysis.local_clustering(graph_fc).mean()

        err = np.abs(c - cmean) / c

        if err > rtol:
            raise RuntimeError("Relative error greater than `rtol`: "
                               "{} > {}".format(err, rtol))

    graph_fc._graph_type = "fixed_clustering"

    return graph_fc


# --------------------- #
# Distance-based models #
# --------------------- #
@@ -1110,6 +1571,7 @@ _di_generator = {
    "newman_watts": newman_watts,
    "price_scale_free": price_scale_free,
    "random_scale_free": random_scale_free,
    "sparse_clustered": sparse_clustered,
    "watts_strogatz": watts_strogatz,
}

diff --git a/nngt/generation/rewiring.py b/nngt/generation/rewiring.py
index ba620cd..921413b 100644
--- a/nngt/generation/rewiring.py
+++ b/nngt/generation/rewiring.py
@@ -236,7 +236,7 @@ def lattice_rewire(g, target_reciprocity=1., node_attr_constraints=None,


def random_rewire(g, constraints=None, node_attr_constraints=None,
                  edge_attr_constraints=None):
                  edge_attr_constraints=None, **kwargs):
    '''
    Generate a new rewired graph from `g`.

@@ -267,6 +267,23 @@ def random_rewire(g, constraints=None, node_attr_constraints=None,
        attributes are randomized by groups (all attributes of a given edge are
        sent to the same new edge). By default, attributes are completely and
        separately randomized.
    **kwargs : optional keyword arguments
        These are optional arguments in the case `constraints` is "clustering".
        In that case, the user can provide both:

        * rtol : float, optional (default: 5%)
          The tolerance on the relative error to the average clustering for the
          rewired graph.
        * connected : bool, optional (default: False)
          Whether the generated graph should be connected (this reduces the
          precision of the final clustering).
        * method : str, optional (default: "star")
          Defines how the initially disconnected components of the generated
          graph should be connected among themselves.
          Available methods are "sequential", where the components are
          connected sequentially, forming a long thread and increasing the
          graph's diameter, "star", where all components are connected to the
          largest one, and "random", where components are connected randomly.
    '''
    directed  = g.is_directed()
    num_nodes = g.node_nb()
@@ -305,8 +322,15 @@ def random_rewire(g, constraints=None, node_attr_constraints=None,
        new_graph = gc.from_degree_list(degrees, constraints,
                                        directed=directed)
    elif constraints == "clustering":
        raise NotImplementedError("Rewiring with constrained clustering is "
                                  "not yet available.")
        rtol = kwargs.get("rtol", 0.05)
        connected = kwargs.get("connected", False)
        method = kwargs.get("method", "star")

        c = nngt.analysis.local_clustering(g).mean()

        new_graph = gc.sparse_clustered(
            c, edges=g.edge_nb(), nodes=num_nodes, connected=connected,
            method=method, exact_edge_nb=True, rtol=rtol)

    rng = nngt._rng

diff --git a/nngt/lib/connect_tools.py b/nngt/lib/connect_tools.py
index db4d483..178e4e1 100755
--- a/nngt/lib/connect_tools.py
+++ b/nngt/lib/connect_tools.py
@@ -409,3 +409,95 @@ def _set_default_edge_attributes(g, attributes, num_edges):
                attributes[k] = [0 for _ in range(num_edges)]
            elif not skip:
                attributes[k] = [None for _ in range(num_edges)]


def _connect_ccs(g, cids, cc1, cc2, cc, c, cmean, edges, bridges):
    '''
    Connect connected components together.

    Parameters
    ----------
    g : graph
    cids : index of the connected component to which each node belongs
    cc1 : first connected component index
    cc2 : second connected component index
    cc : list of local clustering coefficients
    c : target average clustering
    cmean : current average clustering
    edges : set of edges to update
    bridges : set of edges bridging connected components
    '''
    # we connect nodes with the lowest/highest clustering
    select = np.argmin if cmean < c else np.argmax

    # find suitable nodes
    other_nodes = np.where(cids == cc1)[0]
    inodes = np.where(cids == cc2)[0]

    def compare(a, b, cmean, c):
        if cmean < c:
            return a < b

        return a > b

    # get those with desired clustering
    change_nodes = True

    while change_nodes:
        u = other_nodes[select(cc[other_nodes])]
        v = inodes[select(cc[inodes])]

        change_nodes = (u, v) in edges or (u, v) in bridges or u == v

        if not g.is_directed():
            change_nodes += (v, u) in edges or (v, u) in bridges

        if change_nodes:
            # remove the node with the larger or smaller cc based on cmean & c
            if compare(cc[u], cc[v], cmean, c) and len(inodes) > 1:
                inodes = list(set(inodes) - {v})
                v = inodes[select(cc[inodes])]
            elif len(other_nodes) > 1:
                other_nodes = list(set(other_nodes) - {u})
                u = other_nodes[select(cc[other_nodes])]
            else:
                raise RuntimeError("Algorithm did not converge.")

    if (u, v) not in edges:
        g.new_edge(u, v)
        edges.add((u, v))
        bridges.add((u, v))


def _independent_edges(g):
    '''
    Return the independent edges in the graph, i.e. the those that do not
    participate to any triangle.

    Parameters
    ----------
    g : :class:`~nngt.Graph`
        Graph to analyze.
    '''
    # prepare adjacency matrices
    A = g.adjacency_matrix()

    A.setdiag(0)  # remove self-loops for the computation

    A2 = A*A
    A2.data = np.ones(len(A2.data))
    A2.setdiag(1)

    ind_edges = set()

    if g.is_directed():
        [ind_edges.add((i, j)) for i in g.get_nodes()
            for j in np.where(
                A.getrow(i).todense().A1*(1 - A2.getrow(i).todense().A1))[0]]
    else:
        [ind_edges.add((i, j)) for i in g.get_nodes()
            for j in np.where(
                A.getrow(i).todense().A1*(1 - A2.getrow(i).todense().A1))[0]
            if (j, i) not in ind_edges]

    return ind_edges
diff --git a/testing/test_generation2.py b/testing/test_generation2.py
index b3641e2..f8ab9eb 100644
--- a/testing/test_generation2.py
+++ b/testing/test_generation2.py
@@ -524,6 +524,47 @@ def test_circular():
    assert np.isclose(na.reciprocity(gc), recip, 1e-3)


@pytest.mark.mpi_skip
def test_sparse_clustered():
    ccs = np.linspace(0.1, 0.9, 4)

    num_nodes = 500

    degrees = [10, 40]

    methods = ["star", "sequential", "random", "star"]

    for directed in (True, False):
        # check errors
        with pytest.raises(ValueError):
            g = ng.sparse_clustered(0, nodes=num_nodes, avg_deg=10,
                                    directed=directed)

        with pytest.raises(RuntimeError):
            g = ng.sparse_clustered(1, nodes=num_nodes, avg_deg=10,
                                    directed=directed, rtol=1e-10)

        # check graphs
        for i, c in enumerate(ccs):
            for deg in degrees:
                if c*num_nodes > deg:
                    g = ng.sparse_clustered(
                        c, nodes=num_nodes, avg_deg=deg, connected=False,
                        directed=directed, rtol=0.08)

                    g = ng.sparse_clustered(
                        c, nodes=num_nodes, avg_deg=deg, directed=directed,
                        exact_edge_nb=True)

                    assert g.edge_nb() == deg*num_nodes

                    g = ng.sparse_clustered(
                        c, nodes=num_nodes, avg_deg=deg, directed=directed,
                        connected=True, method=methods[i])

                    assert g.is_connected()


if __name__ == "__main__":
    if not nngt.get_config("mpi"):
        test_circular()
@@ -535,6 +576,7 @@ if __name__ == "__main__":
        test_distances()
        test_price()
        test_connect_switch_distance_rule_max_proba()
        test_sparse_clustered()

    if nngt.get_config("mpi"):
        test_mpi_from_degree_list()
diff --git a/testing/test_rewire.py b/testing/test_rewire.py
index 261e74f..924ffae 100644
--- a/testing/test_rewire.py
+++ b/testing/test_rewire.py
@@ -120,6 +120,57 @@ def test_random_rewire():
                    assert set(old_attr) == set(new_attr)


@pytest.mark.mpi_skip
def test_clst_rewire():
    ''' Test rewire preserving clustering '''
    num_nodes = 200
    coord_nb  = 8
    recip     = 0.7
    shortcut  = 0.1

    g = ng.newman_watts(coord_nb, shortcut, reciprocity_circular=recip,
                        nodes=num_nodes)

    num_edges = g.edge_nb()

    # make some node and edge attributes
    g.new_node_attribute("random_int", "int",
                         values=nngt._rng.integers(1, 2000, num_nodes))
    g.new_node_attribute("attr2", "float",
                         values=nngt._rng.uniform(size=num_nodes))

    ww = np.arange(1, num_edges + 1, dtype=float)
    g.set_weights(ww)

    g.new_edge_attribute("my-edge-attr", "int", values=-ww[::-1].astype(int))

    # rewire
    rc = ng.random_rewire(g, constraints="clustering",
                          node_attr_constraints="preserve",
                          edge_attr_constraints="together", rtol=0.1)

    r1 = ng.random_rewire(g)

    c0 = nngt.analysis.local_clustering(g).mean()
    c1 = nngt.analysis.local_clustering(r1).mean()
    cc = nngt.analysis.local_clustering(rc).mean()

    assert c0 - c1 > c0 - cc
    assert np.abs(c0 - cc) / c0 < 0.1

    # check node attributes
    assert g.node_attributes == rc.node_attributes
    assert np.array_equal(g.node_attributes["random_int"],
                          rc.node_attributes["random_int"])
    assert np.all(np.isclose(g.node_attributes["attr2"],
                             rc.node_attributes["attr2"]))
    # check that attributes were moved together
    weights  = rc.get_weights()
    my_eattr = rc.edge_attributes["my-edge-attr"]

    assert np.array_equal(my_eattr, (weights - (num_edges + 1)).astype(int))


@pytest.mark.mpi_skip
def test_complete_lattice_rewire():
    ''' Check lattice rewiring method. '''
@@ -260,6 +311,7 @@ def test_reciprocity_lattice_rewire():
if __name__ == "__main__":
    if not nngt.get_config("mpi"):
        test_random_rewire()
        test_clst_rewire()
        test_complete_lattice_rewire()
        test_incomplete_lattice_rewire()
        test_reciprocity_lattice_rewire()
-- 
2.32.0