#***************************************************************************** # Copyright (C) 2008 Carlo Hamalainen , # # Distributed under the terms of the GNU General Public License (GPL) # # This code is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # The full text of the GPL is available at: # # http://www.gnu.org/licenses/ #***************************************************************************** import sets import copy load "dfs_latin.spyx" load "dancing_links.spyx" load "dancing_links.sage" def tau123(T1, T2): assert is_bitrade(T1, T2) cells_map = filled_cells_map(T1) t1 = tau1(T1, T2, cells_map) t2 = tau2(T1, T2, cells_map) t3 = tau3(T1, T2, cells_map) return (cells_map, t1, t2, t3) def print_cycles(t): for x in t.to_cycles(): print x def apply_isotopism(P, row_perm, col_perm, sym_perm): """ EXAMPLE: sage: B = back_circulant(5) sage: alpha = isotopism((0,1,2,3,4)) sage: beta = isotopism((1,0,2,3,4)) sage: gamma = isotopism((2,1,0,3,4)) sage: print apply_isotopism(B, alpha, beta, gamma) [3 4 2 0 1] [0 2 3 1 4] [1 3 0 4 2] [4 0 1 2 3] [2 1 4 3 0] """ Q = matrix(ZZ, P.nrows(), P.ncols()) for r in range(P.nrows()): for c in range(P.ncols()): #Q[r, c] = sym_perm[P[row_perm[r]-1, col_perm[c]-1]] - 1 try: s2 = sym_perm[P[r, c]] - 1 except IndexError: s2 = P[r, c] # we must be leaving the symbol fixed? #Q[row_perm[r]-1, col_perm[c]-1] = sym_perm[P[r, c]] - 1 Q[row_perm[r]-1, col_perm[c]-1] = s2 return Q #def identity_permutation(n): #return Permutation([1..n]) def isotopism(p): """ Convert a list of tuples or whatever fixme to a Permutation object. """ # Identity isotopism on p points: if type(p) == Integer: return Permutation(range(1, p+1)) if type(p) == sage.groups.perm_gps.permgroup_element.PermutationGroupElement: # fixme Ask the Sage mailing list about the tuple/list issue! return Permutation(list(p.tuple())) if type(p) == list: # We expect a list like [0,3,2,1] which means # that 0 goes to 0, 1 goes to 3, etc. return Permutation(map(lambda x: x+1, p)) if type(p) == tuple: # We have a single cycle: if type(p[0]) == Integer: return Permutation(tuple(map(lambda x: x+1, p))) # We have a tuple of cycles: if type(p[0]) == tuple: x = isotopism(p[0]) for i in range(1, len(p)): x = x * isotopism(p[i]) return x # Not sure what we got! raise NotImplemented def filled_cells_map(P): """ fixme """ cells_map = {} k = 1 for r in range(P.nrows()): for c in range(P.ncols()): e = P[r, c] if e < 0: continue cells_map[ (r,c) ] = k cells_map[k] = (r,c) k += 1 return cells_map def cells_map_as_square(cells_map, n): assert n > 1 L = matrix(ZZ, n, n) for r in range(n): for c in range(n): try: L[r, c] = cells_map[ (r,c) ] except KeyError: # There is no cell (r,c) so skip it L[r, c] = -1 return L def beta1(rce, P, Q): """ fixme """ r = rce[0] c = rce[1] e = rce[2] assert P[r, c] == e assert e >= 0 for x in range(P.nrows()): if Q[x, c] == e: return (x, c, e) raise PairNotBitrade def beta2(rce, P, Q): """ fixme """ r = rce[0] c = rce[1] e = rce[2] assert P[r, c] == e assert e >= 0 for x in range(P.ncols()): if Q[r, x] == e: return (r, x, e) raise PairNotBitrade def beta3(rce, P, Q): """ fixme """ r = rce[0] c = rce[1] e = rce[2] assert P[r, c] == e assert e >= 0 # fixme this range could be iffy if we # work with latin bitrade rectangles... for x in range(P.nrows()): if Q[r, c] == x: return (r, c, x) raise PairNotBitrade def test_beta(): (a, b, c, G) = alternating_group(1) (T1, T2) = bitrade_pair(a, b, c, G) assert is_bitrade(T1, T2) print T1 print "" print T2 for r in range(T1.nrows()): for c in range(T1.ncols()): e = T1[r, c] if e < 0: continue (r2, c2, e2) = beta2( (r,c,e), T1, T2) print (r, c, e), "|-->", (r2, c2, e2) def tau1(T1, T2, cells_map): """ The definition of tau1 is: tau1 : T1 -> T1 tau1 = beta^(-1)_2 beta_3 (composing left to right) where beta_i : T2 -> T1 changes just the i-th coordinate of a triple. """ # The cells_map has both directions, i.e. integer to # cell and cell to integer, so the size of T1 is # just half of len(cells_map). x = (int(len(cells_map)/2) + 1) * [-1] for r in range(T1.nrows()): for c in range(T1.ncols()): e = T1[r, c] if e < 0: continue (r2, c2, e2) = beta2( (r,c,e), T1, T2) (r3, c3, e3) = beta3( (r2,c2,e2), T2, T1) x[ cells_map[(r,c)] ] = cells_map[ (r3,c3) ] x.pop(0) # remove the head of the list since we # have permutations on 1..(something). return Permutation(x) def tau2(T1, T2, cells_map): """ The definition of tau1 is: tau1 : T1 -> T1 tau1 = beta^(-1)_2 beta_3 (composing left to right) where beta_i : T2 -> T1 changes just the i-th coordinate of a triple. """ # The cells_map has both directions, i.e. integer to # cell and cell to integer, so the size of T1 is # just half of len(cells_map). x = (int(len(cells_map)/2) + 1) * [-1] for r in range(T1.nrows()): for c in range(T1.ncols()): e = T1[r, c] if e < 0: continue (r2, c2, e2) = beta3( (r,c,e), T1, T2) (r3, c3, e3) = beta1( (r2,c2,e2), T2, T1) x[ cells_map[(r,c)] ] = cells_map[ (r3,c3) ] x.pop(0) # remove the head of the list since we # have permutations on 1..(something). return Permutation(x) def tau3(T1, T2, cells_map): """ The definition of tau1 is: tau1 : T1 -> T1 tau1 = beta^(-1)_2 beta_3 (composing left to right) where beta_i : T2 -> T1 changes just the i-th coordinate of a triple. """ # The cells_map has both directions, i.e. integer to # cell and cell to integer, so the size of T1 is # just half of len(cells_map). x = (int(len(cells_map)/2) + 1) * [-1] for r in range(T1.nrows()): for c in range(T1.ncols()): e = T1[r, c] if e < 0: continue (r2, c2, e2) = beta1( (r,c,e), T1, T2) (r3, c3, e3) = beta2( (r2,c2,e2), T2, T1) x[ cells_map[(r,c)] ] = cells_map[ (r3,c3) ] x.pop(0) # remove the head of the list since we # have permutations on 1..(something). return Permutation(x) def back_circulant(n): """ The back-circulant latin square of order n is the Cayley table for (Z_n, +), the integers under addition modulo n. INPUT: n -- int; order of the latin square. EXAMPLE: sage: print back_circulant(5) [0 1 2 3 4] [1 2 3 4 0] [2 3 4 0 1] [3 4 0 1 2] [4 0 1 2 3] """ assert n >= 1 L = matrix(ZZ, n, n) for r in range(n): for c in range(n): L[r, c] = (r + c) % n return L def forward_circulant(n): """ fixme """ assert n >= 1 L = matrix(ZZ, n, n) for r in range(n): for c in range(n): L[r, c] = (n-c+r) % n return L def direct_product(L1, L2, L3, L4): """ ----------- | L1 | L2 | ----------- | L3 | L4 | ----------- """ assert L1.nrows() == L2.nrows() == L3.nrows() == L4.nrows() assert L1.ncols() == L2.ncols() == L3.ncols() == L4.ncols() assert L1.nrows() == L1.ncols() n = L1.nrows() D = matrix(ZZ, 2*n, 2*n) for r in range(n): for c in range(n): D[r, c] = L1[r, c] D[r, c+n] = L2[r, c] + n D[r+n, c] = L3[r, c] + n D[r+n, c+n] = L4[r, c] return D def abelian_2group(s): """ Returns the latin square based on the Cayley table for the abelian 2-group of order 2^s. INPUT: s -- int; order of the latin square will be 2^s. EXAMPLE: sage: print abelian_2group(3) [0 1 2 3 4 5 6 7] [1 0 3 2 5 4 7 6] [2 3 0 1 6 7 4 5] [3 2 1 0 7 6 5 4] [4 5 6 7 0 1 2 3] [5 4 7 6 1 0 3 2] [6 7 4 5 2 3 0 1] [7 6 5 4 3 2 1 0] """ assert s > 0 if s == 1: L = matrix(ZZ, 2, 2) L[0, 0] = 0 L[0, 1] = 1 L[1, 0] = 1 L[1, 1] = 0 return L else: L_prev = abelian_2group(s-1) L = matrix(ZZ, 2^s, 2^s) offset = L.nrows()/2 for r in range(L_prev.nrows()): for c in range(L_prev.ncols()): L[r, c] = L_prev[r, c] L[r+offset, c] = L_prev[r, c] + offset L[r, c+offset] = L_prev[r, c] + offset L[r+offset, c+offset] = L_prev[r, c] return L def empty_latin_square(n): """ Returns an empty latin square of order n. Note that L[r, c] is an empty cell if and only if L[r, c] = -1. INPUT: n -- int; order of the latin square. EXAMPLE: sage: print empty_latin_square(5) [-1 -1 -1 -1 -1] [-1 -1 -1 -1 -1] [-1 -1 -1 -1 -1] [-1 -1 -1 -1 -1] [-1 -1 -1 -1 -1] """ assert n >= 1 L = matrix(ZZ, n, n) # fixme: there are so many nested loops like this. Is there a nicer # way to iterate over the cells of a matrix? for r in range(L.nrows()): for c in range(L.ncols()): L[r, c] = -1 return L def top_left_empty_cell(L): """ Returns the least [r, c] such that L[r, c] is an empty cell. INPUT: L -- matrix; a latin square. EXAMPLE: sage: B = back_circulant(5) sage: B[3, 4] = -1 sage: print top_left_empty_cell(B) [3, 4] """ for r in range(L.nrows()): for c in range(L.ncols()): if L[r, c] < 0: return [r, c] return None def permissable_values(L, r, c): """ Find all values that do not appear in row r and column c of the latin square L. If L[r, c] is filled then we return the empty list. INPUT: L -- matrix; a latin square r -- int; row of the latin square c -- int; column of the latin square EXAMPLE: sage: L = back_circulant(5) sage: L[0, 0] = -1 sage: print permissable_values(L, 0, 0) [0] """ if L[r, c] >= 0: return [] assert L.nrows() == L.ncols() n = L.nrows() vals = {} for e in range(n): vals[e] = True for i in range(n): if L[i, c] >= 0: del vals[ L[i, c] ] for j in range(n): if L[r, j] >= 0: try: del vals[ L[r, j] ] except KeyError: # We may have already removed a symbol # in the previous for-loop. pass return vals.keys() def random_empty_cell(P): """ Find an empty cell of P, uniformly at random. INPUT: P -- matrix; a partial latin square RETURNS: [r, c] -- cell such that P[r, c] is empty, or returns None if P is a latin square. """ cells = {} for r in range(P.nrows()): for c in range(P.ncols()): if P[r, c] < 0: cells[ (r,c) ] = True cells = cells.keys() if len(cells) == 0: return None rc = cells[ ZZ.random_element(len(cells)) ] return [rc[0], rc[1]] def dfs_py(P, random_search = False, nr_to_find = 0): """ Depth first search to find a completion of the partial latin square P. INPUT: P -- matrix; partial latin square nr_to_find -- int; optional parameter for the maximum number of completions to find. RETURNS: A list of latin squares such that P is contained in each such square. EXAMPLE: sage: P = empty_latin_square(3) sage: P[0, 0] = 0 sage: all_comps = dfs_py(P) sage: some_comps = dfs_py(P, 2) NOTES: There are much quicker ways of finding partial latin square completions (e.g. Knuth's "Dancing Links" algorithm for solving 0-1 matrix exact cover) so this function should only be used for small examples or testing! """ completions = {} def search(P): if random_search: rc = random_empty_cell(P) else: rc = top_left_empty_cell(P) if rc == None: P.set_immutable() completions[P] = True return r = rc[0] c = rc[1] for e in permissable_values(P, r, c): P2 = copy.deepcopy(P) P2[r, c] = e search(P2) # Make sure that we don't generate too many solutions. if nr_to_find > 0 and len(completions) >= nr_to_find: return search(P) return completions.keys() def is_partial_latin_square(P): """ P is a partial latin square if it is an n by n matrix, and each symbol in [0, 1, ..., n-1] appears at most once in each row, and at most once in each column. """ assert P.nrows() == P.ncols() n = P.nrows() for r in range(n): vals_in_row = {} for c in range(n): e = P[r, c] # Entry out of range 0, 1, ..., n-1: if e >= n: return False # Entry has already appeared in this row: if vals_in_row.has_key(e): return False vals_in_row[e] = True for c in range(n): vals_in_col = {} for r in range(n): e = P[r, c] # Entry out of range 0, 1, ..., n-1: if e >= n: return False # Entry has already appeared in this column: if vals_in_col.has_key(e): return False vals_in_col[e] = True return True def is_latin_square(L): """ L is a partial latin square if it is an n by n matrix, and each symbol in [0, 1, ..., n-1] appears exactly once in each row, and exactly once in each column. """ # We don't allow latin rectangles: if L.nrows() != L.ncols(): return False # Every cell must be filled: if len(filter(lambda x: x >= 0, L.list())) != L.nrows()*L.ncols(): return False # By necessity L must be a partial latin square: if not is_partial_latin_square(L): return False return True def is_uniquely_completable(P): """ Returns True if the partial latin square P has exactly one completion to a latin square. """ return len(dfs_pyrex(P, random_search = True, nr_to_find = 2)) == 1 def is_completable(P): """ Returns True if the partial latin square can be completed to a latin square. EXAMPLE: The following partial latin square has no completion because there is nowhere that we can place the symbol 0 in the third row: sage: B = empty_latin_square(3) sage: B[0, 0] = 0 sage: B[1, 1] = 0 sage: B[2, 2] = 1 sage: print B [ 0 -1 -1] [-1 0 -1] [-1 -1 1] sage: print is_completable(B) False """ return len(dfs_pyrex(P, random_search = True, nr_to_find = 1)) > 0 def gcs(L, has_unique_completion, random_search = False): """ A greedy critical set of a latin square L is found by successively removing elements in a row-wise (bottom-up) manner, checking for unique completion at each step. """ n = L.nrows() G = L.copy() for r in range(n-1, -1, -1): for c in range(n-1, -1, -1): e = G[r, c] G[r, c] = -1 #if not is_uniquely_completable(G): #dlxcpp_has_unique_completion #if not dlxcpp_has_unique_completion(G, random_search = False): if not has_unique_completion(G, random_search = False): G[r, c] = e return G def coin(): """ Simulates a fair coin (returns True or False) using ZZ.random_element(2). """ return ZZ.random_element(2) == 0 def row_containing_sym(L, c, x): """ Given an improper latin square L with L[r1, c] = L[r2, c] = x, return r1 or r2 with equal probability. This is an internal function and should only be used in latin_square_generator(). """ r1 = -1 r2 = -1 for r in range(L.nrows()): if r1 >= 0 and r2 >= 0: break if L[r, c] == x and r1 < 0: r1 = r continue if L[r, c] == x and r2 < 0: r2 = r break assert r1 >= 0 and r2 >= 0 if coin(): return r1 else: return r2 def column_containing_sym(L, r, x): """ Given an improper latin square L with L[r, c1] = L[r, c2] = x, return c1 or c2 with equal probability. This is an internal function and should only be used in latin_square_generator(). """ c1 = -1 c2 = -1 for c in range(L.ncols()): if c1 >= 0 and c2 >= 0: break if L[r, c] == x and c1 < 0: c1 = c continue if L[r, c] == x and c2 < 0: c2 = c break assert c1 >= 0 and c2 >= 0 if coin(): return c1 else: return c2 def next_conjugate(L): """ Permute L[r, c] = e to the conjugate L[c, e] = r We assume that L is an n by n matrix and has values in the range 0, 1, ..., n-1. """ assert L.nrows() == L.ncols() n = L.nrows() C = matrix(ZZ, n, n) for r in range(n): for c in range(n): e = L[r, c] assert e >= 0 and e < n C[c, e] = r return C def latin_square_generator(L_start, check_assertions = False): """ Generator for a sequence of uniformly distributed latin squares, given L_start as the initial latin square. This code implements the Markov chain algorithm of Jacobson and Matthews (1996), see below for the BibTex entry. This generator will never throw the StopIteration exception, so it provides an infinite sequence of latin squares. EXAMPLE: Use the back circulant latin square of order 4 as the initial square and print the next two latin squares given by the Markov chain: sage: g = latin_square_generator(back_circulant(4)) sage: print g.next() sage: print g.next() REFERENCE: @article{MR1410617, AUTHOR = {Jacobson, Mark T. and Matthews, Peter}, TITLE = {Generating uniformly distributed random {L}atin squares}, JOURNAL = {J. Combin. Des.}, FJOURNAL = {Journal of Combinatorial Designs}, VOLUME = {4}, YEAR = {1996}, NUMBER = {6}, PAGES = {405--437}, ISSN = {1063-8539}, MRCLASS = {05B15 (60J10)}, MRNUMBER = {MR1410617 (98b:05021)}, MRREVIEWER = {Lars D{\o}vling Andersen}, } """ if check_assertions: assert is_latin_square(L_start) n = L_start.nrows() r1 = r2 = c1 = c2 = x = y = z = -1 proper = True L = L_start.copy() L_rce = L L_cer = matrix(ZZ, n, n) L_erc = matrix(ZZ, n, n) while True: if proper: if check_assertions: assert is_latin_square(L) ################################# # Update the other two conjugates for r in range(n): for c in range(n): e = L[r, c] L_cer[c, e] = r L_erc[e, r] = c ################################# yield L r1 = ZZ.random_element(n) c1 = ZZ.random_element(n) x = L[r1, c1] y = x while y == x: y = ZZ.random_element(n) # Now find y in row r1 and column c1. if check_assertions: r2 = 0 c2 = 0 while L[r1, c2] != y: c2 += 1 while L[r2, c1] != y: r2 += 1 assert L_erc[y, r1] == c2 assert L_cer[c1, y] == r2 c2 = L_erc[y, r1] r2 = L_cer[c1, y] if check_assertions: assert L[r1, c2] == y if check_assertions: assert L[r2, c1] == y L[r1, c1] = y L[r1, c2] = x L[r2, c1] = x # Now deal with the unknown point. # We want to form z + (y - x) z = L[r2, c2] if z == x: L[r2, c2] = y else: # z and y have positive coefficients # x is the improper term with a negative coefficient proper = False else: # improper square, # L[r2, c2] = y + z - x # y and z are proper while x is the # improper symbol in the cell L[r2, c2]. r1 = row_containing_sym(L, c2, x) c1 = column_containing_sym(L, r2, x) if check_assertions: assert L[r1, c2] == x if check_assertions: assert L[r2, c1] == x # choose one of the proper symbols # uniformly at random (we will use whatever # lands in variable y). if coin(): y, z = z, y # Add/subtract the symbolic difference (y - x) L[r2, c2] = z L[r1, c2] = y L[r2, c1] = y if L[r1, c1] == y: L[r1, c1] = x proper = True else: # got another improper square z = L[r1, c1] x, y = y, x r2 = r1 c2 = c1 # Now we have L[r2, c2] = z+y-x as # usual proper = False # for emphasis def cayley_table(G): """ Construct a latin square on the symbols [0, 1, ..., n-1] for a group with an n by n Cayley table. EXAMPLES: sage: print cayley_table(DihedralGroup(2)) [0 1 2 3] [1 0 3 2] [2 3 0 1] [3 2 1 0] sage: G = gap.Group(PermutationGroupElement((1,2,3))) sage: print cayley_table(G) [0 1 2] [1 2 0] [2 0 1] """ if isinstance(G, sage.interfaces.gap.GapElement): rows = map(lambda x: list(x), list(gap.MultiplicationTable(G))) new_rows = [] for x in rows: new_rows.append(map(lambda x: int(x)-1, x)) return matrix(new_rows) # Otherwise we must have some kind of Sage permutation group object, # such as sage.groups.perm_gps.permgroup.PermutationGroup_generic # or maybe sage.groups.perm_gps.permgroup_named. # We should have a cayley_table() function. try: C = G.cayley_table() except: # We got some other kind of group object? raise NotImplemented L = matrix(ZZ, C.nrows(), C.ncols()) k = 0 entries = {} for r in range(C.nrows()): for c in range(C.ncols()): e_table = C[r, c] if entries.has_key(e_table): L[r, c] = entries[e_table] else: entries[e_table] = k L[r, c] = k k += 1 return L def count_frequencies(n): freq = {} gen = latin_square_generator(back_circulant(n), check_assertions = True) nr_squares = 30000 for _ in range(nr_squares): L = gen.next() M = L.copy() M.set_immutable() try: freq[M] += 1 except KeyError: freq[M] = 1 # fixme: surely there's a built in function to calculate # a frequency distribution in Sage? freq_list = map(lambda x: x/float(nr_squares), freq.values()) return freq_list def alternating_group(m): assert m >= 1 a = tuple(range(1, 2*m+1 + 1)) b = tuple(range(m+1, 0, -1) + range(2*m+2, 3*m+1 + 1)) a = PermutationGroupElement(a) b = PermutationGroupElement(b) c = PermutationGroupElement((a*b)^(-1)) G = PermutationGroup([a, b]) return (a, b, c, G) def pq_group(p, q): assert is_prime(p) assert is_prime(q) assert (q % p) == 1 # beta is a primitive root of the # congruence x^p = 1 mod q F = FiniteField(q) fgen = F.multiplicative_generator() beta = fgen^((q-1)/p) assert beta != 1 assert (beta^p % q) == 1 Q = tuple(range(1, q+1)) P = [] seenValues = {} for i in range(2, q): if seenValues.has_key(i): continue cycle = [] for k in range(p): x = (1 + (i-1)*beta^k) % q if x == 0: x = q seenValues[x] = True cycle.append(x) P.append(tuple(cycle)) G = PermutationGroup([P, Q]) assert G.order() == p*q assert not G.is_abelian() a = PermutationGroupElement(P) b = PermutationGroupElement(Q) c = PermutationGroupElement((a*b)^(-1)) return (a, b, c, PermutationGroup([P, Q])) def p3_group(p): assert is_prime(p) F = gap.new("FreeGroup(3)") a = F.gen(1) b = F.gen(2) c = F.gen(3) rels = [] rels.append( a^p ) rels.append( b^p ) rels.append( c^p ) rels.append( a*b*((b*a*c)^(-1)) ) rels.append( c*a*((a*c)^(-1)) ) rels.append( c*b*((b*c)^(-1)) ) G = F.FactorGroupFpGroupByRels(rels) N = gap.NormalClosure(F, gap.Subgroup(F, rels)) niso = gap.NaturalHomomorphismByNormalSubgroupNC(F, N) a = PermutationGroupElement(gap.Image(niso, a)) b = PermutationGroupElement(gap.Image(niso, b)) c = PermutationGroupElement(gap.Image(niso, c)) G = PermutationGroup([a, b, c]) c = PermutationGroupElement((a*b)^(-1)) return (a, b, c, G) def check_bitrade_properties(a, b, c, G): A = PermutationGroup([a]) B = PermutationGroup([b]) C = PermutationGroup([c]) assert a*b == c^(-1) #assert AB.Size() == AC.Size() == BC.Size() == 1 X = gap.Intersection(gap.Intersection(A, B), C) assert X.Size() == 1 def is_bitrade(T1, T2): if not is_disjoint(T1, T2): return False if not is_same_shape(T1, T2): return False if not is_row_and_col_balanced(T1, T2): return False return True def is_primary_bitrade(a, b, c, G): H = PermutationGroup([a, b, c]) return G.order() == H.order() def coset_representatives(x, G): rcosets = gap.RightCosets(G, PermutationGroup([x])) cosetReps = [] for i in range(1, len(rcosets)+1): cosetReps.append(gap.Representative(rcosets[i])) assert gap.Index(G, PermutationGroup([x])) == len(cosetReps) return cosetReps # The rows of the bitrade are labeled using # the cosets of in G. def bitrade_labels(a, b, c, G): rLabels = coset_representatives(a, G) cLabels = coset_representatives(b, G) eLabels = coset_representatives(c, G) return (rLabels, cLabels, eLabels) # There should be a nicer way to do this, surely? def entry_label(eLabels, C, x): """ Arguments: eLabels --- list of coset representatives of C in G C --- subgroup G of G x --- element of G Returns: integer i such that C*eLabels[i] == C*x """ Cx = gap.RightCoset(C, x) for i in range(len(eLabels)): y = eLabels[i] Cy = gap.RightCoset(C, y) if Cx == Cy: return i assert False def bitrade_pair(a, b, c, G): Tleft = {} Tright = {} (rLabels, cLabels, eLabels) = bitrade_labels(a, b, c, G) #print rLabels #print cLabels #print eLabels T1 = matrix(ZZ, len(rLabels), len(cLabels)) T2 = matrix(ZZ, len(rLabels), len(cLabels)) for i in range(len(rLabels)): for j in range(len(cLabels)): T1[i, j] = -1 T2[i, j] = -1 A = gap.Group([a]) B = gap.Group([b]) C = gap.Group([c]) for i in range(len(rLabels)): for j in range(len(cLabels)): for k in range(len(eLabels)): cosetRow = gap.RightCoset(A, rLabels[i]) cosetCol = gap.RightCoset(B, cLabels[j]) cosetEnt = gap.RightCoset(C, eLabels[k]) #print "---------------------" #print cosetRow #print cosetCol #print cosetEnt g = gap.Intersection(gap.Elements(cosetRow), gap.Elements(cosetCol), gap.Elements(cosetEnt)) #if len(gap.Intersection(cosetRow, cosetCol)) == 0: continue #if len(gap.Intersection(cosetRow, cosetEnt)) == 0: continue #if len(gap.Intersection(cosetCol, cosetEnt)) == 0: continue # fixme: why did the previous 3 lines not work? #print g, i, j, k if len(g) == 0: continue # At this point we know that a, b, c are the right # generators for a group-based bitrade so we can get # the *unique* element g in the intersection of # cosetRow, cosetCol, and cosetEnt by getting the # intersection of just one pair. g = gap.Intersection(cosetRow, cosetCol) assert len(g) == 1 g = g[1] #eIdx = entry_label(eLabels, C, a^(-1) * g) eIdx = entry_label(eLabels, C, a * g) # fixme we're doing things with right instead of # left cosets now, arrrgh! #Tleft[(rLabels[i], cLabels[j])] = eLabels[k] #Tright[(rLabels[i], cLabels[j])] = eLabels[eIdx] T1[i, j] = k T2[i, j] = eIdx #return (Tleft, Tright) return (T1, T2) def is_disjoint(Tleft, Tright): assert type(Tleft) == type(Tright) if is_Matrix(Tleft): return is_disjoint_matrix(Tleft, Tright) if isinstance(Tleft, dict): return is_disjoint_dictionary(Tleft, Tright) def is_disjoint_dictionary(Tleft, Tright): #print Tleft.keys() #print Tright.keys() #print Tleft.keys() == Tright.keys() assert Tleft.keys() == Tright.keys() for ij in Tleft.keys(): if Tleft[ij] == Tright[ij]: return False return True def is_disjoint_matrix(Tleft, Tright): for i in range(Tleft.nrows()): for j in range(Tleft.ncols()): if Tleft[i, j] < 0 and Tright[i, j] < 0: continue if Tleft[i, j] == Tright[i, j]: return False return True def is_same_shape(T1, T2): for i in range(T1.nrows()): for j in range(T1.ncols()): if T1[i, j] < 0 and T2[i, j] < 0: continue if T1[i, j] >= 0 and T2[i, j] >= 0: continue return False return True def is_row_and_col_balanced(T1, T2): for r in range(T1.nrows()): val1 = sets.Set(filter(lambda x: x >= 0, T1.row(r))) val2 = sets.Set(filter(lambda x: x >= 0, T2.row(r))) if val1 != val2: return False for c in range(T1.ncols()): val1 = sets.Set(filter(lambda x: x >= 0, T1.column(c))) val2 = sets.Set(filter(lambda x: x >= 0, T2.column(c))) if val1 != val2: return False return True # examples def example1(): B = back_circulant(6) print B print "Top left entry:", B[0, 0] print "Top-left-most empty cell:", top_left_empty_cell(B) B[3, 4] = -1 print B print "Top-left-most empty cell:", top_left_empty_cell(B) def example2(): A = abelian_2group(3) # abelian 2-group of order 2^3 print A G = gcs(A) print G print "is_uniquely_completable(G) =", is_uniquely_completable(G) print "" G[0, 0] = -1 print G print "is_uniquely_completable(G) =", is_uniquely_completable(G) def example3(): G = gcs(back_circulant(3)) G[0, 0] = -1 print G print "" comps = dfs_pyrex(G, random_search = True, nr_to_find = 2) for L in comps: print "Completion of G:" print L print "" def example4(): B = back_circulant(4) B[0, 1] = 0 B[1, 0] = -1 print B print "is_partial_latin_square(B) =", is_partial_latin_square(B) def example5(): B = empty_latin_square(3) B[0, 0] = 0 B[1, 1] = 0 B[2, 2] = 1 print B print "is_completable(B) =", is_completable(B) print "" B[2, 2] = 0 print B print "is_completable(B) =", is_completable(B) def example6(): L = abelian_2group(3) print L print "Next conjugate:" print next_conjugate(L) # If we do the conjugate three times we should get back # to the original latin square. print "Conjugate^3(L) == L:", next_conjugate(next_conjugate(next_conjugate(L))) == L def example7(): g = latin_square_generator(back_circulant(4)) print g.next() print "" print g.next() print "" def example8(): P = empty_latin_square(8) # Time some non-randomised searches: for nr_to_find in range(1, 1000, 50): print "aaa", nr_to_find, time dfs_pyrex(P, False, nr_to_find) print "bbb", nr_to_find, time dfs_py(P, False, nr_to_find) # example8() # example7() def is_maximal_subgroup(G, H): """ Returns True if H is a maximal subgroup of G. """ pts = G.orbits() # A group should only have one orbit. assert len(pts) == 1 pts = pts[0] return gap.IsPrimitive(H, pts) def is_maximal_cyclic_subgroup(G, H): """ Is H a maximal (cyclic) subgroup of G? """ assert gap.IsCyclic(H) #gens = list(H.gens()) gens = list(gap.GeneratorsOfGroup(H)) # fixme Brute force, must be a better way? for g in gap.Elements(G): if g in H: pass K = gap.Group(list(gens) + [g]) if not gap.IsCyclic(K): continue # We need K and H to be of the same order if gap.Order(K) != gap.Order(H): print gens, g return False return True def gap_to_sage(G): gens = map(lambda x: PermutationGroupElement(x), list(gap.GeneratorsOfGroup(G))) return PermutationGroup(gens) def sage_to_gap(G): return gap.Group(G.gens()) if False: (a, b, c, G) = p3_group(11) #(a, b, c, G) = alternating_group(1) A = PermutationGroup([a]) B = PermutationGroup([b]) C = PermutationGroup([c]) # print is_maximal_subgroup(G, gap.Subgroup(G, [a])) # print is_maximal_subgroup(G, gap.Subgroup(G, [b])) # print is_maximal_subgroup(G, gap.Subgroup(G, [c])) #print is_maximal_cyclic_subgroup(G, gap.Subgroup(G, [a])) #print is_maximal_cyclic_subgroup(G, gap.Subgroup(G, [b])) #print is_maximal_cyclic_subgroup(G, gap.Subgroup(G, [c])) #[(1,2,3,4,5,6,7)] ( 8, 9,10) #g1 = PermutationGroupElement((1,2,3,4,5,6,7)) #g2 = PermutationGroupElement( (8,9,10) ) #G1 = gap.Group( g1 ) #G2 = gap.Group( g1, g2 ) partition = {} partition[A] = True partition[B] = True partition[C] = True if False: for g in list(G): for X in [A, B, C]: H = gap_to_sage(gap.ConjugateGroup(X, g)) if partition.has_key(H): continue partition[H] = True # check something here? for H in partition.keys(): for K in partition.keys(): if H == K: continue if len(sets.Set(list(H)).intersection(sets.Set(list(K)))) > 1: print H print K assert False print "win!" if False: D_list = sets.Set(list(D)) A_list = sets.Set(list(A)) B_list = sets.Set(list(B)) C_list = sets.Set(list(C)) assert len(D_list.intersection(A_list)) == 1 assert len(D_list.intersection(B_list)) == 1 assert len(D_list.intersection(C_list)) == 1 print G.order() - (A.order() + B.order() + C.order() + D.order()) #print D_list == list(A) #print D_list == list(B) #print D_list == list(C) if False: G = sage_to_gap(G) A = sage_to_gap(A) B = sage_to_gap(B) C = sage_to_gap(C) clA = gap.ConjugacyClassSubgroups(G, A) clB = gap.ConjugacyClassSubgroups(G, B) clC = gap.ConjugacyClassSubgroups(G, C) # Check that these classes form a partition for i in range(1, clA.Size()+1): for j in range(1, clB.Size()+1): #print clA[i], clB[j] x = gap.Intersection(clA[i], clB[j]) assert x.Order() == 1 if True: (a, b, c, G) = alternating_group(1) (T1, T2) = bitrade_pair(a, b, c, G) is_bitrade(T1, T2) cells_map = filled_cells_map(T1) t1 = tau1(T1, T2, cells_map) t2 = tau2(T1, T2, cells_map) t3 = tau3(T1, T2, cells_map) #test_beta() if False: B = back_circulant(5) print B print "" x = isotopism( (0,1,2,3,4) ) id = isotopism(5) C = apply_isotopism(B, id, x, id) print C print "" for _ in range(5): B = apply_isotopism(B, id, x, id) print B def blah(L): assert L.nrows() == L.ncols() n = L.nrows() cells_map = filled_cells_map(L) g = graphs.EmptyGraph() # fixme: dummy isolated vertex with label 0 # so that we can verify the isomorphism explicitly. g.add_vertex(0) # n^2 vertices corresponding to entries g.add_vertices([1..(n^2)]) # n row vertices row_vertex = [(n^2 + 1)..(n^2 + n)] col_vertex = [(n^2 + n + 1)..(n^2 + 2*n)] sym_vertex = [(n^2 + 2*n + 1)..(n^2 + 3*n)] # Special R, C, S vertices so that we only check for # isotopic latin squares, not conjugates as well. R = n^2 + 3*n + 1 C = n^2 + 3*n + 2 S = n^2 + 3*n + 3 #print "row, col, sym etc vertex lists:" #print [1..(n^2)] print row_vertex #print col_vertex #print sym_vertex #print R, C, S #print "" assert len(row_vertex) == n assert len(col_vertex) == n assert len(sym_vertex) == n g.add_vertices(row_vertex) g.add_vertices(col_vertex) g.add_vertices(sym_vertex) #print row_vertex #print col_vertex #print sym_vertex #print R, C, S # R is connected to each row_vertex # C is connected to each col_vertex # S is connected to each sym_vertex g.add_vertex(R) g.add_vertex(C) g.add_vertex(S) for x in range(n): g.add_edge(R, row_vertex[x]) g.add_edge(C, col_vertex[x]) g.add_edge(S, sym_vertex[x]) # If (r,c,s) is in L then we create edges: # ( cells_map[(r,c)], row_vertex[r] ) # ( cells_map[(r,c)], col_vertex[c] ) # ( cells_map[(r,c)], sym_vertex[s] ) for r in range(L.nrows()): for c in range(L.ncols()): s = L[r, c] g.add_edge(cells_map[(r,c)], row_vertex[r]) g.add_edge(cells_map[(r,c)], col_vertex[c]) g.add_edge(cells_map[(r,c)], sym_vertex[s]) # k = n^2 + 3*n + 100 k = len(g.vertices()) # fixme we should know the max degree and not have # to calculate it. m = max(g.degree()) # fixme deal with this return g # Make R have degree m+1 while g.degree(R) < m+1: g.add_vertex(k) g.add_edge(R, k) k += 1 # Make C have degree m+2 while g.degree(C) < m+2: g.add_vertex(k) g.add_edge(C, k) k += 1 # Make S have degree m+2 while g.degree(S) < m+3: g.add_vertex(k) g.add_edge(S, k) k += 1 return g if False: n = 4 G = SymmetricGroup(n) #alpha = isotopism(G.random_element()) alpha = isotopism([3,1,2,0]) beta = isotopism(n) gamma = isotopism(n) print "Annoyingly, we now think of rows 1..n, etc:" print alpha #print beta #print gamma print "" L1 = back_circulant(n) L2 = apply_isotopism(L1, alpha, beta, gamma) print L1 print "" print L2 print "" g1 = blah(L1) g2 = blah(L2) print g1.is_isomorphic(g2, certify = True) # Let's construct the relabelling ourselves! cells_map = filled_cells_map(L1) C = cells_map_as_square(cells_map, n) #print C #print "" C2 = apply_isotopism(C, alpha, beta, gamma) #print C2 perm = {} for x in [1..(n^2)]: r, c = cells_map[(x)] perm[x] = C2[r, c] row_vertex = [(n^2 + 1)..(n^2 + n)] sym_vertex = [(n^2 + 2*n + 1)..(n^2 + 3*n)] for r in range(n): perm[row_vertex[r]] = row_vertex[alpha[r]-1] #print "alpha =", alpha if False: h = g1.copy() h.relabel(perm) print h.is_isomorphic(g2, certify = True) print "g1 edges:" for (a,b,_) in g1.edges(): print str(a) + ", " + str(b) print "g2 edges:" for (a,b,_) in g2.edges(): print str(a) + ", " + str(b) #print g1.edges_incident(5) #print g2.edges_incident(13) h = g1.copy() h.relabel(perm) print h.is_isomorphic(g2, certify=True) print g1.edges() print g2.edges() if False: (a, b, c, G) = alternating_group(1) #omega = G.orbits() #assert len(omega) == 1 #omega = omega[0] N = gap.Stabilizer(G, 1) (T1, T2) = bitrade_pair(a, b, c, G) ############ (cells_map, t1, t2, t3) = tau123(T1, T2) U1 = PermutationGroup([t1]) U2 = PermutationGroup([t2]) U3 = PermutationGroup([t3]) def vals_in_row(P, r): n = P.ncols() vals_in_row = {} for c in range(n): e = P[r, c] if e >= 0: vals_in_row[e] = True return vals_in_row def vals_in_col(P, c): n = P.nrows() vals_in_col = {} for r in range(n): e = P[r, c] if e >= 0: vals_in_col[e] = True return vals_in_col def dlx_find_completions(P, nr_to_find = None, random_search = False): global SOLUTIONS assert P.nrows() == P.ncols() n = P.nrows() # We will need 3n^2 columns in total: # # n^2 for the xCy columns # n^2 for the xRy columns # n^2 for the xy columns dlx_rows = [] cmap = {} k = 1 for r in range(n): valsrow = vals_in_row(P, r) for c in range(n): valscol = vals_in_col(P, c) for e in range(n): # These should be constants c_OFFSET = 1 + e + c*n r_OFFSET = 1 + e + r*n + n*n xy_OFFSET = 1 + 2*n*n + r*n + c cmap[(c_OFFSET, r_OFFSET, xy_OFFSET)] = (r,c,e) #print "possibility: ", r, c, e, "offsets:", c_OFFSET, r_OFFSET, xy_OFFSET #if P[r, c] >= 0: continue # We only want the correct value to pop in here if P[r, c] >= 0 and P[r, c] != e: continue if P[r, c] < 0 and valsrow.has_key(e): continue if P[r, c] < 0 and valscol.has_key(e): continue dlx_rows.append([k, [c_OFFSET, r_OFFSET, xy_OFFSET]]) k += 1 #print dlx_rows SOLUTIONS = {} m = dlx.DLXMatrix() m.solve(dlx_rows, callback=solution_accumulator, nr_to_find=nr_to_find, do_we_terminate = do_we_terminate, random_search = True) #print SOLUTIONS comps = [] #for i in range(len(SOLUTIONS)): for i in SOLUTIONS.keys(): soln = list(i) Q = copy.deepcopy(P) for x in soln: [k, [c_OFFSET, r_OFFSET, xy_OFFSET]] = dlx_rows[x-1] (r,c,e) = cmap[(c_OFFSET, r_OFFSET, xy_OFFSET)] if Q[r, c] >= 0: assert Q[r, c] == e else: Q[r, c] = e comps.append(Q) #print Q #print "" #SOLUTIONS[i] = Q return comps # fixme I use the following variable as a global to accumulate # the solutions using the callback function solution_accumulator. # What's the right way (i.e. no global variable) to accumulate # values somewhere? SOLUTIONS = {} def do_we_terminate(nr_to_find): global SOLUTIONS #print "do_we_terminate:", nr_to_find return len(SOLUTIONS) >= nr_to_find def solution_accumulator(x): global SOLUTIONS #print "solution_accumulator:", x # Lists are not immutable so we hash on tuples (bit of a dirty trick). SOLUTIONS[tuple(x)] = True def dlx_has_completion(P, random_search = None): return len(dlx_find_completions(P, nr_to_find = 1, random_search = random_search)) > 0 def py_has_unique_completion(P, random_search = None): return len(dfs_py(P, random_search = False, nr_to_find = 2)) == 1 def pyrex_has_unique_completion(P, random_search = None): return len(dfs_pyrex(P, random_search = False, nr_to_find = 2)) == 1 def dlx_has_unique_completion(P, random_search = None): return len(dlx_find_completions(P, nr_to_find = 2, random_search = random_search)) == 1 def dlxcpp_has_unique_completion(P, random_search = None): return len(dlxcpp_find_completions(P, nr_to_find = 2)) == 1 def dlxcpp_rows_and_map(P): assert P.nrows() == P.ncols() n = P.nrows() # We will need 3n^2 columns in total: # # n^2 for the xCy columns # n^2 for the xRy columns # n^2 for the xy columns dlx_rows = [] cmap = {} for r in range(n): valsrow = vals_in_row(P, r) for c in range(n): valscol = vals_in_col(P, c) for e in range(n): # These should be constants c_OFFSET = e + c*n r_OFFSET = e + r*n + n*n xy_OFFSET = 2*n*n + r*n + c cmap[(c_OFFSET, r_OFFSET, xy_OFFSET)] = (r,c,e) #print "possibility: ", r, c, e, "offsets:", c_OFFSET, r_OFFSET, xy_OFFSET #if P[r, c] >= 0: continue # We only want the correct value to pop in here if P[r, c] >= 0 and P[r, c] != e: continue if P[r, c] < 0 and valsrow.has_key(e): continue if P[r, c] < 0 and valscol.has_key(e): continue dlx_rows.append([c_OFFSET, r_OFFSET, xy_OFFSET]) return dlx_rows, cmap def dlxcpp_row_to_entry(x, dlx_rows, cmap): [c_OFFSET, r_OFFSET, xy_OFFSET] = dlx_rows[x] (r,c,e) = cmap[(c_OFFSET, r_OFFSET, xy_OFFSET)] return (r,c,e) def dlxcpp_find_completions(P, nr_to_find = None): assert P.nrows() == P.ncols() n = P.nrows() dlx_rows, cmap = dlxcpp_rows_and_map(P) SOLUTIONS = {} for x in DLXCPP(dlx_rows): x.sort() SOLUTIONS[tuple(x)] = True if len(SOLUTIONS) >= nr_to_find: break comps = [] for i in SOLUTIONS.keys(): soln = list(i) Q = copy.deepcopy(P) for x in soln: (r, c, e) = dlxcpp_row_to_entry(x, dlx_rows, cmap) if Q[r, c] >= 0: assert Q[r, c] == e else: Q[r, c] = e comps.append(Q) return comps """ P = empty_latin_square(2) P[0,0] = 1 #P[0,1] = 1 #print dlx_has_completion(P) #print dlx_has_unique_completion(P) for x in dlxcpp_find_completions(P, nr_to_find = 2): print x print "" #dlxcpp_has_unique_completion """ """ for n in range(5, 15): print n time gcs(back_circulant(n), dlx_has_unique_completion) time gcs(back_circulant(n), dlxcpp_has_unique_completion) """