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

module Math.Combinatorics.Graph where

import qualified Data.List as L
import Data.Maybe (isJust)
import qualified Data.Map as M
import qualified Data.Set as S
import Control.Arrow ( (&&&) )

import Math.Common.ListSet
import Math.Algebra.Group.PermutationGroup
import Math.Algebra.Group.SchreierSims as SS

-- Main source: Godsil & Royle, Algebraic Graph Theory

-- COMBINATORICS
-- Some functions we'll use

set xs = map head \$ L.group \$ L.sort xs

-- subsets of a set (returned in "binary" order)
powerset [] = [[]]
powerset (x:xs) = let p = powerset xs in p ++ map (x:) p

-- subsets of size k (returned in ascending order)
combinationsOf 0 _ = [[]]
combinationsOf _ [] = []
combinationsOf k (x:xs) = map (x:) (combinationsOf (k-1) xs) ++ combinationsOf k xs

-- GRAPH

data Graph a = G [a] [[a]] deriving (Eq,Ord,Show)

-- we require that vs, es, and each individual e are sorted
isSetSystem xs bs = isListSet xs && isListSet bs && all isListSet bs && all (`isSubset` xs) bs

isGraph vs es = isSetSystem vs es && all ( (==2) . length) es

graph (vs,es) | isGraph vs es = G vs es
-- isValid g = g where g = G vs es

toGraph (vs,es) | isGraph vs' es' = G vs' es' where
vs' = L.sort vs
es' = L.sort \$ map L.sort es
-- note that calling isListSet on a sorted list still checks that there are no duplicates

vertices (G vs _) = vs

edges (G _ es) = es

-- OTHER REPRESENTATIONS

-- incidence matrix of a graph
-- (rows and columns indexed by edges and vertices respectively)
-- (warning: in the literature it is often the other way round)
incidenceMatrix (G vs es) = [ [if v `elem` e then 1 else 0 | v <- vs] | e <- es]

fromIncidenceMatrix m = graph (vs,es) where
n = L.genericLength \$ head m
vs = [1..n]
es = L.sort \$ map edge m
edge row = [v | (1,v) <- zip row vs]

