text: """Social distance attachment network model."""
  # pylint: disable=R0902
  # pylint: disable=arguments-differ,too-many-locals
  from collections import OrderedDict
  from itertools import product
  import numpy as np
  from scipy.sparse import csr_matrix
  from scipy.optimize import toms748
  from sdnet.networks import make_adjacency_matrix, rewire_edges
  from sdnet.utils import random_integer
  
  
  class SDA:
      """Social distance attachment network model.
  
      Attributes
      ----------
      P : (N, N) array_like
          Edge formation probability matrix.
      k : float
          Expected average node degree.
      b : float
          Characteristic length scale (positive value).
          This is the distance value at which edge formation
          probability equals ``1``.
      alpha : float
          Homophily value. It determines rate at which edge formation
          probability decreases/increases with the distance from the characteristic
          length scale `b`. Usually it should be set to a positive value.
      p_rewire : float
          Defualt probability of random edge rewiring when generating
          adjacency matrices. This is used to ensure small worldedness.
      directed : bool
          Should directed networks be generated by default.
      degseq : (N,) array_like, optional
          Default node degree sequence used in the configuration model.
      """
      def __init__(self, P, k, b, alpha, p_rewire=0.01, directed=False):
          """Initialization method."""
          self.P = P
          self.k = k
          self.b = b
          self.alpha = alpha
          self.p_rewire = p_rewire
          self.directed = directed
          self.degseq = None
  
      def __repr__(self):
          nm = self.__class__.__name__
          D = 'D' if self.directed else 'U'
          a = self.alpha
          p = self.p_rewire
          return f"<{nm} {D}{self.N} k={self.k} b={self.b:.2} a={a} p={p}>"
  
      @property
      def N(self):
          return self.P.shape[0]
  
      def sort_degseq(self, degseq):
          """Sort degseq accordingly to centrality of nodes in the social space."""
          w = self.P.sum(axis=1)
          return degseq[np.lexsort((w, w.argsort()[::-1]))]
  
      def set_degseq(self, degseq, sort=True):
          """Set default degree sequence.
  
          Parameters
          ----------
          sort : bool
              Should degree sequence be sorted to align with highest
              average edge formation probabilities per node.
          """
          degseq = degseq.copy()
          if sort and degseq is not None:
              degseq = self.sort_degseq(degseq)
          self.degseq = degseq
  
      @staticmethod
      def prob_measure(D, b, alpha):
          """Compute *SDA model* probability measure.
  
          Parameters
          ----------
          D : (N, N) array_like
              Distance matrix.
          b : float
              Characteristic length scale (positive value).
          alpha :
              Homophily value.
          """
          if b == 0:
              return np.zeros_like(D, dtype=float)
          P = 1 / (1 + (D/b)**alpha)
          np.fill_diagonal(P, 0)
          return P
  
      @classmethod
      def optim_b(cls, D, k, alpha, b_optim_min=0, b_optim_max=None,
                  xtol=10**-12, **kwds):
          """Find optimal value of the characteristic length scale `b`.
  
          Parameters
          ----------
          D : (N, N) array_like
              Distance matrix.
          b : float
              Characteristic length scale (positive value).
          alpha :
              Homophily value.
          b_optim_min : float
              Lower bound of the interval used in optimizing `b`.
          b_optim_max : float, optional
              Uppe bound of the interval used in optimizing `b`.
              If ``None`` then defaults to ``D.max()``.
          **kwds :
              Keyword parameters passed to :py:func:`scipy.optimize.brentq`.
          """
          def b_optim(b):
              P = cls.prob_measure(D, b, alpha)
              Ek = P.sum(axis=1).mean()
              return k - Ek
  
          if b_optim_max is None:
              b_optim_max = D.max()
          b = toms748(b_optim, a=b_optim_min, b=b_optim_max, xtol=xtol, **kwds)
          return b
  
      @classmethod
      def from_dist_matrix(cls, D, k, alpha, p_rewire=0.01, directed=False, **kwds):
          """Constructor method based on a distance matrix.
  
          Parameters
          ----------
          D : (N, N) array_like
              Distance matrix.
          **kwds :
              Keyword parameters passed to `optim_b`.
          """
          b = cls.optim_b(D, k, alpha, **kwds)
          P = cls.prob_measure(D, b, alpha)
          return cls(P, k, b, alpha, p_rewire=p_rewire, directed=directed)
  
      @classmethod
      def from_weighted_dist_matrices(cls, k, alpha, dm, weights=None, p_rewire=0.01,
                                      directed=False, **kwds):
          """Constructor method based on a weighted sequence of distance matrices.
  
          Parameters
          ----------
          dm : (m,) sequence of (N, N) array_like
              Sequence of distance matrices.
          weights : (m,) array_like, optional
              Weights used for combining arrays.
              If ``None`` then defaults to equal weights.
          **kwds :
              Keyword parameters passed to `optim_b`.
          """
          P = None
          sum_w = 0
          loop = product(dm, (1,)) if weights is None else zip(dm, weights)
          for D, w in loop:
              b = cls.optim_b(D, k, alpha, **kwds)
              _P = cls.prob_measure(D, b, alpha)
              if P is None:
                  P = _P*w
              else:
                  P += _P*w
              sum_w += w
          P = P / sum_w
          return cls(P, k, b, alpha, p_rewire=p_rewire, directed=directed)
  
      def adjacency_matrix(self, sparse=True, p_rewire=None, directed=None):
          """Generate an adjacency matrix.
  
          Parameters
          ----------
          sparse : bool
              Should sparse matrix be used.
              If ``True`` then :py:class:`scipy.sparse.csr_matrix`
              is used.
          p_rewire : float, optional
              Random edge rewiring probability.
              If ``None`` then defaults to the class attribute value.
          directed : bool or None
              Should directed graph be generated.
              If ``None`` then class attribute is used.
          """
          if p_rewire is None:
              p_rewire = self.p_rewire
          if directed is None:
              directed = self.directed
          A = make_adjacency_matrix(self.P, directed=directed)
          if p_rewire > 0:
              A = rewire_edges(A, p=p_rewire, directed=directed, copy=False)
          if sparse:
              A = csr_matrix(A)
          return A
  
      def conf_model(self, simplify=True, sparse=True, p_rewire=None, directed=None,
                     p_malformed=10**-9):
          """Generate an adjacency matrix from the configuration model.
  
          Parameters
          ----------
          simplify : bool
              Should resulting adjacency matrix be simplified by removing
              self-loops and multiple edges, which are rare in larger graphs
              but sometimes can not be avoided.
          sparse : bool
              Should sparse matrix be used.
              If ``True`` then :py:class:`scipy.sparse.csr_matrix`
              is used.
          p_rewire : float, optional
              Random edge rewiring probability.
              If ``None`` then defaults to the class attribute value.
          directed : bool or None
              Should directed graph be generated.
              If ``None`` then class attribute is used.
          configuration : bool
              Should configuration model be used.
              If ``False`` then standard SDA network is generated.
          p_malformed : float
              Probability assigned to self-loops, multiple edges and impossible
              edges in the case of hard RGG (infinite `alpha`).
              This should be very small, but non-zero value.
          self_loop_p : float
              Probability of self-loop. This must be nonzero but preferably
              very small.
          multiple_edges_p : float
              Probability of generating multiple egdes. This must be nonzero
              but preferably very small.
          """
          # pylint: disable=too-many-branches
          # pylint: disable=too-many-statements
  
          # Argument checks and setup of control variables
          if p_rewire is None:
              p_rewire = self.p_rewire
          if directed is None:
              directed = self.directed
          degseq = self.degseq
          if degseq is None:
              raise ValueError("'degseq' must be set with 'set_degseq' method")
          if directed:
              out_stubs = degseq[:, 0].copy()
              in_stubs = degseq[:, 1].copy()
          else:
              in_stubs = out_stubs = degseq.copy()
          if not directed and degseq.sum() % 2 != 0:
              raise ValueError("'degseq' must be even for undirect graphs")
          if directed and in_stubs.sum() != out_stubs.sum():
              raise ValueError("sums of in and out degrees must be equal for directed graphs")
          if directed:
              n_stubs = in_stubs.sum() + out_stubs.sum()
          else:
              n_stubs = degseq.sum()
          # Configuration algorithm
          P = self.P.copy()
          np.fill_diagonal(P, p_malformed)
          # Assign non-zero probabilities to edges impossible due to homophily
          # i.e. in the case of hard RGG
          P[P == 0] = p_malformed
          # Assign zero probabilities to edges with no stubs
          idx = np.where(out_stubs == 0)
          P[idx, :] = 0
          if not directed:
              P[:, idx] = 0
          if directed:
              idx = np.where(in_stubs == 0)
              P[:, idx] = 0
          # Main configuration algorithm loop
          W = P.sum(axis=1)
          A = np.zeros(self.P.shape, dtype=int)
          while n_stubs > 0:
              i = random_integer(W/W.sum())
              j = random_integer(P[i, :]/P[i, :].sum())
              # Decrement and check out stubs
              out_stubs[i] -= 1
              if out_stubs[i] == 1:
                  if not directed or in_stubs[i] == 0:
                      W[i] -= P[i, i]
                      P[i, i] = 0
              elif out_stubs[i] == 0:
                  W[i] = 0
                  P[i, :] = 0
                  if not directed:
                      W -= P[:, i]
                      P[:, i] = 0
              # Decrement and check in stubs
              in_stubs[j] -= 1
              if in_stubs[j] == 1:
                  if not directed or out_stubs[j] == 0:
                      W[j] -= P[j, j]
                      P[j, j] = 0
              elif in_stubs[j] == 0:
                  W -= P[:, j]
                  P[:, j] = 0
                  if not directed:
                      W[j] = 0
                      P[j, :] = 0
              # Decrement total stubs and update adjacency matrix
              n_stubs -= 2
              A[i, j] += 1
              if P[i, j] > 0:
                  W[i] += p_malformed - P[i, j]
                  P[i, j] = p_malformed
              if not directed:
                  A[j, i] += 1
                  if P[j, i] > 0:
                      W[j] += p_malformed - P[j, i]
                      P[j, i] = p_malformed
          # Final processing
          if simplify:
              np.fill_diagonal(A, 0)
              A[A > 1] = 1
          if p_rewire > 0:
              A = rewire_edges(A, p=p_rewire, directed=directed, copy=False)
          if sparse:
              A = csr_matrix(A)
          return A
  
      def to_dict(self, P=False, am='normal', ds=False, **kwds):
          """Dump to ordered dictionary.
  
          Parameters
          ----------
          P : bool
              Should edge formation probability matrix be returned.
          am : {'normal', 'configuration'} or bool or None
              Return stadard adjacency matrix if ``normal`` or ``True``.
              Return configuration model adjacency matrix if ``configuration``.
              Return ``None`` if ``None`` or ``False``.
          ds : bool
              Should degree sequence for configuration model be returned.
          **kwds :
              Keyword arguments passed to the adjacency matrix method.
  
          Raises
          ------
          ValueError
              If `A` value is incorrect.
          """
          if not am:
              A = None
          elif am is True or am == 'normal':
              A = self.adjacency_matrix(**kwds)
          elif am == 'configuration':
              A = self.conf_model(**kwds)
          else:
              raise ValueError("incorrect 'am' value")
          _ds = kwds.get('degseq', self.degseq)
          return OrderedDict([
              ('P', self.P if P else None),
              ('A', A),
              ('N', self.N),
              ('k', self.k),
              ('b', self.b),
              ('alpha', self.alpha),
              ('p_rewire', self.p_rewire),
              ('directed', self.directed),
              ('degseq', _ds if ds else None)
          ])
  