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

    module Math.Combinatorics.GraphAuts (isVertexTransitive, isEdgeTransitive,
                                         isArcTransitive, is2ArcTransitive, is3ArcTransitive, isnArcTransitive,
                                         isDistanceTransitive,
                                         graphAuts, incidenceAuts,
                                         graphIsos, incidenceIsos,
                                         isGraphIso, isIncidenceIso) where

    import Data.Either (lefts)
    import qualified Data.List as L
    import qualified Data.Map as M
    import qualified Data.Set as S
    import Data.Maybe

    import Math.Common.ListSet
    import Math.Core.Utils (combinationsOf, pairs)
    import Math.Combinatorics.Graph
    -- import Math.Combinatorics.StronglyRegularGraph
    -- import Math.Combinatorics.Hypergraph -- can't import this, creates circular dependency
    import Math.Algebra.Group.PermutationGroup
    import Math.Algebra.Group.SchreierSims as SS

    -- The code for finding automorphisms - "graphAuts" - follows later on in file

    -- TRANSITIVITY PROPERTIES OF GRAPHS

    -- |A graph is vertex-transitive if its automorphism group acts transitively on the vertices. Thus, given any two distinct vertices, there is an automorphism mapping one to the other.
    isVertexTransitive :: (Ord t) => Graph t -> Bool
    isVertexTransitive (G [] []) = True -- null graph is trivially vertex transitive
    isVertexTransitive g@(G (v:vs) es) = orbitV auts v == v:vs where
        auts = graphAuts g

    -- |A graph is edge-transitive if its automorphism group acts transitively on the edges. Thus, given any two distinct edges, there is an automorphism mapping one to the other.
    isEdgeTransitive :: (Ord t) => Graph t -> Bool
    isEdgeTransitive (G _ []) = True
    isEdgeTransitive g@(G vs (e:es)) = orbitE auts e == e:es where
        auts = graphAuts g

    arc ->^ g = map (.^ g) arc
    -- unlike edges/blocks, arcs are directed, so the action on them does not sort

    -- Godsil & Royle 59-60
    -- |A graph is arc-transitive (or flag-transitive) if its automorphism group acts transitively on arcs. (An arc is an ordered pair of adjacent vertices.)
    isArcTransitive :: (Ord t) => Graph t -> Bool
    isArcTransitive (G _ []) = True -- empty graphs are trivially arc transitive
    isArcTransitive g@(G vs es) = orbit (->^) a auts == a:as where
    -- isArcTransitive g@(G vs es) = closure [a] [ ->^ h | h <- auts] == a:as where
        a:as = L.sort $ es ++ map reverse es
        auts = graphAuts g

    isArcTransitive' g@(G (v:vs) es) =
        orbitP auts v == v:vs && -- isVertexTransitive g
        orbitP stab n == n:ns
        where auts = graphAuts g
              stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
              n:ns = nbrs g v

    -- execution time of both of the above is dominated by the time to calculate the graph auts, so their performance is similar

    -- then k n, kb n n, q n, other platonic solids, petersen graph, heawood graph, pappus graph, desargues graph are all arc-transitive

    -- find arcs of length l from x using dfs - results returned in order
    -- an arc is a sequence of vertices connected by edges, no doubling back, but self-crossings allowed
    findArcs g@(G vs es) x l = map reverse $ dfs [ ([x],0) ] where
        dfs ( (z1:z2:zs,l') : nodes)
            | l == l'   = (z1:z2:zs) : dfs nodes
            | otherwise = dfs $ [(w:z1:z2:zs,l'+1) | w <- nbrs g z1, w /= z2] ++ nodes
        dfs ( ([z],l') : nodes)
            | l == l'   = [z] : dfs nodes
            | otherwise = dfs $ [([w,z],l'+1) | w <- nbrs g z] ++ nodes
        dfs [] = []

    -- note that a graph with triangles can't be 3-arc transitive, etc, because an aut can't map a self-crossing arc to a non-self-crossing arc

    -- |A graph is n-arc-transitive is its automorphism group is transitive on n-arcs. (An n-arc is an ordered sequence (v0,...,vn) of adjacent vertices, with crossings allowed but not doubling back.)
    isnArcTransitive :: (Ord t) => Int -> Graph t -> Bool
    isnArcTransitive _ (G [] []) = True
    isnArcTransitive n g@(G (v:vs) es) =
        orbitP auts v == v:vs && -- isVertexTransitive g
        orbit (->^) a stab == a:as
        -- closure [a] [ ->^ h | h <- stab] == a:as
        where auts = graphAuts g
              stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
              a:as = findArcs g v n

    is2ArcTransitive :: (Ord t) => Graph t -> Bool
    is2ArcTransitive g = isnArcTransitive 2 g

    is3ArcTransitive :: (Ord t) => Graph t -> Bool
    is3ArcTransitive g = isnArcTransitive 3 g

    -- Godsil & Royle 66-7
    -- |A graph is distance transitive if given any two ordered pairs of vertices (u,u') and (v,v') with d(u,u') == d(v,v'),
    -- there is an automorphism of the graph that takes (u,u') to (v,v')
    isDistanceTransitive :: (Ord t) => Graph t -> Bool
    isDistanceTransitive (G [] []) = True
    isDistanceTransitive g@(G (v:vs) es)
        | isConnected g =
            orbitP auts v == v:vs && -- isVertexTransitive g
            length stabOrbits == diameter g + 1 -- the orbits under the stabiliser of v coincide with the distance partition from v
        | otherwise = error "isDistanceTransitive: only defined for connected graphs"
        where auts = graphAuts g
              stab = dropWhile (\p -> v .^ p /= v) auts -- we know that graphAuts are returned in this order
              stabOrbits = let os = orbits stab in os ++ map (:[]) ((v:vs) L.\\ concat os) -- include fixed point orbits

    -- GRAPH AUTOMORPHISMS

    -- !! Note, in the literature the following is just called the intersection of two partitions
    -- !! Refinement actually refers to the process of refining to an equitable partition

    -- refine one partition by another
    refine p1 p2 = filter (not . null) $ refine' p1 p2
    -- Refinement preserves ordering within cells but not between cells
    -- eg the cell [1,2,3,4] could be refined to [2,4],[1,3]

    -- refine, but leaving null cells in
    -- we use this in the graphAuts functions when comparing two refinements to check that they split in the same way
    refine' p1 p2 = concat [ [c1 `intersect` c2 | c2 <- p2] | c1 <- p1]

    isGraphAut (G vs es) h = all (`S.member` es') [e -^ h | e <- es]
        where es' = S.fromList es
    -- this works best on sparse graphs, where p(edge) < 1/2
    -- if p(edge) > 1/2, it would be better to test on the complement of the graph

    -- Calculate a map consisting of neighbour lists for each vertex in the graph
    -- If a vertex has no neighbours then it is left out of the map
    adjLists (G vs es) = adjLists' M.empty es
        where adjLists' nbrs ([u,v]:es) =
                  adjLists' (M.insertWith' (flip (++)) v [u] $ M.insertWith' (flip (++)) u [v] nbrs) es
              adjLists' nbrs [] = nbrs

    -- ALTERNATIVE VERSIONS OF GRAPH AUTS
    -- (showing how we got to the final version)

    -- return all graph automorphisms, using naive depth first search
    graphAuts1 (G vs es) = dfs [] vs vs
        where dfs xys (x:xs) ys =
                  concat [dfs ((x,y):xys) xs (L.delete y ys) | y <- ys, isCompatible (x,y) xys]
              dfs xys [] [] = [fromPairs xys]
              isCompatible (x,y) xys = and [([x',x] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x',y') <- xys]
              es' = S.fromList es

    -- return generators for graph automorphisms
    -- (using Lemma 9.1.1 from Seress p203 to prune the search tree)
    graphAuts2 (G vs es) = graphAuts' [] vs
        where graphAuts' us (v:vs) =
                  let uus = zip us us
                  in concat [take 1 $ dfs ((v,w):uus) vs (v : L.delete w vs) | w <- vs, isCompatible (v,w) uus]
                  ++ graphAuts' (v:us) vs
                  -- stab us == transversal for stab (v:us) ++ stab (v:us)  (generators thereof)
              graphAuts' _ [] = [] -- we're not interested in finding the identity element
              dfs xys (x:xs) ys =
                  concat [dfs ((x,y):xys) xs (L.delete y ys) | y <- ys, isCompatible (x,y) xys]
              dfs xys [] [] = [fromPairs xys]
              isCompatible (x,y) xys = and [([x',x] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x',y') <- xys]
              es' = S.fromList es

    -- Now using distance partitions
    -- Note that because of the use of distance partitions, this is only valid for connected graphs
    graphAuts3 g@(G vs es) = graphAuts' [] [vs] where
        graphAuts' us ((x:ys):pt) =
            let px = refine' (ys : pt) (dps M.! x)
                p y = refine' ((x : L.delete y ys) : pt) (dps M.! y)
                uus = zip us us
                p' = L.sort $ filter (not . null) $ px
            in concat [take 1 $ dfs ((x,y):uus) px (p y) | y <- ys]
            ++ graphAuts' (x:us) p'
        graphAuts' us ([]:pt) = graphAuts' us pt
        graphAuts' _ [] = []
        dfs xys p1 p2
            | map length p1 /= map length p2 = []
            | otherwise =
                let p1' = filter (not . null) p1
                    p2' = filter (not . null) p2
                in if all isSingleton p1'
                   then let xys' = xys ++ zip (concat p1') (concat p2')
                        in if isCompatible xys' then [fromPairs' xys'] else []
                        -- we shortcut the search when we have all singletons, so must check isCompatible to avoid false positives
                   else let (x:xs):p1'' = p1'
                            ys:p2'' = p2'
                        in concat [dfs ((x,y):xys)
                                       (refine' (xs : p1'') (dps M.! x))
                                       (refine' ((L.delete y ys):p2'') (dps M.! y))
                                       | y <- ys]
        isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
        dps = M.fromList [(v, distancePartition g v) | v <- vs]
        es' = S.fromList es

    isSingleton [_] = True
    isSingleton _ = False

    -- Now we try to use generators we've already found at a given level to save us having to look for others
    -- For example, if we have found (1 2)(3 4) and (1 3 2), then we don't need to look for something taking 1 -> 4
    graphAuts4 g@(G vs es) = graphAuts' [] [vs] where
        graphAuts' us p@((x:ys):pt) =
            -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
            let p' = L.sort $ refine (ys:pt) (dps M.! x)
            in level us p x ys []
            ++ graphAuts' (x:us) p'
        graphAuts' us ([]:pt) = graphAuts' us pt
        graphAuts' _ [] = []
        level us p@(ph:pt) x (y:ys) hs =
            let px = refine' (L.delete x ph : pt) (dps M.! x)
                py = refine' (L.delete y ph : pt) (dps M.! y)
                uus = zip us us
            in case dfs ((x,y):uus) px py of
               []  -> level us p x ys hs
               h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
        level _ _ _ [] _ = []
        dfs xys p1 p2
            | map length p1 /= map length p2 = []
            | otherwise =
                let p1' = filter (not . null) p1
                    p2' = filter (not . null) p2
                in if all isSingleton p1'
                   then let xys' = xys ++ zip (concat p1') (concat p2')
                        in if isCompatible xys' then [fromPairs' xys'] else []
                   else let (x:xs):p1'' = p1'
                            ys:p2'' = p2'
                        in concat [dfs ((x,y):xys)
                                       (refine' (xs : p1'') (dps M.! x))
                                       (refine' ((L.delete y ys):p2'') (dps M.! y))
                                       | y <- ys]
        isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
        dps = M.fromList [(v, distancePartition g v) | v <- vs]
        es' = S.fromList es

    -- contrary to first thought, you can't stop when a level is null - eg kb 2 3, the third level is null, but the fourth isn't

    -- an example for equitable partitions
    -- this is a graph whose distance partition (from any vertex) can be refined to an equitable partition
    eqgraph = G vs es where
        vs = [1..14]
        es = L.sort $ [[1,14],[2,13]] ++ [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, v1+1 == v2 || v1+3 == v2 && even v2]

    -- refine a partition to give an equitable partition
    toEquitable g cells = L.sort $ toEquitable' [] cells where
        toEquitable' ls (r:rs) =
            let (lls,lrs) = L.partition isSingleton $ map (splitNumNbrs r) ls
                -- so the lrs split, and the lls didn't
                rs' = concatMap (splitNumNbrs r) rs
            in if isSingleton r -- then we know it won't split further, so can remove it from further processing
               then r : toEquitable' (concat lls) (concat lrs ++ rs')
               else toEquitable' (r : concat lls) (concat lrs ++ rs')
        toEquitable' ls [] = ls
        splitNumNbrs t c = map (map snd) $ L.groupBy (\x y -> fst x == fst y) $ L.sort
                        [ (length ((nbrs_g M.! v) `intersect` t), v) | v <- c]
        nbrs_g = M.fromList [(v, nbrs g v) | v <- vertices g]

    -- try to refine two partitions in parallel, failing if they become mismatched
    toEquitable2 nbrs_g psrc ptrg = unzip $ L.sort $ toEquitable' [] (zip psrc ptrg) where
        toEquitable' ls (r:rs) =
            let ls' = map (splitNumNbrs nbrs_g r) ls
                (lls,lrs) = L.partition isSingleton $ map fromJust ls'
                rs' = map (splitNumNbrs nbrs_g r) rs
            in if any isNothing ls' || any isNothing rs'
               then []
               else
                   {- if (isSingleton . fst) r
                   then r : toEquitable' (concat lls) (concat lrs ++ concatMap fromJust rs')
                   else -} toEquitable' (r : concat lls) (concat lrs ++ concatMap fromJust rs')
        toEquitable' ls [] = ls

    splitNumNbrs nbrs_g (t_src,t_trg) (c_src,c_trg) =
        let src_split = L.groupBy (\x y -> fst x == fst y) $ L.sort
                        [ (length ((nbrs_g M.! v) `intersect` t_src), v) | v <- c_src]
            trg_split = L.groupBy (\x y -> fst x == fst y) $ L.sort
                        [ (length ((nbrs_g M.! v) `intersect` t_trg), v) | v <- c_trg]
        in if map length src_split == map length trg_split
           && map (fst . head) src_split == map (fst . head) trg_split
           then Just $ zip (map (map snd) src_split) (map (map snd) trg_split)
           else Nothing
           -- else error (show (src_split, trg_split)) -- for debugging

    -- Now, every time we intersect two partitions, refine to an equitable partition

    -- |Given a graph g, @graphAuts g@ returns generators for the automorphism group of g.
    -- If g is connected, then the generators will be a strong generating set.
    graphAuts :: (Ord a) => Graph a -> [Permutation a]
    graphAuts g = autsWithinComponents ++ isosBetweenComponents
        where cs = map (inducedSubgraph g) (components g)
              -- autsWithinComponents = concatMap graphAutsCon cs
              autsWithinComponents = concatMap graphAuts4 cs
              isosBetweenComponents = map swapFromIso $ concat [take 1 (graphIsos ci cj) | (ci,cj) <- pairs cs]
              swapFromIso xys = fromPairs (xys ++ map swap xys)
              swap (x,y) = (y,x)
    -- Using graphAuts4 instead of graphAutsCon as latter appears to have a bug, eg
    -- > graphAuts4 $ G [1..3] [[1,2],[2,3]]
    -- [[[1,3]]]
    -- > graphAutsCon $ G [1..3] [[1,2],[2,3]]
    -- []

    -- Automorphisms of a connected graph
    graphAutsCon g@(G vs es)
        | isConnected g = graphAuts' [] (toEquitable g $ valencyPartition g)
        | otherwise = error "graphAutsCon: graph is not connected"
        where graphAuts' us p@((x:ys):pt) =
                  let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
                  in level us p x ys []
                  ++ graphAuts' (x:us) p'
              graphAuts' us ([]:pt) = graphAuts' us pt
              graphAuts' _ [] = []
              level us p@(ph:pt) x (y:ys) hs =
                  let px = refine' (L.delete x ph : pt) (dps M.! x)
                      py = refine' (L.delete y ph : pt) (dps M.! y)
                      uus = zip us us
                  in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
                     []  -> level us p x ys hs
                     h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
              level _ _ _ [] _ = []
              dps = M.fromList [(v, distancePartition g v) | v <- vs]
              es' = S.fromList es
              nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]

    dfsEquitable (dps,es',nbrs_g) xys p1 p2 = dfs xys p1 p2 where
        dfs xys p1 p2
            | map length p1 /= map length p2 = []
            | otherwise =
                let p1' = filter (not . null) p1
                    p2' = filter (not . null) p2
                    (p1e,p2e) = toEquitable2 nbrs_g p1' p2'
                in if null p1e
                   then []
                   else
                       if all isSingleton p1e
                       then let xys' = xys ++ zip (concat p1e) (concat p2e)
                            in if isCompatible xys' then [fromPairs' xys'] else []
                       else let (x:xs):p1'' = p1e
                                ys:p2'' = p2e
                            in concat [dfs ((x,y):xys)
                                           (refine' (xs : p1'') (dps M.! x))
                                           (refine' ((L.delete y ys):p2'') (dps M.! y))
                                           | y <- ys]
        isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']

    -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH

    -- based on graphAuts as applied to the incidence graph, but modified to avoid point-block crossover auts

    -- |Given the incidence graph of an incidence structure between points and blocks
    -- (for example, a set system),
    -- @incidenceAuts g@ returns generators for the automorphism group of the incidence structure.
    -- The generators are represented as permutations of the points.
    -- The incidence graph should be represented with the points on the left and the blocks on the right.
    -- If the incidence graph is connected, then the generators will be a strong generating set.
    incidenceAuts :: (Ord p, Ord b) => Graph (Either p b) -> [Permutation p]
    incidenceAuts g = autsWithinComponents ++ isosBetweenComponents
        where cs = map (inducedSubgraph g) (components g)
              autsWithinComponents = concatMap incidenceAutsCon cs
              isosBetweenComponents = map swapFromIso $ concat [take 1 (incidenceIsos ci cj) | (ci,cj) <- pairs cs]
              swapFromIso xys = fromPairs (xys ++ map swap xys)
              swap (x,y) = (y,x)

    incidenceAutsCon g@(G vs es)
        | isConnected g = map points (incidenceAuts' [] [vs])
        | otherwise = error "incidenceAutsCon: graph is not connected"
        where points h = fromPairs [(x,y) | (Left x, Left y) <- toPairs h] -- filtering out the action on blocks
              incidenceAuts' us p@((x@(Left _):ys):pt) =
                  -- let p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
                  let p' = L.sort $ refine (ys:pt) (dps M.! x)
                  in level us p x ys []
                  ++ incidenceAuts' (x:us) p'
              incidenceAuts' us ([]:pt) = incidenceAuts' us pt
              incidenceAuts' _ (((Right _):_):_) = [] -- if we fix all the points, then the blocks must be fixed too
              incidenceAuts' _ [] = []
              level us p@(ph:pt) x (y@(Left _):ys) hs =
                  let px = refine' (L.delete x ph : pt) (dps M.! x)
                      py = refine' (L.delete y ph : pt) (dps M.! y)
                      uus = zip us us
                  in case dfsEquitable (dps,es',nbrs_g) ((x,y):uus) px py of
                     []  -> level us p x ys hs
                     h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
              level _ _ _ _ _ = [] -- includes the case where y matches Right _, which can only occur on first level, before we've distance partitioned
              dps = M.fromList [(v, distancePartition g v) | v <- vs]
              es' = S.fromList es
              nbrs_g = M.fromList [(v, nbrs g v) | v <- vs]

    -- GRAPH ISOMORPHISMS

    -- !! not yet using equitable partitions, so could probably be more efficient

    -- graphIsos :: (Ord a, Ord b) => Graph a -> Graph b -> [[(a,b)]]
    graphIsos g1 g2
        | length cs1 /= length cs2 = []
        | otherwise = graphIsos' cs1 cs2
        where cs1 = map (inducedSubgraph g1) (components g1)
              cs2 = map (inducedSubgraph g2) (components g2)
              graphIsos' (ci:cis) cjs =
                  [iso ++ iso' | cj <- cjs,
                                 iso <- graphIsosCon ci cj,
                                 let cjs' = L.delete cj cjs,
                                 iso' <- graphIsos' cis cjs']
              graphIsos' [] [] = [[]]

    -- isos between connected graphs
    graphIsosCon g1 g2
        | isConnected g1 && isConnected g2
            = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2)
                     | v1 <- take 1 (vertices g1), v2 <- vertices g2]
                     -- the take 1 handles the case where g1 is the null graph
        | otherwise = error "graphIsosCon: either or both graphs are not connected"
        where dfs xys p1 p2
                  | map length p1 /= map length p2 = []
                  | otherwise =
                      let p1' = filter (not . null) p1
                          p2' = filter (not . null) p2
                      in if all isSingleton p1'
                         then let xys' = xys ++ zip (concat p1') (concat p2')
                              in if isCompatible xys' then [xys'] else []
                         else let (x:xs):p1'' = p1'
                                  ys:p2'' = p2'
                              in concat [dfs ((x,y):xys)
                                             (refine' (xs : p1'') (dps1 M.! x))
                                             (refine' ((L.delete y ys):p2'') (dps2 M.! y))
                                             | y <- ys]
              isCompatible xys = and [([x,x'] `S.member` es1) == (L.sort [y,y'] `S.member` es2) | (x,y) <- xys, (x',y') <- xys, x < x']
              dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
              dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
              es1 = S.fromList $ edges g1
              es2 = S.fromList $ edges g2

    -- |Are the two graphs isomorphic?
    isGraphIso :: (Ord a, Ord b) => Graph a -> Graph b -> Bool
    isGraphIso g1 g2 = (not . null) (graphIsos g1 g2)
    -- !! If we're only interested in seeing whether or not two graphs are iso,
    -- !! then the cost of calculating distancePartitions may not be warranted
    -- !! (see Math.Combinatorics.Poset: orderIsos01 versus orderIsos)

    -- !! deprecate
    isIso g1 g2 = (not . null) (graphIsos g1 g2)

    -- the following differs from graphIsos in only two ways
    -- we avoid Left, Right crossover isos, by insisting that a Left is taken to a Left (first two lines)
    -- we return only the action on the Lefts, and unLeft it
    -- incidenceIsos :: (Ord p1, Ord b1, Ord p2, Ord b2) =>
    --     Graph (Either p1 b1) -> Graph (Either p2 b2) -> [[(p1,p2)]]

    incidenceIsos g1 g2
        | length cs1 /= length cs2 = []
        | otherwise = incidenceIsos' cs1 cs2
        where cs1 = map (inducedSubgraph g1) (filter (not . null . lefts) $ components g1)
              cs2 = map (inducedSubgraph g2) (filter (not . null . lefts) $ components g2)
              incidenceIsos' (ci:cis) cjs =
                  [iso ++ iso' | cj <- cjs,
                                 iso <- incidenceIsosCon ci cj,
                                 let cjs' = L.delete cj cjs,
                                 iso' <- incidenceIsos' cis cjs']
              incidenceIsos' [] [] = [[]]

    incidenceIsosCon g1 g2
        | isConnected g1 && isConnected g2
            = concat [dfs [] (distancePartition g1 v1) (distancePartition g2 v2)
                     | v1@(Left _) <- take 1 (vertices g1), v2@(Left _) <- vertices g2]
                     -- g1 may have no vertices
        | otherwise = error "incidenceIsos: one or both graphs not connected"
        where dfs xys p1 p2
                  | map length p1 /= map length p2 = []
                  | otherwise =
                      let p1' = filter (not . null) p1
                          p2' = filter (not . null) p2
                      in if all isSingleton p1'
                         then let xys' = xys ++ zip (concat p1') (concat p2')
                              in if isCompatible xys' then [[(x,y) | (Left x, Left y) <- xys']] else []
                         else let (x:xs):p1'' = p1'
                                  ys:p2'' = p2'
                              in concat [dfs ((x,y):xys)
                                             (refine' (xs : p1'') (dps1 M.! x))
                                             (refine' ((L.delete y ys):p2'') (dps2 M.! y))
                                             | y <- ys]
              isCompatible xys = and [([x,x'] `S.member` es1) == (L.sort [y,y'] `S.member` es2) | (x,y) <- xys, (x',y') <- xys, x < x']
              dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1]
              dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2]
              es1 = S.fromList $ edges g1
              es2 = S.fromList $ edges g2

    -- |Are the two incidence structures represented by these incidence graphs isomorphic?
    isIncidenceIso :: (Ord p1, Ord b1, Ord p2, Ord b2) =>
         Graph (Either p1 b1) -> Graph (Either p2 b2) -> Bool
    isIncidenceIso g1 g2 = (not . null) (incidenceIsos g1 g2)

    {-
    removeGens x gs = removeGens' [] gs where
        baseOrbit = x .^^ gs
        removeGens' ls (r:rs) =
            if x .^^ (ls++rs) == baseOrbit
            then removeGens' ls rs
            else removeGens' (r:ls) rs
        removeGens' ls [] = reverse ls
    -- !! reverse is probably pointless

    -- !! DON'T THINK THIS IS WORKING PROPERLY
    -- eg graphAutsSGSNew $ toGraph ([1..7],[[1,3],[2,3],[3,4],[4,5],[4,6],[4,7]])
    -- returns [[[1,2]],[[5,6]],[[5,7,6]],[[6,7]]]
    -- whereas [[6,7]] was a Schreier generator, so shouldn't have been listed

    -- Using Schreier generators to seed the next level
    -- At the moment this is slower than the above
    -- (This could be modified to allow us to start the search with a known subgroup)
    graphAutsNew g@(G vs es) = graphAuts' [] [] [vs] where
        graphAuts' us hs p@((x:ys):pt) =
            let ys' = ys L.\\ (x .^^ hs) -- don't need to consider points which can already be reached from Schreier generators
                hs' = level us p x ys' []
                p' = L.sort $ filter (not . null) $ refine' (ys:pt) (dps M.! x)
                reps = cosetRepsGx (hs'++hs) x
                schreierGens = removeGens x $ schreierGeneratorsGx (x,reps) (hs'++hs)
            in hs' ++ graphAuts' (x:us) schreierGens p'
        graphAuts' us hs ([]:pt) = graphAuts' us hs pt
        graphAuts' _ _ [] = []
        level us p@(ph:pt) x (y:ys) hs =
            let px = refine' (L.delete x ph : pt) (dps M.! x)
                py = refine' (L.delete y ph : pt) (dps M.! y)
                uus = zip us us
            in if map length px /= map length py
               then level us p x ys hs
               else case dfs ((x,y):uus) (filter (not . null) px) (filter (not . null) py) of
                    []  -> level us p x ys hs
                    h:_ -> let hs' = h:hs in h : level us p x (ys L.\\ (x .^^ hs')) hs'
                    -- if h1 = (1 2)(3 4), and h2 = (1 3 2), then we can remove 4 too
        level _ _ _ [] _ = []
        dfs xys p1 p2
            | map length p1 /= map length p2 = []
            | otherwise =
                let p1' = filter (not . null) p1
                    p2' = filter (not . null) p2
                in if all isSingleton p1'
                   then let xys' = xys ++ zip (concat p1') (concat p2')
                        in if isCompatible xys' then [fromPairs' xys'] else []
                   else let (x:xs):p1'' = p1'
                            ys:p2'' = p2'
                        in concat [dfs ((x,y):xys)
                                       (refine' (xs : p1'') (dps M.! x))
                                       (refine' ((L.delete y ys):p2'') (dps M.! y))
                                       | y <- ys]
        isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x']
        dps = M.fromList [(v, distancePartition g v) | v <- vs]
        es' = S.fromList es
    -}