Merge branch 'u/dimpase/muzS6' of trac.sagemath.org:sage into iokar

src/sage/graphs/generators/families.py

@@ -11,6 +11,9 @@
 #                              and Emily A. Kirkman
 #           Copyright (C) 2009 Michael C. Yurko <myurko@gmail.com>
 #
+#           Copyright (C) 2016 Rowan Schrecker <rowan.schrecker@hertford.ox.ac.uk>
+#            (Rowan Schrecker supported by UK EPSRC grant EP/K040251/2)
+#
 # Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
 #                         http://www.gnu.org/licenses/
 ###########################################################################
@@ -2680,3 +2683,262 @@ def TuranGraph(n,r):
     g.name('Turan Graph with n: {}, r: {}'.format(n,r))
 
     return g
+
+def MuzychukS6Graph(n, d, Phi='fixed', Sigma='fixed', verbose=False):
+    r"""
+    Return a strongly regular graph of S6 type from [Mu07]_ on `n^d((n^d-1)/(n-1)+1)` vertices
+
+    The construction depends upon a number of parameters, two of then, `n` and `d`, mandatory,
+    and `\Phi` and `\Sigma` are mappings defined in [Mu07]_.
+    Let `m=\frac{n^d-1}{n-1}+1`. Then these graphs have parameters
+    `(mn^d, n^{d-1}(m-1) - 1,\mu - 2,\mu)`, with `\mu=\frac{n^{d-1}-1}{n-1}n^{d-1}`.
+    Let `L` be the complete graph on `M:=\{0,..., m-1\}` with the matching
+    `\{(2i,2i+1) | i=0,...,m/2\}` removed.
+
+
+    INPUT:
+
+    - ``n`` (integer)-- a prime power
+
+    - ``d`` (integer)-- must be odd if `n` is odd
+
+    - ``Phi`` is an optional parameter of the construction; it must be either
+
+        - 'fixed'-- this will generate fixed default `\Phi_i`, for `i \in M`, or
+
+        - 'random'-- `\Phi_i` are generated at random, or
+
+        - A dictionary describing the functions `\Phi_i`; for `i \in M`,
+          Phi[(i, T)] in `M`, for each edge T of `L` on `i`.
+          Also, each `\Phi_i` must be injective.
+
+    - ``Sigma`` is an optional parameter of the construction; it must be either
+
+        - 'fixed'-- this will generate a fixed default `\Sigma`, or
+
+        - 'random'-- `\Sigma` is generated at random.
+
+    - ``verbose`` (Boolean)-- default is False. If True, print progress information
+
+    .. TODO::
+
+        Implement the possibility to explicitly supply the parameter `\Sigma`
+        of the construction.
+
+    EXAMPLES::
+
+        sage: graphs.MuzychukS6Graph(3, 3).is_strongly_regular(parameters=True)
+        (378, 116, 34, 36)
+        sage: phi={(2,(0,2)):0,(1,(1,3)):1,(0,(0,3)):1,(2,(1,2)):1,(1,(1,
+        ....:  2)):0,(0,(0,2)):0,(3,(0,3)):0,(3,(1,3)):1}
+        sage: graphs.MuzychukS6Graph(2,2,Phi=phi).is_strongly_regular(parameters=True)
+        (16, 5, 0, 2)
+
+    TESTS::
+
+        sage: graphs.MuzychukS6Graph(2,2,Phi='random',Sigma='random').is_strongly_regular(parameters=True)
+        (16, 5, 0, 2)
+        sage: graphs.MuzychukS6Graph(3,3,Phi='random',Sigma='random').is_strongly_regular(parameters=True)
+        (378, 116, 34, 36)
+        sage: graphs.MuzychukS6Graph(3,2)
+        Traceback (most recent call last):
+        ...
+        AssertionError: n must be even or d must be odd
+        sage: graphs.MuzychukS6Graph(6,2)
+        Traceback (most recent call last):
+        ...
+        AssertionError: n must be a prime power
+        sage: graphs.MuzychukS6Graph(3,1)
+        Traceback (most recent call last):
+        ...
+        AssertionError: d must be at least 2
+        sage: graphs.MuzychukS6Graph(3,3,Phi=42)
+        Traceback (most recent call last):
+        ...
+        AssertionError: Phi must be a dictionary or 'random' or 'fixed'
+        sage: graphs.MuzychukS6Graph(3,3,Sigma=42)
+        Traceback (most recent call last):
+        ...
+        ValueError: Sigma must be 'random' or 'fixed'
+
+    REFERENCE:
+
+    .. [Mu07] M. Muzychuk.
+       A generalization of Wallis-Fon-Der-Flaass construction of strongly regular graphs.
+       J. Algebraic Combin., 25(2):169–187, 2007.
+    """
+    ### TO DO: optimise
+    ###        add option to return phi, sigma? generate phi, sigma from seed? (int say?)
+
+    from sage.combinat.designs.block_design import ProjectiveGeometryDesign
+    from sage.misc.prandom import randrange
+    from sage.misc.functional import is_even
+    from sage.arith.misc import is_prime_power
+    from sage.graphs.generators.basic import CompleteGraph
+    from sage.rings.finite_rings.finite_field_constructor import GF
+    from sage.matrix.special import ones_matrix
+    from sage.matrix.constructor import matrix
+    from sage.rings.rational_field import QQ
+    from sage.rings.integer_ring import ZZ
+    from time import time
+    from __builtin__ import range # we cannot use xrange here
+
+    assert is_even(n * (d-1)), 'n must be even or d must be odd'
+    assert is_prime_power(n),  'n must be a prime power'
+    assert d > 1,              'd must be at least 2'
+    t = time()
+
+    #build L, L_i and the design
+    m = int((n**d-1)/(n-1) + 1) #from m = p + 1, p = (n^d-1) / (n-1)
+    L = CompleteGraph(m)
+    L.delete_edges([(2*x, 2*x + 1) for x in range(m/2)])
+    L_i = [0]*m
+    for x in range(m):
+        L_i[x] = L.edges_incident(x, labels=False)
+    Design = ProjectiveGeometryDesign(d, d-1, GF(n, 'a'), point_coordinates=False)
+    projBlocks = Design.blocks()
+    atInf = projBlocks[-1]
+    Blocks = [[x for x in block if x not in atInf] for block in projBlocks[:-1]]
+    if verbose:
+        print('finished preamble at %f (+%f)' % (time() - t, time() - t))
+    t1 = time()
+
+    #sort the hyperplanes into parallel classes
+    ParClasses = [Blocks]
+    while ParClasses[0]:
+        nextHyp = ParClasses[0].pop()
+        for C in ParClasses[1:]:
+            listC = sum(C,[])
+            for x in nextHyp:
+                if x in listC:
+                    break
+            else:
+                C.append(nextHyp)
+                break
+        else:
+            ParClasses.append([nextHyp])
+    del ParClasses[0]
+    if verbose:
+        print('finished ParClasses at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+
+    #build E^C_j
+    E = {}
+    v = ZZ(n**d)
+    k = ZZ(n**(d-1))
+    ones = ones_matrix(v)
+    for C in ParClasses:
+        EC = matrix(QQ, v)
+        for line in C:
+            for i in line:
+                for j in line:
+                    EC[i, j] = 1/k
+        EC -= ones/v
+        E[tuple(C[0])] = EC
+    if verbose:
+        print('finished E at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+
+    #handle Phi
+    if Phi == 'random':
+        Phi = {}
+        for x in range(m):
+            temp = range(len(ParClasses))
+            for line in L_i[x]:
+                rand = randrange(0, len(temp))
+                Phi[(x, line)] = temp.pop(rand)
+    elif Phi == 'fixed':
+        Phi = {}
+        for x in range(m):
+            val = 0
+            for line in L_i[x]:
+                Phi[(x, line)] = val
+                val+=1
+    else:
+        assert isinstance(Phi, dict), \
+            "Phi must be a dictionary or 'random' or 'fixed'"
+        assert set(Phi.keys()) == \
+        set([(x, line) for x in range(m) for line in L_i[x]]), \
+        'each Phi_i must have domain L_i'
+        for x in range(m):
+            assert m - 2 == len(set([val
+                for (key, val) in Phi.items() if key[0] == x])), \
+            'each phi_i must be injective'
+        for val in Phi.values():
+            assert val in range(m-1), \
+            'codomain should be {0,..., (n^d - 1)/(n - 1) - 1}'
+    phi = {}
+    for x in range(m):
+        for line in L_i[x]:
+            phi[(x, line)] = ParClasses[Phi[(x, line)]]
+    if verbose:
+        print('finished phi at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+
+    #handle sigma
+    sigma = {}
+    if Sigma == 'random':
+        for x in range(m):
+            for line in L_i[x]:
+                [i, j] = line
+                temp = phi[(j, line)][:]
+                for hyp in phi[(i, line)]:
+                    rand = randrange(0, len(temp))
+                    sigma[(i, j, tuple(hyp))] = temp[rand]
+                    sigma[(j, i, tuple(temp[rand]))] = hyp
+                    del temp[rand]
+    elif Sigma == 'fixed':
+        for x in range(m):
+            for line in L_i[x]:
+                [i, j] = line
+                temp = phi[(j, line)][:]
+                for hyp in phi[(i, line)]:
+                    val = temp.pop()
+                    sigma[(i, j, tuple(hyp))] = val
+                    sigma[(j, i, tuple(val))] = hyp
+    else:
+        raise ValueError("Sigma must be 'random' or 'fixed'")
+    if verbose:
+        print('finished sigma at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+
+    #build V
+    edges = [] ###how many? *m^2*n^2
+    for (i, j) in L.edges(labels=False):
+        for hyp in phi[(i, (i, j))]:
+            for x in hyp:
+                newEdges = [((i, x), (j, y))
+                            for y in sigma[(i, j, tuple(hyp))]]
+                edges.extend(newEdges)
+    if verbose:
+        print('finished edges at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+    V = Graph(edges)
+    if verbose:
+        print('finished V at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+
+    #build D_i, F_i and A_i
+    D_i = [0]*m
+    for x in range(m):
+        D_i[x] = sum([E[tuple(phi[x, line][0])] for line in L_i[x]])
+    F_i = [1 - D_i[x] - ones/v for x in range(m)]
+    #as the sum of (1/v)*J_\Omega_i, D_i, F_i is identity
+    A_i = [0]*m
+    for x in range(m):
+        A_i[x] = ((v-k)/v)*ones - k*F_i[x]
+        #we know A_i = k''*(1/v)*J_\Omega_i + r''*D_i + s''*F_i,
+        #and (k'', s'', r'') = (v - k, 0, -k)
+    if verbose:
+        print('finished D, F and A at %f (+%f)' % (time() - t, time() - t1))
+    t1 = time()
+
+    #add the edges of the graph of B to V
+    for i in range(m):
+        V.add_edges([((i, x), (i, y)) for x in range(v)
+                     for y in range(v) if not A_i[i][(x, y)]])
+
+    V.name('Muzychuk S6 graph with parameters ('+str(n)+','+str(d)+')')
+    if verbose:
+        print('finished at %f (+%f)' % ((time() - t), time() - t1))
+    return V

src/sage/graphs/graph_generators.py

@@ -218,6 +218,7 @@ def __append_to_doc(methods):
      "line_graph_forbidden_subgraphs",
      "MathonPseudocyclicMergingGraph",
      "MathonPseudocyclicStronglyRegularGraph",
+     "MuzychukS6Graph",
      "MycielskiGraph",
      "MycielskiStep",
      "NKStarGraph",
@@ -1999,6 +2000,7 @@ def quadrangulations(self, order, minimum_degree=None, minimum_connectivity=None
     line_graph_forbidden_subgraphs = staticmethod(sage.graphs.generators.families.line_graph_forbidden_subgraphs)
     MathonPseudocyclicMergingGraph = staticmethod(sage.graphs.generators.families.MathonPseudocyclicMergingGraph)
     MathonPseudocyclicStronglyRegularGraph = staticmethod(sage.graphs.generators.families.MathonPseudocyclicStronglyRegularGraph)
+    MuzychukS6Graph        = staticmethod(sage.graphs.generators.families.MuzychukS6Graph)
     MycielskiGraph         = staticmethod(sage.graphs.generators.families.MycielskiGraph)
     MycielskiStep          = staticmethod(sage.graphs.generators.families.MycielskiStep)
     NKStarGraph            = staticmethod(sage.graphs.generators.families.NKStarGraph)

src/sage/graphs/strongly_regular_db.pyx

@@ -137,6 +137,42 @@ def is_mathon_PC_srg(int v,int k,int l,int mu):
                                     MathonPseudocyclicStronglyRegularGraph
                 return (MathonPseudocyclicStronglyRegularGraph,t)
 
+@cached_function
+def is_muzychuk_S6(int v, int k, int l, int mu):
+    r"""
+    Test whether some Muzychuk S6 graph is (v, k, l, mu)-strongly regular.
+
+    INPUT:
+
+    - ``v, k, l, mu`` (integers)
+
+    OUTPUT:
+
+    A tuple ``t`` such that ``t[0](*t[1:])`` builds the required graph if it exists,
+    and ``None`` otherwise.
+
+    EXAMPLES::
+
+        sage: from sage.graphs.strongly_regular_db import is_muzychuk_S6
+        sage: t = is_muzychuk_S6(378, 116, 34, 36)
+        sage: G = t[0](*t[1:]); G
+        Muzychuk S6 graph with parameters (3,3): Graph on 378 vertices
+        sage: G.is_strongly_regular(parameters=True)
+        (378, 116, 34, 36)
+        sage: t = is_muzychuk_S6(5, 5, 5, 5); t
+    """
+    cdef int n, d
+    from sage.rings.integer_ring import ZZ
+    n_list = [n for n in range(l-1) if ZZ(n).is_prime_power()]
+    for n in n_list:
+        d = 2
+        while n**d * ((n**d-1)/(n-1)+1) <= v:
+            if v == n**d * ((n**d-1)/(n-1)+1) and k == n**(d-1)*(n**d-1)/(n-1) - 1\
+            and l == mu - 2 and mu == n**(d-1) * (n**(d-1)-1) / (n-1):
+                from sage.graphs.generators.families import MuzychukS6Graph
+                return (MuzychukS6Graph, n, d)
+            d += 1
+
 @cached_function
 def is_orthogonal_array_block_graph(int v,int k,int l,int mu):
     r"""
@@ -1512,9 +1548,9 @@ def is_taylor_twograph_srg(int v,int k,int l,int mu):
 
 def is_switch_skewhad(int v, int k, int l, int mu):
     r"""
-    Test whether some `switch skewhad^2+*` is `(v,k,\lambda,\mu)`-strongly regular.
+    Test whether some ``switch skewhad^2+*`` is `(v,k,\lambda,\mu)`-strongly regular.
 
-    The `switch skewhad^2+*` graphs appear on `Andries Brouwer's database
+    The ``switch skewhad^2+*`` graphs appear on `Andries Brouwer's database
     <http://www.win.tue.nl/~aeb/graphs/srg/srgtab.html>`__ and are built by
     adding an isolated vertex to the complement of
     :func:`~sage.graphs.graph_generators.GraphGenerators.SquaredSkewHadamardMatrixGraph`,
@@ -2851,6 +2887,7 @@ def strongly_regular_graph(int v,int k,int l,int mu=-1,bint existence=False,bint
                       is_haemers,
                       is_cossidente_penttila,
                       is_mathon_PC_srg,
+                      is_muzychuk_S6,
                       is_switch_skewhad]
 
     # Going through all test functions, for the set of parameters and its