    -- 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
    -}