Mercurial > repos > shellac > sam_consensus_v3
comparison env/lib/python3.9/site-packages/networkx/linalg/algebraicconnectivity.py @ 0:4f3585e2f14b draft default tip
"planemo upload commit 60cee0fc7c0cda8592644e1aad72851dec82c959"
author | shellac |
---|---|
date | Mon, 22 Mar 2021 18:12:50 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:4f3585e2f14b |
---|---|
1 """ | |
2 Algebraic connectivity and Fiedler vectors of undirected graphs. | |
3 """ | |
4 from functools import partial | |
5 import networkx as nx | |
6 from networkx.utils import not_implemented_for | |
7 from networkx.utils import reverse_cuthill_mckee_ordering | |
8 from networkx.utils import random_state | |
9 | |
10 try: | |
11 from numpy import array, asarray, dot, ndarray, ones, sqrt, zeros, atleast_2d | |
12 from numpy.linalg import norm, qr | |
13 from scipy.linalg import eigh, inv | |
14 from scipy.sparse import csc_matrix, spdiags | |
15 from scipy.sparse.linalg import eigsh, lobpcg | |
16 | |
17 __all__ = ["algebraic_connectivity", "fiedler_vector", "spectral_ordering"] | |
18 except ImportError: | |
19 __all__ = [] | |
20 | |
21 try: | |
22 from scipy.linalg.blas import dasum, daxpy, ddot | |
23 except ImportError: | |
24 if __all__: | |
25 # Make sure the imports succeeded. | |
26 # Use minimal replacements if BLAS is unavailable from SciPy. | |
27 dasum = partial(norm, ord=1) | |
28 ddot = dot | |
29 | |
30 def daxpy(x, y, a): | |
31 y += a * x | |
32 return y | |
33 | |
34 | |
35 class _PCGSolver: | |
36 """Preconditioned conjugate gradient method. | |
37 | |
38 To solve Ax = b: | |
39 M = A.diagonal() # or some other preconditioner | |
40 solver = _PCGSolver(lambda x: A * x, lambda x: M * x) | |
41 x = solver.solve(b) | |
42 | |
43 The inputs A and M are functions which compute | |
44 matrix multiplication on the argument. | |
45 A - multiply by the matrix A in Ax=b | |
46 M - multiply by M, the preconditioner surragate for A | |
47 | |
48 Warning: There is no limit on number of iterations. | |
49 """ | |
50 | |
51 def __init__(self, A, M): | |
52 self._A = A | |
53 self._M = M or (lambda x: x.copy()) | |
54 | |
55 def solve(self, B, tol): | |
56 B = asarray(B) | |
57 X = ndarray(B.shape, order="F") | |
58 for j in range(B.shape[1]): | |
59 X[:, j] = self._solve(B[:, j], tol) | |
60 return X | |
61 | |
62 def _solve(self, b, tol): | |
63 A = self._A | |
64 M = self._M | |
65 tol *= dasum(b) | |
66 # Initialize. | |
67 x = zeros(b.shape) | |
68 r = b.copy() | |
69 z = M(r) | |
70 rz = ddot(r, z) | |
71 p = z.copy() | |
72 # Iterate. | |
73 while True: | |
74 Ap = A(p) | |
75 alpha = rz / ddot(p, Ap) | |
76 x = daxpy(p, x, a=alpha) | |
77 r = daxpy(Ap, r, a=-alpha) | |
78 if dasum(r) < tol: | |
79 return x | |
80 z = M(r) | |
81 beta = ddot(r, z) | |
82 beta, rz = beta / rz, beta | |
83 p = daxpy(p, z, a=beta) | |
84 | |
85 | |
86 class _CholeskySolver: | |
87 """Cholesky factorization. | |
88 | |
89 To solve Ax = b: | |
90 solver = _CholeskySolver(A) | |
91 x = solver.solve(b) | |
92 | |
93 optional argument `tol` on solve method is ignored but included | |
94 to match _PCGsolver API. | |
95 """ | |
96 | |
97 def __init__(self, A): | |
98 if not self._cholesky: | |
99 raise nx.NetworkXError("Cholesky solver unavailable.") | |
100 self._chol = self._cholesky(A) | |
101 | |
102 def solve(self, B, tol=None): | |
103 return self._chol(B) | |
104 | |
105 try: | |
106 from scikits.sparse.cholmod import cholesky | |
107 | |
108 _cholesky = cholesky | |
109 except ImportError: | |
110 _cholesky = None | |
111 | |
112 | |
113 class _LUSolver: | |
114 """LU factorization. | |
115 | |
116 To solve Ax = b: | |
117 solver = _LUSolver(A) | |
118 x = solver.solve(b) | |
119 | |
120 optional argument `tol` on solve method is ignored but included | |
121 to match _PCGsolver API. | |
122 """ | |
123 | |
124 def __init__(self, A): | |
125 if not self._splu: | |
126 raise nx.NetworkXError("LU solver unavailable.") | |
127 self._LU = self._splu(A) | |
128 | |
129 def solve(self, B, tol=None): | |
130 B = asarray(B) | |
131 X = ndarray(B.shape, order="F") | |
132 for j in range(B.shape[1]): | |
133 X[:, j] = self._LU.solve(B[:, j]) | |
134 return X | |
135 | |
136 try: | |
137 from scipy.sparse.linalg import splu | |
138 | |
139 _splu = partial( | |
140 splu, | |
141 permc_spec="MMD_AT_PLUS_A", | |
142 diag_pivot_thresh=0.0, | |
143 options={"Equil": True, "SymmetricMode": True}, | |
144 ) | |
145 except ImportError: | |
146 _splu = None | |
147 | |
148 | |
149 def _preprocess_graph(G, weight): | |
150 """Compute edge weights and eliminate zero-weight edges. | |
151 """ | |
152 if G.is_directed(): | |
153 H = nx.MultiGraph() | |
154 H.add_nodes_from(G) | |
155 H.add_weighted_edges_from( | |
156 ((u, v, e.get(weight, 1.0)) for u, v, e in G.edges(data=True) if u != v), | |
157 weight=weight, | |
158 ) | |
159 G = H | |
160 if not G.is_multigraph(): | |
161 edges = ( | |
162 (u, v, abs(e.get(weight, 1.0))) for u, v, e in G.edges(data=True) if u != v | |
163 ) | |
164 else: | |
165 edges = ( | |
166 (u, v, sum(abs(e.get(weight, 1.0)) for e in G[u][v].values())) | |
167 for u, v in G.edges() | |
168 if u != v | |
169 ) | |
170 H = nx.Graph() | |
171 H.add_nodes_from(G) | |
172 H.add_weighted_edges_from((u, v, e) for u, v, e in edges if e != 0) | |
173 return H | |
174 | |
175 | |
176 def _rcm_estimate(G, nodelist): | |
177 """Estimate the Fiedler vector using the reverse Cuthill-McKee ordering. | |
178 """ | |
179 G = G.subgraph(nodelist) | |
180 order = reverse_cuthill_mckee_ordering(G) | |
181 n = len(nodelist) | |
182 index = dict(zip(nodelist, range(n))) | |
183 x = ndarray(n, dtype=float) | |
184 for i, u in enumerate(order): | |
185 x[index[u]] = i | |
186 x -= (n - 1) / 2.0 | |
187 return x | |
188 | |
189 | |
190 def _tracemin_fiedler(L, X, normalized, tol, method): | |
191 """Compute the Fiedler vector of L using the TraceMIN-Fiedler algorithm. | |
192 | |
193 The Fiedler vector of a connected undirected graph is the eigenvector | |
194 corresponding to the second smallest eigenvalue of the Laplacian matrix of | |
195 of the graph. This function starts with the Laplacian L, not the Graph. | |
196 | |
197 Parameters | |
198 ---------- | |
199 L : Laplacian of a possibly weighted or normalized, but undirected graph | |
200 | |
201 X : Initial guess for a solution. Usually a matrix of random numbers. | |
202 This function allows more than one column in X to identify more than | |
203 one eigenvector if desired. | |
204 | |
205 normalized : bool | |
206 Whether the normalized Laplacian matrix is used. | |
207 | |
208 tol : float | |
209 Tolerance of relative residual in eigenvalue computation. | |
210 Warning: There is no limit on number of iterations. | |
211 | |
212 method : string | |
213 Should be 'tracemin_pcg', 'tracemin_chol' or 'tracemin_lu'. | |
214 Otherwise exception is raised. | |
215 | |
216 Returns | |
217 ------- | |
218 sigma, X : Two NumPy arrays of floats. | |
219 The lowest eigenvalues and corresponding eigenvectors of L. | |
220 The size of input X determines the size of these outputs. | |
221 As this is for Fiedler vectors, the zero eigenvalue (and | |
222 constant eigenvector) are avoided. | |
223 """ | |
224 n = X.shape[0] | |
225 | |
226 if normalized: | |
227 # Form the normalized Laplacian matrix and determine the eigenvector of | |
228 # its nullspace. | |
229 e = sqrt(L.diagonal()) | |
230 D = spdiags(1.0 / e, [0], n, n, format="csr") | |
231 L = D * L * D | |
232 e *= 1.0 / norm(e, 2) | |
233 | |
234 if normalized: | |
235 | |
236 def project(X): | |
237 """Make X orthogonal to the nullspace of L. | |
238 """ | |
239 X = asarray(X) | |
240 for j in range(X.shape[1]): | |
241 X[:, j] -= dot(X[:, j], e) * e | |
242 | |
243 else: | |
244 | |
245 def project(X): | |
246 """Make X orthogonal to the nullspace of L. | |
247 """ | |
248 X = asarray(X) | |
249 for j in range(X.shape[1]): | |
250 X[:, j] -= X[:, j].sum() / n | |
251 | |
252 if method == "tracemin_pcg": | |
253 D = L.diagonal().astype(float) | |
254 solver = _PCGSolver(lambda x: L * x, lambda x: D * x) | |
255 elif method == "tracemin_chol" or method == "tracemin_lu": | |
256 # Convert A to CSC to suppress SparseEfficiencyWarning. | |
257 A = csc_matrix(L, dtype=float, copy=True) | |
258 # Force A to be nonsingular. Since A is the Laplacian matrix of a | |
259 # connected graph, its rank deficiency is one, and thus one diagonal | |
260 # element needs to modified. Changing to infinity forces a zero in the | |
261 # corresponding element in the solution. | |
262 i = (A.indptr[1:] - A.indptr[:-1]).argmax() | |
263 A[i, i] = float("inf") | |
264 if method == "tracemin_chol": | |
265 solver = _CholeskySolver(A) | |
266 else: | |
267 solver = _LUSolver(A) | |
268 else: | |
269 raise nx.NetworkXError("Unknown linear system solver: " + method) | |
270 | |
271 # Initialize. | |
272 Lnorm = abs(L).sum(axis=1).flatten().max() | |
273 project(X) | |
274 W = ndarray(X.shape, order="F") | |
275 | |
276 while True: | |
277 # Orthonormalize X. | |
278 X = qr(X)[0] | |
279 # Compute iteration matrix H. | |
280 W[:, :] = L @ X | |
281 H = X.T @ W | |
282 sigma, Y = eigh(H, overwrite_a=True) | |
283 # Compute the Ritz vectors. | |
284 X = X @ Y | |
285 # Test for convergence exploiting the fact that L * X == W * Y. | |
286 res = dasum(W @ Y[:, 0] - sigma[0] * X[:, 0]) / Lnorm | |
287 if res < tol: | |
288 break | |
289 # Compute X = L \ X / (X' * (L \ X)). | |
290 # L \ X can have an arbitrary projection on the nullspace of L, | |
291 # which will be eliminated. | |
292 W[:, :] = solver.solve(X, tol) | |
293 X = (inv(W.T @ X) @ W.T).T # Preserves Fortran storage order. | |
294 project(X) | |
295 | |
296 return sigma, asarray(X) | |
297 | |
298 | |
299 def _get_fiedler_func(method): | |
300 """Returns a function that solves the Fiedler eigenvalue problem. | |
301 """ | |
302 if method == "tracemin": # old style keyword <v2.1 | |
303 method = "tracemin_pcg" | |
304 if method in ("tracemin_pcg", "tracemin_chol", "tracemin_lu"): | |
305 | |
306 def find_fiedler(L, x, normalized, tol, seed): | |
307 q = 1 if method == "tracemin_pcg" else min(4, L.shape[0] - 1) | |
308 X = asarray(seed.normal(size=(q, L.shape[0]))).T | |
309 sigma, X = _tracemin_fiedler(L, X, normalized, tol, method) | |
310 return sigma[0], X[:, 0] | |
311 | |
312 elif method == "lanczos" or method == "lobpcg": | |
313 | |
314 def find_fiedler(L, x, normalized, tol, seed): | |
315 L = csc_matrix(L, dtype=float) | |
316 n = L.shape[0] | |
317 if normalized: | |
318 D = spdiags(1.0 / sqrt(L.diagonal()), [0], n, n, format="csc") | |
319 L = D * L * D | |
320 if method == "lanczos" or n < 10: | |
321 # Avoid LOBPCG when n < 10 due to | |
322 # https://github.com/scipy/scipy/issues/3592 | |
323 # https://github.com/scipy/scipy/pull/3594 | |
324 sigma, X = eigsh(L, 2, which="SM", tol=tol, return_eigenvectors=True) | |
325 return sigma[1], X[:, 1] | |
326 else: | |
327 X = asarray(atleast_2d(x).T) | |
328 M = spdiags(1.0 / L.diagonal(), [0], n, n) | |
329 Y = ones(n) | |
330 if normalized: | |
331 Y /= D.diagonal() | |
332 sigma, X = lobpcg( | |
333 L, X, M=M, Y=atleast_2d(Y).T, tol=tol, maxiter=n, largest=False | |
334 ) | |
335 return sigma[0], X[:, 0] | |
336 | |
337 else: | |
338 raise nx.NetworkXError(f"unknown method '{method}'.") | |
339 | |
340 return find_fiedler | |
341 | |
342 | |
343 @random_state(5) | |
344 @not_implemented_for("directed") | |
345 def algebraic_connectivity( | |
346 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None | |
347 ): | |
348 """Returns the algebraic connectivity of an undirected graph. | |
349 | |
350 The algebraic connectivity of a connected undirected graph is the second | |
351 smallest eigenvalue of its Laplacian matrix. | |
352 | |
353 Parameters | |
354 ---------- | |
355 G : NetworkX graph | |
356 An undirected graph. | |
357 | |
358 weight : object, optional (default: None) | |
359 The data key used to determine the weight of each edge. If None, then | |
360 each edge has unit weight. | |
361 | |
362 normalized : bool, optional (default: False) | |
363 Whether the normalized Laplacian matrix is used. | |
364 | |
365 tol : float, optional (default: 1e-8) | |
366 Tolerance of relative residual in eigenvalue computation. | |
367 | |
368 method : string, optional (default: 'tracemin_pcg') | |
369 Method of eigenvalue computation. It must be one of the tracemin | |
370 options shown below (TraceMIN), 'lanczos' (Lanczos iteration) | |
371 or 'lobpcg' (LOBPCG). | |
372 | |
373 The TraceMIN algorithm uses a linear system solver. The following | |
374 values allow specifying the solver to be used. | |
375 | |
376 =============== ======================================== | |
377 Value Solver | |
378 =============== ======================================== | |
379 'tracemin_pcg' Preconditioned conjugate gradient method | |
380 'tracemin_chol' Cholesky factorization | |
381 'tracemin_lu' LU factorization | |
382 =============== ======================================== | |
383 | |
384 seed : integer, random_state, or None (default) | |
385 Indicator of random number generation state. | |
386 See :ref:`Randomness<randomness>`. | |
387 | |
388 Returns | |
389 ------- | |
390 algebraic_connectivity : float | |
391 Algebraic connectivity. | |
392 | |
393 Raises | |
394 ------ | |
395 NetworkXNotImplemented | |
396 If G is directed. | |
397 | |
398 NetworkXError | |
399 If G has less than two nodes. | |
400 | |
401 Notes | |
402 ----- | |
403 Edge weights are interpreted by their absolute values. For MultiGraph's, | |
404 weights of parallel edges are summed. Zero-weighted edges are ignored. | |
405 | |
406 To use Cholesky factorization in the TraceMIN algorithm, the | |
407 :samp:`scikits.sparse` package must be installed. | |
408 | |
409 See Also | |
410 -------- | |
411 laplacian_matrix | |
412 """ | |
413 if len(G) < 2: | |
414 raise nx.NetworkXError("graph has less than two nodes.") | |
415 G = _preprocess_graph(G, weight) | |
416 if not nx.is_connected(G): | |
417 return 0.0 | |
418 | |
419 L = nx.laplacian_matrix(G) | |
420 if L.shape[0] == 2: | |
421 return 2.0 * L[0, 0] if not normalized else 2.0 | |
422 | |
423 find_fiedler = _get_fiedler_func(method) | |
424 x = None if method != "lobpcg" else _rcm_estimate(G, G) | |
425 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed) | |
426 return sigma | |
427 | |
428 | |
429 @random_state(5) | |
430 @not_implemented_for("directed") | |
431 def fiedler_vector( | |
432 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None | |
433 ): | |
434 """Returns the Fiedler vector of a connected undirected graph. | |
435 | |
436 The Fiedler vector of a connected undirected graph is the eigenvector | |
437 corresponding to the second smallest eigenvalue of the Laplacian matrix of | |
438 of the graph. | |
439 | |
440 Parameters | |
441 ---------- | |
442 G : NetworkX graph | |
443 An undirected graph. | |
444 | |
445 weight : object, optional (default: None) | |
446 The data key used to determine the weight of each edge. If None, then | |
447 each edge has unit weight. | |
448 | |
449 normalized : bool, optional (default: False) | |
450 Whether the normalized Laplacian matrix is used. | |
451 | |
452 tol : float, optional (default: 1e-8) | |
453 Tolerance of relative residual in eigenvalue computation. | |
454 | |
455 method : string, optional (default: 'tracemin_pcg') | |
456 Method of eigenvalue computation. It must be one of the tracemin | |
457 options shown below (TraceMIN), 'lanczos' (Lanczos iteration) | |
458 or 'lobpcg' (LOBPCG). | |
459 | |
460 The TraceMIN algorithm uses a linear system solver. The following | |
461 values allow specifying the solver to be used. | |
462 | |
463 =============== ======================================== | |
464 Value Solver | |
465 =============== ======================================== | |
466 'tracemin_pcg' Preconditioned conjugate gradient method | |
467 'tracemin_chol' Cholesky factorization | |
468 'tracemin_lu' LU factorization | |
469 =============== ======================================== | |
470 | |
471 seed : integer, random_state, or None (default) | |
472 Indicator of random number generation state. | |
473 See :ref:`Randomness<randomness>`. | |
474 | |
475 Returns | |
476 ------- | |
477 fiedler_vector : NumPy array of floats. | |
478 Fiedler vector. | |
479 | |
480 Raises | |
481 ------ | |
482 NetworkXNotImplemented | |
483 If G is directed. | |
484 | |
485 NetworkXError | |
486 If G has less than two nodes or is not connected. | |
487 | |
488 Notes | |
489 ----- | |
490 Edge weights are interpreted by their absolute values. For MultiGraph's, | |
491 weights of parallel edges are summed. Zero-weighted edges are ignored. | |
492 | |
493 To use Cholesky factorization in the TraceMIN algorithm, the | |
494 :samp:`scikits.sparse` package must be installed. | |
495 | |
496 See Also | |
497 -------- | |
498 laplacian_matrix | |
499 """ | |
500 if len(G) < 2: | |
501 raise nx.NetworkXError("graph has less than two nodes.") | |
502 G = _preprocess_graph(G, weight) | |
503 if not nx.is_connected(G): | |
504 raise nx.NetworkXError("graph is not connected.") | |
505 | |
506 if len(G) == 2: | |
507 return array([1.0, -1.0]) | |
508 | |
509 find_fiedler = _get_fiedler_func(method) | |
510 L = nx.laplacian_matrix(G) | |
511 x = None if method != "lobpcg" else _rcm_estimate(G, G) | |
512 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed) | |
513 return fiedler | |
514 | |
515 | |
516 @random_state(5) | |
517 def spectral_ordering( | |
518 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None | |
519 ): | |
520 """Compute the spectral_ordering of a graph. | |
521 | |
522 The spectral ordering of a graph is an ordering of its nodes where nodes | |
523 in the same weakly connected components appear contiguous and ordered by | |
524 their corresponding elements in the Fiedler vector of the component. | |
525 | |
526 Parameters | |
527 ---------- | |
528 G : NetworkX graph | |
529 A graph. | |
530 | |
531 weight : object, optional (default: None) | |
532 The data key used to determine the weight of each edge. If None, then | |
533 each edge has unit weight. | |
534 | |
535 normalized : bool, optional (default: False) | |
536 Whether the normalized Laplacian matrix is used. | |
537 | |
538 tol : float, optional (default: 1e-8) | |
539 Tolerance of relative residual in eigenvalue computation. | |
540 | |
541 method : string, optional (default: 'tracemin_pcg') | |
542 Method of eigenvalue computation. It must be one of the tracemin | |
543 options shown below (TraceMIN), 'lanczos' (Lanczos iteration) | |
544 or 'lobpcg' (LOBPCG). | |
545 | |
546 The TraceMIN algorithm uses a linear system solver. The following | |
547 values allow specifying the solver to be used. | |
548 | |
549 =============== ======================================== | |
550 Value Solver | |
551 =============== ======================================== | |
552 'tracemin_pcg' Preconditioned conjugate gradient method | |
553 'tracemin_chol' Cholesky factorization | |
554 'tracemin_lu' LU factorization | |
555 =============== ======================================== | |
556 | |
557 seed : integer, random_state, or None (default) | |
558 Indicator of random number generation state. | |
559 See :ref:`Randomness<randomness>`. | |
560 | |
561 Returns | |
562 ------- | |
563 spectral_ordering : NumPy array of floats. | |
564 Spectral ordering of nodes. | |
565 | |
566 Raises | |
567 ------ | |
568 NetworkXError | |
569 If G is empty. | |
570 | |
571 Notes | |
572 ----- | |
573 Edge weights are interpreted by their absolute values. For MultiGraph's, | |
574 weights of parallel edges are summed. Zero-weighted edges are ignored. | |
575 | |
576 To use Cholesky factorization in the TraceMIN algorithm, the | |
577 :samp:`scikits.sparse` package must be installed. | |
578 | |
579 See Also | |
580 -------- | |
581 laplacian_matrix | |
582 """ | |
583 if len(G) == 0: | |
584 raise nx.NetworkXError("graph is empty.") | |
585 G = _preprocess_graph(G, weight) | |
586 | |
587 find_fiedler = _get_fiedler_func(method) | |
588 order = [] | |
589 for component in nx.connected_components(G): | |
590 size = len(component) | |
591 if size > 2: | |
592 L = nx.laplacian_matrix(G, component) | |
593 x = None if method != "lobpcg" else _rcm_estimate(G, component) | |
594 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed) | |
595 sort_info = zip(fiedler, range(size), component) | |
596 order.extend(u for x, c, u in sorted(sort_info)) | |
597 else: | |
598 order.extend(component) | |
599 | |
600 return order |