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

    module Math.Algebra.Group.SchreierSims where

    import qualified Data.List as L
    import Data.Maybe (isNothing, isJust)
    import qualified Data.Set as S
    import qualified Data.Map as M
    import Math.Algebra.Group.PermutationGroup hiding (elts, order, orderBSGS, gens, isMember, isSubgp, isNormal, reduceGens, normalClosure, commutatorGp, derivedSubgp)
    import Math.Common.ListSet (toListSet)
    import Math.Core.Utils hiding (elts)

    -- COSET REPRESENTATIVES FOR STABILISER OF A POINT

    -- Given a group G = <gs>, and a point x, find (right) coset representatives for Gx (stabiliser of x) in G
    -- In other words, for each x' in the orbit of x under G, we find a g <- G taking x to x'
    -- The code is similar to the code for calculating orbits, but modified to keep track of the group elements that we used to get there
    cosetRepsGx gs x = cosetRepsGx' gs M.empty (M.singleton x 1) where
        cosetRepsGx' gs interior boundary
            | M.null boundary = interior
            | otherwise =
                let interior' = M.union interior boundary
                    boundary' = M.fromList [(p .^ g, h*g) | g <- gs, (p,h) <- M.toList boundary] M.\\ interior'
                in cosetRepsGx' gs interior' boundary'

    -- SCHREIER GENERATORS

    -- toSet xs = (map head . group . sort) xs

    -- Generators for Gx, the stabiliser of x, given that G is generated by gs, and rs is a set of coset representatives for Gx in G.
    -- Schreier's Lemma states that if H < G = <S>, and R is a set of coset reps for H in G
    -- then H is generated by { rs(rs)*^-1 | r <- R, s <- S } (where * means "the coset representative of")
    -- In particular, with H = Gx, this gives us a way of finding a set of generators for Gx
    schreierGeneratorsGx (x,rs) gs = L.nub $ filter (/= 1) [schreierGenerator r g | r <- M.elems rs,  g <- gs]
        where schreierGenerator r g = let h = r * g
                                          h' = rs M.! (x .^ h)
                                      in h * inverse h'

    -- SCHREIER-SIMS ALGORITHM

    -- Given a list of right transversals for a stabiliser chain, sift a group element through it
    -- Note, this version assumes the base is non-redundant
    sift _ 1 = Nothing
    sift ((b,t):bts) g = case M.lookup (b .^ g) t of
                         Nothing -> Just g
                         -- Nothing -> sift bts g -- if we allow redundant levels
                         Just h -> sift bts (g * inverse h)
    sift [] g = Just g -- g == 1 case already caught above

    -- findBase gs = minimum $ concatMap supp gs
    findBase gs = minimum $ map minsupp gs
    {-
    -- Find base and strong generating set using Schreier-Sims algorithm
    bsgs gs | all (/= 1) gs = map fst $ ss [newLevel gs] []

    newLevel s = let b = findBase s
                     t = cosetRepsGx s b
                 in ((b,t),s)

    ss (bad@((b,t),s):bads) goods =
        let bts = map fst goods
            sgs = schreierGeneratorsGx (b,t) s
            siftees = filter isJust $ map (sift bts) sgs
        in  if null siftees
            then ss bads (bad:goods)
            else let Just h = head siftees
                 in if null goods
                    then ss (newLevel [h] : bad : bads) []
                    else let ((b_,t_),s_) = head goods
                             s' = h:s_
                             t' = cosetRepsGx s' b_
                         in ss (((b_,t'),s') : bad : bads) (tail goods)
    ss [] goods = goods
    -}

    -- |Given generators for a permutation group, return a strong generating set.
    -- The result is calculated using Schreier-Sims algorithm, and is relative to the base implied by the Ord instance
    sgs :: (Ord a, Show a) => [Permutation a] -> [Permutation a]
    sgs gs = toListSet $ concatMap snd $ ss bs gs
        where bs = toListSet $ concatMap supp gs

    -- Find base and strong generating set using Schreier-Sims algorithm
    -- !! This function is poorly named - it actually finds you a base and sets of transversals
    -- This version guarantees to use bases in order
    bsgs gs = bsgs' bs gs
        where bs = toListSet $ concatMap supp gs

    -- This version lets you pass in bases in the order you want them (or [], and it will find its own)
    bsgs' bs gs = map fst $ ss bs gs

    -- For example, bsgs (_A 5) uses [1,2,3] as the bases, but bsgs' [] (_A 5) uses [1,3,2]

    newLevel (b:bs) s = (bs, newLevel' b s)
    newLevel [] s = ([], newLevel' b s) where b = findBase s
    newLevel' b s = ((b,t),s) where t = cosetRepsGx s b

    ss bs gs = ss' bs' [level] []
        where (bs',level) = newLevel bs $ filter (/=1) gs

    ss' bs (bad@((b,t),s):bads) goods =
        let bts = map fst goods
            sgs = schreierGeneratorsGx (b,t) s
            siftees = filter isJust $ map (sift bts) sgs
        in  if null siftees
            then ss' bs bads (bad:goods)
            else let Just h = head siftees
                 in if null goods
                    then let (bs', level) = newLevel bs [h]
                         in ss' bs' (level : bad : bads) []
                    else let ((b_,t_),s_) = head goods
                             s' = h:s_
                             t' = cosetRepsGx s' b_
                         in ss' bs (((b_,t'),s') : bad : bads) (tail goods)
    ss' _ [] goods = goods
    {-
    extendbsgs [] g = bsgs [g]
    extendbsgs (((b,t),s):bts) g = ss (((b,t),g:s):bts) []

    bsgs' gs = map fst $ foldl extendbsgs [] gs
    -}

    -- The above is written for simplicity.
    -- Its efficiency could be improved by incrementally updating the transversals,
    -- and keeping track of Schreier generators we have already tried.
    -- (Remember to add new Schreier generators every time the generating set or transversal is augmented.)

    -- USING THE SCHREIER-SIMS TRANSVERSALS

    isMemberBSGS bts g = isNothing $ sift bts g

    -- By Lagrange's thm, every g <- G can be written uniquely as g = r_m ... r_1 (Seress p56)
    -- Note that we have to reverse the list of coset representatives
    eltsBSGS bts = map (product . reverse) (cartProd ts)
    	where ts = map (M.elems . snd) bts

    cartProd (set:sets) = [x:xs | x <- set, xs <- cartProd sets]
    cartProd [] = [[]]

    orderBSGS bts = product (map (toInteger . M.size . snd) bts)

    -- |Given generators for a group, determine whether a permutation is a member of the group, using Schreier-Sims algorithm
    isMember :: (Ord t, Show t) => [Permutation t] -> Permutation t -> Bool
    isMember gs h = isMemberBSGS (bsgs gs) h

    -- |Given generators for a group, return a (sorted) list of all elements of the group, using Schreier-Sims algorithm
    elts :: (Ord t, Show t) => [Permutation t] -> [Permutation t]
    elts [] = [1]
    elts gs = eltsBSGS $ bsgs gs

    -- |Given generators for a group, return the order of the group (the number of elements), using Schreier-Sims algorithm
    order :: (Ord t, Show t) => [Permutation t] -> Integer
    order [] = 1
    order gs = orderBSGS $ bsgs gs

    isSubgp hs gs = all (isMemberBSGS gs') hs
        where gs' = bsgs gs

    isNormal hs gs = hs `isSubgp` gs
                  && all (isMemberBSGS hs') [h~^g | h <- hs, g <- gs]
        where hs' = bsgs hs

    index gs hs = order gs `div` order hs

    -- given list of generators, try to find a shorter list
    reduceGens gs = fst $ reduceGensBSGS (filter (/=1) gs)

    reduceGensBSGS (g:gs) = reduceGens' ([g],bsgs [g]) gs where
        reduceGens' (gs,bsgsgs) (h:hs) =
            if isMemberBSGS bsgsgs h
            then reduceGens' (gs,bsgsgs) hs
            else reduceGens' (h:gs, bsgs $ h:gs) hs
        reduceGens' (gs,bsgsgs) [] = (reverse gs,bsgsgs)
    reduceGensBSGS [] = ([],[])

    -- normal closure of H in G
    -- for efficiency, should be called with gs and hs already reduced sets of generators
    normalClosure gs hs = reduceGens $ hs ++ [h ~^ g | h <- hs, g <- gs']
        where gs' = gs ++ map inverse gs

    -- commutator gp of H and K
    commutatorGp hs ks = normalClosure (hsks) [h^-1 * k^-1 * h * k | h <- hs', k <- ks']
        where hs' = reduceGens hs
              ks' = reduceGens ks
              hsks = reduceGens (hs' ++ ks')
              -- no point processing more potential generators than we have to

    -- derived subgroup (or commutator subgroup)
    derivedSubgp gs = normalClosure gs' [g^-1 * h^-1 * g * h | g <- gs', h <- gs']
        where gs' = reduceGens gs
        -- == commutatorGp gs gs

    {-
    isPerfect gs = order gs == order (derivedSubgp gs)
    -- We compare orders rather than the generators themselves, because derivedSubgp will usually find different generators

    derivedSeries gs = derivedSeries' (gs, order gs) where
        derivedSeries' ([],1) = [[]]
        derivedSeries' (hs, orderhs) =
            let hs' = derivedSubgp hs
                orderhs' = order hs'
            in if orderhs' == orderhs
               then [hs]
               else hs : derivedSeries' (hs',orderhs')

    lowerCentralSeries gs = lowerCentralSeries' (gs, order gs) where
        lowerCentralSeries' ([],1) = [[]]
        lowerCentralSeries' (hs, orderhs) =
            let hs' = commutatorGp gs hs
                orderhs' = order hs'
            in if orderhs' == orderhs
               then [hs]
               else hs : lowerCentralSeries' (hs',orderhs')
    -}