# 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 """ Integer vectors AUTHORS: *   Mike Hanson (2007) - original module *   Nathann Cohen, David Joyner (2009-2010) - Gale-Ryser stuff """ #***************************************************************************** #       Copyright (C) 2007 Mike Hansen , 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 js[n-k]: break A0 = _slider01(A0,s[n-k],n-k) return A0 def _default_function(l, default, i): """ EXAMPLES:: """ return infinity class IntegerVectors_nk(CombinatorialClass): def __init__(self, n, k): """