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

{-# LANGUAGE NoMonomorphismRestriction, DeriveFunctor #-}

-- |A module providing functions to construct and investigate (small, finite) matroids.
module Math.Combinatorics.Matroid where

-- Source: Oxley, Matroid Theory (second edition)

import Math.Core.Utils
import Math.Core.Field hiding (f7)

import Math.Common.ListSet as LS -- set operations on strictly ascending lists
-- import Math.Algebra.Field.Base hiding (Q, F2, F3, F5, F7, F11, f2, f3, f5, f7, f11)
import Math.Algebra.LinearAlgebra hiding (rank, (*>))

import qualified Math.Combinatorics.Graph as G -- hiding (combinationsOf, restriction, component, isConnected)
import Math.Combinatorics.FiniteGeometry
import Math.Combinatorics.GraphAuts
import Math.Algebra.Group.PermutationGroup hiding (closure)

import Math.Algebras.VectorSpace hiding (dual)
import Math.Algebras.Structures
import Math.Algebras.Commutative

import Data.List as L
import qualified Data.Map as M
import qualified Data.Set as S

implies p q = q || not p
exists = not . null
unique [x] = x

shortlex xs ys = case compare (length xs) (length ys) of
    LT -> LT
    EQ -> compare xs ys
    GT -> GT

isShortlex xs = foldcmpl (\x1 x2 -> shortlex x1 x2 /= GT) xs

toShortlex xs = map snd $ L.sort [(length x, x) | x <- xs]

isClutter ss = and [ (s1 `LS.isSubset` s2) `implies` (s1 == s2) | s1 <- ss, s2 <- ss ]
-- note that this definition would allow duplicates

-- single-element deletions of xs
-- if xs is sorted, then so is reverse (deletions xs)
deletions xs = zipWith (++) (inits xs) (tail $ tails xs)

closedUnderSubsets xss = and [xs' `S.member` xss' | xs <- xss, xs' <- deletions xs]
    where xss' = S.fromList xss


-- |The data structure that we use to store the bases of the matroid
data TrieSet a = TS [(a, TrieSet a)] deriving (Eq,Ord,Functor)
-- Note that in a trie we would normally have a Bool at each node saying whether or not the node is terminal
-- (ie corresponds to a member and not just a prefix of a member).
-- However, since we intend to use the trie to store the bases, and they are all the same length,
-- we can use lack of children to detect that a node is terminal.

-- We could try storing the children in a map rather than a list

-- for debugging
tsshow (TS xts) = "TS [" ++ concatMap (\(x,t) -> "(" ++ show x ++ "," ++ tsshow t ++ ")") xts ++ "]"

instance Show a => Show (TrieSet a) where
    show = show . tstolist

tsempty = TS []

tsinsert (x:xs) (TS ts) =
    case L.lookup x ts of
    Nothing -> let t = tsinsert xs (TS [])
               in TS $ L.insert (x,t) ts
    Just t -> let t' = tsinsert xs t
              in TS $ L.insert (x,t') $ L.delete (x,t) ts
tsinsert [] t = t

tsmember (x:xs) (TS ts) =
    case lookup x ts of
    Nothing -> False
    Just t -> tsmember xs t
tsmember [] (TS []) = True -- the node has no children, and hence is terminal
tsmember [] _ = False

-- xs is a subset of a member of t
tssubmember (x:xs) (TS ts) = or [ case compare x y of
                                  LT -> False
                                  EQ -> tssubmember xs t
                                  GT -> tssubmember (x:xs) t
                                | (y,t) <- ts ]
tssubmember [] _ = True

tstolist (TS []) = [[]]
tstolist (TS xts) = concatMap (\(x,t) -> map (x:) (tstolist t)) xts 

tsfromlist = foldl' (flip tsinsert) tsempty


-- |A datatype to represent a matroid. @M es bs@ is the matroid whose elements are @es@, and whose bases are @bs@.
-- The normal form is for the @es@ to be in order, for each of the @bs@ individually to be in order.
-- (So the TrieSet should have the property that any path from the root to a leaf is strictly increasing).
data Matroid a = M [a] (TrieSet a) deriving (Eq,Show,Functor)

-- |Return the elements over which the matroid is defined.
elements :: Matroid t -> [t]
elements (M es bs) = es

-- |Return all the independent sets of a matroid, in shortlex order.
indeps :: (Ord a) => Matroid a -> [[a]]
indeps m = bfs [ ([],es) ]
    where es = elements m
          bfs ( (ls,rs) : nodes ) =
              let ls' = reverse ls in
              if isIndependent m ls'
              then ls' : bfs ( nodes ++ successors (ls,rs) )
              else bfs nodes
          bfs [] = []
          successors (ls,rs) = [ (r:ls, rs') | r:rs' <- L.tails rs ]

isIndependent :: (Ord a) => Matroid a -> [a] -> Bool
isIndependent (M es bs) xs = xs `tssubmember` bs

isDependent :: (Ord a) => Matroid a -> [a] -> Bool
isDependent m = not . isIndependent m

-- |Are these the independent sets of a matroid? (The sets must individually be ordered.)
isMatroidIndeps :: (Ord a) => [[a]] -> Bool
isMatroidIndeps is =
    [] `elem` is &&
    closedUnderSubsets is &&
    and [ (l1 < l2) `implies` exists [e | e <- i2 LS.\\ i1, L.insert e i1 `elem` is]
        | i1 <- is, let l1 = length i1, i2 <- is, let l2 = length i2 ]

-- |Construct a matroid from its elements and its independent sets.
fromIndeps :: (Ord a) => [a] -> [[a]] -> Matroid a
fromIndeps es is = fromBases es bs
    where bs = dfs [] [([],es)]
          dfs bs ( node@(ls,rs) : nodes ) =
              let succs = successors node
              in if null succs then dfs (ls:bs) nodes else dfs bs (succs ++ nodes)
          dfs ls [] = let r = length $ last ls -- we know that the first one we found is a true base
                      in map reverse $ filter (\b -> length b == r) ls
                      -- we might have had null succs simply because we ran out of vectors to add
          successors (ls,rs) = [ (r:ls, rs') | r:rs' <- L.tails rs, (r:ls) `S.member` is' ]
          is' = S.fromList $ map reverse is

-- this seems to be just slightly slower
fromIndeps1 es is = fromBases es bs
    where b = greedy [] es -- first find any basis
          greedy ls (r:rs) = if (r:ls) `S.member` ris'
                             then greedy (r:ls) rs
                             else greedy ls rs
          greedy ls [] = reverse ls
          ris' = S.fromList $ map reverse is
          bs = closure S.empty (S.singleton b) -- now find all other bases by passing to "neighbouring" bases
          closure interior boundary =
              if S.null boundary
              then S.toList interior
              else let interior' = interior `S.union` boundary
                       boundary' = S.fromList [ b' | b <- S.toList boundary,
                                                     x <- b, y <- es LS.\\ b,
                                                     let b' = L.insert y (L.delete x b),
                                                     b' `S.notMember` interior',
                                                     b' `S.member` is' ]
                   in closure interior' boundary'
          is' = S.fromList is
-- The basis properties imply that the set of all bases is connected under the "neighbour" relation


-- |Given a matrix, represented as a list of rows, number the columns [1..],
-- and construct the matroid whose independent sets correspond to those sets of columns which are linearly independent
-- (or in case there are repetitions, those multisets of columns which are sets, and which are linearly independent).
vectorMatroid :: (Fractional k) => [[k]] -> Matroid Int
vectorMatroid = vectorMatroid' . L.transpose

-- |Given a list of vectors (or rows of a matrix), number the vectors (rows) [1..], and construct the matroid whose independent sets
-- correspond to those sets of vectors (rows) which are linearly independent
-- (or in case there are repetitions, those multisets which are sets, and which are linearly independent).
vectorMatroid' :: (Fractional k) => [[k]] -> Matroid Int
vectorMatroid' vs = fromBases (map fst vs') bs
    where vs' = zip [1..] vs
          bs = dfs [] [([],[],vs')]
          dfs ls (r@(i,ref,es) : rs) =
              let succs = successors r in
              if null succs then dfs (i:ls) (succs ++ rs) else dfs ls (succs ++ rs)
          dfs ls [] = let r = length $ last ls -- we know that the first one we found is a true base
                      in map reverse $ filter (\b -> length b == r) ls
                      -- we might have had null succs simply because we ran out of vectors to add
          successors (i,ref,es) = [(i',ref',es') | (j,e):es' <- L.tails es,
                                                   not (inSpanRE ref e),
                                                   let ref' = rowEchelonForm (ref ++ [e]), -- is this really better than e:ref?
                                                   let i' = j : i ]

-- |Given the edges of an undirected graph, number the edges [1..], and construct the matroid whose independent sets
-- correspond to those sets of edges which contain no cycle. The bases therefore correspond to maximal forests within the graph.
-- The edge set is allowed to contain loops or parallel edges.
cycleMatroid :: (Ord a) => [[a]] -> Matroid Int
cycleMatroid es = fromBases (map fst es') bs
    where es' = zip [1..] es
          bs = dfs [] [([], M.empty, es')]
          dfs ls (r@(i,ref,es) : rs) =
              let succs = successors r
              in if null succs then dfs (i:ls) (succs ++ rs) else dfs ls (succs ++ rs)
          dfs ls [] = let r = length $ last ls -- we know that the first one we found is a true base
                      in map reverse $ filter (\b -> length b == r) ls
                      -- we might have had null succs simply because we ran out of edges to add
          successors (i, reps, (j,[u,v]):es' ) =
              if u == v
              then successors (i, reps, es')
              else case (M.lookup u reps, M.lookup v reps) of
                   (Nothing, Nothing) -> (j:i, M.insert u u $ M.insert v u reps, es') : successors (i, reps, es')
                                         -- neither of these vertices has been seen before, so add this edge as a new tree in the forest
                   (Just u', Nothing) -> (j:i, M.insert v u' reps, es') : successors (i, reps, es')
                                         -- we have seen u before but not v, so add v to the u-tree
                   (Nothing, Just v') -> (j:i, M.insert u v' reps, es') : successors (i, reps, es')
                                         -- we have seen v before but not u, so add u to the v-tree
                   (Just u', Just v') -> if u' == v'
                                         then successors (i,reps,es')
                                              -- u and v are already in the same tree, so adding this edge would make a cycle
                                         else (j:i, M.map (\w -> if w == v' then u' else w) reps, es') : successors (i, reps, es')
                                              -- u and v are in different trees, so join the trees together
          successors (_, _, []) = []

-- A version of cycle matroid that retains the original edges rather than relabeling with integers.
-- Not really valid if there are parallel edges (because they can't be distinguished in the output).
-- Required for "markedfcim" below.
cycleMatroid' es = fmap lookupEdge $ cycleMatroid es
    where table = M.fromList $ zip [1..] es
          lookupEdge = (M.!) table


-- |Given a matroid over an arbitrary type, relabel to obtain a matroid over the integers.
to1n :: (Ord a) => Matroid a -> Matroid Int
to1n m = fmap to1n' m
    where es = elements m
          table = M.fromList $ zip es [1..]
          to1n' = (M.!) table


-- ISOMORPHISMS AND AUTOMORPHISMS

incidenceGraphB m = G.G vs' es'
    where es = elements m
          bs = bases m
          vs' = map Left es ++ map Right bs
          es' = L.sort [ [Left e, Right b] | b <- bs, e <- b ]
-- incidence graph for the matroid considered as an incidence structure between elements and bases

incidenceGraphH m = G.G vs' es'
    where es = elements m
          hs = L.sort $ hyperplanes m
          vs' = map Left es ++ map Right hs
          es' = L.sort [ [Left e, Right h] | h <- hs, e <- h ]
-- incidence graph for the matroid considered as an incidence structure between elements and hyperplanes
-- for "sparse" matroids, there are likely to be fewer hyperplanes than bases (why?)
-- incidenceGraphH may not be connected - eg u 2 4

-- (So a rank 3 or higher matroid, provided it has more than one hyperplane, certainly has a connected incidenceGraphH)

matroidIsos m1 m2 = incidenceIsos (incidenceGraphH m1) (incidenceGraphH m2)

-- |Are the two matroids isomorphic?
isMatroidIso :: (Ord a, Ord b) => Matroid a -> Matroid b -> Bool
isMatroidIso m1 m2 = isIncidenceIso (incidenceGraphH m1) (incidenceGraphH m2)

-- |Return the automorphisms of the matroid.
matroidAuts :: (Ord a) => Matroid a -> [Permutation a]
matroidAuts m
    | G.isConnected hgraph = incidenceAuts hgraph
    | G.isConnected bgraph = incidenceAuts bgraph
    | otherwise = error "matroidAuts: incidence graph is not connected"
    where hgraph = incidenceGraphH m
          bgraph = incidenceGraphB m
-- Note that the results aren't always what one intuitively expects from the geometric representation.
-- This is because geometric representations suggest additional structure beyond matroid structure.
-- For example, for the Vamos matroid v8,
-- it returns auts which are not "geometric" auts of the geometric representation.
-- Matroids are really combinatorial objects, not geometric.
-- In the case of v8, the key is that the only thing the auts have to preserve are the planes, not the lines


-- CIRCUITS

-- |A circuit in a matroid is a minimal dependent set.
isCircuit :: (Ord a) => Matroid a -> [a] -> Bool
isCircuit m c =
    isDependent m c &&
    all (isIndependent m) (deletions c)

-- |Return all circuits for the given matroid, in shortlex order.
circuits :: (Ord a) => Matroid a -> [[a]]
circuits m = toShortlex $ dfs S.empty [L.insert e b | b <- bs, e <- es LS.\\ b]
    where es = elements m
          bs = bases m
          dfs vs (c:cs) | c `S.member` vs = dfs vs cs
                        | otherwise = let cs' = successors c
                                          vs' = S.insert c vs
                                      in if null cs' then c : dfs vs' cs else dfs vs' (cs' ++ cs)
          dfs _ [] = []
          successors c = [c' | c' <- deletions c, isDependent m c' ]

-- Oxley p10
-- |Are the given sets the circuits of some matroid?
isMatroidCircuits :: (Ord a) => [[a]] -> Bool
isMatroidCircuits cs =
    [] `notElem` cs &&
    and [(c1 `LS.isSubset` c2) `implies` (c1 == c2) | c1 <- cs, c2 <- cs] &&
    and [ exists [c3 | c3 <- cs, c3 `LS.isSubset` c12']
        | c1 <- cs, c2 <- cs, c1 /= c2,
          e <- c1 `LS.intersect` c2, let c12' = L.delete e (c1 `LS.union` c2)] 

-- |Reconstruct a matroid from its elements and circuits.
fromCircuits :: (Ord a) => [a] -> [[a]] -> Matroid a
fromCircuits es cs = fromBases es bs
    where b = greedy [] es -- first find any basis
          greedy ls (r:rs) = let ls' = ls ++ [r] in
                             if isIndep ls'
                             then greedy ls' rs
                             else greedy ls rs
          greedy ls [] = ls
          bs = closure S.empty (S.singleton b) -- now find all other bases by passing to "neighbouring" bases
          closure interior boundary =
              if S.null boundary
              then S.toList interior
              else let interior' = interior `S.union` boundary
                       boundary' = S.fromList [ b' | b <- S.toList boundary,
                                                     x <- b, y <- es LS.\\ b,
                                                     let b' = L.insert y (L.delete x b),
                                                     b' `S.notMember` interior',
                                                     isIndep b' ]
                   in closure interior' boundary'
          isIndep xs = not (any (`LS.isSubset` xs) cs)


-- |An element e in a matroid M is a loop if {e} is a circuit of M.
isLoop :: (Ord a) => Matroid a -> a -> Bool
isLoop m e = isCircuit m [e]
-- isLoop (M es is) e = [e] `notElem` is

-- |Elements f and g in a matroid M are parallel if {f,g} is a circuit of M.
isParallel :: (Ord a) => Matroid a -> a -> a -> Bool
isParallel m f g = isCircuit m [f,g]

-- |A matroid is simple if it has no loops or parallel elements
isSimple :: (Ord a) => Matroid a -> Bool
isSimple m = all ( (>2) . length ) (circuits m)


-- BASES

-- |A base or basis in a matroid is a maximal independent set.
isBase :: (Ord a) => Matroid a -> [a] -> Bool
isBase (M es bs) b = b `tsmember` bs

-- |Return all bases for the given matroid
bases :: (Ord a) => Matroid a -> [[a]]
bases (M es bs) = tstolist bs

-- |Are the given sets the bases of some matroid?
isMatroidBases :: (Ord a) => [[a]] -> Bool
isMatroidBases bs =
    (not . null) bs &&
    and [ exists [y | y <- b2 LS.\\ b1, L.insert y (L.delete x b1) `elem` bs]
        | b1 <- bs, b2 <- bs, x <- b1 LS.\\ b2 ]

-- |Reconstruct a matroid from its elements and bases.
fromBases :: (Ord a) => [a] -> [[a]] -> Matroid a
fromBases es bs = M es (tsfromlist bs)
-- The elements are required because a loop does not appear in any basis.

-- Oxley p17
-- |Given a matroid m, a basis b, and an element e, @fundamentalCircuit m b e@ returns the unique circuit contained in b union {e},
-- which is called the fundamental circuit of e with respect to b.
fundamentalCircuit :: (Ord a) => Matroid a -> [a] -> a -> [a]
fundamentalCircuit m b e = unique [c | c <- circuits m, c `LS.isSubset` be]
    where be = L.insert e b

uniformMatroid m n | m <= n = fromBases es bs
    where es = [1..n]
          bs = combinationsOf m es

-- |The uniform matroid U m n is the matroid whose independent sets are all subsets of [1..n] with m or fewer elements.
u :: Int -> Int -> Matroid Int
u = uniformMatroid


-- RANK FUNCTION

-- Oxley p103, 3.1.14
restriction1 m xs = fromBases xs bs'
    where bs = bases m
          is' = toShortlex [b `LS.intersect` xs | b <- bs]
          r = length $ last is'
          bs' = dropWhile ( (< r) . length ) is'

-- |The restriction of a matroid to a subset of its elements
restriction :: (Ord a) => Matroid a -> [a] -> Matroid a
restriction m@(M es bs) xs = M xs bs'
    where (_,bs') = balance $ prune bs
          prune (TS yts) = let (ins, outs) = L.partition (\(y,t) -> y `elem` xs) yts
                               ins' = [(y, prune t) | (y,t) <- ins]
                               outs' = concat [ zts | (y,t) <- outs, let TS zts = prune t ]
                           in TS $ ins' ++ outs'
          balance (TS yts) = let dyt's = [(d',(y,t')) | (y,t) <- yts, let (d',t') = balance t]
                                 d = maximum $ 0 : map fst dyt's
                             in (d+1, TS $ toSet [(y,t') | (d',(y,t')) <- dyt's, d' == d])
-- we have to make the toSet call *after* we balance, otherwise two trees may appear unequal because of undergrowth that is going to be removed

-- !! Need thorough testing to prove that restriction == restriction1


-- |Given a matroid m, @rankfun m@ is the rank function on subsets of its element set
rankfun :: (Ord a) => Matroid a -> [a] -> Int
rankfun m xs = (length . head . bases) (restriction m xs)
-- no danger of head [], because bases must be non-null

-- |The rank of a matroid is the cardinality of a basis
rank :: (Ord a) => Matroid a -> Int
rank m = length $ head $ bases m
-- rank m@(M es bs) = rankfun m es

-- Oxley p23
-- |Reconstruct a matroid from its elements and rank function
fromRankfun :: (Ord a) => [a] -> ([a] -> Int) -> Matroid a
fromRankfun es rkf = fromBases es bs
    where b = greedy 0 [] es -- first find any basis
          greedy rk ls (r:rs) = let ls' = ls ++ [r] in
                                if rkf ls' == rk+1
                                then greedy (rk+1) ls' rs
                                else greedy rk ls rs
          greedy _ ls [] = ls
          rk = rkf b
          isBasis b' = rkf b' == rk
          bs = closure S.empty (S.singleton b) S.empty -- now find all other bases by passing to "neighbouring" bases
          closure interior boundary exterior =
              if S.null boundary
              then S.toList interior
              else let interior' = interior `S.union` boundary
                       candidates = S.fromList [ b' | b <- S.toList boundary,
                                                      x <- b, y <- es LS.\\ b,
                                                      let b' = L.insert y (L.delete x b),
                                                      b' `S.notMember` interior',
                                                      b' `S.notMember` exterior ]
                       (boundary', exterior') = S.partition isBasis candidates
                   in closure interior' boundary' (S.union exterior exterior')
-- The purpose of keeping track of exterior points we have encountered
-- is to avoid making repeat calls to the rank function rkf


-- CLOSURE OPERATOR AND FLATS

-- |Given a matroid m, @closure m@ is the closure operator on subsets of its element set
closure :: (Ord a) => Matroid a -> [a] -> [a]
closure m xs = [x | x <- es, x `elem` xs || rankfun m (L.insert x xs) == rankxs]
    where es = elements m
          rankxs = rankfun m xs
-- The intuition is that closure xs is all elements within the span of xs

-- |Reconstruct a matroid from its elements and closure operator
fromClosure :: (Ord a) => [a] -> ([a] -> [a]) -> Matroid a
fromClosure es cl = fromBases es bs
    where b = greedy (cl []) [] es -- first find any basis
          greedy span ls (r:rs) = let ls' = ls ++ [r] in
                                  if r `notElem` span -- r is independent relative to ls
                                  then greedy (cl ls') ls' rs
                                  else greedy span ls rs
          greedy _ ls [] = ls
          rk = length b
          isBasis b' = cl b' == es
          bs = closure S.empty (S.singleton b) S.empty -- now find all other bases by passing to "neighbouring" bases
          closure interior boundary exterior =
              if S.null boundary
              then S.toList interior
              else let interior' = interior `S.union` boundary
                       candidates = S.fromList [ b' | b <- S.toList boundary,
                                                      x <- b, y <- es LS.\\ b,
                                                      let b' = L.insert y (L.delete x b),
                                                      b' `S.notMember` interior',
                                                      b' `S.notMember` exterior ]
                       (boundary', exterior') = S.partition isBasis candidates
                   in closure interior' boundary' (S.union exterior exterior')
-- The purpose of keeping track of exterior points we have encountered
-- is to avoid making repeat calls to the closure operator cl

{-
-- Not quite sure why the following is so much slower than the above
-- May just be because Data.Set is compiled, and this would be comparable if also compiled
fromClosure2 es cl = M es bs
    where b = greedy (cl []) [] es -- first find any basis
          greedy span ls (r:rs) = let ls' = ls ++ [r] in
                                  if r `notElem` span -- r is independent relative to ls
                                  then greedy (cl ls') ls' rs
                                  else greedy span ls rs
          greedy _ ls [] = ls
          rk = length b
          isBasis b' = cl b' == es
          bs = closure tsempty [b] tsempty -- now find all other bases by passing to "neighbouring" bases
          closure interior boundary exterior =
              if null boundary
              then interior
              else let interior' = foldl' (flip tsinsert) interior boundary
                       candidates = [ b' | b <- boundary,
                                           x <- b, y <- es LS.\\ b,
                                           let b' = L.insert y (L.delete x b),
                                           not (b' `tsmember` interior'),
                                           not (b' `tsmember` exterior) ]
                       (boundary', exterior') = L.partition isBasis candidates
                   in closure interior' boundary' (foldl' (flip tsinsert) exterior exterior')
-}

-- |A flat in a matroid is a closed set, that is a set which is equal to its own closure
isFlat :: (Ord a) => Matroid a -> [a] -> Bool
isFlat m xs = closure m xs == xs

flats1 m = [xs | xs <- powersetbfs es, isFlat m xs]
    where es = elements m
-- first, inefficient, implementation

-- given xs, a flat in m, return the flats which cover xs
-- these have the property that they partition es \\ xs
coveringFlats m xs = coveringFlats' (es LS.\\ xs)
    where es = elements m
          coveringFlats' (y:ys) = let zs = closure m (L.insert y xs)
                                  in zs : coveringFlats' (ys LS.\\ zs)
          coveringFlats' [] = []

-- since we are dealing with finite matroids, the lattice of flats is finite, so it has a minimal element
minimalFlat m = head $ filter (isFlat m) $ powersetbfs $ elements m

-- |The flats of a matroid are its closed sets. They form a lattice under inclusion.
flats :: (Ord a) => Matroid a -> [[a]]
flats m = flats' S.empty [minimalFlat m]
    where flats' ls (r:rs) = if r `S.member` ls
                             then flats' ls rs
                             else flats' (S.insert r ls) (rs ++ coveringFlats m r)
          flats' ls [] = toShortlex $ S.toList ls
          -- this is just breadth-first search

-- isMatroidFlats
-- Oxley p31-32
{-
isMatroidFlats es fs =
    es `elem` fs &&
    [ (f1 `LS.intersect` f2) `elem` fs | f1 <- fs, f2 <- fs ] &&
    -- for all flats f, the minimal flats containing f partition es-f
-}

-- |Reconstruct a matroid from its flats. (The flats must be given in shortlex order.)
fromFlats :: (Ord a) => [[a]] -> Matroid a
fromFlats fs | isShortlex fs = fromFlats' fs
             | otherwise = error "fromFlats: flats must be in shortlex order"

fromFlats' fs = fromClosure es cl
    where es = last fs -- es is a flat, and last in shortlex order
          cl xs = head [f | f <- fs, xs `LS.isSubset` f] -- the first flat is minimal, because of shortlex order
-- !! we can probably do better (efficiency-wise)
-- eg by constructing the lattice as a DAG, and climbing up it

-- |A subset of the elements in a matroid is spanning if its closure is all the elements
isSpanning :: (Ord a) => Matroid a -> [a] -> Bool
isSpanning m xs = closure m xs == es
    where es = elements m

-- |A hyperplane is a flat whose rank is one less than that of the matroid
isHyperplane :: (Ord a) => Matroid a -> [a] -> Bool
isHyperplane m xs = isFlat m xs && rankfun m xs == rank m - 1

hyperplanes1 m = [h | h <- flats m, rankfun m h == rk - 1]
    where rk = rank m

-- Oxley p65: h is a hyperplane iff its complement is a cocircuit
hyperplanes :: (Ord a) => Matroid a -> [[a]]
hyperplanes m = toShortlex $ map complement $ cocircuits m
    where es = elements m
          complement cc = es LS.\\ cc
-- This appears to be faster than hyperplanes1

-- Oxley p70
isMatroidHyperplanes :: (Ord a) => [a] -> [[a]] -> Bool
isMatroidHyperplanes es hs =
    es `notElem` hs &&
    isClutter hs &&
    and [ exists [h3 | h3 <- hs, h12e `LS.isSubset` h3] | (h1,h2) <- pairs hs,
                                                          e <- es LS.\\ (LS.union h1 h2),
                                                          let h12e = L.insert e (LS.intersect h1 h2) ]
-- Note that contrary to what one might initially think
-- it does *not* follow from the third condition that every element is in some hyperplane
-- since there might be only one hyperplane, eg fromBases [1,2] [[2]] has [1] as it's only hyperplane

-- Haven't actually proven that the following is always correct
-- but it seems to work
fromHyperplanes1 es hs = fromFlats $ closure S.empty (S.fromList hs)
    where closure interior boundary =
              if S.null boundary
              then (toShortlex $ S.toList interior) ++ [es]
              else let interior' = S.union interior boundary
                       boundary' = S.fromList [ f1 `LS.intersect` f2 | (f1,f2) <- pairs (S.toList boundary) ]
                                   S.\\ interior'
                   in closure interior' boundary'

-- |Reconstruct a matroid from its elements and hyperplanes
fromHyperplanes :: (Ord a) => [a] -> [[a]] -> Matroid a
fromHyperplanes es hs = fromCocircuits es $ map complement hs
    where fromCocircuits es = dual . fromCircuits es
          complement xs = es LS.\\ xs


-- GEOMETRIC REPRESENTATION

-- |Given a list of points in k^n, number the points [1..], and construct the matroid whose independent sets
-- correspond to those sets of points which are affinely independent.
--
-- A multiset of points in k^n is said to be affinely dependent if it contains two identical points,
-- or three collinear points, or four coplanar points, or ... - and affinely independent otherwise.
affineMatroid :: (Fractional k) => [[k]] -> Matroid Int
affineMatroid vs = vectorMatroid' $ map (1:) vs

-- |fromGeoRep returns a matroid from a geometric representation consisting of dependent flats of various ranks.
-- Given lists of dependent rank 0 flats (loops), rank 1 flats (points), rank 2 flats (lines) and rank 3 flats (planes),
-- @fromGeoRep loops points lines planes@ returns the matroid having these as dependent flats.
-- Note that if all the elements lie in the same plane, then this should still be listed as an argument.
fromGeoRep :: (Ord a) => [[a]] -> [[a]] -> [[a]] -> [[a]] -> Matroid a
fromGeoRep loops points lines planes =
    fromCircuits es $ minimal $
        loops ++
        concatMap (combinationsOf 2) points ++
        concatMap (combinationsOf 3) lines ++
        concatMap (combinationsOf 4) planes ++
        combinationsOf 5 es
    where es = toSet $ concat loops ++ concat points ++ concat lines ++ concat planes

-- Note that we don't check that the inputs are valid

-- xss assumed to be in shortlex order
minimal xss = minimal' [] xss
    where minimal' ls (r:rs) = if any (`LS.isSubset` r) ls
                               then minimal' ls rs
                               else minimal' (r:ls) rs
          minimal' ls [] = reverse ls

-- |A simple matroid has no loops or parallel elements, hence its geometric representation has no loops or dependent points.
-- @simpleFromGeoRep lines planes@ returns the simple matroid having these dependent flats.
simpleFromGeoRep :: (Ord a) => [[a]] -> [[a]] -> Matroid a
simpleFromGeoRep lines planes = fromGeoRep [] []  lines planes

-- Oxley p37-8
isSimpleGeoRep lines planes =
    all ( (<= 1) . length ) [ l1 `LS.intersect` l2 | (l1,l2) <- pairs lines ] &&
    all ( \i -> length i <= 2 || i `elem` lines ) [ p1 `LS.intersect` p2 | (p1,p2) <- pairs planes ] &&
    and [ any (u `LS.isSubset`) planes | (l1,l2) <- pairs lines, length (l1 `LS.intersect` l2) == 1, let u = l1 `LS.union` l2 ] &&
    and [ length i == 1 || i == l | l <- lines, p <- planes, let i = l `LS.intersect` p ]


isCircuitHyperplane m xs = isCircuit m xs && isHyperplane m xs

-- |List the circuit-hyperplanes of a matroid.
circuitHyperplanes :: (Ord a) => Matroid a -> [[a]]
circuitHyperplanes m = [ h | h <- hyperplanes m, isCircuit m h ] 

-- Oxley p39
-- |Given a matroid m, and a set of elements b which is both a circuit and a hyperplane in m,
-- then @relaxation m b@ is the matroid which is obtained by adding b as a new basis.
-- This corresponds to removing b from the geometric representation of m.
relaxation :: (Ord a) => Matroid a -> [a] -> Matroid a
relaxation m b
    | isCircuitHyperplane m b = fromBases es bs
    | otherwise = error "relaxation: not a circuit-hyperplane"
    where es = elements m
          bs = b : bases m 


-- TRANSVERSAL MATROIDS

ex161 = [ [1,2,6], [3,4,5,6], [2,3], [2,4,6] ]

-- the edges of the bipartite graph
transversalGraph as = [(Left x, Right i) | (a,i) <- zip as [1..], x <- a]

partialMatchings es = dfs [(S.empty, [], es)]
    where dfs (node@(vs,ls,rs): nodes) = ls : dfs (successors node ++ nodes)
          dfs [] = []
          successors (vs,ls,rs) =  [ (S.insert u $ S.insert v vs, L.insert r ls, rs')
                                   | r@(u,v):rs' <- L.tails rs, u `S.notMember` vs, v `S.notMember` vs ]

-- |Given a set of elements es, and a sequence as = [a1,...,am] of subsets of es,
-- return the matroid whose independent sets are the partial transversals of the as.
transversalMatroid :: (Ord a) => [a] -> [[a]] -> Matroid a
transversalMatroid es as = fromBases es bs
    where is@(i:_) = reverse $ toShortlex $ toSet $ (map . map) (unLeft . fst) $ partialMatchings (transversalGraph as)
          unLeft (Left x) = x
          l = length i
          bs = reverse $ takeWhile ( (== l) . length ) is
-- In this case, as is called a presentation of the matroid
-- Note that there may be partial transversals even if there are no transversals


-- Not obvious how to efficiently find the bases without finding the independent sets,
-- since neighbouring partial transversals might arise in different ways
-- (Although Oxley p93 seems to be saying that they don't)
{-
transversalMatroid2 es as = fromBases es bs
    where es' = transversalGraph as
          b = greedy [] es' -- first find any basis
          greedy ls (r@(u,v):rs) = if u `notElem` map fst ls && v `notElem` map snd ls
                                   then greedy (r:ls) rs
                                   else greedy ls rs
          greedy ls [] = ls
-- now choose the subset of the as indicated by the Right part of b
-- and seek all full transversals
          bs = closure S.empty (S.singleton b) -- now find all other bases by passing to "neighbouring" bases
          closure interior boundary =
              if S.null boundary
              then S.toList interior
              else let interior' = interior `S.union` boundary
                       boundary' = S.fromList [ b' | b <- S.toList boundary,
                                                     x <- b, y <- es LS.\\ b,
                                                     let b' = L.insert y (L.delete x b),
                                                     b' `S.notMember` interior',
                                                     cl b' == es ] -- b' is spanning, and same size as b, hence must be basis
                                   -- S.\\ interior'
                   in closure interior' boundary'
-}


-- DUALITY

-- |The dual matroid
dual :: (Ord a) => Matroid a -> Matroid a
dual m = fromBases es bs'
    where es = elements m
          bs = bases m
          bs' = L.sort $ map (es LS.\\) bs

isCoindependent m xs = isIndependent (dual m) xs

isCobase m xs = isBase (dual m) xs
-- quicker but less clear to calculate this directly

isCocircuit m xs = isCircuit (dual m) xs

cocircuits :: (Ord a) => Matroid a -> [[a]]
cocircuits m = circuits (dual m)

isColoop m e = isLoop (dual m) e

isCoparallel m f g = isParallel (dual m) f g


-- MINORS

deletion :: (Ord a) => Matroid a -> [a] -> Matroid a
deletion m xs = restriction m (es LS.\\ xs)
    where es = elements m

(\\\) = deletion

contraction :: (Ord a) => Matroid a -> [a] -> Matroid a
contraction m xs = dual (deletion (dual m) xs)

(///) = contraction


-- CONNECTIVITY

-- Oxley p120
-- |A matroid is (2-)connected if, for every pair of distinct elements, there is a circuit containing both
isConnected :: (Ord a) => Matroid a -> Bool
isConnected m = and [any (pair `LS.isSubset`) cs | pair <- combinationsOf 2 es]
    where es = elements m
          cs = circuits m

component m x = closure S.empty (S.singleton x)
    where cs = circuits m
          closure interior boundary =
              if S.null boundary
              then S.toList interior
              else let interior' = S.union interior boundary
                       boundary' = S.fromList (concat [c | c <- cs, (not . null) (LS.intersect c (S.toList boundary)) ])
                                   S.\\ interior'
                   in closure interior' boundary'

-- |The direct sum of two matroids
dsum :: (Ord a, Ord b) => Matroid a -> Matroid b -> Matroid (Either a b)
dsum m1 m2 = fromBases es bs
    where es = map Left (elements m1) ++ map Right (elements m2)
          bs = [map Left b1 ++ map Right b2 | b1 <- bases m1, b2 <- bases m2]


-- REPRESENTABILITY

-- |@matroidPG n fq@ returns the projective geometry PG(n,Fq), where fq is a list of the elements of Fq
matroidPG :: (Fractional a) => Int -> [a] -> Matroid Int
matroidPG n fq = vectorMatroid' $ ptsPG n fq

-- |@matroidAG n fq@ returns the affine geometry AG(n,Fq), where fq is a list of the elements of Fq
matroidAG :: (Fractional a) => Int -> [a] -> Matroid Int
matroidAG n fq = vectorMatroid' $ ptsAG n fq


-- Oxley p182
-- |Given a matroid m, the fundamental-circuit incidence matrix relative to a base b
-- has rows indexed by the elements of b, and columns indexed by the elements not in b.
-- The bi, ej entry is 1 if bi is in the fundamental circuit of ej relative to b, and 0 otherwise.
fundamentalCircuitIncidenceMatrix :: (Ord a, Num k) => Matroid a -> [a] -> [[k]]
fundamentalCircuitIncidenceMatrix m b = L.transpose $ fundamentalCircuitIncidenceMatrix' m b

fundamentalCircuitIncidenceMatrix' m b =
    [ [if e `elem` fundamentalCircuit m b e' then 1 else 0 | e <- b]
    | e' <- elements m LS.\\ b ]

fcim = fundamentalCircuitIncidenceMatrix
fcim' = fundamentalCircuitIncidenceMatrix'

-- Then fcim w4 [1..4] == L.transpose (fcim (dual w4) [5..8])


-- Given a matrix of 0s and 1s, return a matrix of 0s, 1s and *s
-- where 0 -> 0, and 1 -> 1 if it is the first 1 in either a row or column, 1 -> * otherwise
markNonInitialRCs mx = mark (replicate w False) mx
    where w = length $ head mx -- the width
          mark cms (r:rs) = let (cms', r') = mark' False [] cms [] r in r' : mark cms' rs
          mark _ [] = []
          mark' rm cms' (cm:cms) ys (x:xs)
              | x == 0 = mark' rm (cm:cms') cms (Zero:ys) xs
              | x == 1 = if rm && cm
                         then mark' True (True:cms') cms (Star:ys) xs
                         else mark' True (True:cms') cms (One:ys) xs
          mark' _ cms' [] ys [] = (reverse cms', reverse ys)

-- Given a matrix of 0s, 1s and *s, return all distinct matrices that can be obtained
-- by substituting non-zero elements of Fq for the *s
substStars mx fq = substStars' mx
    where fq' = tail fq -- non-zero elts of fq
          substStars' (r:rs) = [r':rs' | r' <- substStars'' r, rs' <- substStars' rs]
          substStars' [] = [[]]
          substStars'' (Zero:xs) = map (0:) $ substStars'' xs
          substStars'' (One:xs) = map (1:) $ substStars'' xs
          substStars'' (Star:xs) = [x':xs' | x' <- fq', xs' <- substStars'' xs]
          substStars'' [] = [[]]

starSubstitutionsV fq' (Zero:xs) = map (0:) $ starSubstitutionsV fq' xs
starSubstitutionsV fq' (One:xs) = map (1:) $ starSubstitutionsV fq' xs
starSubstitutionsV fq' (Star:xs) = [x':xs' | x' <- fq', xs' <- starSubstitutionsV fq' xs]
starSubstitutionsV _ [] = [[]]


-- Oxley p184-5
-- Note that the particular representations you get depend on which basis is used
-- (Perhaps we should let you pass in the basis to use)
representations1 fq m = [ L.transpose d | d <- substStars dhash fq, let mx = ir ++ d,
                                          to1n m == (vectorMatroid' $ map snd $ L.sort $ zip (b ++ b') mx) ]
    where b = head $ bases m
          b' = elements m LS.\\ b
          r = length b -- rank of the matroid
          ir = idMx r -- identity matrix
          dhash = markNonInitialRCs $ fcim' m b

-- edges of the fundamental circuit incidence graph
fcig m b = [ [e,e'] | e <- b, e' <- elements m LS.\\ b, e `elem` fundamentalCircuit m b e' ]

-- the fcim of m relative to b, with 1s and *s marked
markedfcim m b = mark b b' (fcim m b)
    where b' = elements m LS.\\ b -- the elements of b are the row labels, those of b' are the column labels
          entries = fcig m b -- the [row,column] coordinates of the non-zero entries in the fcim
          ones = head $ bases $ cycleMatroid' entries -- a set of entries which we can take to be 1
          stars = entries LS.\\ ones -- the set of entries which we are then still free to assign
          mark (i:is) js (r:rs) = (mark' i js r) : mark is js rs
          mark [] _ [] = []
          mark' i (j:js) (x:xs)
              | x == 0 = Zero : mark' i js xs
              | x == 1 = (if [i,j] `elem` stars then Star else One) : mark' i js xs
          mark' _ [] [] = []
-- The markedfcim will sometimes do better than markNonInitialRCs, and is best possible

-- Oxley p184-5
representations2 fq m = [ L.transpose mx | d <- substStars dhash' fq, let mx = ir ++ d,
                                           m' == (vectorMatroid' $ map snd $ L.sort $ zip (b ++ b') mx) ]
    where m' = to1n m
          es = elements m
          b = head $ bases m
          b' = es LS.\\ b
          r = length b -- rank of the matroid
          ir = idMx r -- identity matrix
          dhash' = L.transpose $ markedfcim m b

-- In the following, we check the representations of the restrictions as we add a column at a time
-- This enables us to early out if a potential representation doesn't work on a restriction

-- |Find representations of the matroid m over fq. Specifically, this function will find one representative
-- of each projective equivalence class of representation.
representations :: (Fractional fq, Ord a) => [fq] -> Matroid a -> [[[fq]]]
representations fq m = map L.transpose $ representations' (reverse $ zip b ir) (zip b' dhash')
    where fq' = tail fq -- fq \ {0}
          b = head $ bases m
          b' = elements m LS.\\ b
          r = length b -- rank of the matroid
          ir = idMx r -- identity matrix
          dhash' = L.transpose $ markedfcim m b
          representations' ls ((i,r):rs) = concat
              [ representations' ((i,r'):ls) rs
              | r' <- starSubstitutionsV fq' r,
                let (is,vs) = unzip $ L.sortBy cmpfst ((i,r'):ls),
                to1n (restriction m is) == (vectorMatroid' vs) ]
          representations' ls [] = [map snd $ reverse ls]

-- |Is the matroid representable over Fq? For example, to find out whether a matroid m is binary, evaluate @isRepresentable f2 m@.
isRepresentable :: (Fractional fq, Ord a) => [fq] -> Matroid a -> Bool
isRepresentable fq m = (not . null) (representations fq m)

-- |A binary matroid is a matroid which is representable over F2
isBinary :: (Ord a) => Matroid a -> Bool
isBinary = isRepresentable f2

-- |A ternary matroid is a matroid which is representable over F3
isTernary :: (Ord a) => Matroid a -> Bool
isTernary = isRepresentable f3


-- CONSTRUCTIONS
-- The next three functions not very thoroughly tested

data LMR a b = L a | Mid | R b deriving (Eq, Ord, Show)

-- Oxley p252
seriesConnection (m1,p1) (m2,p2)
    | not (isLoop m1 p1) && not (isColoop m1 p1) && not (isLoop m2 p2) && not (isColoop m2 p2) =
        fromCircuits es cs
    | otherwise = error "not yet implemented"
        where es = map L (elements m1 LS.\\ [p1]) ++ [Mid] ++ map R (elements m2 LS.\\ [p2])
              cs = (map . map) L (circuits $ m1 \\\ [p1]) ++
                   (map . map) R (circuits $ m2 \\\ [p2]) ++
                   [ map L (L.delete p1 c1) ++ [Mid] ++ map R (L.delete p2 c2)
                   | c1 <- circuits m1, p1 `elem` c1, c2 <- circuits m2, p2 `elem` c2]

parallelConnection (m1,p1) (m2,p2)
    | not (isLoop m1 p1) && not (isColoop m1 p1) && not (isLoop m2 p2) && not (isColoop m2 p2) =
        fromCircuits es cs
    | otherwise = error "not yet implemented"
        where es = map L (elements m1 LS.\\ [p1]) ++ [Mid] ++ map R (elements m2 LS.\\ [p2])
              cs = (map . map) L (circuits $ m1 \\\ [p1]) ++
                   [ map L (L.delete p1 c1) ++ [Mid] | c1 <- circuits m1, p1 `elem` c1 ] ++
                   (map . map) R (circuits $ m2 \\\ [p2]) ++
                   [ [Mid] ++ map R (L.delete p2 c2) | c2 <- circuits m2, p2 `elem` c2 ] ++
                   [ map L (L.delete p1 c1) ++ map R (L.delete p2 c2)
                   | c1 <- circuits m1, p1 `elem` c1, c2 <- circuits m2, p2 `elem` c2 ]

twoSum (m1,p1) (m2,p2)
    | not (isLoop m1 p1) && not (isColoop m1 p1) && not (isLoop m2 p2) && not (isColoop m2 p2) =
        fromCircuits es cs
    | otherwise = error "not yet implemented"
        where es = map L (elements m1 LS.\\ [p1]) ++ [Mid] ++ map R (elements m2 LS.\\ [p2])
              cs = (map . map) L (circuits $ m1 \\\ [p1]) ++
                   (map . map) R (circuits $ m2 \\\ [p2]) ++
                   [ map L (L.delete p1 c1) ++ map R (L.delete p2 c2)
                   | c1 <- circuits m1, p1 `elem` c1, c2 <- circuits m2, p2 `elem` c2]
-- Note: The only different from seriesConnection is that we don't have Mid in the last set

-- Oxley p427
matroidUnion m1 m2 = fromBases es bs
    where es = LS.union (elements m1) (elements m2)
          is = toShortlex $ toSet [ LS.union b1 b2 | b1 <- bases m1, b2 <- bases m2 ]
          r = length $ last is
          bs = dropWhile ( (< r) . length ) is



-- SOME INTERESTING MATROIDS

-- |The Fano plane F7 = PG(2,F2)
f7 :: Matroid Int
f7 = fromGeoRep [] [] [[1,2,3],[1,4,7],[1,5,6],[2,4,6],[2,5,7],[3,4,5],[3,6,7]] [[1..7]]

-- Oxley p36, fig 1.12:
-- cycleMatroid [[1,2],[1,3],[2,3],[3,4],[2,4],[1,4]] == restriction fanoPlane [1..6]

-- |F7-, the relaxation of the Fano plane by removal of a line
f7m :: Matroid Int
f7m = relaxation f7 [2,4,6]

-- Oxley p39
-- |The Pappus configuration from projective geometry
pappus :: Matroid Int
pappus = fromGeoRep [] [] [[1,2,3],[1,5,7],[1,6,8],[2,4,7],[2,6,9],[3,4,8],[3,5,9],[4,5,6],[7,8,9]] [[1..9]]

-- more logical relabeling
-- pappus' = fromGeoRep [] [] [[1,2,3],[1,5,9],[1,6,8],[2,4,9],[2,6,7],[3,4,8],[3,5,7],[4,5,6],[7,8,9]] [[1..9]]

-- |Relaxation of the Pappus configuration by removal of a line
nonPappus :: Matroid Int
nonPappus = relaxation pappus [7,8,9]
-- fromGeoRep [] [] [[1,2,3],[1,5,7],[1,6,8],[2,4,7],[2,6,9],[3,4,8],[3,5,9],[4,5,6]] [[1..9]]

-- |The Desargues configuration
desargues :: Matroid Int
desargues = fromGeoRep [] [] [[1,2,5],[1,3,6],[1,4,7],[2,3,8],[2,4,9],[3,4,10],[5,6,8],[5,7,9],[6,7,10],[8,9,10]]
                              [[1,2,3,5,6,8],[1,2,4,5,7,9],[1,3,4,6,7,10],[2,3,4,8,9,10],[5,6,7,8,9,10]]
-- desargues == cycleMatroid (combinationsOf 2 [1..5])
-- (ie Desargues = M(K5) )
-- interestingly, although these are all the dependent flats that are evident from the diagram,
-- there are also some rank 3 flats consisting of a line and an "antipodal" point
-- eg, in terms of K5, [1,8,9,10] corresponds to [[1,2],[3,4],[3,5],[4,5]]

-- Oxley p71
vamosMatroid1 = fromHyperplanes [1..8] (hs4 ++ hs3)
    where hs4 = [ [1,2,3,4], [1,4,5,6], [2,3,5,6], [1,4,7,8], [2,3,7,8] ]
          hs3 = [ h3 | h3 <- combinationsOf 3 [1..8], null [h4 | h4 <- hs4, h3 `LS.isSubset` h4] ]
vamosMatroid = fromGeoRep [] [] [] [[1,2,3,4],[1,4,5,6],[2,3,5,6],[1,4,7,8],[2,3,7,8]]

-- |The Vamos matroid V8. It is not representable over any field.
v8 :: Matroid Int
v8 = vamosMatroid
-- v8 is self-dual (isomorphic to its own dual)

-- Oxley p188
-- |P8 is a minor-minimal matroid that is not representable over F4, F8, F16, ... .
-- It is Fq-representable if and only if q is not a power of 2.
p8 :: Matroid Int
p8 = vectorMatroid $
    ( [ [1,0,0,0,  0, 1, 1,-1],
        [0,1,0,0,  1, 0, 1, 1],
        [0,0,1,0,  1, 1, 0, 1],
        [0,0,0,1, -1, 1, 1, 0] ] :: [[F3]] )

p8' = fromGeoRep [] [] [] [ [1,2,3,8], [1,2,4,7], [1,3,4,6], [1,4,5,8], [1,5,6,7], [2,3,4,5], [2,3,6,7], [2,5,6,8], [3,5,7,8], [4,6,7,8] ]

-- |P8- is a relaxation of P8. It is Fq-representable if and only if q >= 4.
p8m :: Matroid Int
p8m = relaxation p8 [2,3,6,7]

-- |P8-- is a relaxation of P8-. It is a minor-minimal matroid that is not representable over F4.
-- It is Fq-representable if and only if q >= 5.
p8mm :: Matroid Int
p8mm = relaxation p8m [1,4,5,8]


-- Oxley p317
-- r-spoked wheel graph
wheelGraph r = G.G vs es
    where vs = [0..r]
          es = [ [0,i] | i <- [1..r] ] ++ [ [i,i+1] | i <- [1..r-1] ] ++ [ [1,r] ]
-- for normal form, should L.sort es

mw4 = cycleMatroid $ G.edges $ wheelGraph 4
-- has [5,6,7,8] as unique circuit-hyperplane

w4' = relaxation mw4 $ unique $ circuitHyperplanes mw4 -- [5,6,7,8]

-- Oxley p183
w4 = fromGeoRep [] [] [[1,2,5],[1,4,8],[2,3,6],[3,4,7]] [[1,2,3,5,6],[1,2,4,5,8],[1,3,4,7,8],[2,3,4,6,7]]


-- BINARY MATROIDS

-- Oxley p344
isBinary2 m = all (even . length) [ c `LS.intersect` cc | c <- circuits m, cc <- cocircuits m ]


-- RANK POLYNOMIAL
-- Godsil & Royle p356

[x,y] = map glexVar ["x","y"] :: [GlexPoly Integer String]

-- first naive version
rankPoly1 m = sum [ x^(rm - r a) * y^(rm' - r' a') | a <- powersetdfs es, let a' = es LS.\\ a ]
    where es = elements m
          rm = rank m
          r = rankfun m
          m' = dual m
          rm' = rank m'
          r' = rankfun m'

-- |Given a matroid m over elements es, the rank polynomial is a polynomial r(x,y),
-- which is essentially a generating function for the subsets of es, enumerated by size and rank.
-- It is efficiently calculated using deletion and contraction.
--
-- It has the property that r(0,0) is the number of bases in m, r(1,0) is the number of independent sets,
-- r(0,1) is the number of spanning sets. It can also be used to derive the chromatic polynomial of a graph,
-- the weight enumerator of a linear code, and more.
rankPoly :: (Ord a) => Matroid a -> GlexPoly Integer String
rankPoly m
    | null es = 1
    | isLoop m e = (1+y) * rankPoly (m \\\ [e]) -- deletion
    | isColoop m e = (1+x) * rankPoly (m /// [e]) -- contraction
    | otherwise = rankPoly (m \\\ [e]) + rankPoly (m /// [e])
    where es = elements m
          e = head es

numBases m = unwrap $ rankPoly m `bind` (\v -> case v of "x" -> 0; "y" -> 0)

numIndeps m = unwrap $ rankPoly m `bind` (\v -> case v of "x" -> 1; "y" -> 0)

numSpanning m = unwrap $ rankPoly m `bind` (\v -> case v of "x" -> 0; "y" -> 1)

-- It is then possible to derive the chromatic poly of a graph from the rankPoly of its cycle matroid, etc


-- Oxley p586

-- How many independent sets are there of each size
indepCounts m = map length $ L.groupBy eqfst $ [(length i, i) | i <- indeps m]
-- relying on shortlex order

{-
-- The following tries to calculate the indep counts
-- as the Hilbert series of the Stanley-Reisner ring.
-- However, it's not giving the right answer.
indepCounts2 m = take r $ indepCounts' (circuits m)
    where n = length (elements m) -- the number of "variables"
          r = rank m
          n' = toInteger n
          indepCounts' [] = map (\k -> (k+n'-1) `choose` (n'-1)) [0..]
              -- the number of ways of choosing n-1 separators, so as to leave a product of k powers
          indepCounts' (c:cs) =
              let d = length c -- the "degree" of c
                  cs' = reduce [] $ toShortlex $ toSet $ map (LS.\\ c) cs -- the quotient of cs by c
              in indepCounts' cs <-> (replicate d 0 ++ indepCounts' cs')
          reduce ls (r:rs) = if any (`LS.isSubset` r) ls then reduce ls rs else reduce (r:ls) rs
          reduce ls [] = reverse ls
          (x:xs) <+> (y:ys) = (x+y) : (xs <+> ys)
          xs <+> ys = xs ++ ys -- one of them is null
          xs <-> ys = xs <+> map negate ys 
-}


-- Whitney numbers of the second kind
-- The number of flats of each rank
whitney2nd m = map length $ L.groupBy eqfst $ L.sort [(rankfun m f, f) | f <- flats m]

-- Whitney numbers of the first kind
-- The number of subsets (of the element set) of each rank
whitney1st m = alternatingSign $ map length $ L.groupBy eqfst $ L.sort [(rankfun m x, x) | x <- powersetdfs (elements m)]
    where alternatingSign (x:xs) = x : alternatingSign (map negate xs)
          alternatingSign [] = []



-- TODO

-- 1. Sort out the isomorphism / automorphism code in the case where the incidence graph isn't connected

-- 2. We could generate the geometric representation from a matroid (provided its rank <= 4)
-- geoRep m = filter (isDependent m) (flats m)

-- 3. isMatroidFlats