Ticket #7301: trac_7301-referee.patch

File trac_7301-referee.patch, 10.5 KB (added by wdj, 13 years ago)

seems to apply to 4.3* and test okay. Apply only this patch.

  • sage/combinat/integer_vector.py

    # HG changeset patch
    # User David Joyner <wdjoyner@gmail.com>
    # Date 1262783697 18000
    # Node ID 44f1def28dfe7c24d79ecc69ae814bbac0632078
    # Parent  afe90cf65c0ebc94601de63ab6b1c44e8a29c872
    followed referees suggestions for trac #7301
    
    diff -r afe90cf65c0e -r 44f1def28dfe 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    INPUT:  r, s - integer sequences `r=(r_1,...,r_m)`
     37            and `s=(s_1,...,s_n)`, where
     38            `r_1\geq r_2 ... \geq r_m\geq 0`,
     39            `s_1\geq s_2 ... \geq s_n\geq 0`.
     40
     41    Let `r^*` denote the conjugate of `r`.
     42 
     43    OUTPUT: The boolean value of `s\preceq r^*`, where `\preceq` is the ordering on sequences above.
     44   
     45    We say that a pair of such sequences `(r,s)` is
     46    ''Gale-Ryser'' if  `s\preceq r^*` (in which case
     47    the hypothesis to the Gale-Ryser Theorem is
     48    satisfied).
     49
     50    EXAMPLES::
     51
     52        sage: from sage.combinat.integer_vector import is_gale_ryser
     53        sage: is_gale_ryser([4,2,2],[3,3,1,1])
     54        True
     55        sage: is_gale_ryser([4,2,1,1],[3,3,1,1])
     56        True
     57        sage: is_gale_ryser([3,2,1,1],[3,3,1,1])
     58        False
     59
     60    REMARK: In the literature, what we are calling a
     61    Gale-Ryser sequence sometimes goes by the (rather
     62    generic-sounding) term ''realizable sequence''.
     63    """
     64    from sage.combinat.partition import Partition
     65    r2 = [x for x in r if x>0] # remove extraneous 0's
     66    s2 = [x for x in s if x>0]
     67    rstar = Partition(r2).conjugate()
     68    m = len(rstar)
     69    if m > len(s2):
     70        return False
     71    for i in range(1,m):
     72        #print i, sum(s2[:i]), sum(rstar[:i])
     73        if sum(rstar[:i]) < sum(s2[:i]):
     74            return False
     75    if sum(rstar) <> sum(s2):
     76            return False
     77    return True
     78
     79def _slider01(A, t, k):
     80    r"""
     81    Assumes `A` is a `(0,1)`-matrix. For each of the
     82    `t` rows with highest row sums, this function
     83    returns a matrix `B` which is the same as `A` except that it
     84    has slid exactly one `1` in `A` over to the `k`-th
     85    column.
     86
     87    This is a 'private' function for use in gale_ryser_theorem.
     88
     89    INPUT::
     90        A - an `m\times n` (0,1) matrix
     91        t, k - integers satisfying `0 < t < m`, `0 < k < n`
     92
     93    OUTPUT::
     94        An `m\times n` (0,1) matrix, which is the same as `A` except
     95        that it has exactly one `1` in `A` slid over to the `k`-th
     96        column.
     97
     98    EXAMPLES::
     99
     100        sage: from sage.combinat.integer_vector import _slider01
     101        sage: A = matrix([[1,1,1,0],[1,1,1,0],[1,0,0,0],[1,0,0,0]])
     102        sage: A
     103        [1 1 1 0]
     104        [1 1 1 0]
     105        [1 0 0 0]
     106        [1 0 0 0]
     107        sage: _slider01(A, 1, 3)
     108        [1 1 1 0]
     109        [1 1 0 1]
     110        [1 0 0 0]
     111        [1 0 0 0]
     112        sage: _slider01(A, 3, 3)
     113        [1 1 1 0]
     114        [1 1 1 0]
     115        [1 0 0 0]
     116        [0 0 0 1]
     117
     118    """
     119    from sage.matrix.constructor import matrix
     120    maxrows = []
     121    m = len(A.rows())
     122    n = len(A.columns())
     123    B = matrix(m,n)
     124    rowsums = [sum(A.row(i)) for i in range(m)]
     125    maxrows = [rowsums[i] for i in range(t)] # t rows with max rowsums
     126    maxRows = [] # want these to correspond to elts of maxrows
     127    B = [list(A.row(i)) for i in range(m)]
     128    for i in range(m):
     129        if sum(A.row(i)) in maxrows:
     130            maxRows.append(i)  # >= t rows numbers with max rowsums
     131    for j in range(len(maxrows)):
     132        i = maxRows[-1]
     133        rw = A.row(i) # to make mutable
     134        # now we want to move the rightmost left 1 to the k-th column
     135        rwlist = list(rw)
     136        if not(0 in rwlist):
     137            N = n-1
     138        else:
     139            N = rwlist.index(0)
     140        rwlist[int(N)-1] = 0
     141        rwlist[k] = 1
     142        B[i] = rwlist
     143    return matrix(B)
     144
     145def gale_ryser_theorem(p1, p2, algorithm="ryser"):
     146        r"""
     147        Returns the binary matrix given by the Gale-Ryser theorem.
     148
     149        The Gale Ryser theorem asserts that if `p_1,p_2` are two
     150        partitions of `n` of respective lengths `k_1,k_2`, then there is
     151        a binary `k_1\times k_2` matrix `M` such that `p_1` is the vector
     152        of row sums and `p_2` is the vector of column sums of `M`, if
     153        and only if `p_2` dominates `p_1`.
     154
     155        INPUT:
     156
     157        - ``p1`` -- the first partition of `n` (trailing 0's allowed)
     158        - ``p2`` -- the second partition of `n` (trailing 0's allowed)
     159
     160        OUTPUT:
     161
     162        - A binary matrix
     163
     164        OPTION:
     165
     166            algorithm - "ryser" (default) implements the construction due
     167                         to Ryser [Ryser63]_.
     168                       - "gale" (anything other than "gale" uses "ryser")
     169                         implements the construction due to Gale [Gale57]_.
     170
     171        ALGORITHM:
     172
     173        (Gale [Gale57]_): The Gale-Ryser theorem asserts that the following Linear
     174        Program admits a solution.
     175
     176        .. MATH::
     177            \forall i&\sum_{j=1}^{k_2} b_{i,j}=p_{1,j}\\
     178            \forall i&\sum_{j=1}^{k_1} b_{j,i}=p_{2,j}\\
     179            &b_{i,j}\mbox{ is a binary variable}
     180
     181        Sage uses MixedIntegerLinearProgram to solve this.
     182
     183        (Ryser [Ryser63]_): The construction of an `m\times n` matrix `A=A_{r,s}`,
     184        due to Ryser, is described as follows. The
     185        construction works if and only if have `s\preceq r^*`.
     186
     187        * Construct the `m\times n` matrix `B` from `r` by
     188        defining the `i`-th row of `B` to be the
     189        vector whose first `r_i` entries are `1`, and the
     190        remainder are `0`'s, `1\leq i\leq m`.
     191        This maximal matrix `B` with row sum `r` and ones
     192        left justified has column sum `r^{*}`.
     193
     194        * Shift the last `1` in certain rows of `B` to
     195        column `n` in order to achieve the sum `s_n`.
     196        Call this `B` again.
     197
     198          * The `1`'s in column n are to appear in those
     199            rows in which `A` has the largest row sums, giving
     200            preference to the bottom-most positions in case of ties.
     201          * Note: When this step automatically "fixes" other columns,
     202            one must skip ahead to the first column index
     203            with a wrong sum in the step below.
     204
     205        * Proceed inductively to construct columns `n-1`, ..., `2`, `1`.
     206
     207        * Set `A = B`. Return `A`.
     208
     209        EXAMPLES:
     210
     211        Computing the matrix for `p_1=p_2=2+2+1` ::
     212
     213            sage: from sage.combinat.integer_vector import gale_ryser_theorem
     214            sage: p1 = [2,2,1]
     215            sage: p2 = [2,2,1]
     216            sage: print gale_ryser_theorem(p1, p2, algorithm="gale") # Optional - requires GLPK or CBC
     217            [0 1 1]
     218            [1 1 0]
     219            [1 0 0]       
     220
     221        Or for a non-square matrix with `p_1=3+3+2+1` and `p_2=3+2+2+1+1` ::
     222
     223            sage: p1 = [3,3,2,1]
     224            sage: p2 = [3,2,2,1,1]
     225            sage: print gale_ryser_theorem(p1, p2, algorithm="gale") # Optional - requires GLPK or CBC
     226            [1 0 1 1 0]
     227            [1 0 1 0 1]
     228            [1 1 0 0 0]
     229            [0 1 0 0 0]
     230            sage: p1 = [3,3,1,1]
     231            sage: p2 = [3,3,1,1]
     232            sage: gale_ryser_theorem(p1, p2)
     233            [1 1 1 0]
     234            [1 1 0 1]
     235            [1 0 0 0]
     236            [0 1 0 0]
     237            sage: p1 = [4,2,2]
     238            sage: p2 = [3,3,1,1]
     239            sage: gale_ryser_theorem(p1, p2)
     240            [1 1 1 1]
     241            [1 1 0 0]
     242            [1 1 0 0]       
     243            sage: p1 = [4,2,2,0]
     244            sage: p2 = [3,3,1,1,0,0]
     245            sage: gale_ryser_theorem(p1, p2)
     246            [1 1 1 1 0 0]
     247            [1 1 0 0 0 0]
     248            [1 1 0 0 0 0]
     249            [0 0 0 0 0 0]
     250
     251        REFERENCES:
     252
     253        ..  [Ryser63] H. J. Ryser, Combinatorial Mathematics,
     254                Carus Monographs, MAA, 1963.
     255        ..  [Gale57] D. Gale, A theorem on flows in networks, Pacific J. Math.
     256                7(1957)1073-1082.
     257        """
     258        from sage.combinat.partition import Partition
     259        from sage.matrix.constructor import matrix
     260        if sum(p1) != sum(p2):
     261            raise ValueError("The two lists must sum to the same value.")
     262        if not Partition(p2).conjugate().dominates(Partition(p1)):
     263            raise ValueError("The domination criterion is not fulfilled.")
     264        if algorithm == "gale":
     265          from sage.numerical.mip import MixedIntegerLinearProgram
     266          p1 = Partition(p1)
     267          p2 = Partition(p2)
     268          k1=p1.length()
     269          k2=p2.length()
     270          p = MixedIntegerLinearProgram()
     271          b = p.new_variable(dim=2)
     272          for (i,c) in enumerate(p1):
     273              p.add_constraint(sum([b[i][j] for j in xrange(k2)]),min=c,max=c)
     274          for (i,c) in enumerate(p2):
     275              p.add_constraint(sum([b[j][i] for j in xrange(k1)]),min=c,max=c)
     276          p.set_objective(None)
     277          p.set_binary(b)
     278          p.solve()
     279          b = p.get_values(b)
     280          M = [[0]*k2 for i in xrange(k1)]
     281          for i in xrange(k1):
     282              for j in xrange(k2):
     283                  M[i][j] = int(b[i][j])
     284          return matrix(M)
     285        else: # ryser's algorithm
     286            r = p1; s = p2
     287            if not(is_gale_ryser(r,s)):
     288                return False
     289            rstar = Partition(r).conjugate()
     290            m = len(r)
     291            n = len(s)
     292            rowsA0 = [[0]*n]*m
     293            for j in range(m):
     294                if j<m:
     295                    rowsA0[j] = [1]*r[j]+[0]*(n-r[j])
     296                else:
     297                    rowsA0[j] = [0]*n
     298            A0 = matrix(rowsA0)
     299            for j in range(1,n-1):
     300                # starts for loop, k = n-1, ..., 1
     301                # which finds the 1st column with
     302                # incorrect column sum. For that bad
     303                # column index, apply slider again
     304                for k in range(1,n):
     305                    if sum(A0.column(n-k))<>s[n-k]:
     306                         break
     307                A0 = _slider01(A0,s[n-k],n-k)
     308            return A0   
     309
     310
    30311def _default_function(l, default, i):
    31312    """
    32313    EXAMPLES::
     
    214495        """
    215496        return infinity
    216497
     498
    217499class IntegerVectors_nk(CombinatorialClass):
    218500    def __init__(self, n, k):
    219501        """