-- Copyright (c) David Amos, 2009-2014. All rights reserved. {-# LANGUAGE NoMonomorphismRestriction, TupleSections, DeriveFunctor #-} module Math.Combinatorics.GraphAuts (isVertexTransitive, isEdgeTransitive, isArcTransitive, is2ArcTransitive, is3ArcTransitive, is4ArcTransitive, isnArcTransitive, isDistanceTransitive, graphAuts, incidenceAuts, graphAuts7, graphAuts8, incidenceAuts2, isGraphAut, isIncidenceAut, graphIsos, incidenceIsos, isGraphIso, isIncidenceIso) where import Data.Either (lefts, rights, partitionEithers) import qualified Data.List as L import qualified Data.Map as M import qualified Data.Set as S import Data.Maybe import Data.Ord (comparing) import qualified Data.Foldable as Foldable import qualified Data.Sequence as Seq import Math.Common.ListSet import Math.Core.Utils (combinationsOf, intersectAsc, pairs, picks, (^-)) 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 = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts -- 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 if 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 = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts -- 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 -- The incidence graphs of the projective planes PG(2,Fq) are 4-arc-transitive is4ArcTransitive :: (Ord t) => Graph t -> Bool is4ArcTransitive g = isnArcTransitive 4 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 implemented for connected graphs" where auts = graphAuts g stab = filter (\p -> v .^ p == v) auts -- relies on v being the first base for the SGS returned by graphAuts -- 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 -- |Is the permutation an automorphism of the graph? isGraphAut :: Ord t => Graph t -> Permutation t -> Bool 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 -- |Is the permutation an automorphism of the incidence structure represented by the graph? -- (Note that an incidence graph colours points as Left, blocks as Right, and a permutation -- that swaps points and blocks, even if it is an automorphism of the graph, does not represent -- an automorphism of the incidence structure. Instead, a point-block crossover is called a duality.) isIncidenceAut :: (Ord p, Ord b) => Graph (Either p b) -> Permutation (Either p b) -> Bool isIncidenceAut (G vs es) h = all (`S.member` es') [e ->^ h | e <- es] -- using ->^ instead of -^ excludes dualities, since each edge is of the form [Left p, Right b] where es' = S.fromList es -- 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) data SearchTree a = T Bool a [SearchTree a] deriving (Eq, Ord, Show, Functor) -- The boolean indicates whether or not this is a terminal / solution node leftDepth (T _ _ []) = 1 leftDepth (T _ _ (t:ts)) = 1 + leftDepth t leftWidths (T _ _ []) = [] leftWidths (T _ _ ts@(t:_)) = length ts : leftWidths t graphAutsEdgeSearchTree (G vs es) = dfs [] vs vs where dfs xys (x:xs) yys = T False xys [dfs ((x,y):xys) xs ys | (y,ys) <- picks yys, isCompatible xys (x,y)] dfs xys [] [] = T True xys [] isCompatible xys (x',y') = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys] es' = S.fromList es graphAuts1 = map fromPairs . terminals . graphAutsEdgeSearchTree terminals (T False _ ts) = concatMap terminals ts terminals (T True xys _) = [xys] -- Using Lemma 9.1.1 from Seress p203 to prune the search tree -- Because auts form a group, it is sufficient to expand only each leftmost branch of the tree in full. -- For every other branch, it is sufficient to find a single representative, since the other elements -- can then be obtained by multiplication in the group (using the leftmost elements). -- In effect, we are finding a transversal generating set. -- Note however, that this transversal generating set is relative to whatever base order the tree uses, -- so for clarity, the tree should use natural vertex order. transversalTerminals (T False _ (t:ts)) = concatMap (take 1 . transversalTerminals) ts ++ transversalTerminals t -- transversalTerminals (T False _ (t:ts)) = transversalTerminals t ++ concatMap (take 1 . transversalTerminals) ts transversalTerminals (T True xys _) = [xys] transversalTerminals _ = [] graphAuts2 = filter (/=1) . map fromPairs . transversalTerminals . graphAutsEdgeSearchTree -- init because last is identity isSingleton [_] = True isSingleton _ = False intersectCells p1 p2 = concat [ [c1 `intersectAsc` c2 | c2 <- p2] | c1 <- p1] -- Intersection preserves ordering within cells but not between cells -- eg the cell [1,2,3,4] could be refined to [2,4],[1,3] graphAutsDistancePartitionSearchTree g@(G vs es) = dfs [] ([vs],[vs]) where dfs xys (srcPart,trgPart) | all isSingleton srcPart = let xys' = zip (concat srcPart) (concat trgPart) in T (isCompatible xys') (xys++xys') [] -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible. -- However, we do need to check that the xys' are edge-compatible with each other. | otherwise = let (x:xs):srcCells = srcPart yys :trgCells = trgPart srcPart' = intersectCells (xs : srcCells) (dps M.! x) in T False xys -- the L.sort in the following line is so that we traverse vertices in natural order [dfs ((x,y):xys) ((unzip . L.sort) (zip (filter (not . null) srcPart') (filter (not . null) trgPart'))) | (y,ys) <- picks yys, let trgPart' = intersectCells (ys : trgCells) (dps M.! y), map length srcPart' == map length trgPart'] isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x'] es' = S.fromList es dps = M.fromAscList [(v, distancePartitionS vs es' v) | v <- vs] graphAuts3 = filter (/=1) . map fromPairs . transversalTerminals . graphAutsDistancePartitionSearchTree -- Whereas transversalTerminals produced a transversal generating set, here we produce a strong generating set. -- In particular, if we have already found (3 4), and then we find (1 2 3), -- then there is no need to look for (1 3 ...) or (1 4 ...), since it is clear that such elements exist -- as products of those we have already found. strongTerminals = strongTerminals' [] where strongTerminals' gs (T False xys ts) = case listToMaybe $ reverse $ filter (\(x,y) -> x /= y) xys of -- the first vertex that isn't fixed Nothing -> L.foldl' (\hs t -> strongTerminals' hs t) gs ts Just (x,y) -> if y `elem` (x .^^ gs) then gs -- Since we're not on the leftmost spine, we can stop as soon as we find one new element else find1New gs ts -- else L.foldl' (\hs t -> if hs /= gs then hs else strongTerminals' hs t) gs ts strongTerminals' gs (T True xys []) = fromPairs xys : gs find1New gs (t:ts) = let hs = strongTerminals' gs t in if take 1 gs /= take 1 hs -- we know a new element would be placed at the front then hs else find1New gs ts find1New gs [] = gs -- |Given a graph g, @graphAuts g@ returns a strong generating set for the automorphism group of g. graphAuts :: (Ord a) => Graph a -> [Permutation a] graphAuts = filter (/=1) . strongTerminals . graphAutsDistancePartitionSearchTree -- Using colourings (M.Map vertex colour, M.Map colour [vertex]), in place of partitions ([[vertex]]) -- This turns out to be slower than using partitions. -- Updating the colour partition incrementally seems to be much less efficient than just recalculating it each time -- (Recalculating each time is O(n), incrementally updating is O(n^2)?) graphAutsDistanceColouringSearchTree g@(G vs es) = dfs [] unitCol unitCol where unitCol = (M.fromList $ map (,[]) vs, M.singleton [] vs) -- "unit colouring" dfs xys srcColouring@(srcVmap,srcCmap) trgColouring@(trgVmap,trgCmap) -- ( | M.map length srcCmap /= M.map length trgCmap = T False xys [] ) | all isSingleton (M.elems srcCmap) = -- discrete colouring let xys' = zip (concat $ M.elems srcCmap) (concat $ M.elems trgCmap) in T (isCompatible xys') (reverse xys'++xys) [] -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible. -- However, we do need to check that the xys' are edge-compatible with each other. | otherwise = let (x,c) = M.findMin srcVmap (xVmap,xCmap) = dcs M.! x ys = trgCmap M.! c srcVmap' = M.delete x (intersectColouring srcVmap xVmap) srcCmap' = colourPartition srcVmap' -- srcCmap' = M.fromAscList [(c1++c2, cell) | (c1,srcCell) <- M.assocs srcCmap, (c2,xCell) <- M.assocs xCmap, -- let cell = L.delete x (intersectAsc srcCell xCell), -- (not . null) cell] in T False xys [dfs ((x,y):xys) (srcVmap',srcCmap') (trgVmap',trgCmap') | y <- ys, let (yVmap,yCmap) = dcs M.! y, let trgVmap' = M.delete y (intersectColouring trgVmap yVmap), let trgCmap' = colourPartition trgVmap', -- let trgCmap' = M.fromAscList [(c1++c2, cell) | (c1,trgCell) <- M.assocs trgCmap, (c2,yCell) <- M.assocs yCmap, -- let cell = L.delete y (intersectAsc trgCell yCell), -- (not . null) cell], M.map length srcCmap' == M.map length trgCmap' ] isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x'] es' = S.fromList es dcs = M.fromAscList [(v, distanceColouring v) | v <- vs] distanceColouring u = let dp = distancePartitionS vs es' u vmap = M.fromList [(v,[c]) | (cell,c) <- zip dp [0..], v <- cell] cmap = M.fromList $ zip (map (:[]) [0..]) dp in (vmap, cmap) {- -- If we are going to recalculate the colour partition each time anyway, -- then we don't need to carry it around, and can simplify the code graphAutsDistanceColouringSearchTree g@(G vs es) = dfs [] initCol initCol where initCol = M.fromList $ map (,[]) vs dfs xys srcCol trgCol | M.map length srcPart /= M.map length trgPart = T False xys [] | all isSingleton (M.elems srcPart) = let xys' = zip (concat $ M.elems srcPart) (concat $ M.elems trgPart) in T (isCompatible xys') (reverse xys'++xys) [] -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible. -- However, we do need to check that the xys' are edge-compatible with each other. | otherwise = let (x,c) = M.findMin srcCol ys = trgPart M.! c srcCol' = M.delete x $ intersectColouring srcCol (dcs M.! x) in T False xys [dfs ((x,y):xys) srcCol' trgCol' | y <- ys, let trgCol' = M.delete y (intersectColouring trgCol (dcs M.! y))] where srcPart = colourPartition srcCol trgPart = colourPartition trgCol isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x'] es' = S.fromList es dcs = M.fromAscList [(v, distanceColouring v) | v <- vs] distanceColouring u = M.fromList [(v,[c]) | (cell,c) <- zip (distancePartitionS vs es' u) [0..], v <- cell] -} distanceColouring (G vs es) u = M.fromList [(v,[c]) | (cell,c) <- zip (distancePartitionS vs es' u) [0..], v <- cell] where es' = S.fromList es intersectColouring c1 c2 = M.intersectionWith (++) c1 c2 colourPartition c = L.foldr (\(k,v) m -> M.insertWith (++) v [k] m) M.empty (M.assocs c) -- Based on McKay’s Canonical Graph Labeling Algorithm, by Stephen G. Hartke and A. J. Radcliffe -- (http://www.math.unl.edu/~aradcliffe1/Papers/Canonical.pdf) equitableRefinement g@(G vs es) p = equitableRefinement' (S.fromList es) p equitableRefinement' edgeset partition = go partition where go cells = let splits = L.zip (L.inits cells) (L.tails cells) shatterPairs = [(L.zip ci counts,ls,rs) | (ls,ci:rs) <- splits, cj <- cells, let counts = map (nbrCount cj) ci, isShatter counts] in case shatterPairs of -- by construction, the lexicographic least (i,j) comes first [] -> cells (vcs,ls,rs):_ -> let fragments = shatter vcs in go (ls ++ fragments ++ rs) isShatter (c:cs) = any (/= c) cs shatter vcs = map (map fst) $ L.groupBy (\x y -> snd x == snd y) $ L.sortBy (comparing snd) $ vcs -- Memoizing here results in about 10% speed improvement. Not worth it for loss of generality (ie requiring HasTrie instances) -- nbrCount = memo2 nbrCount' -- How many neighbours in cell does vertex have nbrCount cell vertex = length (filter (isEdge vertex) cell) isEdge u v = L.sort [u,v] `S.member` edgeset equitablePartitionSearchTree g@(G vs es) p = dfs [] p where dfs bs p = let p' = equitableRefinement' es' p in if all isSingleton p' then T True (p',bs) [] else T False (p',bs) [dfs (b:bs) p'' | (b,p'') <- splits [] p'] -- For now, we just split the first non-singleton cell we find splits ls (r:rs) | isSingleton r = splits (r:ls) rs | otherwise = let ls' = reverse ls in [(x, ls' ++ [x]:xs:rs) | (x,xs) <- picks r] es' = S.fromList es {- -- Using Data.Sequence instead of list for the partitions -- Makes no difference to speed (in fact slightly slower) equitableRefinementSeq' edgeset partition = go partition where go cells = let splits = Seq.zip (Seq.inits cells) (Seq.tails cells) shatterPairs = [(L.zip ci counts,ls,rs') | (ls,rs) <- Foldable.toList splits, (not . Seq.null) rs, let ci Seq.:< rs' = Seq.viewl rs, cj <- Foldable.toList cells, let counts = map (nbrCount cj) ci, isShatter counts] in case shatterPairs of -- by construction, the lexicographic least (i,j) comes first [] -> cells (vcs,ls,rs):_ -> let fragments = Seq.fromList (shatter vcs) in go (ls Seq.>< fragments Seq.>< rs) isShatter (c:cs) = any (/= c) cs shatter vcs = map (map fst) $ L.groupBy (\x y -> snd x == snd y) $ L.sortBy (comparing snd) $ vcs -- How many neighbours in cell does vertex have nbrCount cell vertex = length (filter (isEdge vertex) cell) isEdge u v = L.sort [u,v] `S.member` edgeset equitablePartitionSeqSearchTree g@(G vs es) p = dfs [] (Seq.fromList p) where dfs bs p = let p' = equitableRefinementSeq' es' p in if Foldable.all isSingleton p' then T True (Foldable.toList p',bs) [] else T False (Foldable.toList p',bs) [dfs (b:bs) p'' | (b,p'') <- splits p'] -- For now, we just split the first non-singleton cell we find splits cells = case Seq.findIndexL (not . isSingleton) cells of Just i -> let (ls,rs) = Seq.splitAt i cells r Seq.:< rs' = Seq.viewl rs in [(x, ls Seq.>< ([x] Seq.<| xs Seq.<| rs')) | (x,xs) <- picks r] Nothing -> error "Not possible, as we know there are non-singleton cells" es' = S.fromList es -} -- In this version, whenever we have an equitable partition, we separate out all the singleton cells and put them to one side. -- (Since the partition is equitable, singleton cells have already done any work they are going to do in shattering other cells, -- so they will no longer play any part.) -- This seems to result in about 20% speedup. equitablePartitionSearchTree2 g@(G vs es) p = dfs [] ([],p) where dfs bs (ss,cs) = let (ss',cs') = L.partition isSingleton $ equitableRefinement' es' cs ss'' = ss++ss' in case cs' of [] -> T True (ss'',bs) [] -- We just split the first non-singleton cell -- c:cs'' -> T False (ss''++cs',bs) [dfs (x:bs) (ss'',[x]:xs:cs'') | (x,xs) <- picks c] c:cs'' -> T False (cs'++ss'',bs) [dfs (x:bs) (ss'',[x]:xs:cs'') | (x,xs) <- picks c] es' = S.fromList es -- TODO: On the first level, we can use a stronger partitioning function (eg distance partitions, + see nauty manual, vertex invariants) equitableDistancePartitionSearchTree g@(G vs es) p = dfs [] p where dfs bs p = let p' = equitableRefinement' es' p in if all isSingleton p' then T True (p',bs) [] else T False (p',bs) [dfs (b:bs) p'' | (b,p'') <- splits [] p'] -- For now, we just split the first non-singleton cell we find splits ls (r:rs) | isSingleton r = splits (r:ls) rs | otherwise = [(x, p'') | let ls' = reverse ls, (x,xs) <- picks r, let p' = ls' ++ [x]:xs:rs, let p'' = filter (not . null) (intersectCells p' (dps M.! x))] es' = S.fromList es dps = M.fromAscList [(v, distancePartitionS vs es' v) | v <- vs] {- -- This is just fmap (\(p,bs) -> (p,bs,trace p)) t equitablePartitionTracedSearchTree g@(G vs es) trace p = dfs [] p where dfs bs p = let p' = equitableRefinement' es' p in if all isSingleton p' then T True (p',bs,trace p') [] else T False (p',bs,trace p') [dfs (b:bs) p'' | (b,p'') <- splits [] p'] -- For now, we just split the first non-singleton cell we find splits ls (r:rs) | isSingleton r = splits (r:ls) rs | otherwise = let ls' = reverse ls in [(x, ls' ++ [x]:xs:rs) | (x,xs) <- picks r] es' = S.fromList es -} -- Intended as a node invariant trace1 p = map (\xs@(x:_) -> (x, length xs)) $ L.group $ L.sort $ map length p equitablePartitionGraphSearchTree g@(G vs es) = equitablePartitionSearchTree g unitPartition where unitPartition = [vs] -- The incidence graph has vertices that are coloured left (points) or right (blocks). -- We are not interested in dualities (automorphisms that swap points and blocks), so we look for colour-preserving automorphisms equitablePartitionIncidenceSearchTree g@(G vs es) = equitablePartitionSearchTree g lrPartition where (lefts, rights) = partitionEithers vs lrPartition = [map Left lefts, map Right rights] leftLeaf (T False _ (t:ts)) = leftLeaf t leftLeaf (T True (p,bs) []) = (concat p, reverse bs) {- leftSpine (T False x (t:ts)) = x : leftSpine t leftSpine (T True x []) = [x] -} allLeaves (T False _ ts) = concatMap allLeaves ts allLeaves (T True (p,bs) []) = [(concat p, reverse bs)] {- partitionTransversals tree = [fromPairs (zip canonical partition) | partition <- findTransversals tree] where (_,canonical) = leftLeaf tree findTransversals (T False _ (t:ts)) = concatMap (take 1 . findTransversals) ts ++ findTransversals t findTransversals (T True (_,partition) []) = [concat partition] graphAuts5 = partitionTransversals . equitablePartitionGraphSearchTree -} -- NOT WORKING partitionBSGS0 g@(G vs es) t = (bs, findLevels t) where (p1,bs) = leftLeaf t g1 = fromPairs $ zip p1 vs g1' = g1^-1 es1 = S.fromList $ edges $ fmap (.^ g1) g -- the edges of the isomorph corresponding to p1. (S.fromList makes it unnecessary to call nf.) findLevels (T True (partition,_) []) = [] findLevels (T False (partition,_) (t:ts)) = let hs = findLevels t -- TODO: It might be better to use the b that is added in t to find the cell that splits cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split in findLevel v hs (zip vs ts) findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs then findLevel v hs vts else let h = find1New t' in findLevel v (h++hs) vts findLevel _ hs [] = hs find1New (T False _ ts) = take 1 $ concatMap find1New ts -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.) find1New (T True (partition,_) []) = let h = fromPairs $ zip (concat partition) vs g' = fmap (.^ h) g in if all (`S.member` es1) (edges g') then [h*g1'] else [] -- isAut h = all (`S.member` es') [e -^ h | e <- es] -- es' = S.fromList es -- Given a partition search tree, return a base and strong generating set for graph automorphism group. partitionBSGS g@(G vs es) t = (bs, findLevels t) where (canonical,bs) = leftLeaf t findLevels (T True (partition,_) []) = [] findLevels (T False (partition,_) (t:ts)) = let hs = findLevels t -- TODO: It might be better to use the b that is added in t to find the cell that splits cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split in findLevel v hs (zip vs ts) findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs -- TODO: Memoize this orbit then findLevel v hs vts else let h = find1New t' in findLevel v (h++hs) vts findLevel _ hs [] = hs find1New (T False _ ts) = take 1 $ concatMap find1New ts -- Some leaf nodes correspond to different isomorphs of the graph, and hence don't yield automorphisms find1New (T True (partition,_) []) = let h = fromPairs $ zip canonical (concat partition) in filter isAut [h] isAut h = all (`S.member` es') [e -^ h | e <- es] es' = S.fromList es -- The tree for g1 has leaf nodes of two different isomorphs, as does the tree for incidenceGraphPG 2 f8 -- Returns auts as Right, different isomorphs as Left -- (Must be used with the tree which doesn't put singletons to end) partitionBSGS3 g@(G vs es) t = (bs, findLevels t) where (p1,bs) = leftLeaf t findLevels (T True (partition,_) []) = [] findLevels (T False (partition,_) (t:ts)) = let hs = findLevels t -- TODO: It might be better to use the b that is added in t to find the cell that splits cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split in findLevel v hs (zip vs ts) findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ rights hs then findLevel v hs vts else let h = find1New t' in findLevel v (h++hs) vts findLevel _ hs [] = hs find1New (T False _ ts) = take 1 $ concatMap find1New ts -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.) find1New (T True (partition,_) []) = let h = fromPairs $ zip p1 (concat partition) in if isAut h then [Right h] else [Left h] isAut h = all (`S.member` es') [e -^ h | e <- es] es' = S.fromList es -- TODO: I think we are only justified in doing find1New (ie only finding 1) if we *do* find an aut. -- If we don't, we should potentially keep looking in that subtree -- (See section 6 of paper. If we find isomorphic leaves, then the two subtrees of their common parent are isomorphic, -- so no need to continue searching the second.) -- This is using a node invariant to do more pruning. -- However, seems to be much slower on very regular graphs (where perhaps there is no pruning to be done) -- (This suggests that perhaps using fmap is not good - perhaps a space leak?) -- (Or perhaps it's just that calculating and comparing the node invariants is expensive) -- TODO: Perhaps use something simpler, like just the number of cells in the partition partitionBSGS2 g@(G vs es) t = (bs, findLevels t') where t' = fmap (\(p,bs) -> (p,bs,trace1 p)) t trace1 = length -- the number of cells in the partition (canonical,bs) = leftLeaf t findLevels (T True (partition,_,_) []) = [] findLevels (T False (partition,_,_) (t:ts)) = let (T _ (_,_,trace) _) = t hs = findLevels t -- TODO: It might be better to use the b that is added in t to find the cell that splits cell@(v:vs) = head $ filter (not . isSingleton) partition -- the cell that is going to split vts = filter (\(_,T _ (_,_,trace') _) -> trace == trace') $ zip vs ts in findLevel v hs vts findLevel v hs ((v',t'):vts) = if v' `elem` v .^^ hs then findLevel v hs vts else let h = find1New t' in findLevel v (h++hs) vts findLevel _ hs [] = hs find1New (T False _ ts) = take 1 $ concatMap find1New ts -- There is a leaf for every aut, but not necessarily an aut for every leaf, so we must check we have an aut -- (For example, incidenceGraphPG 2 f8 has leaf nodes which do not correspond to auts.) -- (The graph g1, below, shows a simple example where this will happen.) find1New (T True (partition,_,_) []) = let h = fromPairs $ zip canonical (concat partition) in filter isAut [h] isAut h = all (`S.member` es') [e -^ h | e <- es] es' = S.fromList es graphAuts7 g = (partitionBSGS g) (equitablePartitionGraphSearchTree g) -- This is faster on kneser graphs, but slower on incidenceGraphPG graphAuts8 g = (partitionBSGS g) (equitableDistancePartitionSearchTree g [vertices g]) -- This is a graph where the node invariant should cause pruning. -- The initial equitable partition will be [[1..8],[9,10]], because it can do no better than distinguish by degree -- However, vertices 1..4 and vertices 5..8 are in fact different (there is no aut that takes one set to the other), -- so the subtrees starting 1..4 have a different invariant to those starting 5..8 g1 = G [1..10] [[1,2],[1,3],[1,9],[2,4],[2,10],[3,4],[3,9],[4,10],[5,6],[5,8],[5,9],[6,7],[6,10],[7,8],[7,9],[8,10]] g1' = nf $ fmap (\x -> if x <= 4 then x+4 else if x <= 8 then x-4 else x) g1 -- G [1..10] [[1,2],[1,4],[1,9],[2,3],[2,10],[3,4],[3,9],[4,10],[5,6],[5,7],[5,9],[6,8],[6,10],[7,8],[7,9],[8,10]] g2 = G [1..12] [[1,2],[1,4],[1,11],[2,3],[2,12],[3,4],[3,11],[4,12],[5,6],[5,8],[5,11],[6,9],[6,12],[7,8],[7,10],[7,11],[8,12],[9,10],[9,11],[10,12]] -- NOT WORKING: This fails to find the isomorphism between g1 and g1' above. -- Instead of using left leaf, we need to find the canonical isomorph, as described in the paper. -- (In a graph where not all leaves lead to automorphisms, we might happen to end up with non-isomorphic left leaves) maybeGraphIso g1 g2 = let (vs1,_) = (leftLeaf . equitablePartitionGraphSearchTree) g1 (vs2,_) = (leftLeaf . equitablePartitionGraphSearchTree) g2 f = M.fromList (zip vs1 vs2) in if length vs1 == length vs2 && (nf . fmap (f M.!)) g1 == g2 then Just f else Nothing -- AUTS OF INCIDENCE STRUCTURE VIA INCIDENCE GRAPH -- This code is nearly identical to the corresponding graphAuts code, with two exceptions: -- 1. We start by partitioning into lefts and rights. -- This avoids left-right crossover auts, which while they are auts of the graph, -- are not auts of the incidence structure -- 2. When labelling the nodes, we filter out Right blocks, and unLeft the Left points incidenceAutsDistancePartitionSearchTree g@(G vs es) = dfs [] (lrPart, lrPart) where dfs xys (srcPart,trgPart) | all isSingleton srcPart = let xys' = zip (concat srcPart) (concat trgPart) in T (isCompatible xys') (unLeft $ xys++xys') [] -- Since the xys' are distance-compatible with the xys, they are certainly edge-compatible. -- However, we do need to check that the xys' are edge-compatible with each other. | otherwise = let (x:xs):srcCells = srcPart yys :trgCells = trgPart srcPart' = intersectCells (xs : srcCells) (dps M.! x) in T False (unLeft xys) -- the L.sort in the following line is so that we traverse vertices in natural order [dfs ((x,y):xys) ((unzip . L.sort) (zip (filter (not . null) srcPart') (filter (not . null) trgPart'))) | (y,ys) <- picks yys, let trgPart' = intersectCells (ys : trgCells) (dps M.! y), map length srcPart' == map length trgPart'] isCompatible xys = and [([x,x'] `S.member` es') == (L.sort [y,y'] `S.member` es') | (x,y) <- xys, (x',y') <- xys, x < x'] (lefts, rights) = partitionEithers vs lrPart = [map Left lefts, map Right rights] -- Partition the vertices into left and right, to exclude crossover auts unLeft xys = [(x,y) | (Left x, Left y) <- xys] -- also filters out Rights es' = S.fromList es dps = M.fromList [(v, distancePartitionS vs es' v) | v <- vs] -- |Given the incidence graph of an incidence structure between points and blocks -- (for example, a set system), -- @incidenceAuts g@ returns a strong generating set 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. incidenceAuts :: (Ord p, Ord b) => Graph (Either p b) -> [Permutation p] incidenceAuts = filter (/= p []) . strongTerminals . incidenceAutsDistancePartitionSearchTree -- TODO: Filter out rights, map unLeft - to bs and gs incidenceAuts2 g = (partitionBSGS g) (equitablePartitionIncidenceSearchTree g) where unLeft (Left x) = x -- map (\g -> fromPairs . map (\(Left x, Left y) -> (x,y)) . filter (\(x,y) -> isLeft x) . toPairs) gs -- 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') <- picks cjs, iso <- graphIsosCon ci cj, 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) (intersectCells (xs : p1'') (dps1 M.! x)) (intersectCells (ys': p2'') (dps2 M.! y)) | (y,ys') <- picks 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.fromAscList [(v, distancePartitionS vs1 es1 v) | v <- vs1] dps2 = M.fromAscList [(v, distancePartitionS vs2 es2 v) | v <- vs2] -- dps1 = M.fromList [(v, distancePartition g1 v) | v <- vertices g1] -- dps2 = M.fromList [(v, distancePartition g2 v) | v <- vertices g2] vs1 = vertices g1 vs2 = 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) -- 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') <- picks cjs, iso <- incidenceIsosCon ci cj, 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) (intersectCells (xs : p1'') (dps1 M.! x)) (intersectCells (ys': p2'') (dps2 M.! y)) | (y,ys') <- picks 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)