Statistics
| Branch: | Revision:

iof-tools / networkxMiCe / networkx-master / networkx / algorithms / distance_measures.py @ 5cef0f13

History | View | Annotate | Download (18.5 KB)

1
# -*- coding: utf-8 -*-
2
#    Copyright (C) 2004-2019 by
3
#    Aric Hagberg <hagberg@lanl.gov>
4
#    Dan Schult <dschult@colgate.edu>
5
#    Pieter Swart <swart@lanl.gov>
6
#    William Schwartz <wkschwartz@gmail.com>
7
#    All rights reserved.
8
#    BSD license.
9
#
10
# Authors: Aric Hagberg (hagberg@lanl.gov)
11
#          Dan Schult (dschult@colgate.edu)
12
#          Brian Kiefer (bkiefer@asu.edu)
13
"""Graph diameter, radius, eccentricity and other properties."""
14

    
15
import networkx as nx
16
from networkx.utils import not_implemented_for
17

    
18
__all__ = ['extrema_bounding', 'eccentricity', 'diameter',
19
           'radius', 'periphery', 'center', 'barycenter',
20
           'resistance_distance']
21

    
22

    
23
def extrema_bounding(G, compute="diameter"):
24
    """Compute requested extreme distance metric of undirected graph G
25

26
    Computation is based on smart lower and upper bounds, and in practice
27
    linear in the number of nodes, rather than quadratic (except for some
28
    border cases such as complete graphs or circle shaped graphs).
29

30
    Parameters
31
    ----------
32
    G : NetworkX graph
33
       An undirected graph
34

35
    compute : string denoting the requesting metric
36
       "diameter" for the maximal eccentricity value,
37
       "radius" for the minimal eccentricity value,
38
       "periphery" for the set of nodes with eccentricity equal to the diameter
39
       "center" for the set of nodes with eccentricity equal to the radius
40

41
    Returns
42
    -------
43
    value : value of the requested metric
44
       int for "diameter" and "radius" or
45
       list of nodes for "center" and "periphery"
46

47
    Raises
48
    ------
49
    NetworkXError
50
        If the graph consists of multiple components
51

52
    Notes
53
    -----
54
    This algorithm was proposed in the following papers:
55

56
    F.W. Takes and W.A. Kosters, Determining the Diameter of Small World
57
    Networks, in Proceedings of the 20th ACM International Conference on
58
    Information and Knowledge Management (CIKM 2011), pp. 1191-1196, 2011.
59
    doi: https://doi.org/10.1145/2063576.2063748
60

61
    F.W. Takes and W.A. Kosters, Computing the Eccentricity Distribution of
62
    Large Graphs, Algorithms 6(1): 100-118, 2013.
63
    doi: https://doi.org/10.3390/a6010100
64

65
    M. Borassi, P. Crescenzi, M. Habib, W.A. Kosters, A. Marino and F.W. Takes,
66
    Fast Graph Diameter and Radius BFS-Based Computation in (Weakly Connected)
67
    Real-World Graphs, Theoretical Computer Science 586: 59-80, 2015.
68
    doi: https://doi.org/10.1016/j.tcs.2015.02.033
69
    """
70

    
71
    # init variables
72
    degrees = dict(G.degree())  # start with the highest degree node
73
    minlowernode = max(degrees, key=degrees.get)
74
    N = len(degrees)  # number of nodes
75
    # alternate between smallest lower and largest upper bound
76
    high = False
77
    # status variables
78
    ecc_lower = dict.fromkeys(G, 0)
79
    ecc_upper = dict.fromkeys(G, N)
80
    candidates = set(G)
81

    
82
    # (re)set bound extremes
83
    minlower = N
84
    maxlower = 0
85
    minupper = N
86
    maxupper = 0
87

    
88
    # repeat the following until there are no more candidates
89
    while candidates:
90
        if high:
91
            current = maxuppernode  # select node with largest upper bound
92
        else:
93
            current = minlowernode  # select node with smallest lower bound
94
        high = not high
95

    
96
        # get distances from/to current node and derive eccentricity
97
        dist = dict(nx.single_source_shortest_path_length(G, current))
98
        if len(dist) != N:
99
            msg = ('Cannot compute metric because graph is not connected.')
100
            raise nx.NetworkXError(msg)
101
        current_ecc = max(dist.values())
102

    
103
        # print status update
