Source code for pyamg.graph

"""Algorithms related to graphs"""

__docformat__ = "restructuredtext en"

import numpy
import scipy
from scipy import sparse

import amg_core

__all__ = ['maximal_independent_set', 'vertex_coloring', 'bellman_ford', \
           'lloyd_cluster', 'connected_components']


def max_value(datatype):
    try:
        return numpy.iinfo(datatype).max
    except:
        return numpy.finfo(datatype).max


def asgraph(G):
    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 ------- An array S where 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 = numpy.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, scipy.rand(N)) else: raise ValueError('unknown algorithm (%s)' % algo) else: fn = amg_core.maximal_independent_set_k_parallel fn(N, G.indptr, G.indices, k, mis, scipy.rand(N)) 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 ------- 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 = numpy.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, scipy.rand(N) ) elif method == 'LDF': fn = amg_core.vertex_coloring_LDF fn(N, G.indptr, G.indices, coloring, scipy.rand(N) ) else: raise ValueError('unknown method (%s)' % method) return coloring
[docs]def bellman_ford(G, seeds, maxiter=None): """ Bellman-Ford iteration Parameters ---------- Returns ------- Notes ----- References ---------- CLR Examples -------- """ G = asgraph(G) N = G.shape[0] if maxiter is not None and maxiter < 0: raise ValueError('maxiter must be positive') if G.dtype == complex: raise ValueError('Bellman-Ford algorithm only defined for real weights') seeds = numpy.asarray(seeds, dtype='intc') distances = numpy.empty(N, dtype=G.dtype) distances[:] = max_value(G.dtype) distances[seeds] = 0 nearest_seed = numpy.empty(N, dtype='intc') nearest_seed[:] = -1 nearest_seed[seeds] = seeds old_distances = numpy.empty_like(distances) iter = 0 while maxiter is None or iter < maxiter: old_distances[:] = distances amg_core.bellman_ford( N, G.indptr, G.indices, G.data, distances, nearest_seed) if (old_distances == distances).all(): break return (distances,nearest_seed)
[docs]def lloyd_cluster(G, seeds, maxiter=10): """Perform Lloyd clustering on graph with weighted edges Parameters ---------- G : csr_matrix or 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. 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 = numpy.abs(G) #interpret seeds argument if numpy.isscalar(seeds): seeds = numpy.random.permutation(N)[:seeds] seeds = seeds.astype('intc') else: seeds = numpy.array(seeds, dtype='intc') if len(seeds) < 1: raise ValueError('at least one seed is required') if seeds.min() < 0: raise ValueError('invalid seed index (%d)' % seeds.min()) if seeds.max() >= N: raise ValueError('invalid seed index (%d)' % seeds.max()) clusters = numpy.empty(N, dtype='intc') distances = numpy.empty(N, dtype=G.dtype) for i in range(maxiter): 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 return (distances, clusters, seeds)
def breadth_first_search(G, seed): """Breadth First search of a graph Parameters ---------- Returns ------- Notes ----- References ---------- CLR Examples -------- """ #TODO document G = asgraph(G) N = G.shape[0] #Check symmetry? order = numpy.empty(N, G.indptr.dtype) level = numpy.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 -------- >>> 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] #Check symmetry? components = numpy.empty(N, G.indptr.dtype) fn = amg_core.connected_components fn(N, G.indptr, G.indices, components) return components
def symmetric_rcm(A): """ Symmetric Reverse Cutthill-McKee Get a pseudo-peripheral node, then call BFS return a symmetric permuted matrix Example ------- >>> import pylab >>> 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 >>> #pylab.figure() >>> #pylab.subplot(121) >>> #pylab.spy(S,marker='.') >>> #pylab.subplot(122) >>> #pylab.spy(symmetric_rcm(S),marker='.') See Also -------- pseudo_peripheral_node """ n = A.shape[0] root, order, level = pseudo_peripheral_node(A) Perm = sparse.identity(n) p = level.argsort() Perm = Perm[p,:] return Perm*A*Perm.T def pseudo_peripheral_node(A): """ Algorithm in Saad """ import numpy from pyamg.graph import breadth_first_search n = A.shape[0] valence = numpy.diff(A.indptr) # select an initial node x, set delta = 0 x = int(numpy.random.rand() * n) delta = 0 while 1: # 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 = numpy.where(level == maxlevel)[0] lastnodesvalence = valence[lastnodes] minlastnodesvalence = lastnodesvalence.min() y = numpy.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