Mercurial > repos > guerler > springsuite

comparison planemo/lib/python3.7/site-packages/networkx/algorithms/threshold.py @ 1:56ad4e20f292 draft

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