104
#        print ("ecc of " + str(current) + " (" + str(ecc_lower[current]) + "/"
105
#        + str(ecc_upper[current]) + ", deg: " + str(dist[current]) + ") is "
106
#        + str(current_ecc))
107
#        print(ecc_upper)
108

    
109
        # (re)set bound extremes
110
        maxuppernode = None
111
        minlowernode = None
112

    
113
        # update node bounds
114
        for i in candidates:
115
            # update eccentricity bounds
116
            d = dist[i]
117
            ecc_lower[i] = low = max(ecc_lower[i], max(d, (current_ecc - d)))
118
            ecc_upper[i] = upp = min(ecc_upper[i], current_ecc + d)
119

    
120
            # update min/max values of lower and upper bounds
121
            minlower = min(ecc_lower[i], minlower)
122
            maxlower = max(ecc_lower[i], maxlower)
123
            minupper = min(ecc_upper[i], minupper)
124
            maxupper = max(ecc_upper[i], maxupper)
125

    
126
        # update candidate set
127
        if compute == 'diameter':
128
            ruled_out = {i for i in candidates if ecc_upper[i] <= maxlower and
129
                         2 * ecc_lower[i] >= maxupper}
130

    
131
        elif compute == 'radius':
132
            ruled_out = {i for i in candidates if ecc_lower[i] >= minupper and
133
                         ecc_upper[i] + 1 <= 2 * minlower}
134

    
135
        elif compute == 'periphery':
136
            ruled_out = {i for i in candidates if ecc_upper[i] < maxlower and
137
                         (maxlower == maxupper or ecc_lower[i] > maxupper)}
138

    
139
        elif compute == 'center':
140
            ruled_out = {i for i in candidates if ecc_lower[i] > minupper and
141
                         (minlower == minupper or ecc_upper[i] + 1 < 2 * minlower)}
142

    
143
        elif compute == 'eccentricities':
144
            ruled_out = {}
145

    
146
        ruled_out.update(i for i in candidates if ecc_lower[i] == ecc_upper[i])
147
        candidates -= ruled_out
148

    
149
#        for i in ruled_out:
150
#            print("removing %g: ecc_u: %g maxl: %g ecc_l: %g maxu: %g"%
151
#                    (i,ecc_upper[i],maxlower,ecc_lower[i],maxupper))
152
#        print("node %g: ecc_u: %g maxl: %g ecc_l: %g maxu: %g"%
153
#                    (4,ecc_upper[4],maxlower,ecc_lower[4],maxupper))
154
#        print("NODE 4: %g"%(ecc_upper[4] <= maxlower))
155
#        print("NODE 4: %g"%(2 * ecc_lower[4] >= maxupper))
156
#        print("NODE 4: %g"%(ecc_upper[4] <= maxlower
157
#                            and 2 * ecc_lower[4] >= maxupper))
158

    
159
        # updating maxuppernode and minlowernode for selection in next round
160
        for i in candidates:
161
            if minlowernode is None \
162
                    or (ecc_lower[i] == ecc_lower[minlowernode]
163
                        and degrees[i] > degrees[minlowernode]) \
164
                    or (ecc_lower[i] < ecc_lower[minlowernode]):
165
                minlowernode = i
166

    
167
            if maxuppernode is None \
168
                    or (ecc_upper[i] == ecc_upper[maxuppernode]
169
                        and degrees[i] > degrees[maxuppernode]) \
170
                    or (ecc_upper[i] > ecc_upper[maxuppernode]):
171
                maxuppernode = i
172

    
173
        # print status update
174
#        print (" min=" + str(minlower) + "/" + str(minupper) +
175
#        " max=" + str(maxlower) + "/" + str(maxupper) +
176
#        " candidates: " + str(len(candidates)))
177
#        print("cand:",candidates)
178
#        print("ecc_l",ecc_lower)
179
#        print("ecc_u",ecc_upper)
180
#        wait = input("press Enter to continue")
181

    
182
    # return the correct value of the requested metric
183
    if compute == 'diameter':
184
        return maxlower
185
    elif compute == 'radius':
186
        return minupper
187
    elif compute == 'periphery':
188
        p = [v for v in G if ecc_lower[v] == maxlower]
189
        return p
190
    elif compute == 'center':
191
        c = [v for v in G if ecc_upper[v] == minupper]
192
        return c
193
    elif compute == 'eccentricities':
194
        return ecc_lower
195
    return None
