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

module Math.Algebra.Group.RandomSchreierSims where

import System.Random
import Data.List as L
import qualified Data.Map as M
import Data.Maybe

import Control.Monad
import Data.Array.MArray
import Data.Array.IO
import System.IO.Unsafe

import Math.Common.ListSet (toListSet)
import Math.Core.Utils hiding (elts)
import Math.Algebra.Group.PermutationGroup
import Math.Algebra.Group.SchreierSims (sift, cosetRepsGx, ss')


testProdRepl :: IO ()
testProdRepl = do (r :: Int
r,xs :: IOArray Int (Permutation Integer)
xs) <- [Permutation Integer]
-> IO (Int, IOArray Int (Permutation Integer))
forall a.
(Ord a, Show a) =>
[Permutation a] -> IO (Int, IOArray Int (Permutation a))
initProdRepl ([Permutation Integer]
 -> IO (Int, IOArray Int (Permutation Integer)))
-> [Permutation Integer]
-> IO (Int, IOArray Int (Permutation Integer))
forall a b. (a -> b) -> a -> b
$ Integer -> [Permutation Integer]
forall a. Integral a => a -> [Permutation a]
_D 10
                  [Maybe (Permutation Integer)]
hs <- Int
-> IO (Maybe (Permutation Integer))
-> IO [Maybe (Permutation Integer)]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM 20 (IO (Maybe (Permutation Integer))
 -> IO [Maybe (Permutation Integer)])
-> IO (Maybe (Permutation Integer))
-> IO [Maybe (Permutation Integer)]
forall a b. (a -> b) -> a -> b
$ (Int, IOArray Int (Permutation Integer))
-> IO (Maybe (Permutation Integer))
forall a.
(Ord a, Show a) =>
(Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl (Int
r,IOArray Int (Permutation Integer)
xs)
                  (Maybe (Permutation Integer) -> IO ())
-> [Maybe (Permutation Integer)] -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Maybe (Permutation Integer) -> IO ()
forall a. Show a => a -> IO ()
print [Maybe (Permutation Integer)]
hs

-- Holt p69-71
-- Product replacement algorithm for generating uniformly distributed random elts of a black box group

initProdRepl :: (Ord a, Show a) => [Permutation a] -> IO (Int, IOArray Int (Permutation a))
initProdRepl :: [Permutation a] -> IO (Int, IOArray Int (Permutation a))
initProdRepl gs :: [Permutation a]
gs =
    let n :: Int
n = [Permutation a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Permutation a]
gs
        r :: Int
r = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max 10 Int
n
        xs :: [Permutation a]
xs = (1Permutation a -> [Permutation a] -> [Permutation a]
forall a. a -> [a] -> [a]
:) ([Permutation a] -> [Permutation a])
-> [Permutation a] -> [Permutation a]
forall a b. (a -> b) -> a -> b
$ Int -> [Permutation a] -> [Permutation a]
forall a. Int -> [a] -> [a]
take Int
r ([Permutation a] -> [Permutation a])
-> [Permutation a] -> [Permutation a]
forall a b. (a -> b) -> a -> b
$ [[Permutation a]] -> [Permutation a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[Permutation a]] -> [Permutation a])
-> [[Permutation a]] -> [Permutation a]
forall a b. (a -> b) -> a -> b
$ [Permutation a] -> [[Permutation a]]
forall a. a -> [a]
repeat [Permutation a]
gs 
    in do IOArray Int (Permutation a)
xs' <- (Int, Int) -> [Permutation a] -> IO (IOArray Int (Permutation a))
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
(i, i) -> [e] -> m (a i e)
newListArray (0,Int
r) [Permutation a]
xs
          Int -> IO (Maybe (Permutation a)) -> IO ()
forall (m :: * -> *) a. Applicative m => Int -> m a -> m ()
replicateM_ 60 (IO (Maybe (Permutation a)) -> IO ())
-> IO (Maybe (Permutation a)) -> IO ()
forall a b. (a -> b) -> a -> b
$ (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
forall a.
(Ord a, Show a) =>
(Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl (Int
r,IOArray Int (Permutation a)
xs') -- perform initial mixing
          (Int, IOArray Int (Permutation a))
-> IO (Int, IOArray Int (Permutation a))
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
r,IOArray Int (Permutation a)
xs')

nextProdRepl :: (Ord a, Show a) => (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl :: (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl (r :: Int
r,xs :: IOArray Int (Permutation a)
xs) =
    do Int
s <- (Int, Int) -> IO Int
forall a. Random a => (a, a) -> IO a
randomRIO (1,Int
r)
       Int
t <- (Int, Int) -> IO Int
forall a. Random a => (a, a) -> IO a
randomRIO (1,Int
r)
       Int
u <- (Int, Int) -> IO Int
forall a. Random a => (a, a) -> IO a
randomRIO (0,3 :: Int)
       Maybe (Permutation a)
out <- IOArray Int (Permutation a)
-> Int -> Int -> Int -> IO (Maybe (Permutation a))
forall i (m :: * -> *) (a :: * -> * -> *) a b.
(MArray a a m, Ix i, Num i, Integral b, HasInverses a, Num a) =>
a i a -> i -> i -> b -> m (Maybe a)
updateArray IOArray Int (Permutation a)
xs Int
s Int
t Int
u
       Maybe (Permutation a) -> IO (Maybe (Permutation a))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Permutation a)
out

updateArray :: a i a -> i -> i -> b -> m (Maybe a)
updateArray xs :: a i a
xs s :: i
s t :: i
t u :: b
u =
    let (swap :: b
swap,invert :: b
invert) = b -> b -> (b, b)
forall a. Integral a => a -> a -> (a, a)
quotRem b
u 2 in
    if i
s i -> i -> Bool
forall a. Eq a => a -> a -> Bool
== i
t
    then Maybe a -> m (Maybe a)
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe a
forall a. Maybe a
Nothing
    else do
        a
x_0 <- a i a -> i -> m a
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray a i a
xs 0
        a
x_s <- a i a -> i -> m a
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray a i a
xs i
s
        a
x_t <- a i a -> i -> m a
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> m e
readArray a i a
xs i
t
        let x_s' :: a
x_s' = (b, b) -> a -> a -> a
forall a a a.
(HasInverses a, Num a, Num a, Num a, Eq a, Eq a) =>
(a, a) -> a -> a -> a
mult (b
swap,b
invert) a
x_s a
x_t
            x_0' :: a
x_0' = (b, Integer) -> a -> a -> a
forall a a a.
(HasInverses a, Num a, Num a, Num a, Eq a, Eq a) =>
(a, a) -> a -> a -> a
mult (b
swap,0) a
x_0 a
x_s'
        a i a -> i -> a -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray a i a
xs 0 a
x_0'
        a i a -> i -> a -> m ()
forall (a :: * -> * -> *) e (m :: * -> *) i.
(MArray a e m, Ix i) =>
a i e -> i -> e -> m ()
writeArray a i a
xs i
s a
x_s'
        Maybe a -> m (Maybe a)
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> Maybe a
forall a. a -> Maybe a
Just a
x_0')
    where mult :: (a, a) -> a -> a -> a
mult (swap :: a
swap,invert :: a
invert) a :: a
a b :: a
b = case (a
swap,a
invert) of
                                   (0,0) -> a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b
                                   (0,1) -> a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
ba -> Integer -> a
forall a b. (Num a, HasInverses a, Integral b) => a -> b -> a
^-1
                                   (1,0) -> a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
a
                                   (1,1) -> a
ba -> Integer -> a
forall a b. (Num a, HasInverses a, Integral b) => a -> b -> a
^-1 a -> a -> a
forall a. Num a => a -> a -> a
* a
a


-- Holt p97-8
-- Random Schreier-Sims algorithm, for finding strong generating set of permutation group

-- It's possible that the following code can be improved by introducing levels only as we need them?

-- |Given generators for a permutation group, return a strong generating set.
-- The result is calculated using random Schreier-Sims algorithm, so has a small (\<10^-6) chance of being incomplete.
-- The sgs is relative to the base implied by the Ord instance.
sgs :: (Ord a, Show a) => [Permutation a] -> [Permutation a]
sgs :: [Permutation a] -> [Permutation a]
sgs gs :: [Permutation a]
gs = [Permutation a] -> [Permutation a]
forall b. Ord b => [b] -> [b]
toListSet ([Permutation a] -> [Permutation a])
-> [Permutation a] -> [Permutation a]
forall a b. (a -> b) -> a -> b
$ (((a, Map a (Permutation a)), [Permutation a]) -> [Permutation a])
-> [((a, Map a (Permutation a)), [Permutation a])]
-> [Permutation a]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap ((a, Map a (Permutation a)), [Permutation a]) -> [Permutation a]
forall a b. (a, b) -> b
snd ([((a, Map a (Permutation a)), [Permutation a])]
 -> [Permutation a])
-> [((a, Map a (Permutation a)), [Permutation a])]
-> [Permutation a]
forall a b. (a -> b) -> a -> b
$ [Permutation a] -> [((a, Map a (Permutation a)), [Permutation a])]
forall k.
(Ord k, Show k) =>
[Permutation k] -> [((k, Map k (Permutation k)), [Permutation k])]
rss [Permutation a]
gs

rss :: [Permutation k] -> [((k, Map k (Permutation k)), [Permutation k])]
rss gs :: [Permutation k]
gs = IO [((k, Map k (Permutation k)), [Permutation k])]
-> [((k, Map k (Permutation k)), [Permutation k])]
forall a. IO a -> a
unsafePerformIO (IO [((k, Map k (Permutation k)), [Permutation k])]
 -> [((k, Map k (Permutation k)), [Permutation k])])
-> IO [((k, Map k (Permutation k)), [Permutation k])]
-> [((k, Map k (Permutation k)), [Permutation k])]
forall a b. (a -> b) -> a -> b
$
    do (r :: Int
r,xs :: IOArray Int (Permutation k)
xs) <- [Permutation k] -> IO (Int, IOArray Int (Permutation k))
forall a.
(Ord a, Show a) =>
[Permutation a] -> IO (Int, IOArray Int (Permutation a))
initProdRepl [Permutation k]
gs
       (Int, IOArray Int (Permutation k))
-> [((k, Map k (Permutation k)), [Permutation k])]
-> Integer
-> IO [((k, Map k (Permutation k)), [Permutation k])]
forall t a.
(Num t, Show a, Eq t, Ord a) =>
(Int, IOArray Int (Permutation a))
-> [((a, Map a (Permutation a)), [Permutation a])]
-> t
-> IO [((a, Map a (Permutation a)), [Permutation a])]
rss' (Int
r,IOArray Int (Permutation k)
xs) ([Permutation k] -> [((k, Map k (Permutation k)), [Permutation k])]
forall a k (t :: * -> *) a.
(Num a, Ord k, Foldable t) =>
t (Permutation k) -> [((k, Map k a), [a])]
initLevels [Permutation k]
gs) 0

rss' :: (Int, IOArray Int (Permutation a))
-> [((a, Map a (Permutation a)), [Permutation a])]
-> t
-> IO [((a, Map a (Permutation a)), [Permutation a])]
rss' (r :: Int
r,xs :: IOArray Int (Permutation a)
xs) levels :: [((a, Map a (Permutation a)), [Permutation a])]
levels i :: t
i
    | t
i t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== 25 = [((a, Map a (Permutation a)), [Permutation a])]
-> IO [((a, Map a (Permutation a)), [Permutation a])]
forall (m :: * -> *) a. Monad m => a -> m a
return [((a, Map a (Permutation a)), [Permutation a])]
levels -- stop if we've had 25 successful sifts in a row
    | Bool
otherwise = do Maybe (Permutation a)
g <- (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
forall a.
(Ord a, Show a) =>
(Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl (Int
r,IOArray Int (Permutation a)
xs)
                     let (changed :: Bool
changed,levels' :: [((a, Map a (Permutation a)), [Permutation a])]
levels') = [((a, Map a (Permutation a)), [Permutation a])]
-> Maybe (Permutation a)
-> (Bool, [((a, Map a (Permutation a)), [Permutation a])])
forall c.
Ord c =>
[((c, Map c (Permutation c)), [Permutation c])]
-> Maybe (Permutation c)
-> (Bool, [((c, Map c (Permutation c)), [Permutation c])])
updateLevels [((a, Map a (Permutation a)), [Permutation a])]
levels Maybe (Permutation a)
g
                     (Int, IOArray Int (Permutation a))
-> [((a, Map a (Permutation a)), [Permutation a])]
-> t
-> IO [((a, Map a (Permutation a)), [Permutation a])]
rss' (Int
r,IOArray Int (Permutation a)
xs) [((a, Map a (Permutation a)), [Permutation a])]
levels' (if Bool
changed then 0 else t
it -> t -> t
forall a. Num a => a -> a -> a
+1)
-- if we currently have an sgs for a subgroup of the group, then it must have index >= 2
-- so the chance of a random elt sifting to identity is <= 1/2

initLevels :: t (Permutation k) -> [((k, Map k a), [a])]
initLevels gs :: t (Permutation k)
gs = [((k
b,k -> a -> Map k a
forall k a. k -> a -> Map k a
M.singleton k
b 1),[]) | k
b <- [k]
bs]
    where bs :: [k]
bs = [k] -> [k]
forall b. Ord b => [b] -> [b]
toListSet ([k] -> [k]) -> [k] -> [k]
forall a b. (a -> b) -> a -> b
$ (Permutation k -> [k]) -> t (Permutation k) -> [k]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Permutation k -> [k]
forall a. Permutation a -> [a]
supp t (Permutation k)
gs

updateLevels :: [((c, Map c (Permutation c)), [Permutation c])]
-> Maybe (Permutation c)
-> (Bool, [((c, Map c (Permutation c)), [Permutation c])])
updateLevels levels :: [((c, Map c (Permutation c)), [Permutation c])]
levels Nothing = (Bool
False,[((c, Map c (Permutation c)), [Permutation c])]
levels) -- not strictly correct to increment count on a Nothing
updateLevels levels :: [((c, Map c (Permutation c)), [Permutation c])]
levels (Just g :: Permutation c
g) =
    case [(c, Map c (Permutation c))]
-> Permutation c -> Maybe (Permutation c)
forall k.
Ord k =>
[(k, Map k (Permutation k))]
-> Permutation k -> Maybe (Permutation k)
sift ((((c, Map c (Permutation c)), [Permutation c])
 -> (c, Map c (Permutation c)))
-> [((c, Map c (Permutation c)), [Permutation c])]
-> [(c, Map c (Permutation c))]
forall a b. (a -> b) -> [a] -> [b]
map ((c, Map c (Permutation c)), [Permutation c])
-> (c, Map c (Permutation c))
forall a b. (a, b) -> a
fst [((c, Map c (Permutation c)), [Permutation c])]
levels) Permutation c
g of
    Nothing -> (Bool
False, [((c, Map c (Permutation c)), [Permutation c])]
levels)
    -- Just 1 -> error "Just 1"
    Just g' :: Permutation c
g' -> (Bool
True, [((c, Map c (Permutation c)), [Permutation c])]
-> [((c, Map c (Permutation c)), [Permutation c])]
-> Permutation c
-> c
-> [((c, Map c (Permutation c)), [Permutation c])]
forall t.
Ord t =>
[((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
-> Permutation t
-> t
-> [((t, Map t (Permutation t)), [Permutation t])]
updateLevels' [] [((c, Map c (Permutation c)), [Permutation c])]
levels Permutation c
g' (Permutation c -> c
forall c. Permutation c -> c
minsupp Permutation c
g'))

updateLevels' :: [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
-> Permutation t
-> t
-> [((t, Map t (Permutation t)), [Permutation t])]
updateLevels' ls :: [((t, Map t (Permutation t)), [Permutation t])]
ls (r :: ((t, Map t (Permutation t)), [Permutation t])
r@((b :: t
b,t :: Map t (Permutation t)
t),s :: [Permutation t]
s):rs :: [((t, Map t (Permutation t)), [Permutation t])]
rs) h :: Permutation t
h b' :: t
b' =
    if t
b t -> t -> Bool
forall a. Eq a => a -> a -> Bool
== t
b'
    then [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
forall a. [a] -> [a]
reverse [((t, Map t (Permutation t)), [Permutation t])]
ls [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
forall a. [a] -> [a] -> [a]
++ ((t
b, [Permutation t] -> t -> Map t (Permutation t)
forall k. Ord k => [Permutation k] -> k -> Map k (Permutation k)
cosetRepsGx (Permutation t
hPermutation t -> [Permutation t] -> [Permutation t]
forall a. a -> [a] -> [a]
:[Permutation t]
s) t
b), Permutation t
hPermutation t -> [Permutation t] -> [Permutation t]
forall a. a -> [a] -> [a]
:[Permutation t]
s) ((t, Map t (Permutation t)), [Permutation t])
-> [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
forall a. a -> [a] -> [a]
: [((t, Map t (Permutation t)), [Permutation t])]
rs
    else [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
-> Permutation t
-> t
-> [((t, Map t (Permutation t)), [Permutation t])]
updateLevels' (((t, Map t (Permutation t)), [Permutation t])
r((t, Map t (Permutation t)), [Permutation t])
-> [((t, Map t (Permutation t)), [Permutation t])]
-> [((t, Map t (Permutation t)), [Permutation t])]
forall a. a -> [a] -> [a]
:[((t, Map t (Permutation t)), [Permutation t])]
ls) [((t, Map t (Permutation t)), [Permutation t])]
rs Permutation t
h t
b'
-- updateLevels' ls [] h b' = error $ "updateLevels: " ++ show (ls,[],h,b')

-- used the following in debugging
-- orderLevels levels = product $ [if M.null t then 1 else toInteger (M.size t) | ((b,t),s) <- levels]


-- recover the base tranversals from the sgs. gs must be an sgs
-- baseTransversalsSGS gs = [let hs = [h | h <- gs, b <= minsupp h] in (b, cosetRepsGx hs b) | b <- bs]
baseTransversalsSGS :: [Permutation k] -> [(k, Map k (Permutation k))]
baseTransversalsSGS gs :: [Permutation k]
gs = [let hs :: [Permutation k]
hs = (Permutation k -> Bool) -> [Permutation k] -> [Permutation k]
forall a. (a -> Bool) -> [a] -> [a]
filter ( (k
b k -> k -> Bool
forall a. Ord a => a -> a -> Bool
<=) (k -> Bool) -> (Permutation k -> k) -> Permutation k -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Permutation k -> k
forall c. Permutation c -> c
minsupp ) [Permutation k]
gs in (k
b, [Permutation k] -> k -> Map k (Permutation k)
forall k. Ord k => [Permutation k] -> k -> Map k (Permutation k)
cosetRepsGx [Permutation k]
hs k
b) | k
b <- [k]
bs]
    where bs :: [k]
bs = [k] -> [k]
forall b. Ord b => [b] -> [b]
toListSet ([k] -> [k]) -> [k] -> [k]
forall a b. (a -> b) -> a -> b
$ (Permutation k -> k) -> [Permutation k] -> [k]
forall a b. (a -> b) -> [a] -> [b]
map Permutation k -> k
forall c. Permutation c -> c
minsupp [Permutation k]
gs
    -- where bs = toListSet $ concatMap supp gs

-- |Given a strong generating set gs, isMemberSGS gs is a membership test for the group
isMemberSGS :: (Ord a, Show a) => [Permutation a] -> Permutation a -> Bool
isMemberSGS :: [Permutation a] -> Permutation a -> Bool
isMemberSGS gs :: [Permutation a]
gs h :: Permutation a
h = let bts :: [(a, Map a (Permutation a))]
bts = [Permutation a] -> [(a, Map a (Permutation a))]
forall k. Ord k => [Permutation k] -> [(k, Map k (Permutation k))]
baseTransversalsSGS [Permutation a]
gs in Maybe (Permutation a) -> Bool
forall a. Maybe a -> Bool
isNothing (Maybe (Permutation a) -> Bool) -> Maybe (Permutation a) -> Bool
forall a b. (a -> b) -> a -> b
$ [(a, Map a (Permutation a))]
-> Permutation a -> Maybe (Permutation a)
forall k.
Ord k =>
[(k, Map k (Permutation k))]
-> Permutation k -> Maybe (Permutation k)
sift [(a, Map a (Permutation a))]
bts Permutation a
h


{-
-- Alternative where we carry on with Schreier-Sims when we finish Random Schreier-Sims, just to make sure
-- !! Unfortunately, doesn't appear to work - perhaps ss' doesn't like finding empty levels
sgs2 gs = toListSet $ concatMap snd $ rss2 gs

rss2 gs = unsafePerformIO $
    do (r,xs) <- initProdRepl gs
       levels <- rss' (r,xs) (initLevels gs) 0
       return $ ss' bs (reverse levels) []
    where bs = toListSet $ concatMap supp gs
-}