Mercurial > repos > shellac > sam_consensus_v3
comparison env/lib/python3.9/site-packages/networkx/algorithms/threshold.py @ 0:4f3585e2f14b draft default tip
"planemo upload commit 60cee0fc7c0cda8592644e1aad72851dec82c959"
author | shellac |
---|---|
date | Mon, 22 Mar 2021 18:12:50 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:4f3585e2f14b |
---|---|
1 """ | |
2 Threshold Graphs - Creation, manipulation and identification. | |
3 """ | |
4 from math import sqrt | |
5 import networkx as nx | |
6 from networkx.utils import py_random_state | |
7 | |
8 __all__ = ["is_threshold_graph", "find_threshold_graph"] | |
9 | |
10 | |
11 def 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 | |
42 def 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 | |
65 def 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 | |
124 def 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 | |
168 def 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 | |
192 def 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 | |
239 def 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 | |
302 def 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 | |
353 def 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 | |
368 def 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 | |
402 def 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 | |
434 def 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 | |
454 def 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 | |
482 def 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 | |
499 def 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 | |
516 def 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 | |
528 def 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 | |
562 def 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 | |
615 def 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 | |
661 def 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 | |
701 def 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 | |
759 def 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 | |
780 def 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) | |
823 def 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? | |
857 def 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 | |
891 def 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) | |
927 def 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 |