Commits (5)
 ... ... @@ -4,10 +4,20 @@ Designs and Incidence Structures .. toctree:: :maxdepth: 1 ../sage/combinat/designs/block_design ../sage/combinat/designs/covering_design ../sage/combinat/designs/ext_rep ../sage/combinat/designs/incidence_structures ../sage/combinat/designs/steiner_quadruple_systems ../sage/combinat/designs/design_catalog Constructions ------------- .. toctree:: :maxdepth: 1 ../sage/combinat/designs/block_design ../sage/combinat/designs/bibd ../sage/combinat/designs/steiner_quadruple_systems ../sage/combinat/designs/latin_squares ../sage/combinat/designs/orthogonal_arrays
 ... ... @@ -36,3 +36,6 @@ deprecated_callable_import(14499, del deprecated_callable_import import sage.combinat.designs.steiner_quadruple_systems import sage.combinat.designs.orthogonal_arrays import sage.combinat.designs.latin_squares import sage.combinat.designs.bibd
This diff is collapsed.
 ... ... @@ -59,6 +59,7 @@ from sage.rings.integer_ring import ZZ from sage.rings.arith import binomial, integer_floor from sage.combinat.designs.incidence_structures import IncidenceStructure, IncidenceStructureFromMatrix from sage.misc.decorators import rename_keyword from sage.rings.finite_rings.constructor import FiniteField ### utility functions ------------------------------------------------------- ... ... @@ -84,21 +85,31 @@ def tdesign_params(t, v, k, L): def ProjectiveGeometryDesign(n, d, F, algorithm=None): """ Returns a projective geometry design. A projective geometry design of parameters n,d,F has for points the lines of F^{n+1}, and for blocks the d+1-dimensional subspaces of F^{n+1}, each of which contains \frac {|F|^{d+1}-1} {|F|-1} lines. INPUT: - n is the projective dimension - v is the number of points PPn(GF(q)) - d is the dimension of the subspaces of P = PPn(F) which make up the blocks. - d is the dimension of the subspaces of P = PPn(GF(q)) which make up the blocks - F is a finite field. - b is the number of d-dimensional subspaces of P - algorithm -- set to None by default, which results in using Sage's own implementation. In order to use GAP's implementation instead (i.e. its PGPointFlatBlockDesign function) set algorithm="gap". Note that GAP's "design" package must be available in this case. Wraps GAP Design's PGPointFlatBlockDesign. Does *not* require GAP's Design. EXAMPLES: EXAMPLES:: The points of the following design are the \frac {2^{2+1}-1} {2-1}=7 lines of \mathbb{Z}_2^{2+1}. It has 7 blocks, corresponding to each 2-dimensional subspace of \mathbb{Z}_2^{2+1}:: sage: designs.ProjectiveGeometryDesign(2, 1, GF(2)) Incidence structure with 7 points and 7 blocks ... ... @@ -109,7 +120,7 @@ def ProjectiveGeometryDesign(n, d, F, algorithm=None): q = F.order() from sage.interfaces.gap import gap, GapElement from sage.sets.set import Set if algorithm == None: if algorithm is None: V = VectorSpace(F, n+1) points = list(V.subspaces(1)) flats = list(V.subspaces(d+1)) ... ... @@ -132,6 +143,69 @@ def ProjectiveGeometryDesign(n, d, F, algorithm=None): gB.append([x-1 for x in b]) return BlockDesign(v, gB, name="ProjectiveGeometryDesign") def DesarguesianProjectivePlaneDesign(n): r""" Returns a (Desarguesian) projective plane of order n. A finite projective plane is a 2-design with n^2+n+1 lines (or blocks) and n^2+n+1 points. For more information on finite projective planes, see the :wikipedia:Projective_plane#Finite_projective_planes. INPUT: - n -- the finite projective plane's order OUTPUT: A finite projective plane obtained by considering the 1- and 2- dimensional spaces of F_n^3. .. SEEALSO:: :meth:ProjectiveGeometryDesign EXAMPLES:: sage: designs.DesarguesianProjectivePlaneDesign(2) Incidence structure with 7 points and 7 blocks Non-existent ones:: sage: designs.DesarguesianProjectivePlaneDesign(10) Traceback (most recent call last): ... ValueError: No projective plane design of order 10 exists. sage: designs.DesarguesianProjectivePlaneDesign(14) Traceback (most recent call last): ... ValueError: By the Bruck-Ryser-Chowla theorem, no projective plane of order 14 exists. An unknown one:: sage: designs.DesarguesianProjectivePlaneDesign(12) Traceback (most recent call last): ... ValueError: If such a projective plane exists, we do not know how to build it. """ from sage.rings.arith import two_squares try: F = FiniteField(n, 'x') except ValueError: if n == 10: raise ValueError("No projective plane design of order 10 exists.") try: if (n%4) in [1,2]: two_squares(n) except ValueError: raise ValueError("By the Bruck-Ryser-Chowla theorem, no projective" " plane of order "+str(n)+" exists.") raise ValueError("If such a projective plane exists, " "we do not know how to build it.") return ProjectiveGeometryDesign(2,1,F) def AffineGeometryDesign(n, d, F): r""" INPUT: ... ... @@ -159,14 +233,14 @@ def AffineGeometryDesign(n, d, F): sage: BD = designs.AffineGeometryDesign(3, 1, GF(2)) sage: BD.parameters() (2, 8, 2, 2) (2, 8, 2, 1) sage: BD.is_block_design() (True, [2, 8, 2, 2]) (True, [2, 8, 2, 1]) sage: BD = designs.AffineGeometryDesign(3, 2, GF(2)) sage: BD.parameters() (2, 8, 4, 12) (2, 8, 4, 3) sage: BD.is_block_design() (True, [3, 8, 4, 4]) (True, [3, 8, 4, 1]) """ q = F.order() from sage.interfaces.gap import gap, GapElement ... ... @@ -176,7 +250,7 @@ def AffineGeometryDesign(n, d, F): gap.eval("Subs:=AsSet(Subspaces(V,%s));"%d) gap.eval("CP:=Cartesian(points,Subs)") flats = gap.eval("flats:=List(CP,x->Sum(x))") # affine spaces gblcks = eval(gap.eval("AsSortedList(List(flats,f->Filtered([1..Length(points)],i->points[i] in f)));")) gblcks = eval(gap.eval("Set(List(flats,f->Filtered([1..Length(points)],i->points[i] in f)));")) v = q**n gB = [] for b in gblcks: ... ... @@ -251,89 +325,6 @@ def HadamardDesign(n): # A is the incidence matrix of the block design return IncidenceStructureFromMatrix(A,name="HadamardDesign") def steiner_triple_system(n): r""" Returns a Steiner Triple System A Steiner Triple System (STS) of a set \{0,...,n-1\} is a family S of 3-sets such that for any i \not = j there exists exactly one set of S in which they are both contained. It can alternatively be thought of as a factorization of the complete graph K_n with triangles. A Steiner Triple System of a n-set exists if and only if n \equiv 1 \pmod 6 or n \equiv 3 \pmod 6, in which case one can be found through Bose's and Skolem's constructions, respectively [AndHon97]_. INPUT: - n returns a Steiner Triple System of \{0,...,n-1\} EXAMPLE: A Steiner Triple System on 9 elements :: sage: sts = designs.steiner_triple_system(9) sage: sts Incidence structure with 9 points and 12 blocks sage: list(sts) [[0, 1, 5], [0, 2, 4], [0, 3, 6], [0, 7, 8], [1, 2, 3], [1, 4, 7], [1, 6, 8], [2, 5, 8], [2, 6, 7], [3, 4, 8], [3, 5, 7], [4, 5, 6]] As any pair of vertices is covered once, its parameters are :: sage: sts.parameters() (2, 9, 3, 1) An exception is raised for invalid values of n :: sage: designs.steiner_triple_system(10) Traceback (most recent call last): ... ValueError: Steiner triple systems only exist for n = 1 mod 6 or n = 3 mod 6 REFERENCE: .. [AndHon97] A short course in Combinatorial Designs, Ian Anderson, Iiro Honkala, Internet Editions, Spring 1997, http://www.utu.fi/~honkala/designs.ps """ name = "Steiner Triple System on "+str(n)+" elements" if n%6 == 3: t = (n-3)/6 Z = range(2*t+1) T = lambda (x,y) : x + (2*t+1)*y sts = [[(i,0),(i,1),(i,2)] for i in Z] + \ [[(i,k),(j,k),(((t+1)*(i+j)) % (2*t+1),(k+1)%3)] for k in range(3) for i in Z for j in Z if i != j] elif n%6 == 1: t = (n-1)/6 N = range(2*t) T = lambda (x,y) : x+y*t*2 if (x,y) != (-1,-1) else n-1 L1 = lambda i,j : (i+j) % (int((n-1)/3)) L = lambda i,j : L1(i,j)/2 if L1(i,j)%2 == 0 else t+(L1(i,j)-1)/2 sts = [[(i,0),(i,1),(i,2)] for i in range(t)] + \ [[(-1,-1),(i,k),(i-t,(k+1) % 3)] for i in range(t,2*t) for k in [0,1,2]] + \ [[(i,k),(j,k),(L(i,j),(k+1) % 3)] for k in [0,1,2] for i in N for j in N if i < j] else: raise ValueError("Steiner triple systems only exist for n = 1 mod 6 or n = 3 mod 6") from sage.sets.set import Set sts = Set(map(lambda x: Set(map(T,x)),sts)) return BlockDesign(n, sts, name=name) def BlockDesign(max_pt, blks, name=None, test=True): """ Returns an instance of the :class:IncidenceStructure class. ... ... @@ -362,6 +353,9 @@ def BlockDesign(max_pt, blks, name=None, test=True): else: raise TypeError("parameters are not those of a block design.") # Possibly in the future there will be methods which apply to block designs and # not incidence structures. None have been implemented yet though. The class # name BlockDesign_generic is reserved in the name space in case more ... ...
 ... ... @@ -39,11 +39,16 @@ Currently, this module gathers the following designs : :delim: | :meth:~sage.combinat.designs.block_design.ProjectiveGeometryDesign :meth:~sage.combinat.designs.block_design.DesarguesianProjectivePlaneDesign :meth:~sage.combinat.designs.bibd.BalancedIncompleteBlockDesign :meth:~sage.combinat.designs.block_design.AffineGeometryDesign :meth:~sage.combinat.designs.block_design.WittDesign :meth:~sage.combinat.designs.block_design.HadamardDesign :meth:~sage.combinat.designs.block_design.steiner_triple_system :meth:~sage.combinat.designs.block_design.steiner_quadruple_system :meth:~sage.combinat.designs.latin_squares.mutually_orthogonal_latin_squares :meth:~sage.combinat.designs.orthogonal_arrays.transversal_design :meth:~sage.combinat.designs.orthogonal_arrays.orthogonal_array :meth:~sage.combinat.designs.bibd.steiner_triple_system :meth:~sage.combinat.designs.steiner_quadruple_systems.steiner_quadruple_system And the :meth:designs.best_known_covering_design_from_LJCR  function ... ... @@ -59,11 +64,17 @@ REFERENCES: http://www.ccrwest.org/cover.html """ from sage.combinat.designs.block_design import (ProjectiveGeometryDesign, DesarguesianProjectivePlaneDesign, AffineGeometryDesign, WittDesign, HadamardDesign, steiner_triple_system) HadamardDesign) from sage.combinat.designs.steiner_quadruple_systems import steiner_quadruple_system from sage.combinat.designs.covering_design import best_known_covering_design_www as best_known_covering_design_from_LJCR from sage.combinat.designs.latin_squares import mutually_orthogonal_latin_squares from sage.combinat.designs.orthogonal_arrays import transversal_design, orthogonal_array from sage.combinat.designs.bibd import BalancedIncompleteBlockDesign, steiner_triple_system
 ... ... @@ -533,7 +533,7 @@ class IncidenceStructure(object): (True, [5, 12, 6, 1]) sage: BD = designs.AffineGeometryDesign(3, 1, GF(2)) sage: BD.is_block_design() (True, [2, 8, 2, 2]) (True, [2, 8, 2, 1]) """ from sage.combinat.designs.incidence_structures import coordinatewise_product from sage.combinat.combinat import unordered_tuples ... ...
 r""" Latin Squares and Mutually Orthogonal Latin Squares (MOLS) A Latin square is an n\times n array filled with n different symbols, each occurring exactly once in each row and exactly once in each column. Functions --------- """ def mutually_orthogonal_latin_squares(n,k=None, partitions = False): r""" Returns k Mutually Orthogonal n\times n Latin Squares (MOLS). For more information on Latin Squares and MOLS, see :mod:~sage.combinat.designs.latin_squares or the :wikipedia:Latin_square, or even the :wikipedia:Wikipedia entry on MOLS . INPUT: - n (integer) -- size of the latin square. - k (integer) -- returns k MOLS. If set to None (default), returns the maximum number of MOLS that Sage can build. .. WARNING:: This has no reason to be the maximum number of n\times n MOLS, just the best Sage can do ! - partition (boolean) -- a Latin Square can be seen as 3 partitions of the n^2 cells of the array into n sets of size n, respectively : * The partition of rows * The partition of columns * The partition of number (cells numbered with 0, cells numbered with 1, ...) These partitions have the additional property that any two sets from different partitions intersect on exactly one element. When partition is set to True, this function returns a list of k+2 partitions satisfying this intersection property instead of the k+2 MOLS (though the data is exactly the same in both cases). EXAMPLES:: sage: designs.mutually_orthogonal_latin_squares(5) [ [0 1 2 3 4] [0 1 2 3 4] [0 1 2 3 4] [0 1 2 3 4] [3 0 1 4 2] [4 3 0 2 1] [1 2 4 0 3] [2 4 3 1 0] [4 3 0 2 1] [1 2 4 0 3] [2 4 3 1 0] [3 0 1 4 2] [1 2 4 0 3] [2 4 3 1 0] [3 0 1 4 2] [4 3 0 2 1] [2 4 3 1 0], [3 0 1 4 2], [4 3 0 2 1], [1 2 4 0 3] ] sage: designs.mutually_orthogonal_latin_squares(7,3) [ [0 1 2 3 4 5 6] [0 1 2 3 4 5 6] [0 1 2 3 4 5 6] [4 0 3 1 6 2 5] [5 6 0 4 2 1 3] [6 4 1 0 5 3 2] [5 6 0 4 2 1 3] [6 4 1 0 5 3 2] [1 3 5 2 0 6 4] [6 4 1 0 5 3 2] [1 3 5 2 0 6 4] [2 5 4 6 3 0 1] [1 3 5 2 0 6 4] [2 5 4 6 3 0 1] [3 2 6 5 1 4 0] [2 5 4 6 3 0 1] [3 2 6 5 1 4 0] [4 0 3 1 6 2 5] [3 2 6 5 1 4 0], [4 0 3 1 6 2 5], [5 6 0 4 2 1 3] ] sage: designs.mutually_orthogonal_latin_squares(5,2,partitions=True) [[[0, 1, 2, 3, 4], [5, 6, 7, 8, 9], [10, 11, 12, 13, 14], [15, 16, 17, 18, 19], [20, 21, 22, 23, 24]], [[0, 5, 10, 15, 20], [1, 6, 11, 16, 21], [2, 7, 12, 17, 22], [3, 8, 13, 18, 23], [4, 9, 14, 19, 24]], [[0, 6, 12, 18, 24], [1, 7, 14, 15, 23], [2, 9, 13, 16, 20], [3, 5, 11, 19, 22], [4, 8, 10, 17, 21]], [[0, 7, 13, 19, 21], [1, 9, 10, 18, 22], [2, 8, 11, 15, 24], [3, 6, 14, 17, 20], [4, 5, 12, 16, 23]]] """ from sage.rings.finite_rings.constructor import FiniteField from sage.combinat.designs.block_design import AffineGeometryDesign from sage.rings.arith import is_prime_power from sage.matrix.constructor import Matrix from sage.rings.arith import factor if k is not None and k >= n: raise ValueError("This does not exist (cite theorem)") if is_prime_power(n): if k is None: k = n-1 # Section 6.4.1 of [Stinson2004] Fp = FiniteField(n,'x') B = AffineGeometryDesign(2,1,Fp).blocks() parallel_classes = [[] for _ in range(k+2)] for b in B: for p in parallel_classes: if p == [] or all(i not in p[0] for i in b): p.append(b) break if partitions: return parallel_classes coord = {v:i for i,L in enumerate(parallel_classes[0]) for v in L} coord = {v:(coord[v],i) for i,L in enumerate(parallel_classes[1]) for v in L} matrices = [] for P in parallel_classes[2:]: matrices.append(Matrix({coord[v]:i for i,L in enumerate(P) for v in L })) return matrices else: # Theorem 6.33 of [Stinson2004], MacNeish's theorem. subcases = [p**i for p,i in factor(n)] s = min(subcases)-1 if k is None: k = s elif k > s: raise ValueError("I don't know how to build these MOLS.") subcalls = [mutually_orthogonal_latin_squares(p,k) for p in subcases] matrices = [latin_square_product(*[sc[i] for sc in subcalls]) for i in range(k)] if partitions: partitions = [[[i*n+j for j in range(n)] for i in range(n)], [[j*n+i for j in range(n)] for i in range(n)]] for m in matrices: partition = [[] for i in range(n)] for i in range(n): for j in range(n): partition[m[i,j]].append(i*n+j) partitions.append(partition) return partitions else: return matrices def latin_square_product(M,N,*others): r""" Returns the product of two (or more) latin squares. This is Lemma 6.25 of [Stinson2004]_. INPUT: An arbitrary number of latin squares (greater than 2). EXAMPLES:: sage: from sage.combinat.designs.transversal_design import latin_square_product sage: m=designs.mutually_orthogonal_latin_squares(4)[0] sage: latin_square_product(m,m,m) 64 x 64 sparse matrix over Integer Ring """ from sage.matrix.constructor import Matrix m = M.nrows() n = N.nrows() D = {((i,j),(ii,jj)):(M[i,ii],N[j,jj]) for i in range(m) for ii in range(m) for j in range(n) for jj in range(n)} L = lambda (i,j):i*n+j D = {(L(c[0]),L(c[1])): L(v) for c,v in D.iteritems()} P = Matrix(D) if others: return latin_square_product(P, others[0],*others[1:]) else: return P
 """ Orthogonal arrays This module gathers anything related to orthogonal arrays. And, incidentally, to transversal designs. Functions --------- """ def transversal_design(k,n,t=2,check=True): r""" Returns a transversal design of parameters k,n. INPUT: - n,k,t -- integers. - check (boolean) -- whether to check that output is correct before returning it. As this is expected to be useless (but we are cautious guys), you may want to disable it whenever you want speed. Set to True by default. **Definition** A transversal design of parameters n, k is a collection \mathcal S of k-subsets of V=V_1 \sqcup \dots \sqcup V_k (where |V_i|=n) such that : * Any element S\in \mathcal S intersects each set V_i on exactly one element. * Any two elements v_i\in V_i,v_j\in V_j with i\neq j belong to exactly one element of \mathcal S For more information on transversal designs, see http://mathworld.wolfram.com/TransversalDesign.html_. EXAMPLES:: sage: designs.transversal_design(5,5) [[0, 5, 10, 15, 20], [0, 6, 12, 18, 24], [0, 7, 14, 16, 23], [0, 8, 11, 19, 22], [0, 9, 13, 17, 21], [1, 6, 11, 16, 21], [1, 7, 13, 19, 20], [1, 8, 10, 17, 24], [1, 9, 12, 15, 23], [1, 5, 14, 18, 22], [2, 7, 12, 17, 22], [2, 8, 14, 15, 21], [2, 9, 11, 18, 20], [2, 5, 13, 16, 24], [2, 6, 10, 19, 23], [3, 8, 13, 18, 23], [3, 9, 10, 16, 22], [3, 5, 12, 19, 21], [3, 6, 14, 17, 20], [3, 7, 11, 15, 24], [4, 9, 14, 19, 24], [4, 5, 11, 17, 23], [4, 6, 13, 15, 22], [4, 7, 10, 18, 21], [4, 8, 12, 16, 20]] """ # Section 6.6 OA = orthogonal_array(k,n) TD = [[i*n+c for i,c in enumerate(l)] for l in OA] if check and not is_transversal_design(TD,k,n): raise RuntimeError("Sage returns wrong results ! Please report the bug.") return TD def is_transversal_design(B,k,n): r""" Checks that a given set B of blocks is a transversal design. See :mod:sage.combinat.design.orthogonal_arrays for a definition. INPUT: - B -- the list of blocks - k,n -- integers .. NOTE:: The tranversal design must have 0,...,kn-1 as a ground set, partitionned as k sets of size n : \{0,...,k-1\}\cup \{k,...,2k-1\}\cup\{k(n-1),...,kn-1\}. """ from sage.graphs.generators.basic import CompleteGraph from itertools import combinations g = k*CompleteGraph(n) m = g.size() for X in B: if len(X) != k: return False g.add_edges(list(combinations(X,2))) if g.size() != m+(len(X)*(len(X)-1))/2: return False m = g.size() if not g.is_clique(): return False return True def orthogonal_array(k,n,t=2,check=True): r""" Returns an orthogonal array of parameters k,n,t INPUT: - k (integer) -- number of columns - n (integer) -- number of symbols - t (integer) -- only t=2 is available at the moment. - check (boolean) -- whether to check that output is correct before returning it. As this is expected to be useless (but we are cautious guys), you may want to disable it whenever you want speed. Set to True by default. For more information on orthogonal arrays, see :wikipedia:Orthogonal_array. .. NOTE:: This method implements theorems from [Stinson2004]_. See the code's documentation for details. .. TODO:: Implement Wilson's construction. See page 146 of [Stinson2004]_. EXAMPLES:: sage: from sage.combinat.designs.orthogonal_arrays import orthogonal_array sage: orthogonal_array(5,5) [[0, 0, 0, 0, 0], [0, 1, 2, 3, 4], [0, 2, 4, 1, 3], [0, 3, 1, 4, 2], [0, 4, 3, 2, 1], [1, 1, 1, 1, 1], [1, 2, 3, 4, 0], [1, 3, 0, 2, 4], [1, 4, 2, 0, 3], [1, 0, 4, 3, 2], [2, 2, 2, 2, 2], [2, 3, 4, 0, 1], [2, 4, 1, 3, 0], [2, 0, 3, 1, 4], [2, 1, 0, 4, 3], [3, 3, 3, 3, 3], [3, 4, 0, 1, 2], [3, 0, 2, 4, 1], [3, 1, 4, 2, 0], [3, 2, 1, 0, 4], [4, 4, 4, 4, 4], [4, 0, 1, 2, 3], [4, 1, 3, 0, 2], [4, 2, 0, 3, 1], [4, 3, 2, 1, 0]] """ from sage.rings.arith import is_prime_power from sage.rings.finite_rings.constructor import FiniteField OA = None if t != 2: raise NotImplementedError("Only implemented for t == 2") if k <= 1: raise ValueError("Not defined for this value !") if k == t: from itertools import product OA = [list(t) for t in product(*[range(n) for _ in range(k)])] # Theorem 6.39 from [Stinson2004] elif 2 <= k and k <= n and is_prime_power(n): M = [] Fp = FiniteField(n,'x') vv = list(Fp)[:k] relabel = {x:i for i,x in enumerate(Fp)} for i in Fp: for j in Fp: M.append([relabel[i+j*v] for v in vv]) OA = M # Theorem 6.40 from [Stinson2004] elif k == n+1 and is_prime_power(n): M = orthogonal_array(n,n) for i,l in enumerate(M): l.append(i%n) OA = M # Section 6.5.1 from [Stinson2004] if OA is None: from latin_squares import mutually_orthogonal_latin_squares try: mols = mutually_orthogonal_latin_squares(n,k-2) except ValueError: mols = None if mols: OA = [[i,j]+[m[i,j] for m in mols] for i in range(n) for j in range(n)] if OA is None: raise NotImplementedError("I don't know how to build this OA !") if check and not is_orthogonal_array(OA,k,n,t): raise RuntimeError("Sage returns wrong results ! Please report the bug.") return OA def is_orthogonal_array(M,k,n,t): r""" Checks that the integer matrix M is an OA(k,n,t). INPUT: - M -- an integer matrix of size k^t\times n. - k,n,t -- integers. EXAMPLES:: sage: _ = designs.orthogonal_array(5,9,check=True) # indirect doctest """ if t != 2: raise NotImplementedError if not all([len(l) == k for l in M]): return False from itertools import combinations for S in combinations(range(k),2): fs = frozenset([tuple([l[i] for i in S]) for l in M]) if len(fs) != n**2: return False return True