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