"""Algorithms related to graphs."""
from warnings import warn
import numpy as np
from scipy import sparse
from . import amg_core
[docs]
def asgraph(G):
"""Return (square) matrix as sparse."""
if not (sparse.isspmatrix_csr(G) or sparse.isspmatrix_csc(G)):
G = sparse.csr_matrix(G)
if G.shape[0] != G.shape[1]:
raise ValueError('expected square matrix')
return G
[docs]
def maximal_independent_set(G, algo='serial', k=None):
"""Compute a maximal independent vertex set for a graph.
Parameters
----------
G : sparse matrix
Symmetric matrix, preferably in sparse CSR or CSC format
The nonzeros of G represent the edges of an undirected graph.
algo : {'serial', 'parallel'}
Algorithm used to compute the MIS
* serial : greedy serial algorithm
* parallel : variant of Luby's parallel MIS algorithm
Returns
-------
S : array
S[i] = 1 if vertex i is in the MIS
S[i] = 0 otherwise
Notes
-----
Diagonal entries in the G (self loops) will be ignored.
Luby's algorithm is significantly more expensive than the
greedy serial algorithm.
"""
G = asgraph(G)
N = G.shape[0]
mis = np.empty(N, dtype='intc')
mis[:] = -1
if k is None:
if algo == 'serial':
fn = amg_core.maximal_independent_set_serial
fn(N, G.indptr, G.indices, -1, 1, 0, mis)
elif algo == 'parallel':
fn = amg_core.maximal_independent_set_parallel
fn(N, G.indptr, G.indices, -1, 1, 0, mis, np.random.rand(N), -1)
else:
raise ValueError(f'Unknown algorithm ({algo})')
else:
fn = amg_core.maximal_independent_set_k_parallel
fn(N, G.indptr, G.indices, k, mis, np.random.rand(N), -1)
return mis
[docs]
def vertex_coloring(G, method='MIS'):
"""Compute a vertex coloring of a graph.
Parameters
----------
G : sparse matrix
Symmetric matrix, preferably in sparse CSR or CSC format
The nonzeros of G represent the edges of an undirected graph.
method : string
Algorithm used to compute the vertex coloring:
* 'MIS' - Maximal Independent Set
* 'JP' - Jones-Plassmann (parallel)
* 'LDF' - Largest-Degree-First (parallel)
Returns
-------
coloring : array
An array of vertex colors (integers beginning at 0)
Notes
-----
Diagonal entries in the G (self loops) will be ignored.
"""
G = asgraph(G)
N = G.shape[0]
coloring = np.empty(N, dtype='intc')
if method == 'MIS':
fn = amg_core.vertex_coloring_mis
fn(N, G.indptr, G.indices, coloring)
elif method == 'JP':
fn = amg_core.vertex_coloring_jones_plassmann
fn(N, G.indptr, G.indices, coloring, np.random.rand(N))
elif method == 'LDF':
fn = amg_core.vertex_coloring_LDF
fn(N, G.indptr, G.indices, coloring, np.random.rand(N))
else:
raise ValueError(f'Unknown method ({method})')
return coloring
[docs]
def bellman_ford(G, seeds):
"""Bellman-Ford iteration.
Parameters
----------
G : sparse matrix
Directed graph with positive weights.
seeds : list
Starting seeds
Returns
-------
distances : array
Distance of each point to the nearest seed
nearest_seed : array
Index of the nearest seed
Notes
-----
This should be viewed as the transpose of Bellman-Ford in
scipy.sparse.csgraph. Here, bellman_ford is used to find the shortest path
from any point *to* the seeds. In csgraph, bellman_ford is used to find
"the shortest distance from point i to point j". So csgraph.bellman_ford
could be run `for seed in seeds`. Also note that `test_graph.py` tests
against `csgraph.bellman_ford(G.T)`.
See Also
--------
scipy.sparse.csgraph.bellman_ford
"""
G = asgraph(G)
N = G.shape[0]
if G.nnz > 0:
if G.data.min() < 0:
raise ValueError('Bellman-Ford is defined only for positive weights.')
if G.dtype == complex:
raise ValueError('Bellman-Ford is defined only for real weights.')
seeds = np.asarray(seeds, dtype='intc')
distances = np.full(N, np.inf, dtype=G.dtype)
distances[seeds] = 0
nearest_seed = np.full(N, -1, dtype='intc')
nearest_seed[seeds] = seeds
amg_core.bellman_ford(N, G.indptr, G.indices, G.data, distances, nearest_seed)
return (distances, nearest_seed)
[docs]
def lloyd_cluster(G, seeds, maxiter=10):
"""Perform Lloyd clustering on graph with weighted edges.
Parameters
----------
G : csr_matrix, csc_matrix
A sparse NxN matrix where each nonzero entry G[i,j] is the distance
between nodes i and j.
seeds : int array
If seeds is an integer, then its value determines the number of
clusters. Otherwise, seeds is an array of unique integers between 0
and N-1 that will be used as the initial seeds for clustering.
maxiter : int
The maximum number of iterations to perform.
Returns
-------
distances : array
final distances
clusters : int array
id of each cluster of points
seeds : int array
index of each seed
Notes
-----
If G has complex values, abs(G) is used instead.
"""
G = asgraph(G)
N = G.shape[0]
if G.dtype.kind == 'c':
# complex dtype
G = np.abs(G)
# interpret seeds argument
if np.isscalar(seeds):
seeds = np.random.permutation(N)[:seeds]
seeds = seeds.astype('intc')
else:
seeds = np.array(seeds, dtype='intc')
if len(seeds) < 1:
raise ValueError('at least one seed is required')
if seeds.min() < 0:
raise ValueError(f'Invalid seed index ({seeds.min()})')
if seeds.max() >= N:
raise ValueError(f'Invalid seed index ({seeds.max()})')
clusters = np.empty(N, dtype='intc')
distances = np.empty(N, dtype=G.dtype)
for _it in range(1, maxiter+1):
last_seeds = seeds.copy()
amg_core.lloyd_cluster(N, G.indptr, G.indices, G.data,
len(seeds), distances, clusters, seeds)
if (seeds == last_seeds).all():
break
if _it == maxiter:
warn('Lloyd clustering reached maxiter (did not converge)')
return (distances, clusters, seeds)
[docs]
def breadth_first_search(G, seed):
"""Breadth First search of a graph.
Parameters
----------
G : csr_matrix, csc_matrix
A sparse NxN matrix where each nonzero entry G[i,j] is the distance
between nodes i and j.
seed : int
Index of the seed location
Returns
-------
order : int array
Breadth first order
level : int array
Final levels
Examples
--------
0---2
| /
| /
1---4---7---8---9
| /| /
| / | /
3/ 6/
|
|
5
>>> import numpy as np
>>> import pyamg
>>> import scipy.sparse as sparse
>>> edges = np.array([[0,1],[0,2],[1,2],[1,3],[1,4],[3,4],[3,5],
... [4,6], [4,7], [6,7], [7,8], [8,9]])
>>> N = np.max(edges.ravel())+1
>>> data = np.ones((edges.shape[0],))
>>> A = sparse.coo_matrix((data, (edges[:,0], edges[:,1])), shape=(N,N))
>>> c, l = pyamg.graph.breadth_first_search(A, 0)
>>> print(l)
[0 1 1 2 2 3 3 3 4 5]
>>> print(c)
[0 1 2 3 4 5 6 7 8 9]
"""
G = asgraph(G)
N = G.shape[0]
order = np.empty(N, G.indptr.dtype)
level = np.empty(N, G.indptr.dtype)
level[:] = -1
BFS = amg_core.breadth_first_search
BFS(G.indptr, G.indices, int(seed), order, level)
return order, level
[docs]
def connected_components(G):
"""Compute the connected components of a graph.
The connected components of a graph G, which is represented by a
symmetric sparse matrix, are labeled with the integers 0,1,..(K-1) where
K is the number of components.
Parameters
----------
G : symmetric matrix, preferably in sparse CSR or CSC format
The nonzeros of G represent the edges of an undirected graph.
Returns
-------
components : ndarray
An array of component labels for each vertex of the graph.
Notes
-----
If the nonzero structure of G is not symmetric, then the
result is undefined.
Examples
--------
>>> from pyamg.graph import connected_components
>>> print(connected_components( [[0,1,0],[1,0,1],[0,1,0]] ))
[0 0 0]
>>> print(connected_components( [[0,1,0],[1,0,0],[0,0,0]] ))
[0 0 1]
>>> print(connected_components( [[0,0,0],[0,0,0],[0,0,0]] ))
[0 1 2]
>>> print(connected_components( [[0,1,0,0],[1,0,0,0],[0,0,0,1],[0,0,1,0]] ))
[0 0 1 1]
"""
G = asgraph(G)
N = G.shape[0]
components = np.empty(N, G.indptr.dtype)
fn = amg_core.connected_components
fn(N, G.indptr, G.indices, components)
return components
[docs]
def symmetric_rcm(A):
"""Symmetric Reverse Cutthill-McKee.
Parameters
----------
A : sparse matrix
Sparse matrix
Returns
-------
B : sparse matrix
Permuted matrix with reordering
Notes
-----
Get a pseudo-peripheral node, then call BFS
Examples
--------
>>> from pyamg import gallery
>>> from pyamg.graph import symmetric_rcm
>>> n = 200
>>> density = 1.0/n
>>> A = gallery.sprand(n, n, density, format='csr')
>>> S = A + A.T
>>> # try the visualizations
>>> # import matplotlib.pyplot as plt
>>> # plt.figure()
>>> # plt.subplot(121)
>>> # plt.spy(S,marker='.')
>>> # plt.subplot(122)
>>> # plt.spy(symmetric_rcm(S),marker='.')
See Also
--------
pseudo_peripheral_node
"""
dummy_root, order, dummy_level = pseudo_peripheral_node(A)
p = order[::-1]
return A[p, :][:, p]
[docs]
def pseudo_peripheral_node(A):
"""Find a pseudo peripheral node.
Parameters
----------
A : sparse matrix
Sparse matrix
Returns
-------
x : int
Location of the node
order : array
BFS ordering
level : array
BFS levels
Notes
-----
Algorithm in Saad
"""
n = A.shape[0]
valence = np.diff(A.indptr)
# select an initial node x, set delta = 0
x = int(np.random.rand() * n)
delta = 0
while True:
# do a level-set traversal from x
order, level = breadth_first_search(A, x)
# select a node y in the last level with min degree
maxlevel = level.max()
lastnodes = np.where(level == maxlevel)[0]
lastnodesvalence = valence[lastnodes]
minlastnodesvalence = lastnodesvalence.min()
y = np.where(lastnodesvalence == minlastnodesvalence)[0][0]
y = lastnodes[y]
# if d(x,y)>delta, set, and go to bfs above
if level[y] > delta:
x = y
delta = level[y]
else:
return x, order, level