196

    
197

    
198
def eccentricity(G, v=None, sp=None):
199
    """Returns the eccentricity of nodes in G.
200

201
    The eccentricity of a node v is the maximum distance from v to
202
    all other nodes in G.
203

204
    Parameters
205
    ----------
206
    G : NetworkX graph
207
       A graph
208

209
    v : node, optional
210
       Return value of specified node
211

212
    sp : dict of dicts, optional
213
       All pairs shortest path lengths as a dictionary of dictionaries
214

215
    Returns
216
    -------
217
    ecc : dictionary
218
       A dictionary of eccentricity values keyed by node.
219
    """
220
#    if v is None:                # none, use entire graph
221
#        nodes=G.nodes()
222
#    elif v in G:               # is v a single node
223
#        nodes=[v]
224
#    else:                      # assume v is a container of nodes
225
#        nodes=v
226
    order = G.order()
227

    
228
    e = {}
229
    for n in G.nbunch_iter(v):
230
        if sp is None:
231
            length = nx.single_source_shortest_path_length(G, n)
232
            L = len(length)
233
        else:
234
            try:
235
                length = sp[n]
236
                L = len(length)
237
            except TypeError:
238
                raise nx.NetworkXError('Format of "sp" is invalid.')
239
        if L != order:
240
            if G.is_directed():
241
                msg = ('Found infinite path length because the digraph is not'
242
                       ' strongly connected')
243
            else:
244
                msg = ('Found infinite path length because the graph is not'
245
                       ' connected')
246
            raise nx.NetworkXError(msg)
247

    
248
        e[n] = max(length.values())
249

    
250
    if v in G:
251
        return e[v]  # return single value
252
    else:
253
        return e
254

    
255

    
256
def diameter(G, e=None, usebounds=False):
257
    """Returns the diameter of the graph G.
258

259
    The diameter is the maximum eccentricity.
260

261
    Parameters
262
    ----------
263
    G : NetworkX graph
264
       A graph
265

266
    e : eccentricity dictionary, optional
267
      A precomputed dictionary of eccentricities.
268

269
    Returns
270
    -------
271
    d : integer
272
       Diameter of graph
273

274
    See Also
275
    --------
276
    eccentricity
277
    """
278
    if usebounds is True and e is None and not G.is_directed():
279
        return extrema_bounding(G, compute="diameter")
280
    if e is None:
281
        e = eccentricity(G)
282
    return max(e.values())
283

    
284

    
285
def periphery(G, e=None, usebounds=False):
286
    """Returns the periphery of the graph G.
287

288
    The periphery is the set of nodes with eccentricity equal to the diameter.
289

290
    Parameters
291
    ----------
292
    G : NetworkX graph
293
       A graph
294

295
    e : eccentricity dictionary, optional
296
      A precomputed dictionary of eccentricities.
297

298
    Returns
299
    -------
300
    p : list
301
       List of nodes in periphery
302

303
    See Also
304
    --------
305
    barycenter
306
    center
307
    """
308
    if usebounds is True and e is None and not G.is_directed():
309
        return extrema_bounding(G, compute="periphery")
310
    if e is None:
311
        e = eccentricity(G)
312
    diameter = max(e.values())
313
    p = [v for v in e if e[v] == diameter]
314
    return p
315

    
316

    
317
def radius(G, e=None, usebounds=False):
318
    """Returns the radius of the graph G.
319

320
    The radius is the minimum eccentricity.
321

322
    Parameters
323
    ----------
324
    G : NetworkX graph
325
       A graph
326

327
    e : eccentricity dictionary, optional
328
      A precomputed dictionary of eccentricities.
329

330
    Returns
331
    -------
332
    r : integer
333
       Radius of graph
334
    """
335
    if usebounds is True and e is None and not G.is_directed():
336
        return extrema_bounding(G, compute="radius")
337
    if e is None:
338
        e = eccentricity(G)
339
    return min(e.values())
340

    
341

    
342
def center(G, e=None, usebounds=False):
343
    """Returns the center of the graph G.
344

345
    The center is the set of nodes with eccentricity equal to radius.
346

347
    Parameters
348
    ----------
349
    G : NetworkX graph
350
       A graph
351

352
    e : eccentricity dictionary, optional
353
      A precomputed dictionary of eccentricities.
354

355
    Returns
356
    -------
357
    c : list
358
       List of nodes in center
359

360
    See Also
361
    --------
362
    barycenter
363
    periphery
364
    """
365
    if usebounds is True and e is None and not G.is_directed():
366
        return extrema_bounding(G, compute="center")
367
    if e is None:
368
        e = eccentricity(G)
369
    radius = min(e.values())
370
    p = [v for v in e if e[v] == radius]
371
    return p
