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

    {-# LANGUAGE NoMonomorphismRestriction #-}
    -- Because unRight defined point-free

    module Math.Algebra.Group.Subquotients where

    import qualified Data.List as L
    import qualified Data.Map as M

    import Math.Common.ListSet
    import Math.Algebra.Group.PermutationGroup hiding (ptStab, normalClosure)
    import Math.Algebra.Group.SchreierSims (cosetRepsGx)
    import Math.Algebra.Group.RandomSchreierSims

    -- Source: Seress, Permutation Group Algorithms

    isLeft (Left _) = True
    isLeft (Right _) = False

    isRight (Right _) = True
    isRight (Left _) = False

    unRight = fromPairs . map (\(Right a, Right b) -> (a,b)) . toPairs

    restrictLeft g = fromPairs [(a,b) | (Left a, Left b) <- toPairs g]
    -- note that this is doing a filter - taking only the left part of the action - and a map, unLefting

    -- pointwise stabiliser of xs
    ptStab gs delta = map unRight $ dropWhile (isLeft . minsupp) $ sgs gs' where
        gs' = [ (fromPairs . map (\(a,b) -> (lr a, lr b)) . toPairs) g | g <- gs]
        lr x = if x `elem` delta then Left x else Right x

    {-
    -- !! NEXT TWO FUNCTIONS NOT TESTED
    -- Need some meaningful examples of homomorphisms
    -- eg Sn -> Sym(k-subsets of n)
    -- restrict to a transitive constituent
    -- blocks

    -- Given generators gs for a group G, and f : G -> H a homomorphism,
    -- return the "semi-diagonal" subgroup [(f g, g) | g <- gs] of f(G) * G
    homomorphismConstruction :: (Ord a, Ord b) => [Permutation a] -> (Permutation a -> Permutation b) -> [Permutation (Either b a)]
    homomorphismConstruction gs f = [lift g | g <- gs] where
        lift g = fromPairs $ [(Right x, Right y) | (x,y) <- toPairs g] ++ [(Left x', Left y') | (x',y') <- toPairs (f g)]

    ker gs f = ks where
        gbar = homomorphismConstruction gs f
        gs' = sgs gbar
        ks' = dropWhile (\h -> isLeft $ minsupp h) gs' -- !! should filter isRight - sgs might not be in order
        ks = map unRight ks'
        unRight = fromPairs . map (\(Right a, Right b) -> (a,b)) . toPairs
    -}

    isTransitive :: (Ord t) => [Permutation t] -> Bool
    isTransitive gs = length (orbits gs) == 1

    -- TRANSITIVE CONSTITUENTS

    {-
    -- find largest composition factor of a group which is not transitive
    -- we do this by taking the smallest orbit delta,
    -- then constructing the homomorphism G -> Sym(delta)
    -- and returning the kernel and the image
    factorNotTransitive gs = transitiveConstituentHomomorphism' gs delta where
        delta = smallest $ orbits gs
        sizeSorted lists = map snd $ L.sort $ [(length l, l) | l <- lists]
        smallest = head . sizeSorted
    -}

    -- Seress p81
    -- |Given a group gs and a transitive constituent ys, return the kernel and image of the transitive constituent homomorphism.
    -- That is, suppose that gs acts on a set xs, and ys is a subset of xs on which gs acts transitively.
    -- Then the transitive constituent homomorphism is the restriction of the action of gs to an action on the ys.
    transitiveConstituentHomomorphism
      :: (Ord a, Show a) =>
         [Permutation a] -> [a] -> ([Permutation a], [Permutation a])
    transitiveConstituentHomomorphism gs delta
        | delta == closure delta [(.^ g) | g <- gs] -- delta is closed under action of gs, hence a union of orbits
            = transitiveConstituentHomomorphism' gs delta

    transitiveConstituentHomomorphism' gs delta = (ker, im) where
        gs' = sgs $ map (fromPairs . map (\(a,b) -> (lr a, lr b)) . toPairs) gs
        -- as delta is a transitive constituent, we will always have a and b either both Left or both Right
        lr x = if x `elem` delta then Left x else Right x
        ker = map unRight $ dropWhile (isLeft . minsupp) gs' -- pointwise stabiliser of delta
        im = map restrictLeft $ takeWhile (isLeft . minsupp) gs' -- restriction of the action to delta

    -- BLOCKS OF IMPRIMITIVITY

    -- Holt p83ff (and also Seress p107ff)
    -- Find a minimal block containing ys. ys are assumed to be sorted.
    minimalBlock gs ys@(y1:yt) = minimalBlock' p yt gs where
        xs = foldl union [] $ map supp gs
        p = M.fromList $ [(yi,y1) | yi <- ys] ++ [(x,x) | x <- xs \\ ys]
        minimalBlock' p (q:qs) (h:hs) =
            let r = p M.! q         -- representative of class containing q
                k = p M.! (q .^ h)  -- rep of class (q^h)
                l = p M.! (r .^ h)  -- rep of class (r^h)
            in if k /= l -- then we need to merge the classes
               then let p' = M.map (\x -> if x == l then k else x) p
                        qs' = qs ++ [l]
                    in minimalBlock' p' (q:qs') hs
               else minimalBlock' p (q:qs) hs
        minimalBlock' p (q:qs) [] = minimalBlock' p qs gs
        minimalBlock' p [] _ =
            let reps = toListSet $ M.elems p
            in L.sort [ filter (\x -> p M.! x == r) xs | r <- reps ]
    -- Because the support of the permutations is not constrained to be [1..n], we have to use a map instead of an array
    -- This probably affects the complexity, but isn't a problem in practice

    -- |Given a transitive group gs, find all non-trivial block systems. That is, if gs act on xs,
    -- find all the ways that the xs can be divided into blocks, such that the gs also have a permutation action on the blocks
    blockSystems :: (Ord t) => [Permutation t] -> [[[t]]]
    blockSystems gs
        | isTransitive gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- xs ]
        | otherwise = error "blockSystems: not transitive"
        where x:xs = foldl union [] $ map supp gs

    -- |A more efficient version of blockSystems, if we have an sgs
    blockSystemsSGS :: (Ord a) => [Permutation a] -> [[[a]]]
    blockSystemsSGS gs = toListSet $ filter (/= [x:xs]) $ map (minimalBlock gs) [ [x,x'] | x' <- rs ]
        where x:xs = foldl union [] $ map supp gs
              hs = filter (\g -> x < minsupp g) gs -- sgs for stabiliser Gx
              os = orbits hs
              rs = map head os ++ (xs \\ L.sort (concat os)) -- orbit representatives, including singleton cycles
    -- Perhaps we could have a function which just returns orbit reps for stabiliser

    -- eg for D 10, the stabiliser of 1 is [[2,6],[3,5]] - we need to make sure we don't forget 4

    -- If we didn't have an SGS, we could try to randomly generate a few elts of stabiliser Gx, as that would still be better than nothing
    -- see Holt RandomStab function

    -- |A permutation group is primitive if it has no non-trivial block systems
    isPrimitive :: (Ord t) => [Permutation t] -> Bool
    isPrimitive gs = null (blockSystems gs)

    isPrimitiveSGS :: (Ord a) => [Permutation a] -> Bool
    isPrimitiveSGS gs = null (blockSystemsSGS gs)

    -- There are other optimisations we haven't done
    -- see Holt p86

    -- |Given a transitive group gs, and a block system for gs, return the kernel and image of the block homomorphism
    -- (the homomorphism onto the action of gs on the blocks)
    blockHomomorphism
      :: (Ord t, Show t) =>
         [Permutation t] -> [[t]] -> ([Permutation t], [Permutation [t]])
    blockHomomorphism gs bs
        | bs == closure bs [(-^ g) | g <- gs] -- bs is closed under action of gs
            = blockHomomorphism' gs bs

    blockHomomorphism' gs bs = (ker,im) where
        gs' = sgs $ map lr gs
        lr g = fromPairs $ [(Left b, Left $ b -^ g) | b <- bs] ++ [(Right x, Right y) | (x,y) <- toPairs g]
        ker = map unRight $ dropWhile (isLeft . minsupp) gs' -- stabiliser of the blocks
        im = map restrictLeft $ takeWhile (isLeft . minsupp) gs' -- restriction to the action on blocks

    -- Note that there is a slightly more efficient way to calculate block homomorphism,
    -- but requires change of base algorithm which we haven't implemented yet

    -- NORMAL CLOSURE

    -- Seress 115
    -- Given G, H < Sym(Omega) return <H^G> (the normal closure)
    normalClosure gs hs = map unRight $ dropWhile (isLeft . minsupp) $ sgs ks where
        xs = foldl union [] $ map supp $ gs ++ hs
        ds = map diag gs -- {(g,g) | g <- G}
        diag g = fromPairs $ concat [ [(Left x, Left y) , (Right x, Right y)] | (x,y) <- toPairs g]
        hsR = map inR hs -- {(1,h) | h <- H}
        inR h = fromPairs [(Right x, Right y) | (x,y) <- toPairs h]
        ks = ds ++ hsR

    -- Seress 116
    -- Given G, H < Sym(Omega) return <H^G> `intersection` G
    intersectionNormalClosure gs hs = map unRight $ dropWhile (isLeft . minsupp) $ sgs ks where
        xs = foldl union [] $ map supp $ gs ++ hs
        ds = map diag gs -- {(g,g) | g <- G}
        diag g = fromPairs $ concat [ [(Left x, Left y) , (Right x, Right y)] | (x,y) <- toPairs g]
        hsL = map inL hs -- {(h,1) | h <- H}
        inL h = fromPairs [(Left x, Left y) | (x,y) <- toPairs h]
        ks = ds ++ hsL

    -- CENTRALISER IN THE SYMMETRIC GROUP

    -- Centralizer of G in Sym(X) - transitive case
    centralizerSymTrans gs = filter (/= 1) $ centralizerSymTrans' [] fix_g_a where
        xs@(a:_) = foldl union [] $ map supp gs
        ss = sgs gs
        g_a = dropWhile ( (==a) . minsupp ) ss -- pt stabiliser of a
        fix_g_a = xs \\ (foldl union [] $ map supp g_a) -- the pts fixed by stabiliser of a
        reps_a = cosetRepsGx gs a
        -- xs = M.keys reps_a
        centralizingElt b = fromPairs [ let g = reps_a M.! x in (x, b .^ g) | x <- xs ]
        centralizerSymTrans' ls (r:rs) =
            let c = centralizingElt r
            in c : centralizerSymTrans' (c:ls) (rs \\ orbitP (c:ls) a)
        centralizerSymTrans' _ [] = []