RandomSchreierSims.hs

-- 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 = do (r,xs) <- initProdRepl $ _D 10
                  hs <- replicateM 20 $ nextProdRepl (r,xs)
                  mapM_ print 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 gs =
    let n = length gs
        r = max 10 n
        xs = (1:) $ take r $ concat $ repeat gs 
    in do xs' <- newListArray (0,r) xs
          replicateM_ 60 $ nextProdRepl (r,xs') -- perform initial mixing
          return (r,xs')

nextProdRepl :: (Ord a, Show a) => (Int, IOArray Int (Permutation a)) -> IO (Maybe (Permutation a))
nextProdRepl (r,xs) =
    do s <- randomRIO (1,r)
       t <- randomRIO (1,r)
       u <- randomRIO (0,3 :: Int)
       out <- updateArray xs s t u
       return out

updateArray xs s t u =
    let (swap,invert) = quotRem u 2 in
    if s == t
    then return Nothing
    else do
        x_0 <- readArray xs 0
        x_s <- readArray xs s
        x_t <- readArray xs t
        let x_s' = mult (swap,invert) x_s x_t
            x_0' = mult (swap,0) x_0 x_s'
        writeArray xs 0 x_0'
        writeArray xs s x_s'
        return (Just x_0')
    where mult (swap,invert) a b = case (swap,invert) of
                                   (0,0) -> a * b
                                   (0,1) -> a * b^-1
                                   (1,0) -> b * a
                                   (1,1) -> b^-1 * 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 gs = toListSet $ concatMap snd $ rss gs

rss gs = unsafePerformIO $
    do (r,xs) <- initProdRepl gs
       rss' (r,xs) (initLevels gs) 0

rss' (r,xs) levels i
    | i == 25 = return levels -- stop if we've had 25 successful sifts in a row
    | otherwise = do g <- nextProdRepl (r,xs)
                     let (changed,levels') = updateLevels levels g
                     rss' (r,xs) levels' (if changed then 0 else i+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 gs = [((b,M.singleton b 1),[]) | b <- bs]
    where bs = toListSet $ concatMap supp gs

updateLevels levels Nothing = (False,levels) -- not strictly correct to increment count on a Nothing
updateLevels levels (Just g) =
    case sift (map fst levels) g of
    Nothing -> (False, levels)
    -- Just 1 -> error "Just 1"
    Just g' -> (True, updateLevels' [] levels g' (minsupp g'))

updateLevels' ls (r@((b,t),s):rs) h b' =
    if b == b'
    then reverse ls ++ ((b, cosetRepsGx (h:s) b), h:s) : rs
    else updateLevels' (r:ls) rs h 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 gs = [let hs = filter ( (b <=) . minsupp ) gs in (b, cosetRepsGx hs b) | b <- bs]
    where bs = toListSet $ map minsupp 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 gs h = let bts = baseTransversalsSGS gs in isNothing $ sift bts 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
-}