Statistics
| Branch: | Revision:

## iof-tools / networkxMiCe / networkx-master / networkx / algorithms / threshold.py @ 5cef0f13

 1 ```# Copyright (C) 2004-2019 by ``` ```# Aric Hagberg ``` ```# Dan Schult ``` ```# Pieter Swart ``` ```# All rights reserved. ``` ```# BSD license. ``` ```# ``` ```# Authors: Aric Hagberg (hagberg@lanl.gov) ``` ```# Pieter Swart (swart@lanl.gov) ``` ```# Dan Schult (dschult@colgate.edu) ``` ```""" ``` ```Threshold Graphs - Creation, manipulation and identification. ``` ```""" ``` ```from math import sqrt ``` ```import networkx as nx ``` ```from networkx.utils import py_random_state ``` ```__all__ = ['is_threshold_graph', 'find_threshold_graph'] ``` ```def is_threshold_graph(G): ``` ``` """ ``` ``` Returns True if G is a threshold graph. ``` ``` """ ``` ``` return is_threshold_sequence(list(d for n, d in G.degree())) ``` ```def is_threshold_sequence(degree_sequence): ``` ``` """ ``` ``` Returns True if the sequence is a threshold degree seqeunce. ``` ``` ``` ``` Uses the property that a threshold graph must be constructed by ``` ``` adding either dominating or isolated nodes. Thus, it can be ``` ``` deconstructed iteratively by removing a node of degree zero or a ``` ``` node that connects to the remaining nodes. If this deconstruction ``` ``` failes then the sequence is not a threshold sequence. ``` ``` """ ``` ``` ds = degree_sequence[:] # get a copy so we don't destroy original ``` ``` ds.sort() ``` ``` while ds: ``` ``` if ds == 0: # if isolated node ``` ``` ds.pop(0) # remove it ``` ``` continue ``` ``` if ds[-1] != len(ds) - 1: # is the largest degree node dominating? ``` ``` return False # no, not a threshold degree sequence ``` ``` ds.pop() # yes, largest is the dominating node ``` ``` ds = [d - 1 for d in ds] # remove it and decrement all degrees ``` ``` return True ``` ```def creation_sequence(degree_sequence, with_labels=False, compact=False): ``` ``` """ ``` ``` Determines the creation sequence for the given threshold degree sequence. ``` ``` ``` ``` The creation sequence is a list of single characters 'd' ``` ``` or 'i': 'd' for dominating or 'i' for isolated vertices. ``` ``` Dominating vertices are connected to all vertices present when it ``` ``` is added. The first node added is by convention 'd'. ``` ``` This list can be converted to a string if desired using "".join(cs) ``` ``` ``` ``` If with_labels==True: ``` ``` Returns a list of 2-tuples containing the vertex number ``` ``` and a character 'd' or 'i' which describes the type of vertex. ``` ``` ``` ``` If compact==True: ``` ``` Returns the creation sequence in a compact form that is the number ``` ``` of 'i's and 'd's alternating. ``` ``` Examples: ``` ``` [1,2,2,3] represents d,i,i,d,d,i,i,i ``` ``` [3,1,2] represents d,d,d,i,d,d ``` ``` ``` ``` Notice that the first number is the first vertex to be used for ``` ``` construction and so is always 'd'. ``` ``` ``` ``` with_labels and compact cannot both be True. ``` ``` ``` ``` Returns None if the sequence is not a threshold sequence ``` ``` """ ``` ``` if with_labels and compact: ``` ``` raise ValueError("compact sequences cannot be labeled") ``` ``` # make an indexed copy ``` ``` if isinstance(degree_sequence, dict): # labeled degree seqeunce ``` ``` ds = [[degree, label] for (label, degree) in degree_sequence.items()] ``` ``` else: ``` ``` ds = [[d, i] for i, d in enumerate(degree_sequence)] ``` ``` ds.sort() ``` ``` cs = [] # creation sequence ``` ``` while ds: ``` ``` if ds == 0: # isolated node ``` ``` (d, v) = ds.pop(0) ``` ``` if len(ds) > 0: # make sure we start with a d ``` ``` cs.insert(0, (v, 'i')) ``` ``` else: ``` ``` cs.insert(0, (v, 'd')) ``` ``` continue ``` ``` if ds[-1] != len(ds) - 1: # Not dominating node ``` ``` return None # not a threshold degree sequence ``` ``` (d, v) = ds.pop() ``` ``` cs.insert(0, (v, 'd')) ``` ``` ds = [[d - 1, d] for d in ds] # decrement due to removing node ``` ``` if with_labels: ``` ``` return cs ``` ``` if compact: ``` ``` return make_compact(cs) ``` ``` return [v for v in cs] # not labeled ``` ```def make_compact(creation_sequence): ``` ``` """ ``` ``` Returns the creation sequence in a compact form ``` ``` that is the number of 'i's and 'd's alternating. ``` ``` ``` ``` Examples ``` ``` -------- ``` ``` >>> from networkx.algorithms.threshold import make_compact ``` ``` >>> make_compact(['d', 'i', 'i', 'd', 'd', 'i', 'i', 'i']) ``` ``` [1, 2, 2, 3] ``` ``` >>> make_compact(['d', 'd', 'd', 'i', 'd', 'd']) ``` ``` [3, 1, 2] ``` ``` ``` ``` Notice that the first number is the first vertex ``` ``` to be used for construction and so is always 'd'. ``` ``` ``` ``` Labeled creation sequences lose their labels in the ``` ``` compact representation. ``` ``` ``` ``` >>> make_compact([3, 1, 2]) ``` ``` [3, 1, 2] ``` ``` """ ``` ``` first = creation_sequence ``` ``` if isinstance(first, str): # creation sequence ``` ``` cs = creation_sequence[:] ``` ``` elif isinstance(first, tuple): # labeled creation sequence ``` ``` cs = [s for s in creation_sequence] ``` ``` elif isinstance(first, int): # compact creation sequence ``` ``` return creation_sequence ``` ``` else: ``` ``` raise TypeError("Not a valid creation sequence type") ``` ``` ccs = [] ``` ``` count = 1 # count the run lengths of d's or i's. ``` ``` for i in range(1, len(cs)): ``` ``` if cs[i] == cs[i - 1]: ``` ``` count += 1 ``` ``` else: ``` ``` ccs.append(count) ``` ``` count = 1 ``` ``` ccs.append(count) # don't forget the last one ``` ``` return ccs ``` ```def uncompact(creation_sequence): ``` ``` """ ``` ``` Converts a compact creation sequence for a threshold ``` ``` graph to a standard creation sequence (unlabeled). ``` ``` If the creation_sequence is already standard, return it. ``` ``` See creation_sequence. ``` ``` """ ``` ``` first = creation_sequence ``` ``` if isinstance(first, str): # creation sequence ``` ``` return creation_sequence ``` ``` elif isinstance(first, tuple): # labeled creation sequence ``` ``` return creation_sequence ``` ``` elif isinstance(first, int): # compact creation sequence ``` ``` ccscopy = creation_sequence[:] ``` ``` else: ``` ``` raise TypeError("Not a valid creation sequence type") ``` ``` cs = [] ``` ``` while ccscopy: ``` ``` cs.extend(ccscopy.pop(0) * ['d']) ``` ``` if ccscopy: ``` ``` cs.extend(ccscopy.pop(0) * ['i']) ``` ``` return cs ``` ```def creation_sequence_to_weights(creation_sequence): ``` ``` """ ``` ``` Returns a list of node weights which create the threshold ``` ``` graph designated by the creation sequence. The weights ``` ``` are scaled so that the threshold is 1.0. The order of the ``` ``` nodes is the same as that in the creation sequence. ``` ``` """ ``` ``` # Turn input sequence into a labeled creation sequence ``` ``` first = creation_sequence ``` ``` if isinstance(first, str): # creation sequence ``` ``` if isinstance(creation_sequence, list): ``` ``` wseq = creation_sequence[:] ``` ``` else: ``` ``` wseq = list(creation_sequence) # string like 'ddidid' ``` ``` elif isinstance(first, tuple): # labeled creation sequence ``` ``` wseq = [v for v in creation_sequence] ``` ``` elif isinstance(first, int): # compact creation sequence ``` ``` wseq = uncompact(creation_sequence) ``` ``` else: ``` ``` raise TypeError("Not a valid creation sequence type") ``` ``` # pass through twice--first backwards ``` ``` wseq.reverse() ``` ``` w = 0 ``` ``` prev = 'i' ``` ``` for j, s in enumerate(wseq): ``` ``` if s == 'i': ``` ``` wseq[j] = w ``` ``` prev = s ``` ``` elif prev == 'i': ``` ``` prev = s ``` ``` w += 1 ``` ``` wseq.reverse() # now pass through forwards ``` ``` for j, s in enumerate(wseq): ``` ``` if s == 'd': ``` ``` wseq[j] = w ``` ``` prev = s ``` ``` elif prev == 'd': ``` ``` prev = s ``` ``` w += 1 ``` ``` # Now scale weights ``` ``` if prev == 'd': ``` ``` w += 1 ``` ``` wscale = 1. / float(w) ``` ``` return [ww * wscale for ww in wseq] ``` ``` # return wseq ``` ```def weights_to_creation_sequence(weights, threshold=1, with_labels=False, compact=False): ``` ``` """ ``` ``` Returns a creation sequence for a threshold graph ``` ``` determined by the weights and threshold given as input. ``` ``` If the sum of two node weights is greater than the ``` ``` threshold value, an edge is created between these nodes. ``` ``` ``` ``` The creation sequence is a list of single characters 'd' ``` ``` or 'i': 'd' for dominating or 'i' for isolated vertices. ``` ``` Dominating vertices are connected to all vertices present ``` ``` when it is added. The first node added is by convention 'd'. ``` ``` ``` ``` If with_labels==True: ``` ``` Returns a list of 2-tuples containing the vertex number ``` ``` and a character 'd' or 'i' which describes the type of vertex. ``` ``` ``` ``` If compact==True: ``` ``` Returns the creation sequence in a compact form that is the number ``` ``` of 'i's and 'd's alternating. ``` ``` Examples: ``` ``` [1,2,2,3] represents d,i,i,d,d,i,i,i ``` ``` [3,1,2] represents d,d,d,i,d,d ``` ``` ``` ``` Notice that the first number is the first vertex to be used for ``` ``` construction and so is always 'd'. ``` ``` ``` ``` with_labels and compact cannot both be True. ``` ``` """ ``` ``` if with_labels and compact: ``` ``` raise ValueError("compact sequences cannot be labeled") ``` ``` # make an indexed copy ``` ``` if isinstance(weights, dict): # labeled weights ``` ``` wseq = [[w, label] for (label, w) in weights.items()] ``` ``` else: ``` ``` wseq = [[w, i] for i, w in enumerate(weights)] ``` ``` wseq.sort() ``` ``` cs = [] # creation sequence ``` ``` cutoff = threshold - wseq[-1] ``` ``` while wseq: ``` ``` if wseq < cutoff: # isolated node ``` ``` (w, label) = wseq.pop(0) ``` ``` cs.append((label, 'i')) ``` ``` else: ``` ``` (w, label) = wseq.pop() ``` ``` cs.append((label, 'd')) ``` ``` cutoff = threshold - wseq[-1] ``` ``` if len(wseq) == 1: # make sure we start with a d ``` ``` (w, label) = wseq.pop() ``` ``` cs.append((label, 'd')) ``` ``` # put in correct order ``` ``` cs.reverse() ``` ``` if with_labels: ``` ``` return cs ``` ``` if compact: ``` ``` return make_compact(cs) ``` ``` return [v for v in cs] # not labeled ``` ```# Manipulating NetworkX.Graphs in context of threshold graphs ``` ```def threshold_graph(creation_sequence, create_using=None): ``` ``` """ ``` ``` Create a threshold graph from the creation sequence or compact ``` ``` creation_sequence. ``` ``` ``` ``` The input sequence can be a ``` ``` ``` ``` creation sequence (e.g. ['d','i','d','d','d','i']) ``` ``` labeled creation sequence (e.g. [(0,'d'),(2,'d'),(1,'i')]) ``` ``` compact creation sequence (e.g. [2,1,1,2,0]) ``` ``` ``` ``` Use cs=creation_sequence(degree_sequence,labeled=True) ``` ``` to convert a degree sequence to a creation sequence. ``` ``` ``` ``` Returns None if the sequence is not valid ``` ``` """ ``` ``` # Turn input sequence into a labeled creation sequence ``` ``` first = creation_sequence ``` ``` if isinstance(first, str): # creation sequence ``` ``` ci = list(enumerate(creation_sequence)) ``` ``` elif isinstance(first, tuple): # labeled creation sequence ``` ``` ci = creation_sequence[:] ``` ``` elif isinstance(first, int): # compact creation sequence ``` ``` cs = uncompact(creation_sequence) ``` ``` ci = list(enumerate(cs)) ``` ``` else: ``` ``` print("not a valid creation sequence type") ``` ``` return None ``` ``` G = nx.empty_graph(0, create_using) ``` ``` if G.is_directed(): ``` ``` raise nx.NetworkXError("Directed Graph not supported") ``` ``` G.name = "Threshold Graph" ``` ``` # add nodes and edges ``` ``` # if type is 'i' just add nodea ``` ``` # if type is a d connect to everything previous ``` ``` while ci: ``` ``` (v, node_type) = ci.pop(0) ``` ``` if node_type == 'd': # dominating type, connect to all existing nodes ``` ``` # We use `for u in list(G):` instead of ``` ``` # `for u in G:` because we edit the graph `G` in ``` ``` # the loop. Hence using an iterator will result in ``` ``` # `RuntimeError: dictionary changed size during iteration` ``` ``` for u in list(G): ``` ``` G.add_edge(v, u) ``` ``` G.add_node(v) ``` ``` return G ``` ```def find_alternating_4_cycle(G): ``` ``` """ ``` ``` Returns False if there aren't any alternating 4 cycles. ``` ``` Otherwise returns the cycle as [a,b,c,d] where (a,b) ``` ``` and (c,d) are edges and (a,c) and (b,d) are not. ``` ``` """ ``` ``` for (u, v) in G.edges(): ``` ``` for w in G.nodes(): ``` ``` if not G.has_edge(u, w) and u != w: ``` ``` for x in G.neighbors(w): ``` ``` if not G.has_edge(v, x) and v != x: ``` ``` return [u, v, w, x] ``` ``` return False ``` ```def find_threshold_graph(G, create_using=None): ``` ``` """ ``` ``` Return a threshold subgraph that is close to largest in G. ``` ``` The threshold graph will contain the largest degree node in G. ``` ``` ``` ``` """ ``` ``` return threshold_graph(find_creation_sequence(G), create_using) ``` ```def find_creation_sequence(G): ``` ``` """ ``` ``` Find a threshold subgraph that is close to largest in G. ``` ``` Returns the labeled creation sequence of that threshold graph. ``` ``` """ ``` ``` cs = [] ``` ``` # get a local pointer to the working part of the graph ``` ``` H = G ``` ``` while H.order() > 0: ``` ``` # get new degree sequence on subgraph ``` ``` dsdict = dict(H.degree()) ``` ``` ds = [(d, v) for v, d in dsdict.items()] ``` ``` ds.sort() ``` ``` # Update threshold graph nodes ``` ``` if ds[-1] == 0: # all are isolated ``` ``` cs.extend(zip(dsdict, ['i'] * (len(ds) - 1) + ['d'])) ``` ``` break # Done! ``` ``` # pull off isolated nodes ``` ``` while ds == 0: ``` ``` (d, iso) = ds.pop(0) ``` ``` cs.append((iso, 'i')) ``` ``` # find new biggest node ``` ``` (d, bigv) = ds.pop() ``` ``` # add edges of star to t_g ``` ``` cs.append((bigv, 'd')) ``` ``` # form subgraph of neighbors of big node ``` ``` H = H.subgraph(H.neighbors(bigv)) ``` ``` cs.reverse() ``` ``` return cs ``` ```# Properties of Threshold Graphs ``` ```def triangles(creation_sequence): ``` ``` """ ``` ``` Compute number of triangles in the threshold graph with the ``` ``` given creation sequence. ``` ``` """ ``` ``` # shortcut algorithm that doesn't require computing number ``` ``` # of triangles at each node. ``` ``` cs = creation_sequence # alias ``` ``` dr = cs.count("d") # number of d's in sequence ``` ``` ntri = dr * (dr - 1) * (dr - 2) / 6 # number of triangles in clique of nd d's ``` ``` # now add dr choose 2 triangles for every 'i' in sequence where ``` ``` # dr is the number of d's to the right of the current i ``` ``` for i, typ in enumerate(cs): ``` ``` if typ == "i": ``` ``` ntri += dr * (dr - 1) / 2 ``` ``` else: ``` ``` dr -= 1 ``` ``` return ntri ``` ```def triangle_sequence(creation_sequence): ``` ``` """ ``` ``` Return triangle sequence for the given threshold graph creation sequence. ``` ``` ``` ``` """ ``` ``` cs = creation_sequence ``` ``` seq = [] ``` ``` dr = cs.count("d") # number of d's to the right of the current pos ``` ``` dcur = (dr - 1) * (dr - 2) // 2 # number of triangles through a node of clique dr ``` ``` irun = 0 # number of i's in the last run ``` ``` drun = 0 # number of d's in the last run ``` ``` for i, sym in enumerate(cs): ``` ``` if sym == "d": ``` ``` drun += 1 ``` ``` tri = dcur + (dr - 1) * irun # new triangles at this d ``` ``` else: # cs[i]="i": ``` ``` if prevsym == "d": # new string of i's ``` ``` dcur += (dr - 1) * irun # accumulate shared shortest paths ``` ``` irun = 0 # reset i run counter ``` ``` dr -= drun # reduce number of d's to right ``` ``` drun = 0 # reset d run counter ``` ``` irun += 1 ``` ``` tri = dr * (dr - 1) // 2 # new triangles at this i ``` ``` seq.append(tri) ``` ``` prevsym = sym ``` ``` return seq ``` ```def cluster_sequence(creation_sequence): ``` ``` """ ``` ``` Return cluster sequence for the given threshold graph creation sequence. ``` ``` """ ``` ``` triseq = triangle_sequence(creation_sequence) ``` ``` degseq = degree_sequence(creation_sequence) ``` ``` cseq = [] ``` ``` for i, deg in enumerate(degseq): ``` ``` tri = triseq[i] ``` ``` if deg <= 1: # isolated vertex or single pair gets cc 0 ``` ``` cseq.append(0) ``` ``` continue ``` ``` max_size = (deg * (deg - 1)) // 2 ``` ``` cseq.append(float(tri) / float(max_size)) ``` ``` return cseq ``` ```def degree_sequence(creation_sequence): ``` ``` """ ``` ``` Return degree sequence for the threshold graph with the given ``` ``` creation sequence ``` ``` """ ``` ``` cs = creation_sequence # alias ``` ``` seq = [] ``` ``` rd = cs.count("d") # number of d to the right ``` ``` for i, sym in enumerate(cs): ``` ``` if sym == "d": ``` ``` rd -= 1 ``` ``` seq.append(rd + i) ``` ``` else: ``` ``` seq.append(rd) ``` ``` return seq ``` ```def density(creation_sequence): ``` ``` """ ``` ``` Return the density of the graph with this creation_sequence. ``` ``` The density is the fraction of possible edges present. ``` ``` """ ``` ``` N = len(creation_sequence) ``` ``` two_size = sum(degree_sequence(creation_sequence)) ``` ``` two_possible = N * (N - 1) ``` ``` den = two_size / float(two_possible) ``` ``` return den ``` ```def degree_correlation(creation_sequence): ``` ``` """ ``` ``` Return the degree-degree correlation over all edges. ``` ``` """ ``` ``` cs = creation_sequence ``` ``` s1 = 0 # deg_i*deg_j ``` ``` s2 = 0 # deg_i^2+deg_j^2 ``` ``` s3 = 0 # deg_i+deg_j ``` ``` m = 0 # number of edges ``` ``` rd = cs.count("d") # number of d nodes to the right ``` ``` rdi = [i for i, sym in enumerate(cs) if sym == "d"] # index of "d"s ``` ``` ds = degree_sequence(cs) ``` ``` for i, sym in enumerate(cs): ``` ``` if sym == "d": ``` ``` if i != rdi: ``` ``` print("Logic error in degree_correlation", i, rdi) ``` ``` raise ValueError ``` ``` rdi.pop(0) ``` ``` degi = ds[i] ``` ``` for dj in rdi: ``` ``` degj = ds[dj] ``` ``` s1 += degj * degi ``` ``` s2 += degi**2 + degj**2 ``` ``` s3 += degi + degj ``` ``` m += 1 ``` ``` denom = (2 * m * s2 - s3 * s3) ``` ``` numer = (4 * m * s1 - s3 * s3) ``` ``` if denom == 0: ``` ``` if numer == 0: ``` ``` return 1 ``` ``` raise ValueError("Zero Denominator but Numerator is %s" % numer) ``` ``` return numer / float(denom) ``` ```def shortest_path(creation_sequence, u, v): ``` ``` """ ``` ``` Find the shortest path between u and v in a ``` ``` threshold graph G with the given creation_sequence. ``` ``` ``` ``` For an unlabeled creation_sequence, the vertices ``` ``` u and v must be integers in (0,len(sequence)) referring ``` ``` to the position of the desired vertices in the sequence. ``` ``` ``` ``` For a labeled creation_sequence, u and v are labels of veritices. ``` ``` ``` ``` Use cs=creation_sequence(degree_sequence,with_labels=True) ``` ``` to convert a degree sequence to a creation sequence. ``` ``` ``` ``` Returns a list of vertices from u to v. ``` ``` Example: if they are neighbors, it returns [u,v] ``` ``` """ ``` ``` # Turn input sequence into a labeled creation sequence ``` ``` first = creation_sequence ``` ``` if isinstance(first, str): # creation sequence ``` ``` cs = [(i, creation_sequence[i]) for i in range(len(creation_sequence))] ``` ``` elif isinstance(first, tuple): # labeled creation sequence ``` ``` cs = creation_sequence[:] ``` ``` elif isinstance(first, int): # compact creation sequence ``` ``` ci = uncompact(creation_sequence) ``` ``` cs = [(i, ci[i]) for i in range(len(ci))] ``` ``` else: ``` ``` raise TypeError("Not a valid creation sequence type") ``` ``` verts = [s for s in cs] ``` ``` if v not in verts: ``` ``` raise ValueError("Vertex %s not in graph from creation_sequence" % v) ``` ``` if u not in verts: ``` ``` raise ValueError("Vertex %s not in graph from creation_sequence" % u) ``` ``` # Done checking ``` ``` if u == v: ``` ``` return [u] ``` ``` uindex = verts.index(u) ``` ``` vindex = verts.index(v) ``` ``` bigind = max(uindex, vindex) ``` ``` if cs[bigind] == 'd': ``` ``` return [u, v] ``` ``` # must be that cs[bigind]=='i' ``` ``` cs = cs[bigind:] ``` ``` while cs: ``` ``` vert = cs.pop() ``` ``` if vert == 'd': ``` ``` return [u, vert, v] ``` ``` # All after u are type 'i' so no connection ``` ``` return -1 ``` ```def shortest_path_length(creation_sequence, i): ``` ``` """ ``` ``` Return the shortest path length from indicated node to ``` ``` every other node for the threshold graph with the given ``` ``` creation sequence. ``` ``` Node is indicated by index i in creation_sequence unless ``` ``` creation_sequence is labeled in which case, i is taken to ``` ``` be the label of the node. ``` ``` ``` ``` Paths lengths in threshold graphs are at most 2. ``` ``` Length to unreachable nodes is set to -1. ``` ``` """ ``` ``` # Turn input sequence into a labeled creation sequence ``` ``` first = creation_sequence ``` ``` if isinstance(first, str): # creation sequence ``` ``` if isinstance(creation_sequence, list): ``` ``` cs = creation_sequence[:] ``` ``` else: ``` ``` cs = list(creation_sequence) ``` ``` elif isinstance(first, tuple): # labeled creation sequence ``` ``` cs = [v for v in creation_sequence] ``` ``` i = [v for v in creation_sequence].index(i) ``` ``` elif isinstance(first, int): # compact creation sequence ``` ``` cs = uncompact(creation_sequence) ``` ``` else: ``` ``` raise TypeError("Not a valid creation sequence type") ``` ``` # Compute ``` ``` N = len(cs) ``` ``` spl =  * N # length 2 to every node ``` ``` spl[i] = 0 # except self which is 0 ``` ``` # 1 for all d's to the right ``` ``` for j in range(i + 1, N): ``` ``` if cs[j] == "d": ``` ``` spl[j] = 1 ``` ``` if cs[i] == 'd': # 1 for all nodes to the left ``` ``` for j in range(i): ``` ``` spl[j] = 1 ``` ``` # and -1 for any trailing i to indicate unreachable ``` ``` for j in range(N - 1, 0, -1): ``` ``` if cs[j] == "d": ``` ``` break ``` ``` spl[j] = -1 ``` ``` return spl ``` ```def betweenness_sequence(creation_sequence, normalized=True): ``` ``` """ ``` ``` Return betweenness for the threshold graph with the given creation ``` ``` sequence. The result is unscaled. To scale the values ``` ``` to the iterval [0,1] divide by (n-1)*(n-2). ``` ``` """ ``` ``` cs = creation_sequence ``` ``` seq = [] # betweenness ``` ``` lastchar = 'd' # first node is always a 'd' ``` ``` dr = float(cs.count("d")) # number of d's to the right of curren pos ``` ``` irun = 0 # number of i's in the last run ``` ``` drun = 0 # number of d's in the last run ``` ``` dlast = 0.0 # betweenness of last d ``` ``` for i, c in enumerate(cs): ``` ``` if c == 'd': # cs[i]=="d": ``` ``` # betweennees = amt shared with eariler d's and i's ``` ``` # + new isolated nodes covered ``` ``` # + new paths to all previous nodes ``` ``` b = dlast + (irun - 1) * irun / dr + 2 * irun * (i - drun - irun) / dr ``` ``` drun += 1 # update counter ``` ``` else: # cs[i]="i": ``` ``` if lastchar == 'd': # if this is a new run of i's ``` ``` dlast = b # accumulate betweenness ``` ``` dr -= drun # update number of d's to the right ``` ``` drun = 0 # reset d counter ``` ``` irun = 0 # reset i counter ``` ``` b = 0 # isolated nodes have zero betweenness ``` ``` irun += 1 # add another i to the run ``` ``` seq.append(float(b)) ``` ``` lastchar = c ``` ``` # normalize by the number of possible shortest paths ``` ``` if normalized: ``` ``` order = len(cs) ``` ``` scale = 1.0 / ((order - 1) * (order - 2)) ``` ``` seq = [s * scale for s in seq] ``` ``` return seq ``` ```def eigenvectors(creation_sequence): ``` ``` """ ``` ``` Return a 2-tuple of Laplacian eigenvalues and eigenvectors ``` ``` for the threshold network with creation_sequence. ``` ``` The first value is a list of eigenvalues. ``` ``` The second value is a list of eigenvectors. ``` ``` The lists are in the same order so corresponding eigenvectors ``` ``` and eigenvalues are in the same position in the two lists. ``` ``` ``` ``` Notice that the order of the eigenvalues returned by eigenvalues(cs) ``` ``` may not correspond to the order of these eigenvectors. ``` ``` """ ``` ``` ccs = make_compact(creation_sequence) ``` ``` N = sum(ccs) ``` ``` vec =  * N ``` ``` val = vec[:] ``` ``` # get number of type d nodes to the right (all for first node) ``` ``` dr = sum(ccs[::2]) ``` ``` nn = ccs ``` ``` vec = [1. / sqrt(N)] * N ``` ``` val = 0 ``` ``` e = dr ``` ``` dr -= nn ``` ``` type_d = True ``` ``` i = 1 ``` ``` dd = 1 ``` ``` while dd < nn: ``` ``` scale = 1. / sqrt(dd * dd + i) ``` ``` vec[i] = i * [-scale] + [dd * scale] +  * (N - i - 1) ``` ``` val[i] = e ``` ``` i += 1 ``` ``` dd += 1 ``` ``` if len(ccs) == 1: ``` ``` return (val, vec) ``` ``` for nn in ccs[1:]: ``` ``` scale = 1. / sqrt(nn * i * (i + nn)) ``` ``` vec[i] = i * [-nn * scale] + nn * [i * scale] +  * (N - i - nn) ``` ``` # find eigenvalue ``` ``` type_d = not type_d ``` ``` if type_d: ``` ``` e = i + dr ``` ``` dr -= nn ``` ``` else: ``` ``` e = dr ``` ``` val[i] = e ``` ``` st = i ``` ``` i += 1 ``` ``` dd = 1 ``` ``` while dd < nn: ``` ``` scale = 1. / sqrt(i - st + dd * dd) ``` ``` vec[i] =  * st + (i - st) * [-scale] + [dd * scale] +  * (N - i - 1) ``` ``` val[i] = e ``` ``` i += 1 ``` ``` dd += 1 ``` ``` return (val, vec) ``` ```def spectral_projection(u, eigenpairs): ``` ``` """ ``` ``` Returns the coefficients of each eigenvector ``` ``` in a projection of the vector u onto the normalized ``` ``` eigenvectors which are contained in eigenpairs. ``` ``` ``` ``` eigenpairs should be a list of two objects. The ``` ``` first is a list of eigenvalues and the second a list ``` ``` of eigenvectors. The eigenvectors should be lists. ``` ``` ``` ``` There's not a lot of error checking on lengths of ``` ``` arrays, etc. so be careful. ``` ``` """ ``` ``` coeff = [] ``` ``` evect = eigenpairs ``` ``` for ev in evect: ``` ``` c = sum([evv * uv for (evv, uv) in zip(ev, u)]) ``` ``` coeff.append(c) ``` ``` return coeff ``` ```def eigenvalues(creation_sequence): ``` ``` """ ``` ``` Return sequence of eigenvalues of the Laplacian of the threshold ``` ``` graph for the given creation_sequence. ``` ``` ``` ``` Based on the Ferrer's diagram method. The spectrum is integral ``` ``` and is the conjugate of the degree sequence. ``` ``` ``` ``` See:: ``` ``` ``` ``` @Article{degree-merris-1994, ``` ``` author = {Russel Merris}, ``` ``` title = {Degree maximal graphs are Laplacian integral}, ``` ``` journal = {Linear Algebra Appl.}, ``` ``` year = {1994}, ``` ``` volume = {199}, ``` ``` pages = {381--389}, ``` ``` } ``` ``` ``` ``` """ ``` ``` degseq = degree_sequence(creation_sequence) ``` ``` degseq.sort() ``` ``` eiglist = [] # zero is always one eigenvalue ``` ``` eig = 0 ``` ``` row = len(degseq) ``` ``` bigdeg = degseq.pop() ``` ``` while row: ``` ``` if bigdeg < row: ``` ``` eiglist.append(eig) ``` ``` row -= 1 ``` ``` else: ``` ``` eig += 1 ``` ``` if degseq: ``` ``` bigdeg = degseq.pop() ``` ``` else: ``` ``` bigdeg = 0 ``` ``` return eiglist ``` ```# Threshold graph creation routines ``` ```@py_random_state(2) ``` ```def random_threshold_sequence(n, p, seed=None): ``` ``` """ ``` ``` Create a random threshold sequence of size n. ``` ``` A creation sequence is built by randomly choosing d's with ``` ``` probabiliy p and i's with probability 1-p. ``` ``` ``` ``` s=nx.random_threshold_sequence(10,0.5) ``` ``` ``` ``` returns a threshold sequence of length 10 with equal ``` ``` probably of an i or a d at each position. ``` ``` ``` ``` A "random" threshold graph can be built with ``` ``` ``` ``` G=nx.threshold_graph(s) ``` ``` ``` ``` seed : integer, random_state, or None (default) ``` ``` Indicator of random number generation state. ``` ``` See :ref:`Randomness`. ``` ``` """ ``` ``` if not (0 <= p <= 1): ``` ``` raise ValueError("p must be in [0,1]") ``` ``` cs = ['d'] # threshold sequences always start with a d ``` ``` for i in range(1, n): ``` ``` if seed.random() < p: ``` ``` cs.append('d') ``` ``` else: ``` ``` cs.append('i') ``` ``` return cs ``` ```# maybe *_d_threshold_sequence routines should ``` ```# be (or be called from) a single routine with a more descriptive name ``` ```# and a keyword parameter? ``` ```def right_d_threshold_sequence(n, m): ``` ``` """ ``` ``` Create a skewed threshold graph with a given number ``` ``` of vertices (n) and a given number of edges (m). ``` ``` ``` ``` The routine returns an unlabeled creation sequence ``` ``` for the threshold graph. ``` ``` ``` ``` FIXME: describe algorithm ``` ``` ``` ``` """ ``` ``` cs = ['d'] + ['i'] * (n - 1) # create sequence with n insolated nodes ``` ``` # m n * (n - 1) / 2: ``` ``` raise ValueError("Too many edges for this many nodes.") ``` ``` # connected case m >n-1 ``` ``` ind = n - 1 ``` ``` sum = n - 1 ``` ``` while sum < m: ``` ``` cs[ind] = 'd' ``` ``` ind -= 1 ``` ``` sum += ind ``` ``` ind = m - (sum - ind) ``` ``` cs[ind] = 'd' ``` ``` return cs ``` ```def left_d_threshold_sequence(n, m): ``` ``` """ ``` ``` Create a skewed threshold graph with a given number ``` ``` of vertices (n) and a given number of edges (m). ``` ``` ``` ``` The routine returns an unlabeled creation sequence ``` ``` for the threshold graph. ``` ``` ``` ``` FIXME: describe algorithm ``` ``` ``` ``` """ ``` ``` cs = ['d'] + ['i'] * (n - 1) # create sequence with n insolated nodes ``` ``` # m n * (n - 1) / 2: ``` ``` raise ValueError("Too many edges for this many nodes.") ``` ``` # Connected case when M>N-1 ``` ``` cs[n - 1] = 'd' ``` ``` sum = n - 1 ``` ``` ind = 1 ``` ``` while sum < m: ``` ``` cs[ind] = 'd' ``` ``` sum += ind ``` ``` ind += 1 ``` ``` if sum > m: # be sure not to change the first vertex ``` ``` cs[sum - m] = 'i' ``` ``` return cs ``` ```@py_random_state(3) ``` ```def swap_d(cs, p_split=1.0, p_combine=1.0, seed=None): ``` ``` """ ``` ``` Perform a "swap" operation on a threshold sequence. ``` ``` ``` ``` The swap preserves the number of nodes and edges ``` ``` in the graph for the given sequence. ``` ``` The resulting sequence is still a threshold sequence. ``` ``` ``` ``` Perform one split and one combine operation on the ``` ``` 'd's of a creation sequence for a threshold graph. ``` ``` This operation maintains the number of nodes and edges ``` ``` in the graph, but shifts the edges from node to node ``` ``` maintaining the threshold quality of the graph. ``` ``` ``` ``` seed : integer, random_state, or None (default) ``` ``` Indicator of random number generation state. ``` ``` See :ref:`Randomness`. ``` ``` """ ``` ``` # preprocess the creation sequence ``` ``` dlist = [i for (i, node_type) in enumerate(cs[1:-1]) if node_type == 'd'] ``` ``` # split ``` ``` if seed.random() < p_split: ``` ``` choice = seed.choice(dlist) ``` ``` split_to = seed.choice(range(choice)) ``` ``` flip_side = choice - split_to ``` ``` if split_to != flip_side and cs[split_to] == 'i' and cs[flip_side] == 'i': ``` ``` cs[choice] = 'i' ``` ``` cs[split_to] = 'd' ``` ``` cs[flip_side] = 'd' ``` ``` dlist.remove(choice) ``` ``` # don't add or combine may reverse this action ``` ``` # dlist.extend([split_to,flip_side]) ``` ```# print >>sys.stderr,"split at %s to %s and %s"%(choice,split_to,flip_side) ``` ``` # combine ``` ``` if seed.random() < p_combine and dlist: ``` ``` first_choice = seed.choice(dlist) ``` ``` second_choice = seed.choice(dlist) ``` ``` target = first_choice + second_choice ``` ``` if target >= len(cs) or cs[target] == 'd' or first_choice == second_choice: ``` ``` return cs ``` ``` # OK to combine ``` ``` cs[first_choice] = 'i' ``` ``` cs[second_choice] = 'i' ``` ``` cs[target] = 'd' ``` ```# print >>sys.stderr,"combine %s and %s to make %s."%(first_choice,second_choice,target) ``` ``` return cs ```