Package networkx :: Module threshold
[hide private]
[frames] | no frames]

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