372

    
373

    
374
def barycenter(G, weight=None, attr=None, sp=None):
375
    r"""Calculate barycenter of a connected graph, optionally with edge weights.
376

377
    The :dfn:`barycenter` a
378
    :func:`connected <networkx.algorithms.components.is_connected>` graph
379
    :math:`G` is the subgraph induced by the set of its nodes :math:`v`
380
    minimizing the objective function
381

382
    .. math::
383

384
        \sum_{u \in V(G)} d_G(u, v),
385

386
    where :math:`d_G` is the (possibly weighted) :func:`path length
387
    <networkx.algorithms.shortest_paths.generic.shortest_path_length>`.
388
    The barycenter is also called the :dfn:`median`. See [West01]_, p. 78.
389

390
    Parameters
391
    ----------
392
    G : :class:`networkx.Graph`
393
        The connected graph :math:`G`.
394
    weight : :class:`str`, optional
395
        Passed through to
396
        :func:`~networkx.algorithms.shortest_paths.generic.shortest_path_length`.
397
    attr : :class:`str`, optional
398
        If given, write the value of the objective function to each node's
399
        `attr` attribute. Otherwise do not store the value.
400
    sp : dict of dicts, optional
401
       All pairs shortest path lengths as a dictionary of dictionaries
402

403
    Returns
404
    -------
405
    :class:`list`
406
        Nodes of `G` that induce the barycenter of `G`.
407

408
    Raises
409
    ------
410
    :exc:`networkx.NetworkXNoPath`
411
        If `G` is disconnected. `G` may appear disconnected to
412
        :func:`barycenter` if `sp` is given but is missing shortest path
413
        lengths for any pairs.
414
    :exc:`ValueError`
415
        If `sp` and `weight` are both given.
416

417
    See Also
418
    --------
419
    center
420
    periphery
421
    """
422
    if sp is None:
423
        sp = nx.shortest_path_length(G, weight=weight)
424
    else:
425
        sp = sp.items()
426
        if weight is not None:
427
            raise ValueError('Cannot use both sp, weight arguments together')
428
    smallest, barycenter_vertices, n = float('inf'), [], len(G)
429
    for v, dists in sp:
430
        if len(dists) < n:
431
            raise nx.NetworkXNoPath(
432
                ("Input graph %r is disconnected, so every induced subgraph "
433
                 "has infinite barycentricity.") % G)
434
        barycentricity = sum(dists.values())
435
        if attr is not None:
436
            G.nodes[v][attr] = barycentricity
437
        if barycentricity < smallest:
438
            smallest = barycentricity
439
            barycenter_vertices = [v]
440
        elif barycentricity == smallest:
441
            barycenter_vertices.append(v)
442
    return barycenter_vertices
443

    
444

    
445
def _laplacian_submatrix(node, mat, node_list):
446
    """Removes row/col from a sparse matrix and returns the submatrix
447
    """
448
    j = node_list.index(node)
449
    n = list(range(len(node_list)))
450
    n.pop(j)
451

    
452
    if mat.shape[0] != mat.shape[1]:
453
        raise nx.NetworkXError('Matrix must be square')
454
    elif len(node_list) != mat.shape[0]:
455
        msg = "Node list length does not match matrix dimentions"
456
        raise nx.NetworkXError(msg)
457

    
458
    mat = mat.tocsr()
459
    mat = mat[n, :]
460

    
461
    mat = mat.tocsc()
462
    mat = mat[:, n]
463

    
464
    node_list.pop(j)
465

    
466
    return mat, node_list
467

    
468

    
469
def _count_lu_permutations(perm_array):
470
    """Counts the number of permutations in SuperLU perm_c or perm_r
471
    """
472
    perm_cnt = 0
473
    arr = perm_array.tolist()
474
    for i in range(len(arr)):
475
        if i != arr[i]:
476
            perm_cnt += 1
477
            n = arr.index(i)
478
            arr[n] = arr[i]
479
            arr[i] = i
480

    
481
    return perm_cnt
