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
<sage.combinat.designs.covering_design.best_known_covering_design_www>` 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 <Graeco-Latin_square#Mutually_orthogonal_Latin_squares>`.
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