~tfardet/nngt-developers

NNGT: Generation - Add clustered graph v3 APPLIED

~tfardet
~tfardet: 1
 Generation - Add clustered graph

 9 files changed, 790 insertions(+), 81 deletions(-)
#568954 .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/24468/mbox | git am -3
Learn more about email & git
View this thread in the archives

[PATCH NNGT v3] 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 | 505 +++++++++++++++++++++++++-
 nngt/generation/rewiring.py           |  33 +-
 nngt/lib/connect_tools.py             |  99 +++++
 testing/test_generation2.py           |  42 +++
 testing/test_rewire.py                |  54 +++
 9 files changed, 790 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..b2fffdb 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,493 @@ 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-component",
        where all components are connected to the largest one, and "random",
        where components are connected randomly, or "central-node", where one
        node from the largest component is chosen to reconnect all disconnected
        components.
        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-component"

    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-component":
                    idx = np.argmax(hist)
                    idx_nodes = np.where(cids == idx)[0]

                    select_source = lambda x: rng.choice(len(x))

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

                        if directed:
                            u_idx = np.where(idx_nodes == u)[0][0]
                            select_target = lambda x: u_idx

                            _connect_ccs(
                                graph_fc, cids, i, idx, cc, c, cmean, edges,
                                bridges, select_source=select_source,
                                select_target=select_target)
                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))
                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 == "central-node":
                    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)
                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 +1592,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..a9f3035 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,26 @@ 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-component")
          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-component", where all disconnected components
          are connected to random nodes in the largest component,
          "central-node" , where all disconnected components are connected to
          the same node in the largest component, and "random", where
          components are connected randomly.
    '''
    directed  = g.is_directed()
    num_nodes = g.node_nb()
@@ -305,8 +325,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", g.is_connected())
        method = kwargs.get("method", "star-component")

        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..5ad5dd0 100755
--- a/nngt/lib/connect_tools.py
+++ b/nngt/lib/connect_tools.py
@@ -409,3 +409,102 @@ 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,
                 select_source=None, select_target=None):
    '''
    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
    if select_source is None:
        select_source = np.argmin if cmean < c else np.argmax

    if select_target is None:
        select_target = 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_source(cc[other_nodes])]
        v = inodes[select_target(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_target(cc[inodes])]
            elif len(other_nodes) > 1:
                other_nodes = list(set(other_nodes) - {u})
                u = other_nodes[select_source(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))

    return (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..25f5cd7 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-component", "sequential", "random", "central-node"]

    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..7346e9f 100644
--- a/testing/test_rewire.py
+++ b/testing/test_rewire.py
@@ -120,6 +120,59 @@ 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
    rtol = 0.1

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

    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 < rtol

    # 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 +313,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