| 34 | def 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 | |
| 79 | def _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 | |
| 145 | def 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 | |