Ticket #11584: trac_11584.patch

File trac_11584.patch, 20.4 KB (added by ncohen, 9 years ago)
  • doc/en/reference/combinat/index.rst

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1310313042 -7200
    # Node ID 2a88f8db92e6f0ad6a59ced381e36f4a02bf770a
    # Parent  96a9712b779c3581bf63e90df30d1063f354b638
    Trac #11584 -- DegreeSequences class
    
    diff --git a/doc/en/reference/combinat/index.rst b/doc/en/reference/combinat/index.rst
    a b  
    1414   ../sage/combinat/composition
    1515   ../sage/combinat/core
    1616   ../sage/combinat/debruijn_sequence
     17   ../sage/combinat/degree_sequences
    1718   ../sage/combinat/dlx
    1819   ../sage/combinat/matrices/dlxcpp
    1920   ../sage/combinat/dyck_word
  • module_list.py

    diff --git a/module_list.py b/module_list.py
    a b  
    186186
    187187    Extension('sage.structure.list_clone',
    188188              sources=['sage/structure/list_clone.pyx']),
     189
    189190    Extension('sage.structure.list_clone_timings_cy',
    190191              sources=['sage/structure/list_clone_timings_cy.pyx']),
    191192
    192 
    193193    Extension('sage.sets.finite_set_map_cy',
    194194              sources=['sage/sets/finite_set_map_cy.pyx']),
    195195
     
    215215    Extension('sage.combinat.debruijn_sequence',
    216216              sources=['sage/combinat/debruijn_sequence.pyx']),
    217217
     218    Extension('sage.combinat.degree_sequences',
     219              sources = ['sage/combinat/degree_sequences.pyx']),
     220
    218221    Extension('sage.combinat.combinat_cython',
    219222              sources=['sage/combinat/combinat_cython.pyx'],
    220223              libraries=['gmp']),
  • sage/combinat/all.py

    diff --git a/sage/combinat/all.py b/sage/combinat/all.py
    a b  
    112112
    113113from kazhdan_lusztig import KazhdanLusztigPolynomial
    114114
     115from degree_sequences import DegreeSequences
     116
    115117from cyclic_sieving_phenomenon import CyclicSievingPolynomial, CyclicSievingCheck
    116118
    117119from sidon_sets import sidon_sets
  • new file sage/combinat/degree_sequences.pyx

    diff --git a/sage/combinat/degree_sequences.pyx b/sage/combinat/degree_sequences.pyx
    new file mode 100644
    - +  
     1r"""
     2Degree sequences
     3
     4The present module implements the ``DegreeSequences`` class, whose instances
     5represent the integer sequences of length `n`::
     6
     7    sage: DegreeSequences(6)
     8    Degree sequences on 6 elements
     9
     10With the object ``DegreeSequences(n)``, one can :
     11
     12    * Check whether a sequence is indeed a degree sequence ::
     13
     14        sage: DS = DegreeSequences(5)
     15        sage: [4, 3, 3, 3, 3] in DS
     16        True
     17        sage: [4, 4, 0, 0, 0] in DS
     18        False
     19
     20    * List all the possible degree sequences of length `n`::
     21
     22        sage: for seq in DegreeSequences(4):
     23        ...       print seq
     24        [0, 0, 0, 0]
     25        [1, 1, 0, 0]
     26        [2, 1, 1, 0]
     27        [3, 1, 1, 1]
     28        [1, 1, 1, 1]
     29        [2, 2, 1, 1]
     30        [2, 2, 2, 0]
     31        [3, 2, 2, 1]
     32        [2, 2, 2, 2]
     33        [3, 3, 2, 2]
     34        [3, 3, 3, 3]
     35
     36.. NOTE::
     37
     38    Given a degree sequence, one can obtain a graph realizing it by using
     39    :meth:`sage.graphs.graph_generators.graphs.DegreeSequence`. For instance ::
     40
     41        sage: ds = [3, 3, 2, 2, 2, 2, 2, 1, 1, 0]
     42        sage: g = graphs.DegreeSequence(ds)
     43        sage: g.degree_sequence()
     44        [3, 3, 2, 2, 2, 2, 2, 1, 1, 0]
     45
     46Definitions
     47~~~~~~~~~~~
     48
     49A sequence of integers `d_1,...,d_n` is said to be a *degree sequence* (or
     50*graphic* sequence) if there exists a graph in which vertex `i` is of degree
     51`d_i`. It is often required to be *non-increasing*, i.e. that `d_1 \geq ... \geq
     52d_n`.
     53
     54An integer sequence need not necessarily be a degree sequence. Indeed, in a
     55degree sequence of length `n` no integer can be larger than `n-1` -- the degree
     56of a vertex is at most `n-1` -- and the sum of them is at most `n(n-1)`.
     57
     58Degree sequences are completely characterized by a result from Erdos and Gallai:
     59
     60**Erdos and Gallai:** *The sequence of integers* `d_1\geq ... \geq d_n` *is a
     61degree sequence if and only if* `\forall i`
     62
     63.. MATH::
     64    \sum_{j\leq i}d_j \leq j(j-1) + \sum_{j>i}min(d_j,i)
     65
     66Alternatively, a degree sequence can be defined recursively :
     67
     68**Havel and Hakimi:** *The sequence of integers* `d_1\geq ... \geq d_n` *is a
     69degree sequence if and only if* `d_2-1,...,d_{d_1+1}-1, d_{d_1+2}, ...,d_n` *is
     70also a degree sequence.*
     71
     72Or equivalently :
     73
     74**Havel and Hakimi (bis):** *If there is a realization of an integer sequence as
     75a graph (i.e. if the sequence is a degree sequence), then it can be realized in
     76such a way that the vertex of maximum degree* `\Delta` *is adjacent to the*
     77`\Delta` *vertices of highest degree (except itself, of course).*
     78
     79
     80Algorithms
     81~~~~~~~~~~
     82
     83**Checking whether a given sequence is a degree sequence**
     84
     85This is tested using Erdos and Gallai's criterion. It is also checked that the
     86given sequence is non-increasing and has length `n`.
     87
     88**Iterating through the sequences of length** `n`
     89
     90From Havel and Hakimi's recursive definition of a degree sequence, one can build
     91an enumeration algorithm as done in [RCES]_. It consists in trying to **extend**
     92a current degree sequence on `n` elements into a degree sequence on `n+1`
     93elements by adding a vertex of degree larger than those already present in the
     94sequence. This can be seen as **reversing** the reduction operation described in
     95Havel and Hakimi's characterization. This operation can appear in several
     96different ways :
     97
     98    * Extensions of a degree sequence that do **not** change the value of the
     99      maximum element
     100
     101        * If the maximum element of a given degree sequence is `0`, then one can
     102          remove it to reduce the sequence, following Havel and Hakimi's
     103          rule. Conversely, if the maximum element of the (current) sequence is
     104          `0`, then one can always extend it by adding a new element of degree
     105          `0` to the sequence.
     106
     107          .. MATH::
     108              0, 0, 0 \xrightarrow{Extension} {\bf 0}, 0, 0, 0 \xrightarrow{Extension} {\bf 0}, 0, 0, ..., 0, 0, 0 \xrightarrow{Reduction} 0, 0, 0, 0 \xrightarrow{Reduction} 0, 0, 0
     109
     110        * If there are at least `\Delta+1` elements of (maximum) degree `\Delta`
     111          in a given degree sequence, then one can reduce it by removing a
     112          vertex of degree `\Delta` and decreasing the values of `\Delta`
     113          elements of value `\Delta` to `\Delta-1`. Conversely, if the maximum
     114          element of the (current) sequence is `d>0`, then one can add a new
     115          element of degree `d` to the sequence if it can be linked to `d`
     116          elements of (current) degree `d-1`. Those `d` vertices of degree `d-1`
     117          hence become vertices of degree `d`, and so `d` elements of degree
     118          `d-1` are removed from the sequence while `d+1` elements of degree `d`
     119          are added to it.
     120
     121          .. MATH::
     122              3, 2, 2, 2, 1 \xrightarrow{Extension} {\bf 3}, 3, (2+1), (2+1), (2+1), 1 =  {\bf 3}, 3, 3, 3, 3, 1 \xrightarrow{Reduction} 3, 2, 2, 2, 1
     123
     124    * Extension of a degree sequence that changes the value of the maximum
     125      element :
     126
     127        * In the general case, i.e. when the number of elements of value
     128          `\Delta,\Delta-1` is small compared to `\Delta` (i.e. the maximum
     129          element of a given degree sequence), reducing a sequence strictly
     130          decreases the value of the maximum element. According to Havel and
     131          Hakimi's characterization there is only **one** way to reduce a
     132          sequence, but reversing this operation is more complicated than in the
     133          previous cases. Indeed, the following extensions are perfectly valid
     134          according to the reduction rule.
     135
     136          .. MATH::
     137              2,1,1,0,0\xrightarrow{Extension} {\bf 3}, (2+1), (1+1), (1+1), 0, 0 = 3, 3, 2, 2, 0, 0 \xrightarrow{Reduction} 2, 1, 1, 0, 0\\
     138              2,1,1,0,0\xrightarrow{Extension} {\bf 3}, (2+1), (1+1), 1, (0+1), 0 = 3, 3, 2, 1, 1, 0 \xrightarrow{Reduction} 2, 1, 1, 0, 0\\
     139              2,1,1,0,0\xrightarrow{Extension} {\bf 3}, (2+1), 1, 1, (0+1), (0+1) = 3, 3, 1, 1, 1, 1 \xrightarrow{Reduction} 2, 1, 1, 0, 0\\
     140              ...
     141
     142          In order to extend a current degree sequence while strictly increasing
     143          its maximum degree, it is equivalent to pick a set `I` of elements of
     144          the degree sequence with `|I|>\Delta` in such a way that the
     145          `(d_i+1)_{i\in I}` are the `|I|` maximum elements of the sequence
     146          `(d_i+\genfrac{}{}{0pt}{}{1\text{ if }i\in I}{0\text{ if }i\not \in
     147          I})_{1\leq i \leq n}`, and to add to this new sequence an element of
     148          value `|I|`. The non-increasing sequence containing the elements `|I|`
     149          and `(d_i+\genfrac{}{}{0pt}{}{1\text{ if }i\in I}{0\text{ if }i\not
     150          \in I})_{1\leq i \leq n}` can be reduced to `(d_i)_{1\leq i \leq n}`
     151          by Havel and Hakimi's rule.
     152
     153          .. MATH::
     154              ... 1, 1, 2, {\bf 2}, {\bf 2}, 2, 2, 3, 3, \underline{3}, {\bf 3}, {\bf 3}, {\bf 4}, {\bf 6}, ... \xrightarrow{Extension} ... 1, 1, 2, 2, 2, 3, 3, \underline{3}, {\bf 3}, {\bf 3}, {\bf 4}, {\bf 4}, {\bf 5}, {\bf 7}, ...
     155
     156          The number of possible sets `I` having this property (i.e. the number
     157          of possible extensions of a sequence) is smaller than it
     158          seems. Indeed, by definition, if `j\not \in I` then for all `i\in I`
     159          the inequality `d_j\leq d_i+1` holds. Hence, each set `I` is entirely
     160          determined by the largest element `d_k` of the sequence that it does
     161          **not** contain (hence `I` contains `\{1,...,k-1\}`), and by the
     162          cardinalities of `\{i\in I:d_i= d_k\}` and `\{i\in I:d_i= d_k-1\}`.
     163
     164          .. MATH::
     165              I = \{i \in I : d_i= d_k \} \cup \{i \in I : d_i= d_k-1 \} \cup \{i : d_i> d_k \}
     166
     167          The number of possible extensions is hence at most cubic, and is
     168          easily enumerated.
     169
     170About the implementation
     171~~~~~~~~~~~~~~~~~~~~~~~~
     172
     173In the actual implementation of the enumeration algorithm, the degree sequence
     174is stored differently for reasons of efficiency.
     175
     176Indeed, when enumerating all the degree sequences of length `n`, Sage first
     177allocates an array ``seq`` of `n+1` integers where ``seq[i]`` is the number of
     178elements of value ``i`` in the current sequence. Obviously, ``seq[n]=0`` holds
     179in permanence : it is useful to allocate a larger array than necessary to
     180simplify the code. The ``seq`` array is a global variable.
     181
     182The recursive function ``enum(depth, maximum)`` is the one building the list of
     183sequences. It builds the list of degree sequences of length `n` which *extend*
     184the sequence currently stored in ``seq[0]...seq[depth-1]``. When it is called,
     185``maximum`` must be set to the maximum value of an element in the partial
     186sequence ``seq[0]...seq[depth-1]``.
     187
     188If during its run the function ``enum`` heavily works on the content of the
     189``seq`` array, the value of ``seq`` is the **same** before and after the run of
     190``enum``.
     191
     192**Extending the current partial sequence**
     193
     194The two cases for which the maximum degree of the partial sequence does not
     195change are easy to detect. It is (sligthly) harder to enumerate all the sets `I`
     196corresponding to possible extensions of the partial sequence. As said
     197previously, to each set `I` one can associate an integer ``current_box`` such
     198that `I` contains all the `i` satisfying `d_i>current\_box`. The variable
     199``taken`` represents the number of all such elements `i`, so that when
     200enumerating all possible sets `I` in the algorithm we have the equality
     201
     202.. MATH::
     203    I = \text{taken }+\text{ number of elements of value }current\_box+ \text{ number of elements of value }current\_box-1
     204
     205References
     206~~~~~~~~~~
     207
     208  .. [RCES] Alley CATs in search of good homes
     209    Ruskey, R. Cohen, P. Eades, A. Scott
     210    Congressus numerantium, 1994
     211    Pages 97--110
     212
     213
     214Author
     215~~~~~~
     216
     217Nathann Cohen
     218
     219Tests
     220~~~~~
     221
     222The sequences produced by random graphs *are* degree sequences::
     223
     224    sage: n = 30
     225    sage: DS = DegreeSequences(30)
     226    sage: for i in range(10):
     227    ...      g = graphs.RandomGNP(n,.2)
     228    ...      if not g.degree_sequence() in DS:
     229    ...          print "Something is very wrong !"
     230
     231Checking that we indeed enumerate *all* the degree sequences for `n=5`::
     232
     233    sage: ds1 = Set([tuple(g.degree_sequence()) for g in graphs(5)])
     234    sage: ds2 = Set(map(tuple,list(DegreeSequences(5))))
     235    sage: ds1 == ds2
     236    True
     237
     238Checking the consistency of enumeration and test::
     239
     240    sage: DS = DegreeSequences(6)
     241    sage: all(seq in DS for seq in DS)
     242    True
     243
     244.. WARNING::
     245
     246    For the moment, iterating over all degree sequences involves building the
     247    list of them first, then iterate on this list.  This is obviously bad, as it
     248    requires uselessly a **lot** of memory for large values of `n`.
     249
     250    As soon as the ``yield`` keyword is available in Cython this should be
     251    changed. Updating the code does not require more than a couple of minutes.
     252   
     253"""
     254
     255##############################################################################
     256#       Copyright (C) 2011 Nathann Cohen <nathann.cohen@gmail.com>
     257#  Distributed under the terms of the GNU General Public License (GPL)
     258#  The full text of the GPL is available at:
     259#                  http://www.gnu.org/licenses/
     260##############################################################################
     261
     262from sage.libs.gmp.all cimport mpz_t
     263from sage.libs.gmp.all cimport *
     264from sage.rings.integer cimport Integer
     265include '../../../../devel/sage/sage/ext/stdsage.pxi'
     266include '../ext/cdefs.pxi'           
     267include "../ext/interrupt.pxi"
     268
     269
     270cdef unsigned char * seq
     271cdef list sequences
     272
     273class DegreeSequences:
     274
     275    def __init__(self, n):
     276        r"""
     277        Constructor
     278
     279        TEST::
     280
     281            sage: DegreeSequences(6)
     282            Degree sequences on 6 elements
     283        """
     284        self._n = n
     285
     286    def __contains__(self, seq):
     287        """
     288        Checks whether a given integer sequence is the degree sequence
     289        of a graph on `n` elements
     290
     291        EXAMPLE::
     292
     293            sage: [3,3,2,2,2,2,2,2] in DegreeSequences(8)
     294            True
     295        """
     296        cdef int n = self._n
     297        if len(seq)!=n:
     298            return False
     299       
     300        cdef int S = sum(seq)
     301
     302        # Partial represents the left side of Erdos and Gallai's inequality,
     303        # i.e. the sum of the i first integers.
     304        cdef int partial = 0
     305        cdef int i,d,dd, right
     306
     307        # Temporary variable to ensure that the sequence is indeed
     308        # non-increasing
     309        cdef int prev = n-1
     310
     311        for i,d in enumerate(seq):
     312
     313            # Non-increasing ?
     314            if d > prev:
     315                return False
     316            else:
     317                prev = d
     318
     319            # Updating the partial sum
     320            partial += d
     321
     322            # Evaluating the right hand side
     323            right = i*(i+1)
     324            for dd in seq[i+1:]:
     325                right += min(dd,i+1)
     326
     327            # Comparing the two
     328            if partial > right:
     329                return False
     330
     331        return True
     332
     333    def __repr__(self):
     334        """
     335        Representing the element
     336
     337        TEST::
     338
     339            sage: DegreeSequences(6)
     340            Degree sequences on 6 elements
     341        """
     342        return "Degree sequences on "+str(self._n)+" elements"
     343
     344    def __iter__(self):
     345        """
     346        Iterate over all the degree sequences.
     347
     348        TODO: THIS SHOULD BE UPDATED AS SOON AS THE YIELD KEYWORD APPEARS IN
     349        CYTHON. See comment in the class' documentation.
     350
     351        EXAMPLE::
     352
     353            sage: DS = DegreeSequences(6)
     354            sage: all(seq in DS for seq in DS)
     355            True
     356        """
     357
     358        init(self._n)
     359        return iter(sequences)
     360
     361    def __dealloc__():
     362        """
     363        Freeing the memory
     364        """
     365        if seq != NULL:
     366            free(seq)
     367           
     368cdef init(int n):
     369    """
     370    Initializes the memory and starts the enumeration algorithm.
     371    """
     372    global seq
     373    global N
     374    global sequences
     375
     376    if n == 0:
     377        return [[]]
     378    elif n == 1:
     379        return [[0]]
     380
     381    seq = <unsigned char *> malloc((n+1)*sizeof(unsigned char))
     382    memset(seq,0,(n+1)*sizeof(unsigned char))
     383
     384    # We begin with one vertex of degree 0
     385    seq[0] = 1
     386
     387    N = n
     388    sequences = []
     389    enum(1,0)
     390    free(seq)
     391    return sequences
     392
     393cdef inline add_seq():
     394     """
     395     This function is called whenever a sequence is found. 
     396
     397     Build the degree sequence corresponding to the current state of the
     398     algorithm and adds it to the sequences list.
     399     """
     400     global sequences
     401     global N
     402     global seq
     403
     404     cdef list s = []
     405     cdef int i, j
     406
     407     for N > i >= 0:
     408         for 0<= j < seq[i]:
     409             s.append(i)
     410
     411     sequences.append(s)   
     412
     413cdef void enum(int k, int M):
     414    """
     415    Main function. For an explanation of the algorithm please refer to the
     416    class' documentation.
     417
     418    INPUT:
     419
     420    * ``k`` -- depth of the partial degree sequence
     421    * ``M`` -- value of a maximum element in the partial degree sequence
     422    """
     423    cdef int i,j
     424    global seq
     425    cdef int taken = 0
     426    cdef int current_box
     427    cdef int n_current_box
     428    cdef int n_previous_box
     429    cdef int new_vertex
     430
     431    # Have we found a new degree sequence ? End of recursion !
     432    if k == N:
     433        add_seq()
     434        return
     435
     436    _sig_on
     437
     438    #############################################
     439    # Creating vertices of Vertices of degree M #
     440    #############################################
     441
     442    # If 0 is the current maximum degree, we can always extend the degree
     443    # sequence with another 0
     444    if M == 0:
     445
     446        seq[0] += 1
     447        enum(k+1, M)
     448        seq[0] -= 1
     449
     450    # We need not automatically increase the degree at each step. In this case,
     451    # we have no other choice but to link the new vertex of degree M to vertices
     452    # of degree M-1, which will become vertices of degree M too.
     453    elif seq[M-1] >= M:
     454
     455        seq[M]   += M+1
     456        seq[M-1] -= M
     457
     458        enum(k+1, M)
     459
     460        seq[M]   -= M+1
     461        seq[M-1] += M
     462
     463    ###############################################
     464    # Creating vertices of Vertices of degree > M #
     465    ###############################################
     466
     467    for M >= current_box > 0:
     468
     469        # If there is not enough vertices in the boxes available
     470        if taken + (seq[current_box] - 1) + seq[current_box-1] <= M:
     471            taken += seq[current_box]
     472            seq[current_box+1] += seq[current_box]
     473            seq[current_box] = 0
     474            continue
     475
     476        # The degree of the new vertex will be taken + i + j where :
     477        #
     478        # * i is the number of vertices taken in the *current* box
     479        # * j the number of vertices taken in the *previous* one
     480
     481        n_current_box = seq[current_box]
     482        n_previous_box = seq[current_box-1]
     483
     484        # Note to self, and others :
     485        #
     486        # In the following lines, there are many incrementation/decrementation
     487        # that *may* be replaced by only +1 and -1 and save some
     488        # instructions. This would involve adding several "if", and I feared it
     489        # would make the code even uglier. If you are willing to give it a try,
     490        # **please check the results** ! It is trickier that it seems ! Even
     491        # changing the lower bounds in the for loops would require tests
     492        # afterwards.
     493
     494        for max(0,((M+1)-n_previous_box-taken)) <= i < n_current_box:
     495            seq[current_box] -= i
     496            seq[current_box+1] += i
     497
     498            for max(0,((M+1)-taken-i)) <= j <= n_previous_box:
     499                seq[current_box-1] -= j
     500                seq[current_box] += j
     501
     502                new_vertex = taken + i + j
     503                seq[new_vertex] += 1
     504                enum(k+1,new_vertex)
     505                seq[new_vertex] -= 1
     506
     507                seq[current_box-1] += j
     508                seq[current_box] -= j
     509
     510            seq[current_box] += i
     511            seq[current_box+1] -= i
     512
     513        taken += n_current_box
     514        seq[current_box] = 0
     515        seq[current_box+1] += n_current_box
     516
     517    # Corner case
     518    #
     519    # Now current_box = 0. All the vertices of nonzero degree are taken, we just
     520    # want to know how many vertices of degree 0 will be neighbors of the new
     521    # vertex.
     522    for max(0,((M+1)-taken)) <= i <= seq[0]:
     523
     524        seq[1] += i
     525        seq[0] -= i
     526        seq[taken+i] += 1
     527
     528        enum(k+1, taken+i)
     529
     530        seq[taken+i] -= 1
     531        seq[1] -= i
     532        seq[0] += i
     533
     534    # Shift everything back to normal ! ( cell N is always equal to 0)
     535    for 1 <= i < N:
     536        seq[i] = seq[i+1]
     537
     538    _sig_off