trac #16295 : Faster is_orthogonal_array

src/module_list.py

@@ -255,6 +255,9 @@ def uname_specific(name, value, alternative):
     Extension('sage.combinat.crystals.letters',
               sources=['sage/combinat/crystals/letters.pyx']),
 
+    Extension('sage.combinat.designs.designs_pyx',
+              sources=['sage/combinat/designs/designs_pyx.pyx']),
+
     ################################
     ##
     ## sage.crypto

src/sage/combinat/designs/block_design.py

@@ -299,7 +299,7 @@ def projective_plane_to_OA(pplane, pt=None, check=True):
     assert len(OA) == n**2, "pplane is not a projective plane"
 
     if check:
-        from orthogonal_arrays import is_orthogonal_array
+        from designs_pyx import is_orthogonal_array
         is_orthogonal_array(OA,n+1,n,2)
 
     return OA

src/sage/combinat/designs/designs_pyx.pyx

@@ -0,0 +1,118 @@
+r"""
+Cython functions for combinatorial designs
+
+This module implements the design methods that need to be somewhat efficient.
+
+Functions
+---------
+"""
+include "sage/misc/bitset.pxi"
+
+def is_orthogonal_array(OA, int k, int n, int t=2, verbose=False, terminology="OA"):
+    r"""
+    Check that the integer matrix `M` is an `OA(k,n,t)`.
+
+    See :func:`~sage.combinat.designs.orthogonal_arrays.orthogonal_array`
+    for a definition.
+
+    INPUT:
+
+    - ``OA`` -- the Orthogonal Array to be tested
+
+    - ``k,n,t`` integers.
+
+    - ``verbose`` (boolean) -- whether to display some information when ``OA``
+      is not an orthogona array `OA(k,n)`.
+
+    - ``terminology`` (string) -- how to phrase the information when ``verbose =
+      True``. Possible values are `"OA"`, `"MOLS"`.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.designs.designs_pyx import is_orthogonal_array
+        sage: OA = designs.orthogonal_array(8,9)
+        sage: is_orthogonal_array(OA,8,9)
+        True
+        sage: is_orthogonal_array(OA,8,10)
+        False
+        sage: OA[4][3] = 1
+        sage: is_orthogonal_array(OA,8,9)
+        False
+        sage: is_orthogonal_array(OA,8,9,verbose=True)
+        Rows 0 and 3 are not orthogonal
+        False
+    """
+    cdef int n2 = n*n
+    cdef int x
+    if any(len(R) != k for R in OA):
+        if verbose:
+            print {"OA"   : "Some row does not have length "+str(k),
+                   "MOLS" : "The number of matrices is not "+str(k)}[terminology]
+        return False
+
+    if len(OA) != n2:
+        if verbose:
+            print {"OA"   : "The number of rows is {} instead of {}^2={}".format(len(OA),n,n2),
+                   "MOLS" : "All matrices do not have dimension n^2={}^2".format(n)}[terminology]
+        return False
+
+    cdef int i,j,l
+
+    # A copy of OA
+    cdef unsigned short * OAc = <unsigned short *> sage_malloc(k*n2*sizeof(unsigned short))
+
+    # A cache to multiply by n
+    cdef unsigned int * times_n = <unsigned int *> sage_malloc((n+1)*sizeof(unsigned int))
+
+    cdef unsigned short * C1
+    cdef unsigned short * C2
+
+    # failed malloc ?
+    if OAc is NULL or times_n is NULL:
+        if OAc is not NULL:
+            sage_free(OAc)
+        if times_n is not NULL:
+            sage_free(times_n)
+        raise MemoryError
+
+    # Filling OAc
+    for i,R in enumerate(OA):
+        for j,x in enumerate(R):
+            OAc[j*n2+i] = x
+            if x < 0 or x >= n:
+                if verbose:
+                    print {"OA"   : "{} is not in the interval [0..{}]".format(x,n-1),
+                           "MOLS" : "Entry {} was expected to be in the interbal [0..{}]".format(x,n-1)}[terminology]
+                sage_free(OAc)
+                sage_free(times_n)
+                return False
+
+    # Filling times_n
+    for i in range(n+1):
+        times_n[i] = i*n
+
+    # A bitset to keep trac of pairs of values
+    cdef bitset_t seen
+    bitset_init(seen, n2)
+
+    for i in range(k): # For any column C1
+        C1 = OAc+i*n2
+        for j in range(i+1,k): # For any column C2 > C1
+            C2 = OAc+j*n2
+            bitset_set_first_n(seen, 0) # No pair has been seen yet
+            for l in range(n2):
+                bitset_add(seen,times_n[C1[l]]+C2[l])
+
+            if bitset_len(seen) != n2: # Have we seen all pairs ?
+                sage_free(OAc)
+                sage_free(times_n)
+                bitset_free(seen)
+                if verbose:
+                    print {"OA"   : "Rows {} and {} are not orthogonal".format(i,j),
+                           "MOLS" : "Matrices {} and {} are not orthogonal".format(i,j)}[terminology]
+                return False
+
+    sage_free(OAc)
+    sage_free(times_n)
+    bitset_free(seen)
+    return True

src/sage/combinat/designs/latin_squares.py

@@ -92,7 +92,7 @@ def are_mutually_orthogonal_latin_squares(l, verbose=False):
         sage: are_mutually_orthogonal_latin_squares([m2,m3])
         True
         sage: are_mutually_orthogonal_latin_squares([m1,m2,m3], verbose=True)
