-- Copyright (c) 2008-2015, David Amos. All rights reserved.

-- |A module defining various strongly regular graphs, including the Clebsch, Hoffman-Singleton, Higman-Sims, and McLaughlin graphs.
--
-- A strongly regular graph with parameters (n,k,lambda,mu) is a (simple) graph with n vertices,
-- in which the number of common neighbours of x and y is k, lambda or mu according as whether
-- x and y are equal, adjacent, or non-adjacent. (In particular, it is a k-regular graph.)
--
-- Strongly regular graphs are highly symmetric, and have large automorphism groups.
module Math.Combinatorics.StronglyRegularGraph where

import qualified Data.List as L
import Data.Maybe (isJust)
import qualified Data.Map as M
import qualified Data.Set as S

import Math.Common.ListSet
import Math.Core.Utils (combinationsOf)
import Math.Algebra.Group.PermutationGroup hiding (P)
import Math.Algebra.Group.SchreierSims as SS
import Math.Combinatorics.Graph as G hiding (G)
import Math.Combinatorics.GraphAuts
import Math.Combinatorics.Design as D
import Math.Algebra.LinearAlgebra -- hiding (t)
import Math.Algebra.Field.Base -- for F2
import Math.Combinatorics.FiniteGeometry

-- Sources
-- Godsil & Royle, Algebraic Graph Theory
-- Cameron & van Lint, Designs, Graphs, Codes and their Links
-- van Lint & Wilson, A Course in Combinatorics, 2nd ed


-- STRONGLY REGULAR GRAPHS

-- strongly regular graphs
srgParams g
    | null es = error "srgParams: not defined for null graph"
    | null es' = error "srgParams: not defined for complete graph"
    | otherwise =
        if all (==k) ks && all (==lambda) ls && all (==mu) ms
        then Just (n,k,lambda,mu)
        else Nothing
    where vs = vertices g
          n = length vs
          es = edges g
          es' = combinationsOf 2 vs \\ es -- the non-edges
          k:ks = map (valency g) vs
          lambda:ls = map (length . commonNbrs) es  -- common neighbours of adjacent vertices
          mu:ms = map (length . commonNbrs) es' -- common neighbours of non-adjacent vertices
          commonNbrs [v1,v2] = (nbrs_g M.! v1) `intersect` (nbrs_g M.! v2)
          nbrs_g = M.fromList [ (v, nbrs g v) | v <- vs ]

isSRG g = isJust $ srgParams g


-- SIMPLE EXAMPLES

-- Triangular graph - van Lint & Wilson p262
-- http://mathworld.wolfram.com/TriangularGraph.html
t' m = G.to1n $ t m

