Ticket #7301: trac_7301-referee.patch

sage/combinat/integer_vector.py

@@ -1 +1 @@
 """
 Integer vectors
+
+AUTHORS:
+ *   Mike Hanson (2007) - original module
+ *   Nathann Cohen, David Joyner (2009-2010) - Gale-Ryser stuff
 """
 #*****************************************************************************
 #       Copyright (C) 2007 Mike Hansen <mhansen@gmail.com>, 
@@ -27 +31 @@
 import functools
 
 
+def is_gale_ryser(r,s):
+    r"""
+    INPUT:  r, s - integer sequences `r=(r_1,...,r_m)` 
+            and `s=(s_1,...,s_n)`, where 
+            `r_1\geq r_2 ... \geq r_m\geq 0`,
+            `s_1\geq s_2 ... \geq s_n\geq 0`. 
+
+    Let `r^*` denote the conjugate of `r`.
+ 
+    OUTPUT: The boolean value of `s\preceq r^*`, where `\preceq` is the ordering on sequences above.
+    
+    We say that a pair of such sequences `(r,s)` is
+    ''Gale-Ryser'' if  `s\preceq r^*` (in which case
+    the hypothesis to the Gale-Ryser Theorem is 
+    satisfied).
+
+    EXAMPLES::
+
+        sage: from sage.combinat.integer_vector import is_gale_ryser
+        sage: is_gale_ryser([4,2,2],[3,3,1,1])
+        True
+        sage: is_gale_ryser([4,2,1,1],[3,3,1,1])
+        True
+        sage: is_gale_ryser([3,2,1,1],[3,3,1,1])
+        False
+
+    REMARK: In the literature, what we are calling a 
+    Gale-Ryser sequence sometimes goes by the (rather
+    generic-sounding) term ''realizable sequence''. 
+    """
+    from sage.combinat.partition import Partition
+    r2 = [x for x in r if x>0] # remove extraneous 0's
+    s2 = [x for x in s if x>0]
+    rstar = Partition(r2).conjugate()
+    m = len(rstar)
+    if m > len(s2):
+        return False
+    for i in range(1,m):
+        #print i, sum(s2[:i]), sum(rstar[:i])
+        if sum(rstar[:i]) < sum(s2[:i]):
+            return False
+    if sum(rstar) <> sum(s2):
+            return False
+    return True
+
+def _slider01(A, t, k):
+    r"""
+    Assumes `A` is a `(0,1)`-matrix. For each of the
+    `t` rows with highest row sums, this function
+    returns a matrix `B` which is the same as `A` except that it
+    has slid exactly one `1` in `A` over to the `k`-th 
+    column.
+
+    This is a 'private' function for use in gale_ryser_theorem.
+
+    INPUT::
+        A - an `m\times n` (0,1) matrix
+        t, k - integers satisfying `0 < t < m`, `0 < k < n`
+
+    OUTPUT::
+        An `m\times n` (0,1) matrix, which is the same as `A` except 
+        that it has exactly one `1` in `A` slid over to the `k`-th 
+        column.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.integer_vector import _slider01
+        sage: A = matrix([[1,1,1,0],[1,1,1,0],[1,0,0,0],[1,0,0,0]])
+        sage: A
+        [1 1 1 0]
+        [1 1 1 0]
+        [1 0 0 0]
+        [1 0 0 0]
+        sage: _slider01(A, 1, 3)
+        [1 1 1 0]
+        [1 1 0 1]
+        [1 0 0 0]
+        [1 0 0 0]
+        sage: _slider01(A, 3, 3)
+        [1 1 1 0]
+        [1 1 1 0]
+        [1 0 0 0]
+        [0 0 0 1]
+
+    """
+    from sage.matrix.constructor import matrix
+    maxrows = []
+    m = len(A.rows())
+    n = len(A.columns())
+    B = matrix(m,n)
+    rowsums = [sum(A.row(i)) for i in range(m)]
+    maxrows = [rowsums[i] for i in range(t)] # t rows with max rowsums
+    maxRows = [] # want these to correspond to elts of maxrows
+    B = [list(A.row(i)) for i in range(m)]
+    for i in range(m):
+        if sum(A.row(i)) in maxrows:
+            maxRows.append(i)  # >= t rows numbers with max rowsums
+    for j in range(len(maxrows)):
+        i = maxRows[-1]
+        rw = A.row(i) # to make mutable
+        # now we want to move the rightmost left 1 to the k-th column
+        rwlist = list(rw)
+        if not(0 in rwlist):
+            N = n-1
+        else:
+            N = rwlist.index(0)
+        rwlist[int(N)-1] = 0
+        rwlist[k] = 1
+        B[i] = rwlist
+    return matrix(B)
+
+def gale_ryser_theorem(p1, p2, algorithm="ryser"):
+        r"""
+        Returns the binary matrix given by the Gale-Ryser theorem.
+
+        The Gale Ryser theorem asserts that if `p_1,p_2` are two
+        partitions of `n` of respective lengths `k_1,k_2`, then there is
+        a binary `k_1\times k_2` matrix `M` such that `p_1` is the vector
+        of row sums and `p_2` is the vector of column sums of `M`, if
+        and only if `p_2` dominates `p_1`.
+
+        INPUT:
+
+        - ``p1`` -- the first partition of `n` (trailing 0's allowed)
+        - ``p2`` -- the second partition of `n` (trailing 0's allowed)
+
+        OUTPUT:
+
+        - A binary matrix
+
+        OPTION:
+
+            algorithm - "ryser" (default) implements the construction due 
+                         to Ryser [Ryser63]_.
+                       - "gale" (anything other than "gale" uses "ryser") 
+                         implements the construction due to Gale [Gale57]_.
+
+        ALGORITHM:
+
+        (Gale [Gale57]_): The Gale-Ryser theorem asserts that the following Linear 
+        Program admits a solution.
+
+        .. MATH:: 
+            \forall i&\sum_{j=1}^{k_2} b_{i,j}=p_{1,j}\\
+            \forall i&\sum_{j=1}^{k_1} b_{j,i}=p_{2,j}\\
+            &b_{i,j}\mbox{ is a binary variable} 
+
+        Sage uses MixedIntegerLinearProgram to solve this.
+
+        (Ryser [Ryser63]_): The construction of an `m\times n` matrix `A=A_{r,s}`,
+        due to Ryser, is described as follows. The
+        construction works if and only if have `s\preceq r^*`.
+
+        * Construct the `m\times n` matrix `B` from `r` by 
+        defining the `i`-th row of `B` to be the 
+        vector whose first `r_i` entries are `1`, and the
+        remainder are `0`'s, `1\leq i\leq m`.
+        This maximal matrix `B` with row sum `r` and ones 
+        left justified has column sum `r^{*}`.
+
+        * Shift the last `1` in certain rows of `B` to
+        column `n` in order to achieve the sum `s_n`.
+        Call this `B` again.
+
+          * The `1`'s in column n are to appear in those 
+            rows in which `A` has the largest row sums, giving 
+            preference to the bottom-most positions in case of ties.
+          * Note: When this step automatically "fixes" other columns,
+            one must skip ahead to the first column index
+            with a wrong sum in the step below.
+
+        * Proceed inductively to construct columns `n-1`, ..., `2`, `1`.
+
+        * Set `A = B`. Return `A`.
+
+        EXAMPLES:
+
+        Computing the matrix for `p_1=p_2=2+2+1` ::
+
+            sage: from sage.combinat.integer_vector import gale_ryser_theorem
+            sage: p1 = [2,2,1]
+            sage: p2 = [2,2,1]
+            sage: print gale_ryser_theorem(p1, p2, algorithm="gale") # Optional - requires GLPK or CBC
+            [0 1 1]
+            [1 1 0]
+            [1 0 0]       
+
+        Or for a non-square matrix with `p_1=3+3+2+1` and `p_2=3+2+2+1+1` ::
+
+            sage: p1 = [3,3,2,1]
+            sage: p2 = [3,2,2,1,1]
+            sage: print gale_ryser_theorem(p1, p2, algorithm="gale") # Optional - requires GLPK or CBC
+            [1 0 1 1 0]
+            [1 0 1 0 1]
+            [1 1 0 0 0]
+            [0 1 0 0 0]
+            sage: p1 = [3,3,1,1]
+            sage: p2 = [3,3,1,1]
+            sage: gale_ryser_theorem(p1, p2)
+            [1 1 1 0]
+            [1 1 0 1]
+            [1 0 0 0]
+            [0 1 0 0]
+            sage: p1 = [4,2,2]
+            sage: p2 = [3,3,1,1]
+            sage: gale_ryser_theorem(p1, p2)
+            [1 1 1 1]
+            [1 1 0 0]
+            [1 1 0 0]        
+            sage: p1 = [4,2,2,0]
+            sage: p2 = [3,3,1,1,0,0]
+            sage: gale_ryser_theorem(p1, p2)
+            [1 1 1 1 0 0]
+            [1 1 0 0 0 0]
+            [1 1 0 0 0 0]
+            [0 0 0 0 0 0]
+
+        REFERENCES:
+
+        ..  [Ryser63] H. J. Ryser, Combinatorial Mathematics, 
+                Carus Monographs, MAA, 1963. 
+        ..  [Gale57] D. Gale, A theorem on flows in networks, Pacific J. Math.
+                7(1957)1073-1082.
+        """
+        from sage.combinat.partition import Partition
+        from sage.matrix.constructor import matrix
+        if sum(p1) != sum(p2):
+            raise ValueError("The two lists must sum to the same value.")
+        if not Partition(p2).conjugate().dominates(Partition(p1)):
+            raise ValueError("The domination criterion is not fulfilled.")
+        if algorithm == "gale":
+          from sage.numerical.mip import MixedIntegerLinearProgram
+          p1 = Partition(p1)
+          p2 = Partition(p2)
+          k1=p1.length()
+          k2=p2.length()
+          p = MixedIntegerLinearProgram()
+          b = p.new_variable(dim=2)
+          for (i,c) in enumerate(p1):
+              p.add_constraint(sum([b[i][j] for j in xrange(k2)]),min=c,max=c)
+          for (i,c) in enumerate(p2):
+              p.add_constraint(sum([b[j][i] for j in xrange(k1)]),min=c,max=c)
+          p.set_objective(None)
+          p.set_binary(b)
+          p.solve()
+          b = p.get_values(b)
+          M = [[0]*k2 for i in xrange(k1)]
+          for i in xrange(k1):
+              for j in xrange(k2):
+                  M[i][j] = int(b[i][j])
+          return matrix(M)
+        else: # ryser's algorithm
+            r = p1; s = p2
+            if not(is_gale_ryser(r,s)):
+                return False
+            rstar = Partition(r).conjugate()
+            m = len(r)
+            n = len(s)
+            rowsA0 = [[0]*n]*m
+            for j in range(m):
+                if j<m:
+                    rowsA0[j] = [1]*r[j]+[0]*(n-r[j])
+                else:
+                    rowsA0[j] = [0]*n
+            A0 = matrix(rowsA0) 
+            for j in range(1,n-1):
+                # starts for loop, k = n-1, ..., 1
+                # which finds the 1st column with 
+                # incorrect column sum. For that bad
+                # column index, apply slider again
+                for k in range(1,n):
+                    if sum(A0.column(n-k))<>s[n-k]:
+                         break
+                A0 = _slider01(A0,s[n-k],n-k)
+            return A0    
+
+
 def _default_function(l, default, i):
     """
     EXAMPLES::
@@ -214 +495 @@
         """
         return infinity
 
+
 class IntegerVectors_nk(CombinatorialClass):
     def __init__(self, n, k):
         """