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