HaskellForMaths-0.4.2: Math/Combinatorics/Matroid.hs
-- 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 :: (Eq k, 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' :: (Eq k, 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
incidenceGraphC m = G.G vs' es'
where es = elements m
cs = L.sort $ circuits m
vs' = map Left es ++ map Right cs
es' = L.sort [ [Left e, Right c] | c <- cs, e <- c ]
-- incidence graph for the matroid considered as an incidence structure between elements and circuits
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 = incidenceAuts $ incidenceGraphH 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 :: (Eq k, 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 :: (Eq a, 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 :: (Eq a, 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 :: (Eq fq, 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 :: (Eq fq, 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 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