Source Code for Module networkx.threshold

  1  """ 
  2  Threshold Graphs - Creation, manipulation and identification. 
  3  """ 
  4  __author__ = """Aric Hagberg (hagberg@lanl.gov)\nPieter Swart (swart@lanl.gov)\nDan Schult (dschult@colgate.edu)""" 
  5  __date__ = "$Date: 2005-06-17 08:06:22 -0600 (Fri, 17 Jun 2005) $" 
  6  __credits__ = """""" 
  7  __version__ = "$Revision: 1049 $" 
  8  # $Header$ 
  9   
 10  #    Copyright (C) 2004,2005 by  
 11  #    Aric Hagberg <hagberg@lanl.gov> 
 12  #    Dan Schult <dschult@colgate.edu> 
 13  #    Pieter Swart <swart@lanl.gov> 
 14  #    Distributed under the terms of the GNU Lesser General Public License 
 15  #    http://www.gnu.org/copyleft/lesser.html 
 16  # 
 17   
 18  import random  # for swap_d  
 19  from math import sqrt 
 20  import networkx 
 21   
 22 -def is_threshold_graph(G): 
 23      """ 
 24      Returns True if G is a threshold graph. 
 25      """ 
 26      return is_threshold_sequence(G.degree()) 
 27       
 28 -def is_threshold_sequence(degree_sequence): 
 29      """ 
 30      Returns True if the sequence is a threshold degree seqeunce. 
 31       
 32      Uses the property that a threshold graph must be constructed by 
 33      adding either dominating or isolated nodes. Thus, it can be 
 34      deconstructed iteratively by removing a node of degree zero or a 
 35      node that connects to the remaining nodes.  If this deconstruction 
 36      failes then the sequence is not a threshold sequence. 
 37      """ 
 38      ds=degree_sequence[:] # get a copy so we don't destroy original  
 39      ds.sort() 
 40      while ds:      
 41          if ds[0]==0:      # if isolated node 
 42              ds.pop(0)     # remove it  
 43              continue 
 44          if ds[-1]!=len(ds)-1:  # is the largest degree node dominating? 
 45              return False       # no, not a threshold degree sequence 
 46          ds.pop()               # yes, largest is the dominating node 
 47          ds=[ d-1 for d in ds ] # remove it and decrement all degrees 
 48      return True 
 49   
 50   
 51 -def creation_sequence(degree_sequence,with_labels=False,compact=False): 
 52      """ 
 53      Determines the creation sequence for the given threshold degree sequence. 
 54   
 55      The creation sequence is a list of single characters 'd' 
 56      or 'i': 'd' for dominating or 'i' for isolated vertices. 
 57      Dominating vertices are connected to all vertices present when it 
 58      is added.  The first node added is by convention 'd'. 
 59      This list can be converted to a string if desired using "".join(cs) 
 60   
 61      If with_labels==True: 
 62      Returns a list of 2-tuples containing the vertex number  
 63      and a character 'd' or 'i' which describes the type of vertex. 
 64   
 65      If compact==True: 
 66      Returns the creation sequence in a compact form that is the number 
 67      of 'i's and 'd's alternating. 
 68      Examples: 
 69      [1,2,2,3] represents d,i,i,d,d,i,i,i 
 70      [3,1,2] represents d,d,d,i,d,d 
 71   
 72      Notice that the first number is the first vertex to be used for 
 73      construction and so is always 'd'. 
 74   
 75      with_labels and compact cannot both be True. 
 76   
 77      Returns None if the sequence is not a threshold sequence 
 78      """ 
 79      if with_labels and compact: 
 80          raise ValueError,"compact sequences cannot be labeled" 
 81   
 82      # make an indexed copy 
 83      if isinstance(degree_sequence,dict):   # labeled degree seqeunce 
 84          ds = [ [degree,label] for (label,degree) in degree_sequence.items() ] 
 85      else: 
 86          ds=[ [d,i] for i,d in enumerate(degree_sequence) ]  
 87      ds.sort() 
 88      cs=[]  # creation sequence 
 89      while ds: 
 90          if ds[0][0]==0:     # isolated node 
 91              (d,v)=ds.pop(0) 
 92              if len(ds)>0:    # make sure we start with a d 
 93                  cs.insert(0,(v,'i')) 
 94              else: 
 95                  cs.insert(0,(v,'d')) 
 96              continue 
 97          if ds[-1][0]!=len(ds)-1:     # Not dominating node 
 98              return None  # not a threshold degree sequence 
 99          (d,v)=ds.pop() 