-        matrices 0 and 2 are not orthogonal
+        Matrices 0 and 2 are not orthogonal
         False
 
         sage: m = designs.mutually_orthogonal_latin_squares(8,7)
@@ -102,32 +102,15 @@ def are_mutually_orthogonal_latin_squares(l, verbose=False):
     if not l:
         raise ValueError("the list must be non empty")
 
-    n = l[0].nrows()
-    if any(M.nrows() != n and M.ncols() != n for M in l):
+    n = l[0].ncols()
+    k = len(l)
+    if any(M.ncols() != n or M.nrows() != n for M in l):
         if verbose:
-            print "some matrix has wrong dimension"
+            print "Not all matrices are square matrices of the same dimensions"
         return False
 
-    # check that the matrices in l are actually latin
-    for i,M in enumerate(l):
-        if (any(sorted(r) != range(n) for r in M.rows()) or
-            any(sorted(c) != range(n) for c in M.columns())):
-            if verbose:
-                print "matrix %d is not latin"%i
-            return False
-
-    # check orthogonality of each pair
-    for k1 in xrange(len(l)):
-        M1 = l[k1]
-        for k2 in xrange(k1):
-            M2 = l[k2]
-            L = [(M1[i,j],M2[i,j]) for i in xrange(n) for j in xrange(n)]
-            if len(set(L)) != len(L):
-                if verbose:
-                    print "matrices %d and %d are not orthogonal"%(k2,k1)
-                return False
-
-    return True
+    from designs_pyx import is_orthogonal_array
+    return is_orthogonal_array(zip(*[[x for R in M for x in R] for M in l]),k,n, verbose=verbose, terminology="MOLS")
 
 def mutually_orthogonal_latin_squares(n,k, partitions = False, check = True, existence=False, who_asked=tuple()):
     r"""

src/sage/combinat/designs/orthogonal_arrays.py

@@ -359,28 +359,8 @@ def is_transversal_design(B,k,n, verbose=False):
         sage: is_transversal_design(TD, 4, 4)
         False
     """
-    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:
-            if verbose:
-                print "A set has wrong size"
-            return False
-        g.add_edges(list(combinations(X,2)))
-        if g.size() != m+(len(X)*(len(X)-1))/2:
-            if verbose:
-                print "A pair appears twice"
-            return False
-        m = g.size()
-
-    if not g.is_clique():
-        if verbose:
-            print "A pair did not appear"
-        return False
-
-    return True
+    from designs_pyx import is_orthogonal_array
+    return is_orthogonal_array([[x%n for x in R] for R in B],k,n,verbose=verbose)
 
 @cached_function
 def find_wilson_decomposition(k,n):
@@ -801,12 +781,12 @@ def orthogonal_array(k,n,t=2,check=True,existence=False,who_asked=tuple()):
             if existence:
                 return projective_plane(n, existence=True)
             p = projective_plane(n, check=False)
-            OA = projective_plane_to_OA(p)
+            OA = projective_plane_to_OA(p, check=False)
         else:
             if existence:
                 return True
             p = projective_plane(n, check=False)
-            OA = [l[:k] for l in projective_plane_to_OA(p)]
+            OA = [l[:k] for l in projective_plane_to_OA(p, check=False)]
 
     # Constructions from the database
     elif n in OA_constructions and k <= OA_constructions[n][0]:
@@ -857,6 +837,7 @@ def orthogonal_array(k,n,t=2,check=True,existence=False,who_asked=tuple()):
         raise NotImplementedError("I don't know how to build this orthogonal array!")
 
     if check:
+        from designs_pyx import is_orthogonal_array
         assert is_orthogonal_array(OA,k,n,t)
 
     return OA
@@ -1056,51 +1037,3 @@ def OA_from_wider_OA(OA,k):
     if len(OA[0]) == k:
         return OA
     return [L[:k] for L in OA]
-
-def is_orthogonal_array(M,k,n,t,verbose=False):
-    r"""
-    Check that the integer matrix `M` is an `OA(k,n,t)`.
-
-    See :func:`~sage.combinat.designs.orthogonal_arrays.orthogonal_array`
-    for a definition.
-
-    INPUT:
-
-    - ``M`` -- an integer matrix of size `k^t \times n`
-
-    - ``k, n, t`` -- integers
-
-    - ``verbose`` -- boolean, if ``True`` provide an information on where ``M``
-      fails to be an `OA(k,n,t)`.
-
-    EXAMPLES::
-
-        sage: OA = designs.orthogonal_array(5,9,check=True) # indirect doctest
-        sage: from sage.combinat.designs.orthogonal_arrays import is_orthogonal_array
-        sage: is_orthogonal_array(OA, 5, 9, 2)
-        True
-        sage: is_orthogonal_array(OA, 4, 5, 2)
-        False
-    """
-    if t != 2:
-        raise NotImplementedError("only implemented for t=2")
-
-    if len(M) != n**2:
-        if verbose:
-            print "wrong number of rows"
-        return False
-
-    if any(len(l) != k for l in M):
-        if verbose:
-            print "a row has the wrong size"
-        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:
-            if verbose:
-                print "for the choice %s of columns we do not get all tuples"%(S,)
-            return False
-
-    return True