```-- 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 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
-- A transitive constituent homomorphism is the restriction of G <= Sym(omega) to an orbit delta <= omega
-- This function returns the kernel and the image
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

-- Find all block systems
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

-- More efficient version if we have an sgs
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

isPrimitive gs = null (blockSystems gs)

isPrimitiveSGS gs = null (blockSystemsSGS gs)

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

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' _ [] = []

```