ioftools / networkxMiCe / networkxmaster / networkx / algorithms / distance_measures.py @ 5cef0f13
History  View  Annotate  Download (18.5 KB)
1 
# * coding: utf8 *


2 
# Copyright (C) 20042019 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. 11911196, 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): 100118, 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 BFSBased Computation in (Weakly Connected)

67 
RealWorld Graphs, Theoretical Computer Science 586: 5980, 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 psuedoinverse

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
