Statistics
| Branch: | Revision:

iof-tools / networkxMiCe / networkx-master / networkx / algorithms / threshold.py @ 5cef0f13

History | View | Annotate | Download (29.6 KB)

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