[ [if L.sort [i,j] `S.member` es' then 1 else 0 | j <- vs] | i <- vs]
where es' = S.fromList es

fromAdjacencyMatrix m = graph (vs,es) where
n = L.genericLength m
vs = [1..n]
es = es' 1 m
es' i (r:rs) = [ [i,j] | (j,1) <- drop i (zip vs r)] ++ es' (i+1) rs
es' _ [] = []

-- SOME SIMPLE FAMILIES OF GRAPHS

nullGraph :: Graph Int -- type signature needed
nullGraph = G [] []

-- cyclic graph
c n = graph (vs,es) where
vs = [1..n]
es = L.insert [1,n] [[i,i+1] | i <- [1..n-1]]
-- automorphism group is D2n

-- complete graph
k n = graph (vs,es) where
vs = [1..n]
es = [[i,j] | i <- [1..n-1], j <- [i+1..n]] -- == combinationsOf 2 [1..n]
-- automorphism group is Sn

-- complete bipartite graph
kb m n = to1n \$ kb' m n

kb' m n = graph (vs,es) where
vs = map Left [1..m] ++ map Right [1..n]
es = [ [Left i, Right j] | i <- [1..m], j <- [1..n] ]
-- automorphism group is Sm*Sn (plus a flip if m==n)

-- k-cube
q k = let vs = zip [0..] (powerset [1..k])
es = [ [i,j] | (i,iset) <- vs, (j,jset) <- vs, i < j, length (iset `symDiff` jset) == 1 ]
in graph (map fst vs,es)

q' k = let us = powerset \$ map (2^) [0..k-1]
vs = [0..2^k-1] -- == L.sort \$ map sum us
es = L.sort [ L.sort [sum u, sum v] | [u,v] <- combinationsOf 2 us, length (u `symDiff` v) == 1 ]
in graph (vs, es)

tetrahedron = k 4

cube = q 3

octahedron = graph (vs,es) where
vs = [1..6]
es = combinationsOf 2 vs L.\\ [[1,6],[2,5],[3,4]]

dodecahedron = toGraph (vs,es) where
vs = [1..20]
es = [ [1,2],[2,3],[3,4],[4,5],[5,1],
[6,7],[7,8],[8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15],[15,6],
[16,17],[17,18],[18,19],[19,20],[20,16],
[1,6],[2,8],[3,10],[4,12],[5,14],
[7,16],[9,17],[11,18],[13,19],[15,20] ]

icosahedron = toGraph (vs,es) where
vs = [1..12]
es = [ [1,2],[1,3],[1,4],[1,5],[1,6],
[2,3],[3,4],[4,5],[5,6],[6,2],
[7,12],[8,12],[9,12],[10,12],[11,12],
[7,8],[8,9],[9,10],[10,11],[11,7],
[2,7],[7,3],[3,8],[8,4],[4,9],[9,5],[5,10],[10,6],[6,11],[11,2] ]

-- convert a graph to have [1..n] as vertices
to1n (G vs es) = graph (vs',es') where
mapping = M.fromList \$ zip vs [1..] -- the mapping from vs to [1..n]
vs' = M.elems mapping
es' = [map (mapping M.!) e | e <- es] -- the edges will already be sorted correctly by construction

-- NEW GRAPHS FROM OLD

complement (G vs es) = graph (vs,es') where es' = combinationsOf 2 vs \\ es
-- es' = [e | e <- combinationsOf 2 vs, e `notElem` es]

lineGraph g = to1n \$ lineGraph' g

lineGraph' (G vs es) = graph (es, [ [ei,ej] | ei <- es, ej <- dropWhile (<= ei) es, ei `intersect` ej /= [] ])

petersen = complement \$ lineGraph \$ k 5

-- SIMPLE PROPERTIES OF GRAPHS

order g = length (vertices g)

size g = length (edges g)

-- also called degree
valency (G vs es) v = length \$ filter (v `elem`) es

valencies g@(G vs es) = map (head &&& length) \$ L.group \$ L.sort \$ map (valency g) vs

regularParam g =
case valencies g of
[(v,_)] -> Just v
_       -> Nothing

isRegular g = isJust \$ regularParam g

isCubic g = regularParam g == Just 3

nbrs (G vs es) v = [u | [u,v'] <- es, v == v']
++ [w | [v',w] <- es, v == v']
-- if the graph is valid, then the neighbours will be returned in ascending order

-- find paths from x to y using bfs
-- by definition, a path is a subgraph isomorphic to a "line" - it can't have self-crossings
-- (a walk allows self-crossings, a trail allows self-crossings but no edge reuse)
findPaths g@(G vs es) x y = map reverse \$ bfs [ [x] ] where
bfs ((z:zs) : nodes)
| z == y    = (z:zs) : bfs nodes
| otherwise = bfs (nodes ++ [(w:z:zs) | w <- nbrs g z, w `notElem` zs])
bfs [] = []

-- length of the shortest path from x to y
distance g x y =
case findPaths g x y of
[] -> -1 -- infinite
p:ps -> length p - 1

-- diameter of a graph is maximum distance between two distinct vertices
diameter g@(G vs es)
| isConnected g = maximum \$ map maxDistance vs
| otherwise = -1
where maxDistance v = length (distancePartition g v) - 1

-- find cycles starting at x
-- by definition, a cycle is a subgraph isomorphic to a cyclic graph - it can't have self-crossings
-- (a circuit allows self-crossings but not edge reuse)
findCycles g@(G vs es) x = [reverse (x:z:zs) | z:zs <- bfs [ [x] ], z `elem` nbrsx, length zs > 1] where
nbrsx = nbrs g x
bfs ((z:zs) : nodes) = (z:zs) : bfs (nodes ++ [ w:z:zs | w <- nbrs g z, w `notElem` zs])
bfs [] = []

-- girth of a graph is the size of the smallest cycle it contains
-- Note: If graph contains no cycles, we return -1, representing infinity
girth g@(G vs es) = minimum' \$ map minCycle vs where
minimum' xs = let (zs,nzs) = L.partition (==0) xs in if null nzs then -1 else minimum nzs
minCycle v = case findCycles g v of
[] -> 0
c:cs -> length c - 1 -- because v occurs twice in c, as startpoint and endpoint

-- circumference = max cycle - Bollobas p104

distancePartition g v = distancePartition' S.empty (S.singleton v) where
distancePartition' interior boundary
| S.null boundary = []
| otherwise = let interior' = S.union interior boundary
boundary' = foldl S.union S.empty [S.fromList (nbrs g x) | x <- S.toList boundary] S.\\ interior'
in S.toList boundary : distancePartition' interior' boundary'

-- the connected component to which v belongs
component g v = concat \$ distancePartition g v

isConnected g@(G (v:vs) es) = length (component g v) == length (v:vs)
isConnected (G [] []) = True

-- MORE GRAPHS

-- Generalized Johnson graph, Godsil & Royle p9
j v k i | v >= k && k >= i
= graph (vs,es) where
vs = combinationsOf k [1..v]
es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, length (v1 `intersect` v2) == i ]
-- j v k i is isomorphic to j v (v-k) (v-2k+i), so may as well have v >= 2k

kneser v k | v >= 2*k = j v k 0

johnson v k | v >= 2*k = j v k (k-1)

petersen1 = to1n \$ j 5 2 0

```