Ticket #7301: trac_7301.2.patch

File trac_7301.2.patch, 10.0 KB (added by wdj, 13 years ago)

apply this patch only to 4.3.

  • sage/combinat/integer_vector.py

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