t m | m >= 4 = graph (vs,es) where
    vs = combinationsOf 2 [1..m]
    es = [ [v,v'] | v <- vs, v' <- dropWhile (<= v) vs, not (disjoint v v')]
-- This is just lineGraph (k m), by another name


-- Lattice graph - van Lint & Wilson p262
-- http://mathworld.wolfram.com/LatticeGraph.html
l2' m = G.to1n $ l2 m

l2 m | m >= 2 = graph (vs,es) where
    vs = [ (i,j) | i <- [1..m], j <- [1..m] ]
    es = [ [v,v'] | v@(i,j) <- vs, v'@(i',j') <- dropWhile (<= v) vs, i == i' || j == j']
-- This is lineGraph (kb m m)
-- Automorphism group is Sm * Sm * C2
-- via i -> ig, j -> jg, i <-> j

paleyGraph fq | length fq `mod` 4 == 1 = graph (vs,es) where
    vs = fq
    qs = set [x^2 | x <- vs] \\ [0] -- the non-zero squares in Fq
    es = [ [x,y] | x <- vs, y <- vs, x < y, (x-y) `elem` qs]


-- CLEBSCH GRAPH
-- van Lint & Wilson, p263

clebsch' = G.to1n clebsch

clebsch = graph (vs,es) where
    vs = L.sort $ filter (even . length) $ powerset [1..5]
    es = [ [v,v'] | v <- vs, v' <- dropWhile (<= v) vs, length (symDiff v v') == 4]

-- Alternative construction from Cameron & van Lint p106
clebsch2 = graph (vs,es) where
    D xs bs = pairDesign 5
    vs = [C] ++ [P x | x <- xs] ++ [B b | b <- bs]
    es = L.sort $ [ [B a, B b]  | a <- bs, b <- dropWhile (<=a) bs, disjoint a b]
               ++ [ [P p, B b] | b <- bs, p <- b]
               ++ [ [C, P p]   | p <- xs ]


-- HOFFMAN-SINGLETON GRAPH
-- Cameron, Permutation Groups, p79ff
-- Godsil & Royle, p92ff
-- Aut group is U3(5).2 (Atlas p34)

triples = combinationsOf 3 [1..7]

heptads = [ [a,b,c,d,e,f,g] | a <- triples,
                              b <- triples, a < b, meetOne b a,
                              c <- triples, b < c, all (meetOne c) [a,b],
                              d <- triples, c < d, all (meetOne d) [a,b,c],
                              e <- triples, d < e, all (meetOne e) [a,b,c,d],
                              f <- triples, e < f, all (meetOne f) [a,b,c,d,e],
                              g <- triples, f < g, all (meetOne g) [a,b,c,d,e,f],
                              foldl intersect [1..7] [a,b,c,d,e,f,g] == [] ]
    where meetOne x y = length (intersect x y) == 1
    -- each pair of triples meet in exactly one point, and there is no point in all of them - Godsil & Royle p69
    -- (so these are the projective planes over 7 points)

plane +^ g = L.sort [line -^ g | line <- plane]

plane +^^ gs = orbit (+^) plane gs
-- plane +^^ gs = closure [plane] [ +^ g | g <- gs ]

hoffmanSingleton' = G.to1n hoffmanSingleton

hoffmanSingleton = graph (vs,es) where
    h = head heptads
    hs = h +^^ _A 7 -- an A7 orbit of a heptad
    vs = map Left hs ++ map Right triples
    es = [ [Left h, Right t] | h <- hs, t <- triples, t `elem` h]
      ++ [ [Right t, Right t'] | t <- triples, t' <- dropWhile (<= t) triples, t `disjoint` t']

-- induced action of A7 on Hoffman-Singleton graph
inducedA7 g = fromPairs [(v, v ~^ g) | v <- vs] where
    vs = vertices hoffmanSingleton
    (Left h) ~^ g = Left (h +^ g)
    (Right t) ~^ g = Right (t -^ g)

hsA7 = toSn $ map inducedA7 $ _A 7


-- GEWIRTZ GRAPH
-- van Lint & Wilson p266-7
-- (also called Sims-Gewirtz graph)

gewirtz' = G.to1n gewirtz

gewirtz = graph (vs,es) where
    vs = [xs | xs <- blocks s_3_6_22, 22 `notElem` xs]
    -- The 21 blocks of S(3,6,22) which contain 22 are the lines of PG(2,4) (projective plane over F4)
    -- The 56 blocks which don't are hyperovals in this plane. They form a 2-(21,6,4) design.
    es = [ [v,v'] | v <- vs, v' <- dropWhile (<= v) vs, length (v `intersect` v') == 0]


-- HIGMAN-SIMS GRAPH
-- Aut group is HS.2, where HS is the Higman-Sims sporadic simple group

data DesignVertex = C | P Integer | B [Integer] deriving (Eq,Ord,Show)

higmanSimsGraph' = G.to1n higmanSimsGraph

-- Cameron & van Lint, p107
higmanSimsGraph = graph (vs,es) where
    D xs bs = s_3_6_22
    vs = [C] ++ [P x | x <- xs] ++ [B b | b <- bs]
    es = L.sort $ [ [B a, B b]  | a <- bs, b <- dropWhile (<=a) bs, disjoint a b]
               ++ [ [P p, B b] | b <- bs, p <- b]
               ++ [ [C, P p]   | p <- xs ]
    -- s_3_6_22' = blocks s_3_6_22

-- There is an induced action of M22 on Higman Sims graph

-- induced action of M22 on Higman-Sims graph
inducedM22 g = fromPairs [(v, v ~^ g) | v <- vs] where
    -- G vs _ = higmanSimsGraph'
    vs = vertices higmanSimsGraph
    (B b) ~^ g = B (b -^ g)
    (P p) ~^ g = P (p .^ g)
    C ~^ _ = C

higmanSimsM22 = toSn $ map inducedM22 $ m22sgs
-- all (isGraphAut higmanSimsGraph) higmanSimsM22

-- M22 is one point stabilizer (of C)

-- HS.2, where HS is Higman-Sims sporadic group
_HS2 = SS.reduceGens $ graphAuts higmanSimsGraph
-- (It will actually find 11 strong generators, but the first 4 are sufficient to generate the group)

_HS = SS.derivedSubgp _HS2



-- SYMPLECTIC GRAPHS

-- Godsil & Royle p242
sp2 r = graph (vs,es) where
    vs = tail $ ptsAG (2*r) f2 -- all non-zero pts in F2^2r
    es = [ [u,v] | [u,v] <- combinationsOf 2 vs, u <*>> n <.> v == 1] -- uT N v == 1, ie vectors adjacent if non-orthogonal
    n = fMatrix (2*r) (\i j -> if abs (i-j) == 1 && even (max i j) then 1 else 0) -- matrix defining a symplectic form

sp n | even n = sp2 (n `div` 2)


-- TWO GRAPHS AND SWITCHING

-- SCHLAFLI GRAPH
-- An srg(27,16,10,8)
-- Has geometric interpretation in terms of 27 lines on general cubic surface in projective 3-space
-- Aut group is G.2 where G = U4(2) = S4(3) (Atlas p26)
-- (G.2 is also the Weyl group of E6 - don't know if there's any connection)

-- Godsil & Royle p254ff
switch g us | us `D.isSubset` vs = graph (vs, L.sort switchedes) where
    vs = vertices g
    us' = vs L.\\ us -- complement of us in vs
    es = edges g
    es' = S.fromList es
    switchedes = [e | e@[v1,v2] <- es, (v1 `elem` us) == (v2 `elem` us)]
                 -- edges within us or its complement are unchanged
              ++ [ L.sort [v1,v2] | v1 <- us, v2 <- us', L.sort [v1,v2] `S.notMember` es']
                 -- edges between us and its complement are switched

-- Godsil & Royle p259
schlafli' = G.to1n schlafli

schlafli = graph (vs,es') where
    g = lineGraph $ k 8
    v:vs = vertices g
    es = edges g
    gswitched = switch g (nbrs g v) -- switch off the vertex v
    es' = edges gswitched


-- MCLAUGHLIN GRAPH
-- Aut group is McL.2, where McL is the McLaughlin sporadic simple group
-- http://people.csse.uwa.edu.au/gordon/constructions/mclaughlin/
-- http://mathworld.wolfram.com/McLaughlinGraph.html

mcLaughlin' = G.to1n mcLaughlin

mcLaughlin = graph (vs',es') where
    D xs bs = s_4_7_23
    vs = map P xs ++ map B bs
    es = [ [P x, B b] | x <- xs, b <- bs, x `notElem` b]
      ++ [ [B b1, B b2] | b1 <- bs, b2 <- bs, b1 < b2, length (b1 `intersect` b2) == 1]
    g276 = graph (vs,es)
    g276switched = switch g276 (nbrs g276 (P 0))
    P 0 : vs' = vs -- drop P 0 as it's now not connected
    es' = edges g276switched

_McL2 = SS.reduceGens $ graphAuts mcLaughlin
-- finds 14 auts - but takes half an hour (interpreted) to do so
-- in fact just the first 2 are sufficient to generate the group

_McL = SS.derivedSubgp $ _McL2

{-
-- TWO GRAPH ON 276 VERTICES
-- Has Conway's .3 as automorphism group

-- Godsil & Royle p260ff

twoGraph276 =
    let nt = D.incidenceMatrix s_4_7_23
        n = L.transpose nt -- Godsil & Royle do incidence matrix the other way round to us
        s = L.transpose $
            (j 23 23 <<->> i 23)      +|+ (j 23 253 <<->> 2 *>> n)
                                      ++
            (j 253 23 <<->> 2 *>> nt) +|+ (nt <<*>> n <<->> 5 *>> i 253 <<->> 2 *>> j 253 253)
        a = (map . map) (`div` 2) (j 276 276 <<->> i 276 <<->> s)
    in fromAdjacencyMatrix a
    where j r c = replicate r (replicate c 1)
          i = idMx
          (+|+) = zipWith (++)

-- Its automorphism group *as a two-graph* is .3 (Co3)
-- But its aut group as a graph is only M23

twoGraph276' = graph (vs,es) where
    D xs bs = s_4_7_23
    vs = map P xs ++ map B bs
    es = [ [P x, B b] | x <- xs, b <- bs, x `notElem` b]
      ++ [ [B b1, B b2] | b1 <- bs, b2 <- bs, b1 < b2, length (b1 `intersect` b2) == 1]
-- !! This isn't isomorphic to twoGraph276
-- (Perhaps it is in the same switching class though)
-- We can obtain McLaughlin graph from this by switching in neighbourhood of P 0
-}