Ticket #7301: trac_7301.patch

File trac_7301.patch, 12.9 KB (added by ncohen, 13 years ago)
  • sage/combinat/integer_vector.py

    # HG changeset patch
    # User David Joyner <wdjoyner@gmail.com>
    # Date 1262783697 18000
    # Node ID ff7588ec0fdf194c4fb3254f5fddb53cdff8dacd
    # Parent  a5ee50fe0c1efa6e7b70aaad21314bfc3bcca629
    gale_ryser in integer_vector
    
    diff -r a5ee50fe0c1e -r ff7588ec0fdf sage/combinat/integer_vector.py
    a b  
    11"""
    22Integer vectors
     3
     4AUTHORS:
     5 *   Mike Hanson (2007) - original module
     6 *   Nathann Cohen, David Joyner (2009-2010) - Gale-Ryser stuff
    37"""
    48#*****************************************************************************
    59#       Copyright (C) 2007 Mike Hansen <mhansen@gmail.com>,
     
    2731import functools
    2832
    2933
     34def is_gale_ryser(r,s):
     35    r"""
     36    Tests whether the given sequences satisfy the condition
     37    of the Gale-Ryser theorem.
     38
     39    Given a binary matrix `B` of dimension `n\times m`, the
     40    vector of row sums is defined as the vector whose
     41    `i^{\mbox{th}}` component is equal to the sum of the `i^{\mbox{th}}`
     42    row in `A`. The vector of column sums is defined similarly.
     43
     44    If, given a binary matrix, these two vectors are easy to compute,
     45    the Gale-Ryser theorem lets us decide whether, given two
     46    non-negative vectors `r,s`, there exists a binary matrix
     47    whose row/colum sums vectors are `r` and `s`.
     48
     49    This functions answers accordingly.
     50
     51    INPUT: 
     52
     53    - ``r``, ``s`` -- lists of non-negative integers.
     54
     55    ALGORITHM:
     56
     57    Without loss of generality, we can assume that :
     58   
     59    - The two given sequences do not contain any `0` ( which
     60    would correspond to an empty column/row )
     61
     62    - The two given sequences are ordered in decreasing order
     63    (reordering the sequence of row (resp. column) sums amounts
     64    to reordering the rows (resp. columns) themselves in the
     65    matrix, which does not alter the columns (resp. rows) sums.
     66
     67    We can then assume that `r` and `s` are partitions
     68    (see the corresponding class ``Partition``)
     69
     70    If `r^*` denote the conjugate of `r`, the Gale-Ryser theorem
     71    asserts that a binary Matrix satisfying the constraints exists
     72    if and only if `s\preceq r^*`, where `\preceq` denotes
     73    the domination order on partitions.
     74
     75    EXAMPLES::
     76
     77        sage: from sage.combinat.integer_vector import is_gale_ryser
     78        sage: is_gale_ryser([4,2,2],[3,3,1,1])
     79        True
     80        sage: is_gale_ryser([4,2,1,1],[3,3,1,1])
     81        True
     82        sage: is_gale_ryser([3,2,1,1],[3,3,1,1])
     83        False
     84
     85    REMARK: In the literature, what we are calling a
     86    Gale-Ryser sequence sometimes goes by the (rather
     87    generic-sounding) term ''realizable sequence''.
     88    """
     89
     90    # The sequences only contan non-negative integers
     91    if [x for x in r if x<0] or [x for x in s if x<0]:
     92        return False
     93
     94    # builds the corresponding partitions, i.e.
     95    # removes the 0 and sorts the sequences
     96    from sage.combinat.partition import Partition
     97    r2 = Partition(sorted([x for x in r if x>0], reverse=True))
     98    s2 = Partition(sorted([x for x in s if x>0], reverse=True))
     99    rstar = Partition(r2).conjugate()
     100
     101    #                                same number of 1s           domination
     102    return len(rstar) <= len(s2) and  sum(r2) == sum(s2) and rstar.dominates(s)
     103
     104def _slider01(A, t, k):
     105    r"""
     106    Assumes `A` is a `(0,1)`-matrix. For each of the
     107    `t` rows with highest row sums, this function
     108    returns a matrix `B` which is the same as `A` except that it
     109    has slid exactly one `1` in `A` over to the `k`-th
     110    column.
     111
     112    This is a 'private' function for use in gale_ryser_theorem.
     113
     114    INPUT:
     115
     116    - ``A`` -- an `m\times n` (0,1) matrix
     117    - ``t``, ``k`` -- integers satisfying `0 < t < m`, `0 < k < n`
     118
     119    OUTPUT:
     120
     121    An `m\times n` (0,1) matrix, which is the same as `A` except
     122    that it has exactly one `1` in `A` slid over to the `k`-th
     123    column.
     124
     125    EXAMPLES::
     126
     127        sage: from sage.combinat.integer_vector import _slider01
     128        sage: A = matrix([[1,1,1,0],[1,1,1,0],[1,0,0,0],[1,0,0,0]])
     129        sage: A
     130        [1 1 1 0]
     131        [1 1 1 0]
     132        [1 0 0 0]
     133        [1 0 0 0]
     134        sage: _slider01(A, 1, 3)
     135        [1 1 1 0]
     136        [1 1 0 1]
     137        [1 0 0 0]
     138        [1 0 0 0]
     139        sage: _slider01(A, 3, 3)
     140        [1 1 1 0]
     141        [1 1 1 0]
     142        [1 0 0 0]
     143        [0 0 0 1]
     144
     145    """
     146    from sage.matrix.constructor import matrix
     147    maxrows = []
     148    m = len(A.rows())
     149    n = len(A.columns())
     150    B = matrix(m,n)
     151    rowsums = [sum(A.row(i)) for i in range(m)]
     152    maxrows = [rowsums[i] for i in range(t)] # t rows with max rowsums
     153    maxRows = [] # want these to correspond to elts of maxrows
     154    B = [list(A.row(i)) for i in range(m)]
     155    for i in range(m):
     156        if sum(A.row(i)) in maxrows:
     157            maxRows.append(i)  # >= t rows numbers with max rowsums
     158    for j in range(len(maxrows)):
     159        i = maxRows[-1]
     160        rw = A.row(i) # to make mutable
     161        # now we want to move the rightmost left 1 to the k-th column
     162        rwlist = list(rw)
     163        if not(0 in rwlist):
     164            N = n-1
     165        else:
     166            N = rwlist.index(0)
     167        rwlist[int(N)-1] = 0
     168        rwlist[k] = 1
     169        B[i] = rwlist
     170    return matrix(B)
     171
     172def gale_ryser_theorem(p1, p2, algorithm="ryser"):
     173        r"""
     174        Returns the binary matrix given by the Gale-Ryser theorem.
     175
     176        The Gale Ryser theorem asserts that if `p_1,p_2` are two
     177        partitions of `n` of respective lengths `k_1,k_2`, then there is
     178        a binary `k_1\times k_2` matrix `M` such that `p_1` is the vector
     179        of row sums and `p_2` is the vector of column sums of `M`, if
     180        and only if `p_2` dominates `p_1`.
     181
     182        INPUT:
     183
     184        - ``p1`` -- the first partition of `n` (trailing 0's allowed)
     185        - ``p2`` -- the second partition of `n` (trailing 0's allowed)
     186
     187        - ``algorithm`` -- two possible string values :
     188
     189            - ``"ryser"`` (default) implements the construction due
     190              to Ryser [Ryser63]_.
     191
     192            - ``"gale"`` implements the construction due to Gale [Gale57]_.
     193
     194
     195        OUTPUT:
     196
     197        - A binary matrix if it exists, ``None`` otherwise.
     198
     199
     200        Gale's Algorithm:
     201
     202        (Gale [Gale57]_): A matrix satisfying the constraints of its
     203        sums can be defined as the solution of the following
     204        Linear Program, which Sage knows how to solve (requires
     205        packages GLPK or CBC).
     206
     207        .. MATH::
     208            \forall i&\sum_{j=1}^{k_2} b_{i,j}=p_{1,j}\\
     209            \forall i&\sum_{j=1}^{k_1} b_{j,i}=p_{2,j}\\
     210            &b_{i,j}\mbox{ is a binary variable}
     211
     212
     213        Ryser's Algorithm:
     214
     215        (Ryser [Ryser63]_): The construction of an `m\times n` matrix `A=A_{r,s}`,
     216        due to Ryser, is described as follows. The
     217        construction works if and only if have `s\preceq r^*`.
     218
     219        * Construct the `m\times n` matrix `B` from `r` by
     220        defining the `i`-th row of `B` to be the
     221        vector whose first `r_i` entries are `1`, and the
     222        remainder are `0`'s, `1\leq i\leq m`.
     223        This maximal matrix `B` with row sum `r` and ones
     224        left justified has column sum `r^{*}`.
     225
     226        * Shift the last `1` in certain rows of `B` to
     227        column `n` in order to achieve the sum `s_n`.
     228        Call this `B` again.
     229
     230          * The `1`'s in column n are to appear in those
     231            rows in which `A` has the largest row sums, giving
     232            preference to the bottom-most positions in case of ties.
     233          * Note: When this step automatically "fixes" other columns,
     234            one must skip ahead to the first column index
     235            with a wrong sum in the step below.
     236
     237        * Proceed inductively to construct columns `n-1`, ..., `2`, `1`.
     238
     239        * Set `A = B`. Return `A`.
     240
     241        EXAMPLES:
     242
     243        Computing the matrix for `p_1=p_2=2+2+1` ::
     244
     245            sage: from sage.combinat.integer_vector import gale_ryser_theorem
     246            sage: p1 = [2,2,1]
     247            sage: p2 = [2,2,1]
     248            sage: print gale_ryser_theorem(p1, p2, algorithm="gale") # Optional - requires GLPK or CBC
     249            [0 1 1]
     250            [1 1 0]
     251            [1 0 0]       
     252
     253        Or for a non-square matrix with `p_1=3+3+2+1` and `p_2=3+2+2+1+1` ::
     254
     255            sage: from sage.combinat.integer_vector import gale_ryser_theorem
     256            sage: p1 = [3,3,1,1]
     257            sage: p2 = [3,3,1,1]
     258            sage: gale_ryser_theorem(p1, p2)
     259            [1 1 1 0]
     260            [1 1 0 1]
     261            [1 0 0 0]
     262            [0 1 0 0]
     263            sage: p1 = [4,2,2]
     264            sage: p2 = [3,3,1,1]
     265            sage: gale_ryser_theorem(p1, p2)
     266            [1 1 1 1]
     267            [1 1 0 0]
     268            [1 1 0 0]       
     269            sage: p1 = [4,2,2,0]
     270            sage: p2 = [3,3,1,1,0,0]
     271            sage: gale_ryser_theorem(p1, p2)
     272            [1 1 1 1 0 0]
     273            [1 1 0 0 0 0]
     274            [1 1 0 0 0 0]
     275            [0 0 0 0 0 0]
     276            sage: p1 = [3,3,2,1]
     277            sage: p2 = [3,2,2,1,1]
     278            sage: print gale_ryser_theorem(p1, p2, algorithm="gale") # Optional - requires GLPK or CBC
     279            [1 0 1 1 0]
     280            [1 0 1 0 1]
     281            [1 1 0 0 0]
     282            [0 1 0 0 0]
     283
     284        With `0` in the sequences, and with unordered inputs ::
     285
     286            sage: from sage.combinat.integer_vector import gale_ryser_theorem
     287            sage: gale_ryser_theorem([3,3,0,1,1,0], [3,1,3,1,0])     
     288            [1 1 1 0 0]
     289            [1 0 1 1 0]
     290            [0 0 0 0 0]
     291            [1 0 0 0 0]
     292            [0 0 1 0 0]
     293            [0 0 0 0 0]
     294
     295        REFERENCES:
     296
     297        ..  [Ryser63] H. J. Ryser, Combinatorial Mathematics,
     298                Carus Monographs, MAA, 1963.
     299        ..  [Gale57] D. Gale, A theorem on flows in networks, Pacific J. Math.
     300                7(1957)1073-1082.
     301        """
     302        from sage.combinat.partition import Partition
     303        from sage.matrix.constructor import matrix
     304
     305        if not(is_gale_ryser(p1,p2)):
     306            return False
     307
     308
     309        if algorithm=="ryser": # ryser's algorithm
     310            from sage.combinat.permutation import Permutation
     311
     312            # Sorts the sequences if they are not, and remembers the permutation
     313            # applied
     314            tmp = sorted(enumerate(p1), reverse=True, key=lambda x:x[1])
     315            r = [x[1] for x in tmp if x[1]>0]
     316            r_permutation = [x-1 for x in Permutation([x[0]+1 for x in tmp]).inverse()]
     317            m = len(r)
     318
     319            tmp = sorted(enumerate(p2), reverse=True, key=lambda x:x[1])
     320            s = [x[1] for x in tmp if x[1]>0]
     321            s_permutation = [x-1 for x in Permutation([x[0]+1 for x in tmp]).inverse()]
     322            n = len(s)
     323
     324            rowsA0 = [[0]*n]*m
     325            for j in range(m):
     326                if j<m:
     327                    rowsA0[j] = [1]*r[j]+[0]*(n-r[j])
     328                else:
     329                    rowsA0[j] = [0]*n
     330            A0 = matrix(rowsA0)
     331
     332            for j in range(1,n-1):
     333                # starts for loop, k = n-1, ..., 1
     334                # which finds the 1st column with
     335                # incorrect column sum. For that bad
     336                # column index, apply slider again
     337                for k in range(1,n):
     338                    if sum(A0.column(n-k))<>s[n-k]:
     339                         break
     340                A0 = _slider01(A0,s[n-k],n-k)
     341
     342            # If we need to add empty rows/columns
     343            if len(p1)!=m:
     344                A0 = A0.stack(matrix([[0]*n]*(len(p1)-m)))
     345
     346            if len(p2)!=n:
     347                A0 = A0.transpose().stack(matrix([[0]*len(p1)]*(len(p2)-n))).transpose()
     348               
     349            # Applying the permutations to get a matrix satisfying the
     350            # order given by the input
     351            A0 = A0.matrix_from_rows_and_columns(r_permutation, s_permutation)
     352            return A0   
     353
     354        elif algorithm == "gale":
     355          from sage.numerical.mip import MixedIntegerLinearProgram
     356          k1, k2=len(p1), len(p2)
     357          p = MixedIntegerLinearProgram()
     358          b = p.new_variable(dim=2)
     359          for (i,c) in enumerate(p1):
     360              p.add_constraint(sum([b[i][j] for j in xrange(k2)]),min=c,max=c)
     361          for (i,c) in enumerate(p2):
     362              p.add_constraint(sum([b[j][i] for j in xrange(k1)]),min=c,max=c)
     363          p.set_objective(None)
     364          p.set_binary(b)
     365          p.solve()
     366          b = p.get_values(b)
     367          M = [[0]*k2 for i in xrange(k1)]
     368          for i in xrange(k1):
     369              for j in xrange(k2):
     370                  M[i][j] = int(b[i][j])
     371          return matrix(M)
     372
     373        else:
     374            raise ValueError("The only two algorithms available are \"gale\" and \"ryser\"")
     375
     376
    30377def _default_function(l, default, i):
    31378    """
    32379    EXAMPLES::
     
    214561        """
    215562        return infinity
    216563
     564
    217565class IntegerVectors_nk(CombinatorialClass):
    218566    def __init__(self, n, k):
    219567        """