-- 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