100          cs.insert(0,(v,'d')) 
101          ds=[ [d[0]-1,d[1]] for d in ds ]   # decrement due to removing node 
102   
103      if with_labels: return cs 
104      if compact: return make_compact(cs) 
105      return [ v[1] for v in cs ]   # not labeled 
106   
107 -def make_compact(creation_sequence): 
108      """ 
109      Returns the creation sequence in a compact form 
110      that is the number of 'i's and 'd's alternating. 
111      Examples: 
112      [1,2,2,3] represents d,i,i,d,d,i,i,i. 
113      [3,1,2] represents d,d,d,i,d,d. 
114      Notice that the first number is the first vertex 
115      to be used for construction and so is always 'd'. 
116   
117      Labeled creation sequences lose their labels in the  
118      compact representation. 
119      """ 
120      first=creation_sequence[0] 
121      if isinstance(first,str):    # creation sequence 
122          cs = creation_sequence[:] 
123      elif isinstance(first,tuple):   # labeled creation sequence 
124          cs = [ s[1] for s in creation_sequence ] 
125      elif isinstance(first,int):   # compact creation sequence 
126          return creation_sequence    
127      else: 
128          raise TypeError, "Not a valid creation sequence type" 
129   
130      ccs=[] 
131      count=1 # count the run lengths of d's or i's. 
132      for i in range(1,len(cs)): 
133          if cs[i]==cs[i-1]:  
134              count+=1 
135          else: 
136              ccs.append(count) 
137              count=1 
138      ccs.append(count) # don't forget the last one 
139      return ccs 
140       
141 -def uncompact(creation_sequence): 
142      """ 
143      Converts a compact creation sequence for a threshold 
144      graph to a standard creation sequence (unlabeled). 
145      If the creation_sequence is already standard, return it. 
146      See creation_sequence. 
147      """ 
148      first=creation_sequence[0] 
149      if isinstance(first,str):    # creation sequence 
150          return creation_sequence 
151      elif isinstance(first,tuple):   # labeled creation sequence 
152          return creation_sequence 
153      elif isinstance(first,int):   # compact creation sequence 
154          ccscopy=creation_sequence[:] 
155      else: 
156          raise TypeError, "Not a valid creation sequence type" 
157      cs = [] 
158      while ccscopy: 
159          cs.extend(ccscopy.pop(0)*['d']) 
160          if ccscopy: 
161              cs.extend(ccscopy.pop(0)*['i']) 
162      return cs 
163   
164 -def creation_sequence_to_weights(creation_sequence): 
165      """ 
166      Returns a list of node weights which create the threshold 
167      graph designated by the creation sequence.  The weights 
168      are scaled so that the threshold is 1.0.  The order of the 
169      nodes is the same as that in the creation sequence. 
170      """ 
171      # Turn input sequence into a labeled creation sequence 
172      first=creation_sequence[0] 
173      if isinstance(first,str):    # creation sequence 
174          if isinstance(creation_sequence,list): 
175              wseq = creation_sequence[:] 
176          else: 
177              wseq = list(creation_sequence) # string like 'ddidid' 
178      elif isinstance(first,tuple):   # labeled creation sequence 
179          wseq = [ v[1] for v in creation_sequence] 
180      elif isinstance(first,int):  # compact creation sequence 
181          wseq = uncompact(creation_sequence) 
182      else: 
183          raise TypeError, "Not a valid creation sequence type" 
184      # pass through twice--first backwards 
185      wseq.reverse() 
186      w=0 
187      prev='i' 
188      for j,s in enumerate(wseq): 
189          if s=='i': 
190              wseq[j]=w 
191              prev=s 
192          elif prev=='i': 
193              prev=s 
194              w+=1 
195      wseq.reverse()  # now pass through forwards 
196      for j,s in enumerate(wseq): 
197          if s=='d': 
198              wseq[j]=w 
199              prev=s 
200          elif prev=='d': 
201              prev=s 
202              w+=1 
203      # Now scale weights 
204      if prev=='d': w+=1 
205      wscale=1./float(w) 
206      return [ ww*wscale for ww in wseq] 
207      #return wseq 
208   
209 -def weights_to_creation_sequence(weights,threshold=1,with_labels=False,compact=False): 
210      """ 
211      Returns a creation sequence for a threshold graph  
212      determined by the weights and threshold given as input. 
213      If the sum of two node weights is greater than the 
214      threshold value, an edge is created between these nodes. 
215   
216      The creation sequence is a list of single characters 'd' 
217      or 'i': 'd' for dominating or 'i' for isolated vertices. 
218      Dominating vertices are connected to all vertices present  
219      when it is added.  The first node added is by convention 'd'. 
220   
221      If with_labels==True: 
222      Returns a list of 2-tuples containing the vertex number  
223      and a character 'd' or 'i' which describes the type of vertex. 
224   
225      If compact==True: 
226      Returns the creation sequence in a compact form that is the number 
227      of 'i's and 'd's alternating. 
228      Examples: 
229      [1,2,2,3] represents d,i,i,d,d,i,i,i 
230      [3,1,2] represents d,d,d,i,d,d 
231   
232      Notice that the first number is the first vertex to be used for 
233      construction and so is always 'd'. 
234   
235      with_labels and compact cannot both be True. 
236      """ 
237      if with_labels and compact: 
238          raise ValueError,"compact sequences cannot be labeled" 
239   
240      # make an indexed copy 
241      if isinstance(weights,dict):   # labeled weights 
242          wseq = [ [w,label] for (label,w) in weights.items() ] 
243      else: 
244          wseq = [ [w,i] for i,w in enumerate(weights) ]  
245      wseq.sort() 
246      cs=[]  # creation sequence 
247      cutoff=threshold-wseq[-1][0] 
248      while wseq: 
249          if wseq[0][0]<cutoff:     # isolated node 
250              (w,label)=wseq.pop(0) 
251              cs.append((label,'i')) 
252          else: 
253              (w,label)=wseq.pop() 
254              cs.append((label,'d')) 
255              cutoff=threshold-wseq[-1][0] 
256          if len(wseq)==1:     # make sure we start with a d 
257              (w,label)=wseq.pop() 
258              cs.append((label,'d')) 
259      # put in correct order 
260      cs.reverse() 
261   
262      if with_labels: return cs 
263      if compact: return make_compact(cs) 
264      return [ v[1] for v in cs ]   # not labeled 
265   
266   
267  # Manipulating NetworkX.Graphs in context of threshold graphs 
268 -def threshold_graph(creation_sequence): 
269      """ 
270      Create a threshold graph from the creation sequence or compact 
271      creation_sequence. 
272   
273      The input sequence can be a  
274   
275      creation sequence (e.g. ['d','i','d','d','d','i']) 
276      labeled creation sequence (e.g. [(0,'d'),(2,'d'),(1,'i')]) 
277      compact creation sequence (e.g. [2,1,1,2,0]) 
278       
279      Use cs=creation_sequence(degree_sequence,labeled=True)  
280      to convert a degree sequence to a creation sequence. 
281   
282      Returns None if the sequence is not valid 
283      """ 
284      # Turn input sequence into a labeled creation sequence 
285      first=creation_sequence[0] 
286      if isinstance(first,str):    # creation sequence 
287          ci = list(enumerate(creation_sequence)) 
288      elif isinstance(first,tuple):   # labeled creation sequence 
289          ci = creation_sequence[:] 
290      elif isinstance(first,int):  # compact creation sequence 
291          cs = uncompact(creation_sequence) 
292          ci = list(enumerate(cs)) 
293      else: 
294          print "not a valid creation sequence type" 
295          return None 
296           
297      G=networkx.Graph() 
298      G.name="Threshold Graph" 
299   
300      # add nodes and edges 
301      # if type is 'i' just add nodea 
302      # if type is a d connect to everything previous 
303      while ci: 
304          (v,node_type)=ci.pop(0) 
305          G.add_node(v) 
306          if node_type=='d': # dominating type, connect to all existing nodes 
307              for u in G.nodes():  
308                  G.add_edge(v,u) # will silently ignore self loop 
309      return G 
310   
311   
312   
313 -def find_alternating_4_cycle(G): 
314      """ 
315      Returns False if there aren't any alternating 4 cycles. 
316      Otherwise returns the cycle as [a,b,c,d] where (a,b) 
317      and (c,d) are edges and (a,c) and (b,d) are not. 
318      """ 
319      for (u,v) in G.edges(): 
320          for w in G.nodes(): 
321              if not G.has_edge(u,w) and u!=w: 
322                  for x in G.neighbors(w): 
323                      if not G.has_edge(v,x) and v!=x: 
324                          return [u,v,w,x] 
325      return False 
326   
327   
328   
329 -def find_threshold_graph(G): 
330      """ 
331      Return a threshold subgraph that is close to largest in G. 
332      The threshold graph will contain the largest degree node in G. 
333       
334      """ 
335      return threshold_graph(find_creation_sequence(G)) 
336   
337   
338 -def find_creation_sequence(G): 
339      """ 
340      Find a threshold subgraph that is close to largest in G. 
341      Returns the labeled creation sequence of that threshold graph.   
342      """ 
343      cs=[] 
344      # get a local pointer to the working part of the graph 
345      H=G 
346      while H.order()>0: 
347          # get new degree sequence on subgraph 
348          dsdict=H.degree(with_labels=True) 
349          ds=[ [d,v] for v,d in dsdict.iteritems() ] 
350          ds.sort() 
351          # Update threshold graph nodes 
352          if ds[-1][0]==0: # all are isolated 
353              cs.extend( zip( dsdict, ['i']*(len(ds)-1)+['d']) ) 
354              break   # Done! 
355          # pull off isolated nodes 
356          while ds[0][0]==0: 
357              (d,iso)=ds.pop(0) 
358              cs.append((iso,'i')) 
359          # find new biggest node 
360          (d,bigv)=ds.pop() 
361          # add edges of star to t_g 
362          cs.append((bigv,'d')) 
363          # form subgraph of neighbors of big node 
364          H=H.subgraph(H.neighbors(bigv)) 
365      cs.reverse() 
366      return cs 
367       
368   
369       
370  ### Properties of Threshold Graphs 
371 -def triangles(creation_sequence): 
372      """ 
373      Compute number of triangles in the threshold graph with the 
374      given creation sequence. 
375      """ 
376      # shortcut algoritm that doesn't require computing number 
377      # of triangles at each node. 
378      cs=creation_sequence    # alias 
379      dr=cs.count("d")        # number of d's in sequence 
380      ntri=dr*(dr-1)*(dr-2)/6 # number of triangles in clique of nd d's 
381      # now add dr choose 2 triangles for every 'i' in sequence where 
382      # dr is the number of d's to the right of the current i 
383      for i,typ in enumerate(cs):  
384          if typ=="i": 
385              ntri+=dr*(dr-1)/2 
386          else: 
387              dr-=1 
388      return ntri 
389   
390   
391 -def triangle_sequence(creation_sequence): 
392      """ 
393      Return triangle sequence for the given threshold graph creation sequence. 
394   
395      """ 
396      cs=creation_sequence 
397      seq=[] 
398      dr=cs.count("d")     # number of d's to the right of the current pos 
399      dcur=(dr-1)*(dr-2)/2 # number of triangles through a node of clique dr 
400      irun=0               # number of i's in the last run 
401      drun=0               # number of d's in the last run 
402      for i,sym in enumerate(cs): 
403          if sym=="d": 
404              drun+=1 
405              tri=dcur+(dr-1)*irun    # new triangles at this d 
406          else: # cs[i]="i": 
407              if prevsym=="d":        # new string of i's 
408                  dcur+=(dr-1)*irun   # accumulate shared shortest paths 
409                  irun=0              # reset i run counter 
410                  dr-=drun            # reduce number of d's to right 
411                  drun=0              # reset d run counter 
412              irun+=1 
413              tri=dr*(dr-1)/2      # new triangles at this i 
414          seq.append(tri) 
415          prevsym=sym 
416      return seq 
417   
418 -def cluster_sequence(creation_sequence): 
419      """ 
420      Return cluster sequence for the given threshold graph creation sequence. 
421      """ 
422      triseq=triangle_sequence(creation_sequence) 
423      degseq=degree_sequence(creation_sequence) 
424      cseq=[] 
425      for i,deg in enumerate(degseq): 
426          tri=triseq[i] 
427          if deg <= 1:    # isolated vertex or single pair gets cc 0 
428              cseq.append(0) 
429              continue 
430          max_size=(deg*(deg-1))/2 
431          cseq.append(float(tri)/float(max_size)) 
432      return cseq 
433   
434   
435 -def degree_sequence(creation_sequence): 
436      """ 
437      Return degree sequence for the threshold graph with the given 
438      creation sequence 
439      """ 
440      cs=creation_sequence  # alias 
441      seq=[] 
442      rd=cs.count("d") # number of d to the right 
443      for i,sym in enumerate(cs): 
444          if sym=="d": 
445              rd-=1 
446              seq.append(rd+i) 
447          else: 
448              seq.append(rd) 
449      return seq             
450   
451 -def density(creation_sequence): 
452      """ 
453      Return the density of the graph with this creation_sequence. 
454      The density is the fraction of possible edges present. 
455      """ 
456      N=len(creation_sequence) 
457      two_size=sum(degree_sequence(creation_sequence)) 
458      two_possible=N*(N-1) 
459      den=two_size/float(two_possible) 
460      return den 
461   
462 -def degree_correlation(creation_sequence): 
463      """ 
464      Return the degree-degree correlation over all edges. 
465      """ 
466      cs=creation_sequence 
467      s1=0  # deg_i*deg_j 
468      s2=0  # deg_i^2+deg_j^2 
469      s3=0  # deg_i+deg_j 
470      m=0   # number of edges 
471      rd=cs.count("d") # number of d nodes to the right 
472      rdi=[ i for i,sym in enumerate(cs) if sym=="d"] # index of "d"s 
473      ds=degree_sequence(cs) 
474      for i,sym in enumerate(cs): 
475          if sym=="d":  
476              if i!=rdi[0]: 
477                  print "Logic error in degree_correlation",i,rdi 
478                  raise ValueError 
479              rdi.pop(0) 
480          degi=ds[i] 
481          for dj in rdi: 
482              degj=ds[dj] 
483              s1+=degj*degi 
484              s2+=degi**2+degj**2 
485              s3+=degi+degj 
486              m+=1 
487      denom=(2*m*s2-s3*s3) 
488      numer=(4*m*s1-s3*s3) 
489      if denom==0: 
490          if numer==0:  
491              return 1 
492          raise ValueError,"Zero Denominator but Numerator is %s"%numer 
493      return numer/float(denom) 
494   
495   
496 -def shortest_path(creation_sequence,u,v): 
497      """ 
498      Find the shortest path between u and v in a  
499      threshold graph G with the given creation_sequence. 
500   
501      For an unlabeled creation_sequence, the vertices  
502      u and v must be integers in (0,len(sequence)) refering  
503      to the position of the desired vertices in the sequence. 
504   
505      For a labeled creation_sequence, u and v are labels of veritices. 
506   
507      Use cs=creation_sequence(degree_sequence,with_labels=True)  
508      to convert a degree sequence to a creation sequence. 
509   
510      Returns a list of vertices from u to v. 
511      Example: if they are neighbors, it returns [u,v] 
512      """ 
513      # Turn input sequence into a labeled creation sequence 
514      first=creation_sequence[0] 
515      if isinstance(first,str):    # creation sequence 
516          cs = [(i,creation_sequence[i]) for i in xrange(len(creation_sequence))] 
517      elif isinstance(first,tuple):   # labeled creation sequence 
518          cs = creation_sequence[:] 
519      elif isinstance(first,int):  # compact creation sequence 
520          ci = uncompact(creation_sequence) 
521          cs = [(i,ci[i]) for i in xrange(len(ci))] 
522      else: 
523          raise TypeError, "Not a valid creation sequence type" 
524           
525      verts=[ s[0] for s in cs ] 
526      if v not in verts: 
527          raise ValueError,"Vertex %s not in graph from creation_sequence"%v 
528      if u not in verts: 
529          raise ValueError,"Vertex %s not in graph from creation_sequence"%u 
530      # Done checking 
531      if u==v: return [u] 
532   
533      uindex=verts.index(u) 
534      vindex=verts.index(v) 
535      bigind=max(uindex,vindex) 
536      if cs[bigind][1]=='d': 
537          return [u,v] 
538      # must be that cs[bigind][1]=='i' 
539      cs=cs[bigind:] 
540      while cs: 
541          vert=cs.pop() 
542          if vert[1]=='d': 
543              return [u,vert[0],v] 
544      # All after u are type 'i' so no connection 
545      return -1 
546   
547 -def shortest_path_length(creation_sequence,i): 
548      """ 
549      Return the shortest path length from indicated node to 
550      every other node for the threshold graph with the given 
551      creation sequence. 
552      Node is indicated by index i in creation_sequence unless 
553      creation_sequence is labeled in which case, i is taken to 
554      be the label of the node. 
555   
556      Paths lengths in threshold graphs are at most 2. 
557      Length to unreachable nodes is set to -1. 
558      """ 
559      # Turn input sequence into a labeled creation sequence 
560      first=creation_sequence[0] 
561      if isinstance(first,str):    # creation sequence 
562          if isinstance(creation_sequence,list): 
563              cs = creation_sequence[:] 
564          else: 
565              cs = list(creation_sequence) 
566      elif isinstance(first,tuple):   # labeled creation sequence 
567          cs = [ v[1] for v in creation_sequence] 
568          i = [v[0] for v in creation_sequence].index(i) 
569      elif isinstance(first,int):  # compact creation sequence 
570          cs = uncompact(creation_sequence) 
571      else: 
572          raise TypeError, "Not a valid creation sequence type" 
573       
574      # Compute 
575      N=len(cs) 
576      spl=[2]*N       # length 2 to every node 
577      spl[i]=0        # except self which is 0 
578      # 1 for all d's to the right 
579      for j in range(i+1,N):    
580          if cs[j]=="d": 
581              spl[j]=1 
582      if cs[i]=='d': # 1 for all nodes to the left 
583          for j in range(i): 
584              spl[j]=1 
585      # and -1 for any trailing i to indicate unreachable 
586      for j in range(N-1,0,-1): 
587          if cs[j]=="d": 
588              break 
589          spl[j]=-1 
590      return spl 
591   
592   
593 -def betweenness_sequence(creation_sequence,normalized=True): 
594      """ 
595      Return betweenness for the threshold graph with the given creation 
596      sequence.  The result is unscaled.  To scale the values 
597      to the iterval [0,1] divide by (n-1)*(n-2). 
598      """ 
599      cs=creation_sequence 
600      seq=[]               # betweenness 
601      lastchar='d'         # first node is always a 'd'         
602      dr=float(cs.count("d"))  # number of d's to the right of curren pos 
603      irun=0               # number of i's in the last run 
604      drun=0               # number of d's in the last run 
605      dlast=0.0              # betweenness of last d 
606      for i,c in enumerate(cs): 
607          if c=='d':    #cs[i]=="d": 
608              # betweennees = amt shared with eariler d's and i's 
609              #             + new isolated nodes covered 
610              #             + new paths to all previous nodes 
611              b=dlast + (irun-1)*irun/dr + 2*irun*(i-drun-irun)/dr 
612              drun+=1           # update counter 
613          else:      # cs[i]="i": 
614              if lastchar=='d': # if this is a new run of i's 
615                  dlast=b       # accumulate betweenness 
616                  dr-=drun      # update number of d's to the right 
617                  drun=0        # reset d counter 
618                  irun=0        # reset i counter 
619              b=0      # isolated nodes have zero betweenness  
620              irun+=1  # add another i to the run    
621          seq.append(float(b)) 
622          lastchar=c 
623   
624      # normalize by the number of possible shortest paths 
625      if normalized: 
626          order=len(cs) 
627          scale=1.0/((order-1)*(order-2)) 
628          seq=[ s*scale for s in seq ] 
629   
630      return seq 
631   
632   
633 -def eigenvectors(creation_sequence): 
634      """ 
635      Return a 2-tuple of Laplacian eigenvalues and eigenvectors  
636      for the threshold network with creation_sequence. 
637      The first value is a list of eigenvalues. 
638      The second value is a list of eigenvectors. 
639      The lists are in the same order so corresponding eigenvectors 
640      and eigenvalues are in the same position in the two lists. 
641   
642      Notice that the order of the eigenvalues returned by eigenvalues(cs) 
643      may not correspond to the order of these eigenvectors. 
644      """ 
645      ccs=make_compact(creation_sequence) 
646      N=sum(ccs) 
647      vec=[0]*N 
648      val=vec[:] 
649      # get number of type d nodes to the right (all for first node) 
650      dr=sum(ccs[::2]) 
651   
652      nn=ccs[0] 
653      vec[0]=[1./sqrt(N)]*N 
654      val[0]=0 
655      e=dr 
656      dr-=nn 
657      type_d=True 
658      i=1 
659      dd=1 
660      while dd<nn: 
661          scale=1./sqrt(dd*dd+i) 
662          vec[i]=i*[-scale]+[dd*scale]+[0]*(N-i-1) 
663          val[i]=e 
664          i+=1 
665          dd+=1 
666      if len(ccs)==1: return (val,vec) 
667      for nn in ccs[1:]: 
668          scale=1./sqrt(nn*i*(i+nn)) 
669          vec[i]=i*[-nn*scale]+nn*[i*scale]+[0]*(N-i-nn) 
670          # find eigenvalue 
671          type_d=not type_d 
672          if type_d: 
673              e=i+dr 
674              dr-=nn 
675          else: 
676              e=dr 
677          val[i]=e 
678          st=i 
679          i+=1 
680          dd=1 
681          while dd<nn: 
682              scale=1./sqrt(i-st+dd*dd) 
683              vec[i]=[0]*st+(i-st)*[-scale]+[dd*scale]+[0]*(N-i-1) 
684              val[i]=e 
685              i+=1 
686              dd+=1 
687      return (val,vec) 
688   
689 -def spectral_projection(u,eigenpairs): 
690      """ 
691      Returns the coefficients of each eigenvector 
692      in a projection of the vector u onto the normalized 
693      eigenvectors which are contained in eigenpairs. 
694   
695      eigenpairs should be a list of two objects.  The 
696      first is a list of eigenvalues and the second a list 
697      of eigenvectors.  The eigenvectors should be lists.  
698   
699      There's not a lot of error checking on lengths of  
700      arrays, etc. so be careful. 
701      """ 
702      coeff=[] 
703      evect=eigenpairs[1] 
704      for ev in evect: 
705          c=sum([ evv*uv for (evv,uv) in zip(ev,u)]) 
706          coeff.append(c) 
707      return coeff 
708           
709   
710   
711 -def eigenvalues(creation_sequence): 
712      """ 
713      Return sequence of eigenvalues of the Laplacian of the threshold 
714      graph for the given creation_sequence. 
715   
716      Based on the Ferrer's diagram method.  The spectrum is integral 
717      and is the conjugate of the degree sequence. 
718   
719      See::  
720   
721        @Article{degree-merris-1994, 
722         author =          {Russel Merris}, 
723         title =   {Degree maximal graphs are Laplacian integral}, 
724         journal =         {Linear Algebra Appl.}, 
725         year =    {1994}, 
726         volume =          {199}, 
727         pages =   {381--389}, 
728        } 
729   
730      """ 
731      degseq=degree_sequence(creation_sequence) 
732      degseq.sort() 
733      eiglist=[]   # zero is always one eigenvalue 
734      eig=0 
735      row=len(degseq) 
736      bigdeg=degseq.pop() 
737      while row: 
738          if bigdeg<row: 
739              eiglist.append(eig) 
740              row-=1 
741          else: 
742              eig+=1 
743              if degseq: 
744                  bigdeg=degseq.pop() 
745              else: 
746                  bigdeg=0 
747      return eiglist 
748   
749   
750  ### Threshold graph creation routines 
751   
752 -def random_threshold_sequence(n,p,seed=None): 
753      """ 
754      Create a random threshold sequence of size n. 
755      A creation sequence is built by randomly choosing d's with 
756      probabiliy p and i's with probability 1-p. 
757   
758      >>> s=random_threshold_sequence(10,0.5) 
759   
760      returns a threshold sequence of length 10 with equal 
761      probably of an i or a d at each position. 
762   
763      A "random" threshold graph can be built with 
764   
765      >>> G=threshold_graph(random_threshold_sequence(10,0.5)) 
766   
767      """ 
768      if not seed is None: 
769          random.seed(seed) 
770   
771      if not (p<=1 and p>=0): 
772          raise ValueError,"p must be in [0,1]" 
773   
774      cs=['d']  # threshold sequences always start with a d 
775      for i in range(1,n): 
776          if random.random() < p:  
777              cs.append('d') 
778          else: 
779              cs.append('i') 
780      return cs 
781   
782   
783   
784   
785   
786  # maybe *_d_threshold_sequence routines should 
787  # be (or be called from) a single routine with a more descriptive name 
788  # and a keyword parameter? 
789 -def right_d_threshold_sequence(n,m): 
790      """ 
791      Create a skewed threshold graph with a given number 
792      of vertices (n) and a given number of edges (m). 
793       
794      The routine returns an unlabeled creation sequence  
795      for the threshold graph. 
796   
797      FIXME: describe algorithm 
798   
799      """ 
800      cs=['d']+['i']*(n-1) # create sequence with n insolated nodes 
801   
802      #  m <n : not enough edges, make disconnected  
803      if m < n:  
804          cs[m]='d' 
805          return cs 
806   
807      # too many edges  
808      if m > n*(n-1)/2: 
809          raise ValueError,"Too many edges for this many nodes." 
810   
811      # connected case m >n-1 
812      ind=n-1 
813      sum=n-1 
814      while sum<m: 
815          cs[ind]='d' 
816          ind -= 1 
817          sum += ind 
818      ind=m-(sum-ind) 
819      cs[ind]='d' 
820      return cs 
821   
822 -def left_d_threshold_sequence(n,m): 
823      """ 
824      Create a skewed threshold graph with a given number 
825      of vertices (n) and a given number of edges (m). 
826   
827      The routine returns an unlabeled creation sequence  
828      for the threshold graph. 
829   
830      FIXME: describe algorithm 
831   
832      """ 
833      cs=['d']+['i']*(n-1) # create sequence with n insolated nodes 
834   
835      #  m <n : not enough edges, make disconnected  
836      if m < n:  
837          cs[m]='d' 
838          return cs 
839   
840      # too many edges  
841      if m > n*(n-1)/2: 
842          raise ValueError,"Too many edges for this many nodes." 
843   
844      # Connected case when M>N-1 
845      cs[n-1]='d' 
846      sum=n-1 
847      ind=1 
848      while sum<m: 
849          cs[ind]='d' 
850          sum += ind 
851          ind += 1 
852      if sum>m:    # be sure not to change the first vertex 
853          cs[sum-m]='i' 
854      return cs 
855   
856 -def swap_d(cs,p_split=1.0,p_combine=1.0,seed=None): 
857      """ 
858      Perform a "swap" operation on a threshold sequence. 
859   
860      The swap preserves the number of nodes and edges 
861      in the graph for the given sequence. 
862      The resulting sequence is still a threshold sequence. 
863   
864      Perform one split and one combine operation on the 
865      'd's of a creation sequence for a threshold graph. 
866      This operation maintains the number of nodes and edges 
867      in the graph, but shifts the edges from node to node 
868      maintaining the threshold quality of the graph. 
869      """ 
870      if not seed is None: 
871          random.seed(seed) 
872   
873      # preprocess the creation sequence 
874      dlist= [ i for (i,node_type) in enumerate(cs[1:-1]) if node_type=='d' ] 
875      # split 
876      if random.random()<p_split: 
877          choice=random.choice(dlist) 
878          split_to=random.choice(range(choice)) 
879          flip_side=choice-split_to 
880          if split_to!=flip_side and cs[split_to]=='i' and cs[flip_side]=='i': 
881              cs[choice]='i' 
882              cs[split_to]='d' 
883              cs[flip_side]='d' 
884              dlist.remove(choice) 
885              # don't add or combine may reverse this action             
886              # dlist.extend([split_to,flip_side]) 
887  #            print >>sys.stderr,"split at %s to %s and %s"%(choice,split_to,flip_side) 
888      # combine 
889      if random.random()<p_combine and dlist: 
890          first_choice= random.choice(dlist) 
891          second_choice=random.choice(dlist) 
892          target=first_choice+second_choice 
893          if target >= len(cs) or cs[target]=='d' or first_choice==second_choice: 
894              return cs 
895          # OK to combine 
896          cs[first_choice]='i' 
897          cs[second_choice]='i' 
898          cs[target]='d' 
899  #        print >>sys.stderr,"combine %s and %s to make %s."%(first_choice,second_choice,target) 
900   
901      return cs 
902           
903   
904   
905           
906 -def _test_suite(): 
907      import doctest 
908      suite = doctest.DocFileSuite('tests/threshold.txt', package='networkx') 
909      return suite 
910   
911   
912  if __name__ == "__main__": 
913      import os 
914      import sys 
915      import unittest 
916      if sys.version_info[:2] < (2, 4): 
917          print "Python version 2.4 or later required for tests (%d.%d detected)." %  sys.version_info[:2] 
918          sys.exit(-1) 
919      # directory of networkx package (relative to this) 
920      nxbase=sys.path[0]+os.sep+os.pardir 
921      sys.path.insert(0,nxbase) # prepend to search path 
922      unittest.TextTestRunner().run(_test_suite()) 
923