482

    
483

    
484
@not_implemented_for('directed')
485
def resistance_distance(G, nodeA, nodeB, weight=None, invert_weight=True):
486
    """Returns the resistance distance between node A and node B on graph G.
487

488
    The resistance distance between two nodes of a graph is akin to treating
489
    the graph as a grid of resistorses with a resistance equal to the provided
490
    weight.
491

492
    If weight is not provided, then a weight of 1 is used for all edges.
493

494
    Parameters
495
    ----------
496
    G : NetworkX graph
497
       A graph
498

499
    nodeA : node
500
      A node within graph G.
501

502
    nodeB : node
503
      A node within graph G, exclusive of Node A.
504

505
    weight : string or None, optional (default=None)
506
       The edge data key used to compute the resistance distance.
507
       If None, then each edge has weight 1.
508

509
    invert_weight : boolean (default=True)
510
        Proper calculation of resistance distance requires building the
511
        Laplacian matrix with the reciprocal of the weight. Not required
512
        if the weight is already inverted. Weight cannot be zero.
513

514
    Returns
515
    -------
516
    rd : float
517
       Value of effective resistance distance
518

519
    Notes
520
    -----
521
    Overview discussion:
522
    * https://en.wikipedia.org/wiki/Resistance_distance
523
    * http://mathworld.wolfram.com/ResistanceDistance.html
524

525
    Additional details:
526
    Vaya Sapobi Samui Vos, “Methods for determining the effective resistance,” M.S.,
527
    Mathematisch Instituut, Universiteit Leiden, Leiden, Netherlands, 2016
528
    Available: `Link to thesis <https://www.universiteitleiden.nl/binaries/content/assets/science/mi/scripties/master/vos_vaya_master.pdf>`_
529
    """
530
    import numpy as np
531
    import scipy.sparse
532

    
533
    if not nx.is_connected(G):
534
        msg = ('Graph G must be strongly connected.')
535
        raise nx.NetworkXError(msg)
536
    elif nodeA not in G:
537
        msg = ('Node A is not in graph G.')
538
        raise nx.NetworkXError(msg)
539
    elif nodeB not in G:
540
        msg = ('Node B is not in graph G.')
541
        raise nx.NetworkXError(msg)
542
    elif nodeA == nodeB:
543
        msg = ('Node A and Node B cannot be the same.')
544
        raise nx.NetworkXError(msg)
545

    
546
    G = G.copy()
547
    node_list = list(G)
548

    
549
    if invert_weight and weight is not None:
550
        if G.is_multigraph():
551
            for (u, v, k, d) in G.edges(keys=True, data=True):
552
                d[weight] = 1/d[weight]
553
        else:
554
           for (u, v, d) in G.edges(data=True):
555
               d[weight] = 1/d[weight]
556
    # Replace with collapsing topology or approximated zero?
557

    
558
    # Using determinants to compute the effective resistance is more memory
559
    # efficent than directly calculating the psuedo-inverse
560
    L = nx.laplacian_matrix(G, node_list, weight=weight)
561

    
562
    Lsub_a, node_list_a = _laplacian_submatrix(nodeA, L.copy(),
563
                                               node_list[:])
564
    Lsub_ab, node_list_ab = _laplacian_submatrix(nodeB, Lsub_a.copy(),
565
                                                 node_list_a[:])
566

    
567
    # Factorize Laplacian submatrixes and extract diagonals
568
    # Order the diagonals to minimize the likelihood over overflows
569
    # during computing the determinant
570
    lu_a = scipy.sparse.linalg.splu(Lsub_a, options=dict(SymmetricMode=True))
571
    LdiagA = lu_a.U.diagonal()
572
    LdiagA_s = np.product(np.sign(LdiagA)) * np.product(lu_a.L.diagonal())
573
    LdiagA_s *= (-1)**_count_lu_permutations(lu_a.perm_r)
574
    LdiagA_s *= (-1)**_count_lu_permutations(lu_a.perm_c)
575
    LdiagA = np.absolute(LdiagA)
576
    LdiagA = np.sort(LdiagA)
577

    
578
    lu_ab = scipy.sparse.linalg.splu(Lsub_ab, options=dict(SymmetricMode=True))
579
    LdiagAB = lu_ab.U.diagonal()
580
    LdiagAB_s = np.product(np.sign(LdiagAB)) * np.product(lu_ab.L.diagonal())
581
    LdiagAB_s *= (-1)**_count_lu_permutations(lu_ab.perm_r)
582
    LdiagAB_s *= (-1)**_count_lu_permutations(lu_ab.perm_c)
583
    LdiagAB = np.absolute(LdiagAB)
584
    LdiagAB = np.sort(LdiagAB)
585

    
586
    # Calculate the ratio of determinant, rd = det(Lsub_ab)/det(Lsub_a)
587
    Ldet = np.product(np.divide(np.append(LdiagAB, [1]), LdiagA))
588
    rd = Ldet * LdiagAB_s / LdiagA_s
589

    
590
    return rd