diff --git a/HaskellForMaths.cabal b/HaskellForMaths.cabal --- a/HaskellForMaths.cabal +++ b/HaskellForMaths.cabal @@ -1,5 +1,5 @@ Name: HaskellForMaths - Version: 0.3.2 + Version: 0.3.3 Category: Math Description: A library of maths code in the areas of combinatorics, group theory, commutative algebra, and non-commutative algebra. The library is mainly intended for educational purposes, but does have efficient implementations of several fundamental algorithms. Synopsis: Combinatorics, group theory, commutative algebra, non-commutative algebra @@ -26,7 +26,11 @@ Math/Test/TAlgebras/TTensorProduct.hs Math/Test/TAlgebras/TStructures.hs Math/Test/TAlgebras/TQuaternions.hs + Math/Test/TAlgebras/TMatrix.hs Math/Test/TAlgebras/TGroupAlgebra.hs + Math/Test/TCombinatorics/TPoset.hs + Math/Test/TCombinatorics/TDigraph.hs + Math/Test/TCombinatorics/TIncidenceAlgebra.hs Library Build-Depends: base >= 2 && < 5, containers, array, random, QuickCheck @@ -45,7 +49,8 @@ Math.Algebras.TensorProduct, Math.Algebras.VectorSpace, Math.Combinatorics.Graph, Math.Combinatorics.GraphAuts, Math.Combinatorics.StronglyRegularGraph, Math.Combinatorics.Design, Math.Combinatorics.FiniteGeometry, Math.Combinatorics.Hypergraph, - Math.Combinatorics.LatinSquares, + Math.Combinatorics.LatinSquares, Math.Combinatorics.Poset, Math.Combinatorics.IncidenceAlgebra, + Math.Combinatorics.Digraph, Math.Common.IntegerAsType, Math.Common.ListSet, Math.Projects.RootSystem, Math.Projects.Rubik, Math.Projects.MiniquaternionGeometry, diff --git a/Math/Algebra/Field/Base.hs b/Math/Algebra/Field/Base.hs --- a/Math/Algebra/Field/Base.hs +++ b/Math/Algebra/Field/Base.hs @@ -45,6 +45,8 @@ negate (Fp x) = Fp $ p - x where p = value (undefined :: n) Fp x * Fp y = Fp $ (x*y) `mod` p where p = value (undefined :: n) fromInteger m = Fp $ m `mod` p where p = value (undefined :: n) + abs _ = error "Prelude.Num.abs: inappropriate abstraction" + signum _ = error "Prelude.Num.signum: inappropriate abstraction" -- n must be prime - could perhaps use a type to guarantee this instance IntegerAsType n => Fractional (Fp n) where diff --git a/Math/Algebra/Field/Extension.hs b/Math/Algebra/Field/Extension.hs --- a/Math/Algebra/Field/Extension.hs +++ b/Math/Algebra/Field/Extension.hs @@ -40,11 +40,13 @@ | i > 1 = v ++ "^" ++ show i -- "x^" ++ show i instance Num a => Num (UPoly a) where - UP as + UP bs = toUPoly $ as <+> bs - negate (UP as) = UP $ map negate as - UP as * UP bs = toUPoly $ as <*> bs - fromInteger 0 = UP [] - fromInteger a = UP [fromInteger a] + UP as + UP bs = toUPoly $ as <+> bs + negate (UP as) = UP $ map negate as + UP as * UP bs = toUPoly $ as <*> bs + fromInteger 0 = UP [] + fromInteger a = UP [fromInteger a] + abs _ = error "Prelude.Num.abs: inappropriate abstraction" + signum _ = error "Prelude.Num.signum: inappropriate abstraction" toUPoly as = UP (reverse (dropWhile (== 0) (reverse as))) @@ -113,6 +115,8 @@ Ext x * Ext y = Ext $ (x*y) `modUP` pvalue (undefined :: (k,poly)) negate (Ext x) = Ext $ negate x fromInteger x = Ext $ fromInteger x + abs _ = error "Prelude.Num.abs: inappropriate abstraction" + signum _ = error "Prelude.Num.signum: inappropriate abstraction" instance (Num k, Fractional k, PolynomialAsType k poly) => Fractional (ExtensionField k poly) where recip 0 = error "ExtensionField.recip 0" diff --git a/Math/Algebra/NonCommutative/TensorAlgebra.hs b/Math/Algebra/NonCommutative/TensorAlgebra.hs --- a/Math/Algebra/NonCommutative/TensorAlgebra.hs +++ b/Math/Algebra/NonCommutative/TensorAlgebra.hs @@ -1,5 +1,8 @@ --- Copyright (c) David Amos, 2008. All rights reserved. +-- Copyright (c) 2008, David Amos. All rights reserved. +-- |A module defining the tensor, symmetric, and exterior algebras. +-- This module has been partially superceded by Math.Algebras.TensorAlgebra, which should be used in preference. +-- This module is likely to be removed at some point. module Math.Algebra.NonCommutative.TensorAlgebra where import Math.Algebra.Field.Base diff --git a/Math/Algebras/Commutative.hs b/Math/Algebras/Commutative.hs --- a/Math/Algebras/Commutative.hs +++ b/Math/Algebras/Commutative.hs @@ -19,7 +19,7 @@ -- type GlexMonomialS = GlexMonomial String instance Ord v => Ord (GlexMonomial v) where - compare (Glex si xis) (Glex sj yjs) = compare (-si, xis) (-sj, yjs) + compare (Glex si xis) (Glex sj yjs) = compare (-si, [(x,-i) | (x,i) <- xis]) (-sj, [(y,-j) | (y,j) <- yjs]) instance Show v => Show (GlexMonomial v) where show (Glex _ []) = "1" @@ -48,19 +48,18 @@ -- This is the monoid algebra for commutative monomials (which are the free commutative monoid) instance (Num k, Ord v) => Algebra k (GlexMonomial v) where - unit 0 = zero -- V [] - unit x = V [(munit,x)] where munit = Glex 0 [] - mult (V ts) = nf $ fmap (\(a,b) -> a `mmult` b) (V ts) + unit x = x *> return munit + where munit = Glex 0 [] + mult xy = nf $ fmap (\(a,b) -> a `mmult` b) xy where mmult (Glex si xis) (Glex sj yjs) = Glex (si+sj) $ addmerge xis yjs -{- --- This is just the Set Coalgebra, so better to use a generic instance --- Also, not used anywhere. Hence commented out + +-- GlexPoly can be given the set coalgebra structure, which is compatible with the monoid algebra structure instance Num k => Coalgebra k (GlexMonomial v) where - counit (V ts) = sum [x | (m,x) <- ts] -- trace - comult = fmap (\m -> T m m) -- diagonal - -- comult (V ts) = V [(T m m, x) | (m, x) <- ts] --} + counit = unwrap . nf . fmap (\m -> () ) -- trace + -- counit (V ts) = sum [x | (m,x) <- ts] -- trace + comult = fmap (\m -> (m,m) ) -- diagonal + type GlexPoly k v = Vect k (GlexMonomial v) @@ -77,8 +76,8 @@ -- secondly because Haskell doesn't support type functions. bind :: (Monomial m, Num k, Ord b, Show b, Algebra k b) => Vect k (m v) -> (v -> Vect k b) -> Vect k b -V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] --- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts] +V ts `bind` f = sum [c *> product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] +-- flipbind f = linear (\m -> product [f x ^ i | (x,i) <- powers m]) instance Monomial GlexMonomial where diff --git a/Math/Algebras/Matrix.hs b/Math/Algebras/Matrix.hs --- a/Math/Algebras/Matrix.hs +++ b/Math/Algebras/Matrix.hs @@ -13,18 +13,19 @@ -- Mat2 +{- +-- defined in Math.Algebras.TensorProduct delta i j | i == j = 1 | otherwise = 0 +-} data Mat2 = E2 Int Int deriving (Eq,Ord,Show) -- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere instance Num k => Algebra k Mat2 where - -- unit 0 = zero -- V [] - unit x = x `smultL` V [(E2 i i, 1) | i <- [1..2] ] - -- mult ab = nf $ ab >>= mult' where + unit x = x *> V [(E2 i i, 1) | i <- [1..2] ] mult = linear mult' where - mult' (E2 i j, E2 k l) = delta j k `smultL` return (E2 i l) + mult' (E2 i j, E2 k l) = delta j k *> return (E2 i l) -- In other words -- unit x = x (1 0) @@ -51,7 +52,7 @@ data Mat2' = E2' Int Int deriving (Eq,Ord,Show) --- E' i j represents the dual basis element corresponding to E i j +-- E2' i j represents the dual basis element corresponding to E i j -- Kassel p42 instance Num k => Coalgebra k Mat2' where @@ -69,17 +70,15 @@ --- !! Now do the quickchecks - - - data M3 = E3 Int Int deriving (Eq,Ord,Show) -- E i j represents the elementary matrix with a 1 at the (i,j) position, and 0s elsewhere instance Num k => Algebra k M3 where unit 0 = zero -- V [] unit x = V [(E3 i i, x) | i <- [1..3] ] - mult (V ts) = nf $ V $ map (\((E3 i j, E3 k l), x) -> (E3 i l, delta j k * x)) ts + -- mult (V ts) = nf $ V $ map (\((E3 i j, E3 k l), x) -> (E3 i l, delta j k * x)) ts + mult = linear mult' where + mult' (E3 i j, E3 k l) = delta j k *> return (E3 i l) {- -- Kassel p42 diff --git a/Math/Algebras/NonCommutative.hs b/Math/Algebras/NonCommutative.hs --- a/Math/Algebras/NonCommutative.hs +++ b/Math/Algebras/NonCommutative.hs @@ -56,9 +56,10 @@ class Monomial m where var :: v -> Vect Q (m v) powers :: Eq v => m v -> [(v,Int)] +-- why do we need "powers"?? -V ts `bind` f = sum [c `smultL` product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] --- V ts `bind` f = sum [product [f x ^ i | (x,i) <- powers m] * unit c | (m, c) <- ts] +V ts `bind` f = sum [c *> product [f x ^ i | (x,i) <- powers m] | (m, c) <- ts] +-- flipbind f = linear (\m -> product [f x ^ i | (x,i) <- powers m]) instance Monomial NonComMonomial where var v = V [(NCM 1 [v],1)] diff --git a/Math/Algebras/Quaternions.hs b/Math/Algebras/Quaternions.hs --- a/Math/Algebras/Quaternions.hs +++ b/Math/Algebras/Quaternions.hs @@ -24,35 +24,72 @@ show K = "k" instance (Num k) => Algebra k HBasis where - unit 0 = zero -- V [] - unit x = V [(One,x)] - -- mult x = nf (x >>= m) - mult = linear m - where m (One,b) = return b - m (b,One) = return b - m (I,I) = unit (-1) - m (J,J) = unit (-1) - m (K,K) = unit (-1) - m (I,J) = return K - m (J,I) = -1 *> return K - m (J,K) = return I - m (K,J) = -1 *> return I - m (K,I) = return J - m (I,K) = -1 *> return J + unit x = x *> return One + mult = linear mult' + where mult' (One,b) = return b + mult' (b,One) = return b + mult' (I,I) = unit (-1) + mult' (J,J) = unit (-1) + mult' (K,K) = unit (-1) + mult' (I,J) = return K + mult' (J,I) = -1 *> return K + mult' (J,K) = return I + mult' (K,J) = -1 *> return I + mult' (K,I) = return J + mult' (I,K) = -1 *> return J i,j,k :: Num k => Quaternion k i = return I j = return J k = return K + +one',i',j',k' :: Num k => Vect k (Dual HBasis) +one' = return (Dual One) +i' = return (Dual I) +j' = return (Dual J) +k' = return (Dual K) + +-- Coalgebra structure on the dual vector space to the quaternions +-- The comult is the transpose of mult +instance Num k => Coalgebra k (Dual HBasis) where + counit = unwrap . linear counit' + where counit' (Dual One) = return () + counit' _ = zero + comult = linear comult' + where comult' (Dual One) = return (Dual One, Dual One) <+> + (-1) *> ( return (Dual I, Dual I) <+> return (Dual J, Dual J) <+> return (Dual K, Dual K) ) + comult' (Dual I) = return (Dual One, Dual I) <+> return (Dual I, Dual One) <+> + return (Dual J, Dual K) <+> (-1) *> return (Dual K, Dual J) + comult' (Dual J) = return (Dual One, Dual J) <+> return (Dual J, Dual One) <+> + return (Dual K, Dual I) <+> (-1) *> return (Dual I, Dual K) + comult' (Dual K) = return (Dual One, Dual K) <+> return (Dual K, Dual One) <+> + return (Dual I, Dual J) <+> (-1) *> return (Dual J, Dual I) + {- +-- Of course, we can define this coalgebra structure on the quaternions themselves +-- However, it is not compatible with the algebra structure: we don't get a bialgebra +instance Num k => Coalgebra k HBasis where + counit = unwrap . linear counit' + where counit' One = return () + counit' _ = zero + comult = linear comult' + where comult' One = return (One,One) <+> (-1) *> ( return (I,I) <+> return (J,J) <+> return (K,K) ) + comult' I = return (One,I) <+> return (I,One) <+> return (J,K) <+> (-1) *> return (K,J) + comult' J = return (One,J) <+> return (J,One) <+> return (K,I) <+> (-1) *> return (I,K) + comult' K = return (One,K) <+> return (K,One) <+> return (I,J) <+> (-1) *> return (J,I) +-} + +{- -- Set coalgebra instance instance Num k => Coalgebra k HBasis where counit (V ts) = sum [x | (m,x) <- ts] -- trace comult = fmap (\m -> T m m) -- diagonal -} +{- instance Num k => Coalgebra k HBasis where counit (V ts) = sum [x | (One,x) <- ts] comult = linear cm where cm m = if m == One then return (m,m) else return (m,One) <+> return (One,m) +-} diff --git a/Math/Algebras/Structures.hs b/Math/Algebras/Structures.hs --- a/Math/Algebras/Structures.hs +++ b/Math/Algebras/Structures.hs @@ -35,7 +35,8 @@ comult :: Vect k b -> Vect k (Tensor b b) --- |A bialgebra is an algebra which is also a coalgebra, subject to some compatibility conditions +-- |A bialgebra is an algebra which is also a coalgebra, subject to the compatibility conditions +-- that counit and comult must be algebra morphisms (or equivalently, that unit and mult must be coalgebra morphisms) class (Algebra k b, Coalgebra k b) => Bialgebra k b where {} class Bialgebra k b => HopfAlgebra k b where @@ -43,8 +44,9 @@ instance (Num k, Eq b, Ord b, Show b, Algebra k b) => Num (Vect k b) where - x+y = add x y - negate (V ts) = V $ map (\(b,x) -> (b, negate x)) ts + x+y = x <+> y + negate x = neg x + -- negate (V ts) = V $ map (\(b,x) -> (b, negate x)) ts x*y = mult (x `te` y) fromInteger n = unit (fromInteger n) abs _ = error "Prelude.Num.abs: inappropriate abstraction" @@ -65,37 +67,66 @@ instance Num k => Algebra k () where - unit 0 = zero -- V [] - unit x = V [( (),x)] + unit = wrap + -- unit 0 = zero -- V [] + -- unit x = V [( (),x)] mult (V [( ((),()), x)]) = V [( (),x)] instance Num k => Coalgebra k () where - counit (V []) = 0 - counit (V [( (),x)]) = x + counit = unwrap + -- counit (V []) = 0 + -- counit (V [( (),x)]) = x comult (V [( (),x)]) = V [( ((),()), x)] --- |Trivial k is the field k considered as a k-vector space. In maths, we would not normally make a distinction here, --- but in the code, we need this if we want to be able to put k as one side of a tensor product. -type Trivial k = Vect k () - unit' :: (Num k, Algebra k b) => Trivial k -> Vect k b -unit' = unit . unwrap where unwrap = counit :: Num k => Trivial k -> k +unit' = unit . unwrap -- where unwrap = counit :: Num k => Trivial k -> k counit' :: (Num k, Coalgebra k b) => Vect k b -> Trivial k -counit' = wrap . counit where wrap = unit :: Num k => k -> Trivial k +counit' = wrap . counit -- where wrap = unit :: Num k => k -> Trivial k -- unit' and counit' enable us to form tensors of these functions +-- Kassel p4 +-- |The direct sum of k-algebras can itself be given the structure of a k-algebra. +-- This is the product object in the category of k-algebras. +instance (Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (DSum a b) where + unit k = i1 (unit k) <+> i2 (unit k) + -- unit == (i1 . unit) <<+>> (i2 . unit) + mult = linear mult' + where mult' (Left a1, Left a2) = i1 $ mult $ return (a1,a2) + mult' (Right b1, Right b2) = i2 $ mult $ return (b1,b2) + mult' _ = zero +-- This is the product algebra, which is the product in the category of algebras +-- 1 = (1,1) +-- (a1,b1) * (a2,b2) = (a1*a2, b1*b2) +-- It's not a coproduct, because i1, i2 aren't algebra morphisms (they violate Unit axiom) + +-- |The direct sum of k-coalgebras can itself be given the structure of a k-coalgebra. +-- This is the coproduct object in the category of k-coalgebras. +instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (DSum a b) where + counit = unwrap . linear counit' + where counit' (Left a) = (wrap . counit) (return a) + counit' (Right b) = (wrap . counit) (return b) + -- counit = counit . p1 <<+>> counit . p2 + comult = linear comult' where + comult' (Left a) = fmap (\(a1,a2) -> (Left a1, Left a2)) $ comult $ return a + comult' (Right b) = fmap (\(b1,b2) -> (Right b1, Right b2)) $ comult $ return b + -- comult = ( (i1 `tf` i1) . comult . p1 ) <<+>> ( (i2 `tf` i2) . comult . p2 ) + + + + -- Kassel p32 +-- |The tensor product of k-algebras can itself be given the structure of a k-algebra instance (Num k, Ord a, Ord b, Algebra k a, Algebra k b) => Algebra k (Tensor a b) where - unit 0 = V [] - unit x = x `smultL` (unit 1 `te` unit 1) - -- mult x = nf $ x >>= m where + -- unit 0 = V [] + unit x = x *> (unit 1 `te` unit 1) mult = linear m where m ((a,b),(a',b')) = (mult $ return (a,a')) `te` (mult $ return (b,b')) -- Kassel p42 +-- |The tensor product of k-coalgebras can itself be given the structure of a k-coalgebra instance (Num k, Ord a, Ord b, Coalgebra k a, Coalgebra k b) => Coalgebra k (Tensor a b) where counit = counit . (counit' `tf` counit') -- counit = counit . linear (\(T x y) -> counit' (return x) * counit' (return y)) @@ -103,11 +134,16 @@ . (id `tf` assocL) . assocR . (comult `tf` comult) +-- The set coalgebra - can be defined on any set +instance Num k => Coalgebra k EBasis where + counit (V ts) = sum [x | (ei,x) <- ts] -- trace + comult = fmap ( \ei -> (ei,ei) ) -- diagonal + newtype SetCoalgebra b = SC b deriving (Eq,Ord,Show) instance Num k => Coalgebra k (SetCoalgebra b) where - counit (V ts) = sum [x | (m,x) <- ts] -- trace - comult = fmap ( \m -> (m,m) ) -- diagonal + counit (V ts) = sum [x | (m,x) <- ts] -- trace + comult = fmap ( \m -> (m,m) ) -- diagonal newtype MonoidCoalgebra m = MC m deriving (Eq,Ord,Show) diff --git a/Math/Algebras/TensorAlgebra.hs b/Math/Algebras/TensorAlgebra.hs --- a/Math/Algebras/TensorAlgebra.hs +++ b/Math/Algebras/TensorAlgebra.hs @@ -1,8 +1,8 @@ --- Copyright (c) 2010, David Amos. All rights reserved. +-- Copyright (c) 2010-2011, David Amos. All rights reserved. -{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-} --- |A module defining the tensor algebra, symmetric algebra, and exterior (or alternating) algebra +-- |A module defining the tensor algebra, symmetric algebra, exterior (or alternating) algebra, and tensor coalgebra module Math.Algebras.TensorAlgebra where import qualified Data.List as L @@ -14,35 +14,185 @@ import Math.Algebra.Field.Base -data TensorAlgebra a = TA Int [a] deriving (Eq,Ord,Show) +-- TENSOR ALGEBRA +-- |A data type representing basis elements of the tensor algebra over a set\/type. +-- Elements of the tensor algebra are linear combinations of iterated tensor products of elements of the set\/type. +-- If V = Vect k a is the free vector space over a, then the tensor algebra T(V) = Vect k (TensorAlgebra a) is isomorphic +-- to the infinite direct sum: +-- +-- T(V) = k ⊕ V ⊕ V⊗V ⊕ V⊗V⊗V ⊕ ... +-- where 1 is the unit vector space k =~ Vect k () +data TensorAlgebra a = TA Int [a] deriving (Eq,Ord) + +instance Show a => Show (TensorAlgebra a) where + show (TA _ []) = "1" + show (TA _ xs) = filter (/= '"') $ concat $ L.intersperse "*" $ map show xs + -- show (TA _ xs) = filter (/= '"') $ concat $ L.intersperse "\x2297" $ map show xs + + instance Mon (TensorAlgebra a) where munit = TA 0 [] mmult (TA i xs) (TA j ys) = TA (i+j) (xs++ys) instance (Num k, Ord a) => Algebra k (TensorAlgebra a) where - unit 0 = zero -- V [] - unit x = V [(munit,x)] + unit x = x *> return munit mult = nf . fmap (\(a,b) -> a `mmult` b) +-- The tensor algebra is the free algebra. It has the following universal property: +-- Given f :: a -> Vect k b, where Vect k b is an algebra +-- (which induces a vector space morphism, linear f :: Vect k a -> Vect k b) +-- then we can lift to an algebra morphism, (liftTA f) :: Vect k (TensorAlgebra a) -> Vect k b +-- with (liftTA f) . linear injectTA = linear f -data SymmetricAlgebra a = Sym Int [a] deriving (Eq,Ord,Show) +-- |Inject an element of the free vector space V = Vect k a into the tensor algebra T(V) = Vect k (TensorAlgebra a) +injectTA :: Num k => Vect k a -> Vect k (TensorAlgebra a) +injectTA = fmap (\a -> TA 1 [a]) +-- The Num k context is not strictly necessary +-- |Inject an element of the set\/type A\/a into the tensor algebra T(A) = Vect k (TensorAlgebra a). +injectTA' :: Num k => a -> Vect k (TensorAlgebra a) +injectTA' = injectTA . return +-- injectTA' a = return (TA 1 [a]) + +-- |Given vector spaces A = Vect k a, B = Vect k b, where B is also an algebra, +-- lift a linear map f: A -> B to an algebra morphism f': T(A) -> B, +-- where T(A) is the tensor algebra Vect k (TensorAlgebra a). +-- f' will agree with f on A itself (considered as a subspace of T(A)). +-- In other words, f = f' . injectTA +liftTA :: (Num k, Ord b, Show b, Algebra k b) => + (Vect k a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k b +liftTA f = linear (\(TA _ xs) -> product [f (return x) | x <- xs]) +-- The Show b constraint is required because we use product (and Num requires Show)!! + +-- |Given a set\/type A\/a, and a vector space B = Vect k b, where B is also an algebra, +-- lift a function f: A -> B to an algebra morphism f': T(A) -> B. +-- f' will agree with f on A itself. In other words, f = f' . injectTA' +liftTA' :: (Num k, Ord b, Show b, Algebra k b) => + (a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k b +liftTA' = liftTA . linear +-- liftTA' f = linear (\(TA _ xs) -> product [f x | x <- xs]) +-- The second version might be more efficient + + +-- |Tensor algebra is a functor from k-Vect to k-Alg. +-- The action on objects is Vect k a -> Vect k (TensorAlgebra a). +-- The action on arrows is f -> fmapTA f. +fmapTA :: (Num k, Ord b, Show b) => + (Vect k a -> Vect k b) -> Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b) +fmapTA f = liftTA (injectTA . f) +-- fmapTA f = linear (\(TA _ xs) -> product [injectTA (f (return x)) | x <- xs]) + +-- |If we compose the free vector space functor Set -> k-Vect with the tensor algebra functor k-Vect -> k-Alg, +-- we obtain a functor Set -> k-Alg, the free algebra functor. +-- The action on objects is a -> Vect k (TensorAlgebra a). +-- The action on arrows is f -> fmapTA' f. +fmapTA' :: (Num k, Ord b, Show b) => + (a -> b) -> Vect k (TensorAlgebra a) -> Vect k (TensorAlgebra b) +fmapTA' = fmapTA . fmap +-- fmapTA' f = liftTA' (injectTA' . f) +-- fmapTA' f = linear (\(TA _ xs) -> product [injectTA' (f x) | x <- xs]) + + +bindTA :: (Num k, Ord b, Show b) => + Vect k (TensorAlgebra a) -> (Vect k a -> Vect k (TensorAlgebra b)) -> Vect k (TensorAlgebra b) +bindTA = flip liftTA + +bindTA' :: (Num k, Ord b, Show b) => + Vect k (TensorAlgebra a) -> (a -> Vect k (TensorAlgebra b)) -> Vect k (TensorAlgebra b) +bindTA' = flip liftTA' +-- Another way to think about this is variable substitution + +-- "The algebra is free until we bind it" + + +-- SYMMETRIC ALGEBRA + +-- |A data type representing basis elements of the symmetric algebra over a set\/type. +-- The symmetric algebra is the quotient of the tensor algebra by +-- the ideal generated by all +-- differences of products u⊗v - v⊗u. +data SymmetricAlgebra a = Sym Int [a] deriving (Eq,Ord) + +instance Show a => Show (SymmetricAlgebra a) where + show (Sym _ []) = "1" + show (Sym _ xs) = filter (/= '"') $ concat $ L.intersperse "." $ map show xs + instance Ord a => Mon (SymmetricAlgebra a) where munit = Sym 0 [] mmult (Sym i xs) (Sym j ys) = Sym (i+j) $ L.sort (xs++ys) instance (Num k, Ord a) => Algebra k (SymmetricAlgebra a) where - unit 0 = zero -- V [] - unit x = V [(munit,x)] + unit x = x *> return munit mult = nf . fmap (\(a,b) -> a `mmult` b) +-- |Algebra morphism from tensor algebra to symmetric algebra. +-- The kernel of the morphism is the ideal generated by all +-- differences of products u⊗v - v⊗u. +toSym :: (Num k, Ord a) => + Vect k (TensorAlgebra a) -> Vect k (SymmetricAlgebra a) +toSym = linear toSym' + where toSym' (TA i xs) = return $ Sym i (L.sort xs) -data ExteriorAlgebra a = Ext Int [a] deriving (Eq,Ord,Show) +-- The symmetric algebra is the free commutative algebra. It has the following universal property: +-- Given f :: a -> Vect k b, where Vect k b is a commutative algebra +-- (which induces a vector space morphism, linear f :: Vect k a -> Vect k b) +-- then we can lift to a commutative algebra morphism, (liftSym f) :: Vect k (SymmetricAlgebra a) -> Vect k b +-- with (liftSym f) . injectSym = f + +injectSym :: Num k => Vect k a -> Vect k (SymmetricAlgebra a) +injectSym = fmap (\a -> Sym 1 [a]) + +injectSym' :: Num k => a -> Vect k (SymmetricAlgebra a) +injectSym' = injectSym . return +-- injectSym' a = return (Sym 1 [a]) + +liftSym :: (Num k, Ord b, Show b, Algebra k b) => + (Vect k a -> Vect k b) -> Vect k (SymmetricAlgebra a) -> Vect k b +liftSym f = linear (\(Sym _ xs) -> product [f (return x) | x <- xs]) + +liftSym' :: (Num k, Ord b, Show b, Algebra k b) => + (a -> Vect k b) -> Vect k (SymmetricAlgebra a) -> Vect k b +liftSym' = liftSym . linear +-- liftSym' f = linear (\(Sym _ xs) -> product [f x | x <- xs]) + +fmapSym :: (Num k, Ord b, Show b) => + (Vect k a -> Vect k b) -> Vect k (SymmetricAlgebra a) -> Vect k (SymmetricAlgebra b) +fmapSym f = liftSym (injectSym . f) +-- fmapSym f = linear (\(Sym _ xs) -> product [injectSym (f (return x)) | x <- xs]) + +fmapSym' :: (Num k, Ord b, Show b) => + (a -> b) -> Vect k (SymmetricAlgebra a) -> Vect k (SymmetricAlgebra b) +fmapSym' = fmapSym . fmap +-- fmapSym' f = liftSym' (injectSym' . f) +-- fmapSym' f = linear (\(Sym _ xs) -> product [injectSym' (f x) | x <- xs]) + +bindSym :: (Num k, Ord b, Show b) => + Vect k (SymmetricAlgebra a) -> (Vect k a -> Vect k (SymmetricAlgebra b)) -> Vect k (SymmetricAlgebra b) +bindSym = flip liftSym + +bindSym' :: (Num k, Ord b, Show b) => + Vect k (SymmetricAlgebra a) -> (a -> Vect k (SymmetricAlgebra b)) -> Vect k (SymmetricAlgebra b) +bindSym' = flip liftSym' +-- Another way to think about this is variable substitution + + +-- EXTERIOR ALGEBRA + +-- |A data type representing basis elements of the exterior algebra over a set\/type. +-- The exterior algebra is the quotient of the tensor algebra by +-- the ideal generated by all +-- self-products u⊗u and sums of products u⊗v + v⊗u +data ExteriorAlgebra a = Ext Int [a] deriving (Eq,Ord) + +instance Show a => Show (ExteriorAlgebra a) where + show (Ext _ []) = "1" + show (Ext _ xs) = filter (/= '"') $ concat $ L.intersperse "^" $ map show xs + + instance (Num k, Ord a) => Algebra k (ExteriorAlgebra a) where - unit 0 = zero -- V [] - unit x = V [(Ext 0 [],x)] + unit x = x *> return (Ext 0 []) mult xy = nf $ xy >>= (\(Ext i xs, Ext j ys) -> signedMerge 1 (0,[]) (i,xs) (j,ys)) where signedMerge s (k,zs) (i,x:xs) (j,y:ys) = case compare x y of @@ -52,3 +202,140 @@ in signedMerge s' (k+1,y:zs) (i,x:xs) (j-1,ys) signedMerge s (k,zs) (i,xs) (0,[]) = s *> (return $ Ext (k+i) $ reverse zs ++ xs) signedMerge s (k,zs) (0,[]) (j,ys) = s *> (return $ Ext (k+j) $ reverse zs ++ ys) + + +-- |Algebra morphism from tensor algebra to exterior algebra. +-- The kernel of the morphism is the ideal generated by all +-- self-products u⊗u and sums of products u⊗v + v⊗u +toExt :: (Num k, Ord a) => + Vect k (TensorAlgebra a) -> Vect k (ExteriorAlgebra a) +toExt = linear toExt' + where toExt' (TA i xs) = let (sign,xs') = signedSort 1 True [] xs + in fromInteger sign *> return (Ext i xs') + +signedSort sign done ls (r1:r2:rs) = + case compare r1 r2 of + EQ -> (0,[]) + LT -> signedSort sign done (r1:ls) (r2:rs) + GT -> signedSort (-sign) False (r2:ls) (r1:rs) +signedSort sign done ls rs = + if done then (sign,reverse ls ++ rs) else signedSort sign True [] (reverse ls ++ rs) + +-- !! The above code seems a bit clumsy - can we do better + + +injectExt :: Num k => Vect k a -> Vect k (ExteriorAlgebra a) +injectExt = fmap (\a -> Ext 1 [a]) + +injectExt' :: Num k => a -> Vect k (ExteriorAlgebra a) +injectExt' = injectExt . return +-- injectExt' a = return (Ext 1 [a]) + +liftExt :: (Num k, Ord b, Show b, Algebra k b) => + (Vect k a -> Vect k b) -> Vect k (ExteriorAlgebra a) -> Vect k b +liftExt f = linear (\(Ext _ xs) -> product [f (return x) | x <- xs]) + +liftExt' :: (Num k, Ord b, Show b, Algebra k b) => + (a -> Vect k b) -> Vect k (ExteriorAlgebra a) -> Vect k b +liftExt' = liftExt . linear +-- liftExt' f = linear (\(Ext _ xs) -> product [f x | x <- xs]) + +fmapExt :: (Num k, Ord b, Show b) => + (Vect k a -> Vect k b) -> Vect k (ExteriorAlgebra a) -> Vect k (ExteriorAlgebra b) +fmapExt f = liftExt (injectExt . f) +-- fmapExt f = linear (\(Ext _ xs) -> product [injectExt (f (return x)) | x <- xs]) + +fmapExt' :: (Num k, Ord b, Show b) => + (a -> b) -> Vect k (ExteriorAlgebra a) -> Vect k (ExteriorAlgebra b) +fmapExt' = fmapExt . fmap +-- fmapExt' f = liftExt' (injectExt' . f) +-- fmapExt' f = linear (\(Ext _ xs) -> product [injectExt' (f x) | x <- xs]) + +bindExt :: (Num k, Ord b, Show b) => + Vect k (ExteriorAlgebra a) -> (Vect k a -> Vect k (ExteriorAlgebra b)) -> Vect k (ExteriorAlgebra b) +bindExt = flip liftExt + +bindExt' :: (Num k, Ord b, Show b) => + Vect k (ExteriorAlgebra a) -> (a -> Vect k (ExteriorAlgebra b)) -> Vect k (ExteriorAlgebra b) +bindExt' = flip liftExt' +-- Another way to think about this is variable substitution + + +-- TENSOR COALGEBRA + +-- Kassel p67 +data TensorCoalgebra c = TC Int [c] deriving (Eq,Ord,Show) + +instance (Num k, Ord c) => Coalgebra k (TensorCoalgebra c) where + counit = unwrap . linear counit' + where counit' (TC 0 []) = return () -- 1 + counit' _ = zerov + comult = linear comult' + where comult' (TC d xs) = sumv [return (TC i ls, TC (d-i) rs) | (i,ls,rs) <- L.zip3 [0..] (L.inits xs) (L.tails xs)] + + +-- Now show that the tensor coalgebra is the cofree coalgebra +-- ie that it has the required universal property: +-- coliftTC f is a coalgebra morphism, and f == projectTC . coliftTC f + +-- projection onto the underlying vector space +projectTC :: (Num k, Ord b) => Vect k (TensorCoalgebra b) -> Vect k b +projectTC = linear projectTC' where projectTC' (TC 1 [b]) = return b; projectTC' _ = zerov +-- projectTC t = V [(b,c) | (TC 1 [b], c) <- terms t] + + +-- lift a vector space morphism C -> D to a coalgebra morphism C -> T'(D) +-- this function returns an approximation, valid only up to second order terms +coliftTC :: (Num k, Coalgebra k c, Ord d) => + (Vect k c -> Vect k d) -> Vect k c -> Vect k (TensorCoalgebra d) +coliftTC f = sumf [coliftTC' i f | i <- [0..2] ] + +coliftTC' 0 f = linear f0' + where f0' c = counit (return c) *> return (TC 0 []) +coliftTC' 1 f = linear f1' + where f1' c = fmap (\d -> TC 1 [d]) (f $ return c) +coliftTC' n f = linear fn' + where f1' = coliftTC' 1 f + fn1' = coliftTC' (n-1) f + fn' c = fmap (\(TC 1 [x], TC _ xs) -> TC n (x:xs)) $ ( (f1' `tf` fn1') . comult) (return c) + + +cobindTC :: (Num k, Ord c, Ord d) => + (Vect k (TensorCoalgebra c) -> Vect k d) -> Vect k (TensorCoalgebra c) -> Vect k (TensorCoalgebra d) +cobindTC = coliftTC + +-- So we have a comonad: +-- projectTC is extract :: w a -> a +-- cobindTC is extend :: (w a -> b) -> w a -> w b + +{- +Derivation of coliftTC: +Write f' = f0' + f1' + f2' + ..., +where fn' is the part of f' whose range is the nth iterated tensor product in TC. +Then we can deduce f0' from counit . f' == counit +If f': c -> sum ai*di + terms of other order +then counit c = sum ai*counit di +We can deduce f1' from f == projectTC . f' +We can deduce the rest recursively from comult +Write comult (on TC) = comult00 + (comult01+comult10) + (comult02+comult11+comult20) + ..., +where comultij is that part that operates on the i+j'th tensor product to produce i'th `te` jth +Then comult . f' = (f' `tf` f') . comult +can be expanded as +(comult00 + comult01+comult10 + ...) . (f0' + f1' + ...) = (f0' `tf` f0' + f0' `tf` f1' + f1' `tf` f0' + ...) . comult +Looking at the 1,n-1 term, we see that +comult1,n-1 . fn' = (f1' `tf` fn-1') . comult +-} + +-- For example +{- +> let f = linear (\x -> case x of Dual One -> e1; Dual I -> e2; Dual J -> e3; Dual K -> e 4) +> let f' = sumf [coliftTC' i f | i <- [0..3] ] + +-- then the following agree up to level three (inclusive) +> (comult . f') one' +> ((f' `tf` f') . comult) one' +-} + + + + diff --git a/Math/Algebras/TensorProduct.hs b/Math/Algebras/TensorProduct.hs --- a/Math/Algebras/TensorProduct.hs +++ b/Math/Algebras/TensorProduct.hs @@ -7,6 +7,10 @@ import Math.Algebras.VectorSpace +infix 7 `te`, `tf` +infix 6 `dsume`, `dsumf` + + -- DIRECT SUM -- |A type for constructing a basis for the direct sum of vector spaces. @@ -85,20 +89,25 @@ -- tensor isomorphisms -- in fact, this definition works for any Functor f, not just (Vect k) -assocL :: Vect k (Tensor u (Tensor v w)) -> Vect k (Tensor (Tensor u v) w) +assocL :: Vect k (Tensor a (Tensor b c)) -> Vect k (Tensor (Tensor a b) c) assocL = fmap ( \(a,(b,c)) -> ((a,b),c) ) -assocR :: Vect k (Tensor (Tensor u v) w) -> Vect k (Tensor u (Tensor v w)) +assocR :: Vect k (Tensor (Tensor a b) c) -> Vect k (Tensor a (Tensor b c)) assocR = fmap ( \((a,b),c) -> (a,(b,c)) ) +unitInL :: Vect k a -> Vect k (Tensor () a) unitInL = fmap ( \a -> ((),a) ) +unitOutL :: Vect k (Tensor () a) -> Vect k a unitOutL = fmap ( \((),a) -> a ) +unitInR :: Vect k a -> Vect k (Tensor a ()) unitInR = fmap ( \a -> (a,()) ) +unitOutR :: Vect k (Tensor a ()) -> Vect k a unitOutR = fmap ( \(a,()) -> a ) +twist :: (Num k, Ord a, Ord b) => Vect k (Tensor a b) -> Vect k (Tensor b a) twist v = nf $ fmap ( \(a,b) -> (b,a) ) v -- note the nf call, as f is not order-preserving @@ -121,3 +130,13 @@ -- For example: -- > distrL (e1 `te` i1 e2) :: Vect Q (DSum (Tensor EBasis EBasis) (Tensor EBasis EBasis)) -- Left (e1,e2) + + +ev :: (Num k, Ord b) => Vect k (Tensor (Dual b) b) -> k +ev = unwrap . linear (\(Dual bi, bj) -> delta bi bj *> return ()) +-- slightly cheating, as delta i j is meant to compare indices, not the basis elements themselves + +delta i j = if i == j then 1 else 0 + +reify :: (Num k, Ord b) => Vect k (Dual b) -> (Vect k b -> k) +reify f x = ev (f `te` x) diff --git a/Math/Algebras/VectorSpace.hs b/Math/Algebras/VectorSpace.hs --- a/Math/Algebras/VectorSpace.hs +++ b/Math/Algebras/VectorSpace.hs @@ -12,7 +12,7 @@ infixr 7 *> infixl 7 <* -infixl 6 <+> +infixl 6 <+>, <-> -- |Given a field type k (ie a Fractional instance), Vect k b is the type of the free k-vector space over the basis type b. @@ -25,10 +25,10 @@ where showTerm (b,x) | show b == "1" = show x | show x == "1" = show b | show x == "-1" = "-" ++ show b - | otherwise = (if isAtomic (show x) then show x else "(" ++ show x ++ ")") - -- (let (c:cs) = show x in - -- if any (`elem` "+-") cs then "(" ++ show x ++ ")" else show x) - ++ show b + | otherwise = (if isAtomic (show x) then show x else "(" ++ show x ++ ")") ++ + show b + -- (if ' ' `notElem` show b then show b else "(" ++ show b ++ ")") + -- if we put this here we miss the two cases above concatWithPlus (t1:t2:ts) = if head t2 == '-' then t1 ++ concatWithPlus (t2:ts) else t1 ++ '+' : concatWithPlus (t2:ts) @@ -40,10 +40,17 @@ isAtomic' (c:cs) = isAtomic' cs isAtomic' [] = True --- |The zero vector -zero :: Vect k b +terms (V ts) = ts + +coeff b v = sum [k | (b',k) <- terms v, b' == b] + +-- Deprecated zero = V [] +-- |The zero vector +zerov :: Vect k b +zerov = V [] + -- |Addition of vectors add :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b add (V ts) (V us) = V $ addmerge ts us @@ -60,10 +67,18 @@ addmerge ts [] = ts addmerge [] us = us +-- |Sum of a list of vectors +sumv :: (Ord b, Num k) => [Vect k b] -> Vect k b +sumv = foldl (<+>) zerov + -- |Negation of vector neg :: (Num k) => Vect k b -> Vect k b neg (V ts) = V $ map (\(b,x) -> (b,-x)) ts +-- |Subtraction of vectors +(<->) :: (Ord b, Num k) => Vect k b -> Vect k b -> Vect k b +(<->) u v = u <+> neg v + -- |Scalar multiplication (on the left) smultL :: (Num k) => k -> Vect k b -> Vect k b smultL 0 _ = zero -- V [] @@ -122,9 +137,44 @@ instance Show EBasis where show (E i) = "e" ++ show i -e i = return (E i) +e i = return $ E i e1 = e 1 e2 = e 2 e3 = e 3 -- dual (E i) = E (-i) + + +-- |Trivial k is the field k considered as a k-vector space. In maths, we would not normally make a distinction here, +-- but in the code, we need this if we want to be able to put k as one side of a tensor product. +type Trivial k = Vect k () + +wrap :: Num k => k -> Vect k () +wrap 0 = zero +wrap x = V [( (),x)] + +unwrap :: Num k => Vect k () -> k +unwrap (V []) = 0 +unwrap (V [( (),x)]) = x + +-- |Given a finite vector space basis b, Dual b represents a basis for the dual vector space. (If b is infinite, then Dual b is only a sub-basis.) +newtype Dual b = Dual b deriving (Eq,Ord) + +instance Show basis => Show (Dual basis) where + show (Dual b) = show b ++ "'" + + +e' i = return $ Dual $ E i +e1' = e' 1 +e2' = e' 2 +e3' = e' 3 + +dual :: Vect k b -> Vect k (Dual b) +dual = fmap Dual + + +(f <<+>> g) v = f v <+> g v + +zerof v = zerov + +sumf fs = foldl (<<+>>) zerof fs diff --git a/Math/Combinatorics/Digraph.hs b/Math/Combinatorics/Digraph.hs new file mode 100644 --- /dev/null +++ b/Math/Combinatorics/Digraph.hs @@ -0,0 +1,288 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +{-# LANGUAGE NoMonomorphismRestriction #-} + +-- |A module for working with directed graphs (digraphs). +-- Some of the functions are specifically for working with directed acyclic graphs (DAGs), +-- that is, directed graphs containing no cycles. +module Math.Combinatorics.Digraph where + +import Data.List as L +import qualified Data.Map as M +import qualified Data.Set as S + +toSet = S.toList . S.fromList + +-- |A digraph is represented as DG vs es, where vs is the list of vertices, and es is the list of edges. +-- Edges are directed: an edge (u,v) means an edge from u to v. +-- A digraph is considered to be in normal form if both es and vs are in ascending order. +-- This is the preferred form, and some functions will only work for digraphs in normal form. +data Digraph v = DG [v] [(v,v)] deriving (Eq,Ord,Show) + +instance Functor Digraph where + -- |If f is not order-preserving, then you should call nf afterwards + fmap f (DG vs es) = DG (map f vs) (map (\(u,v)->(f u, f v)) es) + +nf (DG vs es) = DG (L.sort vs) (L.sort es) + +vertices (DG vs _) = vs + +edges (DG _ es) = es + + +-- Is it valid to call them predecessors / successors in the case when the digraph contains cycles? + +predecessors (DG _ es) v = [u | (u,v') <- es, v' == v] + +successors (DG _ es) u = [v | (u',v) <- es, u' == u] + +-- Calculate maps of predecessor and successor lists for each vertex in a digraph. +-- If a vertex has no predecessors (respectively successors), then it is left out of the relevant map +adjLists (DG vs es) = adjLists' (M.empty, M.empty) es + where adjLists' (preds,succs) ((u,v):es) = + adjLists' (M.insertWith' (flip (++)) v [u] preds, M.insertWith' (flip (++)) u [v] succs) es + adjLists' (preds,succs) [] = (preds, succs) + + +digraphIsos1 (DG vsa esa) (DG vsb esb) + | length vsa /= length vsb = [] + | length esa /= length esb = [] + | otherwise = digraphIsos' [] vsa vsb + where digraphIsos' xys [] [] = [xys] + digraphIsos' xys (x:xs) ys = + concat [ digraphIsos' ((x,y):xys) xs (L.delete y ys) + | y <- ys, isCompatible (x,y) xys] + isCompatible (x,y) xys = and [ ((x,x') `elem` esa) == ((y,y') `elem` esb) + && ((x',x) `elem` esa) == ((y',y) `elem` esb) + | (x',y') <- xys ] + +digraphIsos2 a b + | length (vertices a) /= length (vertices b) = [] + | L.sort (M.elems indega) /= L.sort (M.elems indegb) = [] + | L.sort (M.elems outdega) /= L.sort (M.elems outdegb) = [] + | otherwise = dfs [] (vertices a) (vertices b) + where (preda,succa) = adjLists a + (predb,succb) = adjLists b + indega = M.map length preda + indegb = M.map length predb + outdega = M.map length succa + outdegb = M.map length succb + isCompatible (x,y) xys = (M.findWithDefault 0 x indega) == (M.findWithDefault 0 y indegb) + && (M.findWithDefault 0 x outdega) == (M.findWithDefault 0 y outdegb) + && and [ (x' `elem` predx) == (y' `elem` predy) + && (x' `elem` succx) == (y' `elem` succy) + | let predx = M.findWithDefault [] x preda, let predy = M.findWithDefault [] y predb, + let succx = M.findWithDefault [] x succa, let succy = M.findWithDefault [] y succb, + (x',y') <- xys] + dfs xys [] [] = [xys] + dfs xys (x:xs) ys = + concat [ dfs ((x,y):xys) xs (L.delete y ys) + | y <- ys, isCompatible (x,y) xys] + +-- For DAGs, can almost certainly do better than the above by using the height partition +-- However see remarks in Poset on orderIsos: +-- What is most efficient will depend on whether you want to list all of them, or just find out whether there are any or not +-- Could also try refining the height partition by (indegree,outdegree) + + +-- doesn't check whether input is a dag +-- if not, then the output will not contain all the vs +heightPartitionDAG dag@(DG vs es) = heightPartition' S.empty [v | v <- vs, v `M.notMember` preds] -- ie vertices with no predecessors + where (preds,succs) = adjLists dag + heightPartition' interior boundary + | null boundary = [] + | otherwise = let interior' = S.union interior $ S.fromList boundary + boundary' = toSet [v | u <- boundary, v <- M.findWithDefault [] u succs, + all (`S.member` interior') (preds M.! v) ] + in boundary : heightPartition' interior' boundary' + +isDAG dag@(DG vs _) = length vs == length (concat (heightPartitionDAG dag)) + +-- Only valid for DAGs, not for digraphs in general +dagIsos dagA@(DG vsA esA) dagB@(DG vsB esB) + | length vsA /= length (concat heightPartA) = error "dagIsos: dagA is not a DAG" + | length vsB /= length (concat heightPartB) = error "dagIsos: dagB is not a DAG" + | map length heightPartA /= map length heightPartB = [] + | otherwise = dfs [] heightPartA heightPartB + where heightPartA = heightPartitionDAG dagA + heightPartB = heightPartitionDAG dagB + (predsA,_) = adjLists dagA + (predsB,_) = adjLists dagB + dfs xys [] [] = [xys] + dfs xys ([]:las) ([]:lbs) = dfs xys las lbs + dfs xys ((x:xs):las) (ys:lbs) = + concat [ dfs ((x,y):xys) (xs:las) (L.delete y ys : lbs) + | y <- ys, isCompatible (x,y) xys] + isCompatible (x,y) xys = + let preds_x = M.findWithDefault [] x predsA + preds_y = M.findWithDefault [] y predsB + in and [ (x' `elem` preds_x) == (y' `elem` preds_y) | (x',y') <- xys] + -- and [ ((x',x) `elem` esA) == ((y',y) `elem` esB) + -- | (x',y') <- xys ] + -- we only need to check predecessors, not successors, because we proceeding by height ordering + +-- can probably do better by intersecting the height partition with the (indegree,outdegree) partition +-- (although on very symmetrical posets such as B n, this won't help at all) + +-- |Are the two DAGs isomorphic? +isDagIso :: (Ord a, Ord b) => Digraph a -> Digraph b -> Bool +isDagIso dagA dagB = (not . null) (dagIsos dagA dagB) + + +perms [] = [[]] +perms (x:xs) = [ls ++ [x] ++ rs | ps <- perms xs, (ls,rs) <- zip (inits ps) (tails ps)] +-- or use L.permutations + +{- +-- orderings compatible with the height partition +heightOrderingsDAG dag@(DG vs es) = heightOrderings' [[]] (heightPartitionDAG dag) + where heightOrderings' initsegs (level:levels) = + let addsegs = perms level + initsegs' = [init ++ add | init <- initsegs, add <- addsegs] + in heightOrderings' initsegs' levels + heightOrderings' segs [] = segs +-} + +isoRepDAG1 dag@(DG vs es) = isoRepDAG' [M.empty] 1 (heightPartitionDAG dag) + where isoRepDAG' initmaps j (level:levels) = + let j' = j + length level + addmaps = [M.fromList (zip ps [j..]) | ps <- perms level] + initmaps' = [init +++ add | init <- initmaps, add <- addmaps] + in isoRepDAG' initmaps' j' levels + isoRepDAG' maps _ [] = DG [1..length vs] (minimum [L.sort (map (\(u,v) -> (m M.! u, m M.! v)) es) | m <- maps]) + initmap +++ addmap = M.union initmap addmap + +-- For example +-- > isoRepDAG1 (DG ['a'..'e'] [('a','c'),('a','d'),('b','d'),('b','e'),('d','e')]) +-- ([1,2,3,4,5],[(1,3),(1,4),(2,3),(2,5),(3,5)]) +-- > isoRepDAG1 (DG ['a'..'e'] [('a','d'),('a','e'),('b','c'),('b','d'),('d','e')]) +-- ([1,2,3,4,5],[(1,3),(1,4),(2,3),(2,5),(3,5)]) + + +-- Find the minimum height-preserving numberings of the vertices, using dfs +isoRepDAG2 dag@(DG vs es) = minimum $ dfs [] srclevels trglevels + where -- (preds,succs) = adjLists dag + srclevels = heightPartitionDAG dag + trglevels = reverse $ fst $ foldl + (\(tls,is) sl -> let (js,ks) = splitAt (length sl) is in (js:tls,ks)) + ([],[1..]) srclevels + dfs xys [] [] = [xys] + dfs xys ([]:sls) ([]:tls) = dfs xys sls tls + dfs xys ((x:xs):sls) (ys:tls) = + concat [ dfs ((x,y):xys) (xs:sls) (L.delete y ys : tls) | y <- ys] + -- not applying any compatibility condition yet + + +-- Find the height-respecting numbering of the vertices which leads to the minimal numbering of the edges +-- So this is calculating the same function as isoRepDAG1, but more efficiently +-- Uses dfs with pruning, rather than exhaustive search +isoRepDAG3 dag@(DG vs es) = dfs root [root] + where n = length vs + root = ([],(1,0),M.empty,(srclevels,trglevels)) -- root of the search tree + (preds,succs) = adjLists dag + srclevels = heightPartitionDAG dag + trglevels = reverse $ fst $ foldl + (\(tls,is) sl -> let (js,ks) = splitAt (length sl) is in (js:tls,ks)) + ([],[1..]) srclevels + dfs best (node:stack) = + -- node : -- for debugging + case cmpPartial best node of + LT -> dfs best stack -- ie prune the search tree at this node + GT -> dfs node (successors node ++ stack) -- ie replace best with this node + EQ -> dfs best (successors node ++ stack) + -- dfs best [] = [best] -- !! for debugging + dfs best@(es',_,_,_) [] = DG [1..n] es' + successors (es,_,_,([],[])) = [] + successors (es,(i,j),m,([]:sls,[]:tls)) = successors (es,(i,j),m,(sls,tls)) + successors (es,(i,j),m,(xs:sls,(y:ys):tls)) = + [ (es', (i',y), m', (L.delete x xs : sls, ys : tls)) + | x <- xs, + let m' = M.insert x y m, + let es' = L.sort $ es ++ [(m M.! u, y) | u <- M.findWithDefault [] x preds], + let i' = nextunfinished m' i ] + -- a vertex is considered finished when all its successors have assignments in the map + nextunfinished m i = + case [v | (v,i') <- M.assocs m, i' == i] of + [] -> i + [u] -> if all (`M.member` m) (M.findWithDefault [] u succs) + then nextunfinished m (i+1) -- i is finished: all successors already have assignments in the map + else i + cmpPartial (es,_,_,_) (es',(i',j'),_,_) = + cmpPartial' (i',j') es es' + -- where j' = maximum $ 0 : map snd es' + cmpPartial' (i',j') ((u,v):es) ((u',v'):es') = + -- Any new e' that can be added to es' must be greater than (i',j') + -- (we don't care about possible extensions of es, because we're not extending them) + case compare (u,v) (u',v') of + EQ -> cmpPartial' (i',j') es es' + LT -> if (u,v) <= (i',j') then LT else EQ + GT -> GT -- always replace best if you beat it + -- (even if it could improve, it's not going to as we're not progressing it) + cmpPartial' (i',j') ((u,v):es) [] = if (u,v) <= (i',j') then LT else EQ + cmpPartial' _ [] ((u',v'):es') = GT -- always extend an existing partial best + cmpPartial' _ [] [] = EQ + + +-- Now we seek a numbering of the vertices which respects height-ordering, +-- and within each height level respects (indegree,outdegree) ordering. +-- We seek the numbering which minimises the resulting edge list. + + +-- |Given a directed acyclic graph (DAG), return a canonical representative for its isomorphism class. +-- @isoRepDAG dag@ is isomorphic to @dag@. It follows that if @isoRepDAG dagA == isoRepDAG dagB@ then @dagA@ is isomorphic to @dagB@. +-- Conversely, @isoRepDAG dag@ is the minimal element in the isomorphism class, subject to some constraints. +-- It follows that if @dagA@ is isomorphic to @dagB@, then @isoRepDAG dagA == isoRepDAG dagB@. +-- +-- The algorithm of course is faster on some DAGs than others: roughly speaking, +-- it prefers \"tall\" DAGs (long chains) to \"wide\" DAGs (long antichains), +-- and it prefers asymmetric DAGs (ie those with smaller automorphism groups). +isoRepDAG :: (Ord a) => Digraph a -> Digraph Int +isoRepDAG dag@(DG vs es) = dfs root [root] + where n = length vs + root = ([],(1,0),M.empty,(srclevels,trglevels)) -- root of the search tree + (preds,succs) = adjLists dag + indegs = M.map length preds + outdegs = M.map length succs + byDegree vs = (map . map) snd $ L.groupBy (\(du,u) (dv,v) -> du == dv) $ L.sort + [( (M.findWithDefault 0 v indegs, M.findWithDefault 0 v outdegs), v) | v <- vs] + srclevels = concatMap byDegree $ heightPartitionDAG dag + trglevels = reverse $ fst $ foldl + (\(tls,is) sl -> let (js,ks) = splitAt (length sl) is in (js:tls,ks)) + ([],[1..]) srclevels + dfs best (node:stack) = + -- node : -- for debugging + case cmpPartial best node of + LT -> dfs best stack -- ie prune the search tree at this node + GT -> dfs node (successors node ++ stack) -- ie replace best with this node + EQ -> dfs best (successors node ++ stack) + -- dfs best [] = [best] -- !! for debugging + dfs best@(es',_,_,_) [] = DG [1..n] es' + successors (es,_,_,([],[])) = [] + successors (es,(i,j),m,([]:sls,[]:tls)) = successors (es,(i,j),m,(sls,tls)) + successors (es,(i,j),m,(xs:sls,(y:ys):tls)) = + [ (es', (i',y), m', (L.delete x xs : sls, ys : tls)) + | x <- xs, + let m' = M.insert x y m, + let es' = L.sort $ es ++ [(m M.! u, y) | u <- M.findWithDefault [] x preds], + let i' = nextunfinished m' i ] + -- a vertex is considered finished when all its successors have assignments in the map + nextunfinished m i = + case [v | (v,i') <- M.assocs m, i' == i] of + [] -> i + [u] -> if all (`M.member` m) (M.findWithDefault [] u succs) + then nextunfinished m (i+1) -- i is finished: all successors already have assignments in the map + else i + cmpPartial (es,_,_,_) (es',(i',j'),_,_) = + cmpPartial' (i',j') es es' + -- where j' = maximum $ 0 : map snd es' + cmpPartial' (i',j') ((u,v):es) ((u',v'):es') = + -- Any new e' that can be added to es' must be greater than (i',j') + -- (we don't care about possible extensions of es, because we're not extending them) + case compare (u,v) (u',v') of + EQ -> cmpPartial' (i',j') es es' + LT -> if (u,v) <= (i',j') then LT else EQ + GT -> GT -- always replace best if you beat it + -- (even if it could improve, it's not going to as we're not progressing it) + cmpPartial' (i',j') ((u,v):es) [] = if (u,v) <= (i',j') then LT else EQ + cmpPartial' _ [] ((u',v'):es') = GT -- always extend an existing partial best + cmpPartial' _ [] [] = EQ diff --git a/Math/Combinatorics/Graph.hs b/Math/Combinatorics/Graph.hs --- a/Math/Combinatorics/Graph.hs +++ b/Math/Combinatorics/Graph.hs @@ -1,4 +1,4 @@ --- Copyright (c) David Amos, 2008. All rights reserved. +-- Copyright (c) 2008-2011, David Amos. All rights reserved. -- |A module defining a polymorphic data type for (simple, undirected) graphs, -- together with constructions of some common families of graphs, @@ -38,10 +38,22 @@ -- GRAPH -- |Datatype for graphs, represented as a list of vertices and a list of edges. --- Both the list of vertices and the list of edges, and also the 2-element lists representing the edges, --- are required to be in ascending order, without duplicates. +-- For most purposes, graphs are required to be in normal form. +-- A graph G vs es is in normal form if (i) vs is in ascending order without duplicates, +-- (ii) es is in ascending order without duplicates, (iii) each e in es is a 2-element list [x,y], x Graph a -> Graph a +nf (G vs es) = G vs' es' where + vs' = L.sort vs + es' = L.sort (map L.sort es) + -- we require that vs, es, and each individual e are sorted isSetSystem xs bs = isListSet xs && isListSet bs && all isListSet bs && all (`isSubset` xs) bs @@ -89,9 +101,14 @@ -- SOME SIMPLE FAMILIES OF GRAPHS -nullGraph :: Graph Int -- type signature needed -nullGraph = G [] [] +-- |The null graph on n vertices is the graph with no edges +nullGraph :: (Integral t) => t -> Graph t +nullGraph n = G [1..n] [] +-- |The null graph, with no vertices or edges +nullGraph' :: Graph Int -- type signature needed +nullGraph' = G [] [] + -- |c n is the cyclic graph on n vertices c :: (Integral t) => t -> Graph t c n = graph (vs,es) where @@ -124,14 +141,8 @@ -- can probably type-coerce this to be Graph [F2] if required q k = fromBinary $ q' k -{- --- note, this definition only in versions >0.1.3 -q'' k = gmap (\v -> v <.> pows2) (q' k) where - pows2 = reverse $ take k $ iterate (*2) 1 - u <.> v = sum $ zipWith (*) u v - gmap f (G vs es) = G (map f vs) ((map . map) f es) --} + tetrahedron = k 4 cube = q 3 @@ -168,20 +179,24 @@ -- return a graph with vertices which are the numbers obtained by interpreting these as digits, eg 123. -- The caller is responsible for ensuring that this makes sense (eg that the small integers are all < 10) fromDigits :: Integral a => Graph [a] -> Graph a +fromDigits = fmap fromDigits' +{- fromDigits (G vs es) = graph (vs',es') where vs' = map fromDigits' vs es' = (map . map) fromDigits' es +-} -- |Given a graph with vertices which are lists of 0s and 1s, -- return a graph with vertices which are the numbers obtained by interpreting these as binary digits. -- For example, [1,1,0] -> 6. fromBinary :: Integral a => Graph [a] -> Graph a +fromBinary = fmap fromBinary' +{- fromBinary (G vs es) = graph (vs',es') where vs' = map fromBinary' vs es' = (map . map) fromBinary' es - +-} --- this definition only in versions >0.1.3 petersen = graph (vs,es) where vs = combinationsOf 2 [1..5] es = [ [v1,v2] | [v1,v2] <- combinationsOf 2 vs, disjoint v1 v2] @@ -194,6 +209,9 @@ complement (G vs es) = graph (vs,es') where es' = combinationsOf 2 vs \\ es -- es' = [e | e <- combinationsOf 2 vs, e `notElem` es] + +inducedSubgraph g@(G vs es) us = G us (es `restrict` us) + where es `restrict` us = [e | e@[i,j] <- es, i `elem` us, j `elem` us] lineGraph g = to1n $ lineGraph' g diff --git a/Math/Combinatorics/GraphAuts.hs b/Math/Combinatorics/GraphAuts.hs --- a/Math/Combinatorics/GraphAuts.hs +++ b/Math/Combinatorics/GraphAuts.hs @@ -110,6 +110,16 @@ -- if p(edge) > 1/2, it would be better to test on the complement of the graph + + +-- Calculate a map consisting of neighbour lists for each vertex in the graph +-- If a vertex has no neighbours then it is left out of the map +adjLists (G vs es) = adjLists' M.empty es + where adjLists' nbrs ([u,v]:es) = + adjLists' (M.insertWith' (flip (++)) v [u] $ M.insertWith' (flip (++)) u [v] nbrs) es + adjLists' nbrs [] = nbrs + + -- ALTERNATIVE VERSIONS OF GRAPH AUTS -- (showing how we got to the final version) @@ -352,6 +362,9 @@ es1 = S.fromList $ edges g1 es2 = S.fromList $ edges g2 +-- !! If we're only interested in seeing whether or not two graphs are iso, +-- !! then the cost of calculating distancePartitions may not be warranted +-- !! (see Math.Combinatorics.Poset: orderIsos01 versus orderIsos) isIso g1 g2 = (not . null) (graphIsos g1 g2) diff --git a/Math/Combinatorics/IncidenceAlgebra.hs b/Math/Combinatorics/IncidenceAlgebra.hs new file mode 100644 --- /dev/null +++ b/Math/Combinatorics/IncidenceAlgebra.hs @@ -0,0 +1,292 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction #-} + + +module Math.Combinatorics.IncidenceAlgebra where + +import Math.Combinatorics.Digraph +import Math.Combinatorics.Poset + +import Math.Algebra.Field.Base +import Math.Algebras.VectorSpace +import Math.Algebras.TensorProduct +import Math.Algebras.Structures + +import Data.List as L +import qualified Data.Map as M +import qualified Data.Set as S + + +-- INTERVALS IN A POSET + +-- |A type to represent an interval in a poset. The (closed) interval [x,y] is the set {z | x <= z <= y} within the poset. +-- Note that the \"empty interval\" is not an interval - that is, the interval [x,y] is only defined for x <= y. +-- The (closed) intervals within a poset form a basis for the incidence algebra as a k-vector space. +data Interval a = Iv (Poset a) (a,a) + +instance Eq a => Eq (Interval a) where + Iv _ (a,b) == Iv _ (a',b') = (a,b) == (a',b') +-- we don't bother to check that they are from the same poset + +instance Ord a => Ord (Interval a) where + compare (Iv _ (a,b)) (Iv _ (a',b')) = compare (a,b) (a',b') + +instance Show a => Show (Interval a) where + show (Iv _ (a,b)) = "Iv (" ++ show a ++ "," ++ show b ++ ")" + +{- +-- !! This should probably be called heightPartition not rank +-- rank is only well-defined if we don't have cover edges jumping levels +rankPartition (Iv poset@(Poset (set,po)) (a,b)) = rankPartition' S.empty [a] (L.delete a iv) + where rankPartition' _ level [] = [level] + rankPartition' interior boundary exterior = + let interior' = S.union interior (S.fromList boundary) + boundary' = toSet [v | (u,v) <- es, u `elem` boundary, all (`S.member` interior') (predecessors es v)] + exterior' = exterior \\ boundary' + in boundary : rankPartition' interior' boundary' exterior' + iv = interval poset (a,b) + (_,es) = coverGraph (Poset (iv,po)) + predecessors es v = [u | (u,v') <- es, v' == v] +-- !! Can be written more efficiently, eg by memoising predecessors and successors, culling covers as we use them, etc. + +-- The point of rankPartition function is to enable a slightly faster isomorphism test +-- Could do even better by refining with (indegree, outdegree) +-} + +-- The sub-poset defined by an interval +ivPoset (Iv poset@(Poset (_,po)) (x,y)) = Poset (interval poset (x,y), po) + +intervalIsos iv1 iv2 = orderIsos (ivPoset iv1) (ivPoset iv2) + +isIntervalIso iv1 iv2 = isOrderIso (ivPoset iv1) (ivPoset iv2) +-- we're only really interested in comparing intervals in the same poset + +{- +intervalIsoMap1 poset = intervalIsoMap' M.empty [Iv poset xy | xy <- L.sort (intervals poset)] + where intervalIsoMap' m (iv:ivs) = + let reps = [iv' | iv' <- M.keys m, m M.! iv' == Nothing, iv `isIntervalIso` iv'] + in if null reps + then intervalIsoMap' (M.insert iv Nothing m) ivs + else let [iv'] = reps in intervalIsoMap' (M.insert iv (Just iv') m) ivs + intervalIsoMap' m [] = m +-} + +-- A poset on n vertices has at most n(n+1)/2 intervals +-- In the worst case, we might have to compare each interval to all earlier intervals +-- Hence this is O(n^4) +intervalIsoMap poset = isoMap + where ivs = [Iv poset xy | xy <- intervals poset] + isoMap = M.fromList [(iv, isoMap' iv) | iv <- ivs] + isoMap' iv = let reps = [iv' | iv' <- ivs, iv' < iv, isoMap M.! iv' == Nothing, iv `isIntervalIso` iv'] + in if null reps then Nothing else let [rep] = reps in Just rep +-- Once an interval is identified as a representative, it is likely to take part in many isomorphism tests +-- Whereas most intervals take part in only one +-- So perhaps we could make this more efficient by having an isomorphism test which uses a height partition +-- for the LHS but not for the RHS? + +-- |List representatives of the order isomorphism classes of intervals in a poset +intervalIsoClasses :: (Ord a) => Poset a -> [Interval a] +intervalIsoClasses poset = [iv | iv <- M.keys isoMap, isoMap M.! iv == Nothing] + where isoMap = intervalIsoMap poset + + +-- INCIDENCE ALGEBRA + +-- |The incidence algebra of a poset is the free k-vector space having as its basis the set of intervals in the poset, +-- with multiplication defined by concatenation of intervals. +-- The incidence algebra can also be thought of as the vector space of functions from intervals to k, with multiplication +-- defined by the convolution (f*g)(x,y) = sum [ f(x,z) g(z,y) | x <= z <= y ]. +instance (Num k, Ord a) => Algebra k (Interval a) where + -- |Note that we are not able to give a generic definition of unit for the incidence algebra, + -- because it depends on which poset we are working in, + -- and that information is encoded at the value level rather than the type level. See unitIA. + unit 0 = zero -- so that sum works + -- unit x = x *> sumv [return (Iv (a,a)) | a <- poset] -- the delta function + -- but we can't know from the types alone which poset we are working in + mult = linear mult' + where mult' (Iv poset (a,b), Iv _ (c,d)) = if b == c then return (Iv poset (a,d)) else zero + +-- So multiplication in the incidence algebra is about composition of intervals + + +-- |The unit of the incidence algebra of a poset +unitIA :: (Num k, Ord t) => Poset t -> Vect k (Interval t) +unitIA poset@(Poset (set,_)) = sumv [return (Iv poset (x,x)) | x <- set] + +basisIA :: Num k => Poset t -> [Vect k (Interval t)] +basisIA poset = [return (Iv poset xy) | xy <- intervals poset] + +-- |The zeta function of a poset +zetaIA :: (Num k, Ord t) => Poset t -> Vect k (Interval t) +zetaIA poset = sumv $ basisIA poset + +-- Then for example, zeta^2 counts the number of points in each interval +-- See Stanley, Enumerative Combinatorics I, p115ff, for more similar + +-- calculate the mobius function of a poset: naive implementation +muIA1 poset@(Poset (set,po)) = sum [mu (x,y) *> return (Iv poset (x,y)) | x <- set, y <- set] + where mu (x,y) | x == y = 1 + | po x y = negate $ sum [mu (x,z) | z <- set, po x z, po z y, z /= y] + | otherwise = 0 + +-- calculate the mobius function of a poset, with memoization +-- |The Mobius function of a poset +muIA :: (Num k, Ord t) => Poset t -> Vect k (Interval t) +muIA poset@(Poset (set,po)) = sumv [mus M.! (x,y) *> return (Iv poset (x,y)) | x <- set, y <- set] + where mu (x,y) | x == y = 1 + | po x y = negate $ sum [mus M.! (x,z) | z <- set, po x z, po z y, z /= y] + | otherwise = 0 + mus = M.fromList [((x,y), mu (x,y)) | x <- set, y <- set] + +-- calculate the inverse of a function in the incidence algebra: naive implementation +invIA1 f | f == zerov = error "invIA 0" + | any (==0) [f' (x,x) | x <- set] = error "invIA: not invertible" + | otherwise = g + where (Iv poset@(Poset (set,po)) _,_) = head $ terms f + f' (x,y) = coeff (Iv poset (x,y)) f + g = sumv [g' xy *> return (Iv poset xy) | xy <- intervals poset] + g' (x,y) | x == y = 1 / f' (x,x) + | otherwise = (-1 / f' (x,x)) * sum [f' (x,z) * g' (z,y) | z <- interval poset (x,y), x /= z] + +-- Stanley, Enumerative Combinatorics I, p144 +-- |The inverse of an element in the incidence algebra of a poset. +-- This is only defined for elements which are non-zero on all intervals (x,x) +invIA :: (Fractional k, Ord t) => Vect k (Interval t) -> Maybe (Vect k (Interval t)) +invIA f | f == zerov = Nothing -- error "invIA 0" + | any (==0) [f' (x,x) | x <- set] = Nothing -- error "invIA: not invertible" + | otherwise = Just g + where (Iv poset@(Poset (set,po)) _,_) = head $ terms f + f' (x,y) = coeff (Iv poset (x,y)) f + g = sumv [g' xy *> return (Iv poset xy) | xy <- intervals poset] + g' (x,y) | x == y = 1 / f' (x,x) + | otherwise = (-1 / f' (x,x)) * sum [f' (x,z) * (g's M.! (z,y)) | z <- interval poset (x,y), x /= z] + g's = M.fromList [(xy, g' xy) | xy <- intervals poset] + +invIA' f = case invIA f of + Just g -> g + Nothing -> error "invIA': not invertible" + +-- Then for example we can count multichains or chains using the incidence algebra - see Stanley + +-- |A function (ie element of the incidence algebra) that counts the total number of chains in each interval +numChainsIA :: (Ord a) => Poset a -> Vect Q (Interval a) +numChainsIA poset = invIA' (2 *> unitIA poset <-> zetaIA poset) + +-- The eta function on intervals (x,y) is 1 if x -< y (y covers x), 0 otherwise +etaIA poset = let DG vs es = hasseDigraph poset + in sumv [return (Iv poset (x,y)) | (x,y) <- es] + +-- |A function (ie element of the incidence algebra) that counts the number of maximal chains in each interval +numMaximalChainsIA :: (Ord a) => Poset a -> Vect Q (Interval a) +numMaximalChainsIA poset = invIA' (unitIA poset <-> etaIA poset) + + +-- In order to quickCheck this, we would need +-- (i) Custom Arbitrary instance - which uses only valid intervals for the poset (ie elts of the basis) +-- (ii) Custom quickCheck property, which uses the correct unit + + +-- SOME KNOWN MOBIUS FUNCTIONS + +muC n = sum [mu' (a,b) *> return (Iv poset (a,b)) | (a,b) <- intervals poset] + where mu' (a,b) | a == b = 1 + | a+1 == b = -1 + | otherwise = 0 + poset = chainN n + +muB n = sumv [(-1)^(length b - length a) *> return (Iv poset (a,b)) | (a,b) <- intervals poset] + where poset = posetB n +-- van Lint & Wilson p335 + +muL n fq = sumv [ ( (-1)^k * q^(k*(k-1) `div` 2) ) *> return (Iv poset (a,b)) | + (a,b) <- intervals poset, + let k = length b - length a ] -- the difference in dimensions + where q = length fq + poset = posetL n fq +-- van Lint & Wilson p335 + + +-- INCIDENCE COALGEBRA +-- Schmitt, Incidence Hopf Algebras + +instance (Num k, Ord a) => Coalgebra k (Interval a) where + counit = unwrap . linear counit' + where counit' (Iv _ (x,y)) = (if x == y then 1 else 0) *> return () + comult = linear comult' + where comult' (Iv poset (x,z)) = sumv [return (Iv poset (x,y), Iv poset (y,z)) | y <- interval poset (x,z)] + +-- So comultiplication in the incidence coalgebra is about decomposition of intervals into subintervals + + +-- But for incidence Hopf algebras, Schmitt wants the basis elts to be isomorphism classes of intervals, not intervals themselves +-- (ie unlabelled intervals) + +-- |@toIsoClasses@ is the linear map from the incidence Hopf algebra of a poset to itself, +-- in which each interval is mapped to (the minimal representative of) its isomorphism class. +-- Thus the result can be considered as a linear combination of isomorphism classes of intervals, +-- rather than of intervals themselves. +-- Note that if this operation is to be performed repeatedly for the same poset, +-- then it is more efficient to use @toIsoClasses' poset@, which memoizes the isomorphism class lookup table. +toIsoClasses :: (Num k, Ord a) => Vect k (Interval a) -> Vect k (Interval a) +toIsoClasses v + | v == zerov = zerov + | otherwise = toIsoClasses' poset v + where (Iv poset _, _) = head $ terms v + +-- |Given a poset, @toIsoClasses' poset@ is the linear map from the incidence Hopf algebra of the poset to itself, +-- in which each interval is mapped to (the minimal representative of) its isomorphism class. +toIsoClasses' :: (Num k, Ord a) => Poset a -> Vect k (Interval a) -> Vect k (Interval a) +toIsoClasses' poset = linear isoRep + where isoRep iv = case isoMap M.! iv of + Nothing -> return iv + Just iv' -> return iv' + isoMap = intervalIsoMap poset + + +{- +-- for example: + +> toIsoClasses $ zetaIA $ posetP 4 +15Iv ([[1],[2],[3],[4]],[[1],[2],[3],[4]])+31Iv ([[1],[2],[3],[4]],[[1],[2],[3,4]])+10Iv ([[1],[2],[3],[4]],[[1],[2,3,4]])+3Iv ([[1],[2],[3],[4]],[[1,2],[3,4]])+Iv ([[1],[2],[3],[4]],[[1,2,3,4]]) + +-- Can we use this to solve "counting squares" problems + +> let b3 = comult $ return $ Iv (posetB 3) ([],[1,2,3]) +> let isoB3 = toIsoClasses' $ posetB 3 +> (isoB3 `tf` isoB3) b3 +(Iv ([],[]),Iv ([],[1,2,3]))+3(Iv ([],[1]),Iv ([],[1,2]))+3(Iv ([],[1,2]),Iv ([],[1]))+(Iv ([],[1,2,3]),Iv ([],[])) + +-- The incidence coalgebra of the binomial poset is isomorphic to the binomial coalgebra + +-- if we just want to get the coefficients, we don't need to use comult: + +> let poset@(Poset (set,po)) = posetB 3 in toIsoClasses $ sumv [return (Iv poset ([],x)) | x <- set] +Iv ([],[])+3Iv ([],[1])+3Iv ([],[1,2])+Iv ([],[1,2,3]) + +> let n = 4; p = comult $ return $ Iv (posetP n) ([[i] | i<- [1..n]],[[1..n]]); iso = toIsoClasses' (posetP n) in (iso `tf` iso) p +(Iv ([[1],[2],[3],[4]],[[1],[2],[3],[4]]),Iv ([[1],[2],[3],[4]],[[1,2,3,4]]))+ +6(Iv ([[1],[2],[3],[4]],[[1],[2],[3,4]]),Iv ([[1],[2],[3],[4]],[[1],[2,3,4]]))+ +4(Iv ([[1],[2],[3],[4]],[[1],[2,3,4]]),Iv ([[1],[2],[3],[4]],[[1],[2],[3,4]]))+ +3(Iv ([[1],[2],[3],[4]],[[1,2],[3,4]]),Iv ([[1],[2],[3],[4]],[[1],[2],[3,4]]))+ +(Iv ([[1],[2],[3],[4]],[[1,2,3,4]]),Iv ([[1],[2],[3],[4]],[[1],[2],[3],[4]])) + +-- These are multinomial coefficients, OEIS A178867: 1; 1,1; 1,3,1; 1,6,4,3,1; 1,10,10,15,5,10,1; ... +-- Although A036040, which is the same up to ordering, seems a better match. (Our order is fairly arbitrary) + +> let n = 4; p = comult $ return $ Iv (posetL n f2) ([],[[1 :: F2,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]); iso = toIsoClasses' (posetL n f2) in (iso `tf` iso) p +(Iv ([],[]),Iv ([],[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]))+ +15(Iv ([],[[0,0,0,1]]),Iv ([],[[0,1,0,0],[0,0,1,0],[0,0,0,1]]))+ +35(Iv ([],[[0,0,1,0],[0,0,0,1]]),Iv ([],[[0,0,1,0],[0,0,0,1]]))+ +15(Iv ([],[[0,1,0,0],[0,0,1,0],[0,0,0,1]]),Iv ([],[[0,0,0,1]]))+ +(Iv ([],[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]),Iv ([],[])) + +-- With L n fq, we get the q-binomial coefficients, eg OEIS A022166: +1; 1, 1; 1, 3, 1; 1, 7, 7, 1; 1, 15, 35, 15, 1 +-} + + +-- This still isn't quite what Schmitt wants +-- Schmitt, IHA, p6 +-- The incidence Hopf algebra should have as its basis isomorphism classes of intervals, not intervals +-- The mult is defined as direct product of posets diff --git a/Math/Combinatorics/Poset.hs b/Math/Combinatorics/Poset.hs new file mode 100644 --- /dev/null +++ b/Math/Combinatorics/Poset.hs @@ -0,0 +1,282 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +{-# LANGUAGE NoMonomorphismRestriction, TupleSections #-} + + +module Math.Combinatorics.Poset where + +import Math.Algebra.Field.Base +import Math.Combinatorics.FiniteGeometry +import Math.Algebra.LinearAlgebra + +import Math.Combinatorics.Digraph + +import Data.List as L +import qualified Data.Map as M +-- import qualified Data.Set as S + + +-- |A poset is represented as a pair (set,po), where set is the underlying set of the poset, and po is the partial order relation +newtype Poset t = Poset ([t], t -> t -> Bool) + +instance Eq t => Eq (Poset t) where + Poset (set,po) == Poset (set',po') = + set == set' && and [po x y == po' x y | x <- set, y <- set] +-- There may be ways to avoid comparing every pair +-- If we could calculate the coverGraph without comparing every pair, +-- then it would be sufficient to test whether their cover graphs are equal + +instance Show t => Show (Poset t) where + show (Poset (set,po)) = "Poset " ++ show set + +implies p q = q || not p + +isReflexive (set,po) = and [x `po` x | x <- set] +isAntisymmetric (set,po) = and [((x `po` y) && (y `po` x)) `implies` (x == y) | x <- set, y <- set] +isTransitive (set,po) = and [((x `po` y) && (y `po` z)) `implies` (x `po` z) | x <- set, y <- set, z <- set] + +isPoset poset = isReflexive poset && isAntisymmetric poset && isTransitive poset + +poset (set,po) + | isPoset (set,po) = Poset (set,po) + | otherwise = error "poset: Not a partial order" + + +-- Most of the posets we will deal with are in fact lattices, meaning that any two elements +-- have a meet (greatest lower bound) and join (least upper bound) + +intervals (Poset (set,po)) = [(a,b) | a <- set, b <- set, a `po` b] + +interval (Poset (set,po)) (x,z) = [y | y <- set, x `po` y, y `po` z] + + +-- LINEAR ORDER POSET +-- This is of course a lattice, with meet = min, join = max + +-- |A chain is a poset in which every pair of elements is comparable (ie either x <= y or y <= x). +-- It is therefore a linear or total order. +-- chainN n is the poset consisting of the numbers [1..n] ordered by (<=) +chainN :: Int -> Poset Int +chainN n = Poset ( [1..n], (<=) ) + +-- hasseN n = DG [1..n] [(i,i+1) | i <- [1..n-1]] + + +-- |An antichain is a poset in which distinct elements are incomparable. +-- antichainN n is the poset consisting of [1..n], with x <= y only when x == y. +antichainN :: Int -> Poset Int +antichainN n = Poset ( [1..n], (==) ) + + +-- LATTICE OF (POSITIVE) DIVISORS OF N + +divides a b = b `mod` a == 0 + +divisors n | n >= 1 = [a | a <- [1..n], a `divides` n] + +-- |posetD n is the lattice of (positive) divisors of n +posetD :: Int -> Poset Int +posetD n | n >= 1 = Poset ( divisors n, divides ) + + +-- LATTICE OF SUBSETS OF [1..N] ORDERED BY INCLUSION +-- (Boolean lattice) + +powerset [] = [[]] +powerset (x:xs) = let p = powerset xs in p ++ map (x:) p + +-- subset test for sorted lists +isSubset (x:xs) (y:ys) = + case compare x y of + LT -> False + EQ -> isSubset xs ys + GT -> isSubset (x:xs) ys +isSubset [] _ = True +isSubset _ [] = False + +-- |posetB n is the lattice of subsets of [1..n] ordered by inclusion +posetB :: Int -> Poset [Int] +posetB n = Poset ( powerset [1..n], isSubset ) + + +-- LATTICE OF PARTITIONS OF [1..N] ORDERED BY REFINEMENT + +partitions [] = [[]] +partitions [x] = [[[x]]] +partitions (x:xs) = let ps = partitions xs in + map ([x]:) ps ++ [ (x:cell):(L.delete cell p) | p <- ps, cell <- p] +-- if the input is sorted, then so is the output + +isRefinement a b = and [or [acell `isSubset` bcell | bcell <- b] | acell <- a] +-- if we know that a and b are appropriately sorted, then this can probably be done more efficiently + +-- |posetP n is the lattice of partitions of [1..n] ordered by refinement +posetP :: Int -> Poset [[Int]] +posetP n = Poset ( partitions [1..n], isRefinement ) + +-- muP n = ... +-- see van Lint and Wilson p336 + + +-- LATTICE OF INTERVAL PARTITIONS OF [1..N] ORDERED BY REFINEMENT + +intervalPartitions xs = filter (all isInterval) (partitions xs) + +isInterval (x1:x2:xs) = x1+1 == x2 && isInterval (x2:xs) +isInterval _ = True + +intervalPartitions2 [] = [[]] +intervalPartitions2 [x] = [[[x]]] +intervalPartitions2 (x:xs) = let ips = intervalPartitions xs in + map ([x]:) ips ++ [ (x:head):tail | (head:tail) <- ips] +-- we're guaranteed that x+1 is at the head of the head + + +-- LATTICE OF SUBSPACES OF Fq^n + +subspaces fq n = [] : concatMap (flatsPG (n-1) fq) [0..n-1] +-- the PG(n-1,Fq) is the set of subspaces of fq^n + +isZero v = all (==0) v + +-- inSpanRE m v returns whether the vector v is in the span of the matrix m, where m is required to be in row echelon form +inSpanRE ((1:xs):bs) (y:ys) = inSpanRE (map tail bs) (if y == 0 then ys else ys <-> y *> xs) +inSpanRE ((0:xs):bs) (y:ys) = if y == 0 then inSpanRE (xs : map tail bs) ys else False +inSpanRE _ ys = isZero ys + +isSubspace s1 s2 = all (inSpanRE s2) s1 + +-- This is the projective geometry PG(n,q) +-- |posetL n fq is the lattice of subspaces of the vector space Fq^n, ordered by inclusion. +-- Subspaces are represented by their reduced row echelon form. +posetL :: FiniteField fq => Int -> [fq] -> Poset [[fq]] +posetL n fq = Poset ( subspaces fq n, isSubspace ) + + +-- choose n k = product [n-k+1..n] `div` product [1..k] + + +-- NEW FROM OLD CONSTRUCTIONS + + +-- |The subposet of a poset satisfying a predicate +subposet :: Poset a -> (a -> Bool) -> Poset a +subposet (Poset (set,po)) p = Poset (filter p set, po) + +-- |The direct sum of two posets +dsum :: Poset a -> Poset b -> Poset (Either a b) +dsum (Poset (setA,poA)) (Poset (setB,poB)) = Poset (set,po) + where set = map Left setA ++ map Right setB + po (Left a1) (Left a2) = poA a1 a2 + po (Right b1) (Right b2) = poB b1 b2 + po _ _ = False + +-- |The direct product of two posets +dprod :: Poset a -> Poset b -> Poset (a,b) +dprod (Poset (setA,poA)) (Poset (setB,poB)) = + Poset ( [(a,b) | a <- setA, b <- setB], \(a1,b1) (a2,b2) -> (a1 `poA` a2) && (b1 `poB` b2) ) + +-- |The dual of a poset +dual :: Poset a -> Poset a +dual (Poset (set, po)) = Poset (set, po') + where po' x y = po y x + + +-- ANALYSIS OF POSETS + +-- |Given a poset (X,<=), we say that y covers x, written x -< y, if x < y and there is no z in X with x < z < y. +-- The Hasse digraph of a poset is the digraph whose vertices are the elements of the poset, +-- with an edge between every pair (x,y) with x -< y. +-- The Hasse digraph can be represented diagrammatically as a Hasse diagram, by drawing x below y whenever x -< y. +hasseDigraph :: (Eq a) => Poset a -> Digraph a +hasseDigraph (Poset (set,po)) = DG set [(x,y) | x <- set, y <- set, x -< y] + where x -< y = x /= y && x `po` y && null [z | z <- set, x `po` z, x /= z, z `po` y, z /= y] +-- The partial order can be recovered as the transitive closure of the covers relation +-- !! Can we construct the cover graph without having to compare every pair ?? + +-- If we know in advance the poset we're interested in, +-- then we're probably better off constructing the Hasse digraph directly + +-- (In effect, the transitive closure of the edge relation) +-- |Given a DAG (directed acyclic graph), return the poset consisting of the vertices of the DAG, ordered by reachability. +-- This can be used to recover a poset from its Hasse digraph. +reachabilityPoset :: (Ord a) => Digraph a -> Poset a +reachabilityPoset (DG vs es) = Poset (vs,tc') -- \u v -> tc M.! (u,v) + where tc = M.fromList [ ((u,v), tc' u v) | u <- vs, v <- vs] + tc' u v | u == v = True + | otherwise = or [tc M.! (w,v) | w <- successors u] + successors u = [v | (u',v) <- es, u' == u] +-- !! looks like we could memoise more than we are doing + + +{- +-- For example: +> let poset = posetB 3 in poset == reachabilityPoset (hasseDigraph poset) +True +-} + + +isOrderPreserving :: (a -> b) -> Poset a -> Poset b -> Bool +isOrderPreserving f (Poset (seta,poa)) (Poset (setb,pob)) = + and [ x `poa` y == f x `pob` f y | x <- seta, y <- seta ] + +-- Find all order isomorphisms between two posets +-- This algorithm is faster to find out whether or not there are any +orderIsos01 (Poset (seta,poa)) (Poset (setb,pob)) + | length seta /= length setb = [] + | otherwise = orderIsos' [] seta setb + where orderIsos' xys [] [] = [xys] + orderIsos' xys (x:xs) ys = + concat [ orderIsos' ((x,y):xys) xs (L.delete y ys) + | y <- ys, and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ] + +-- |Are the two posets order-isomorphic? +isOrderIso :: (Eq a, Eq b) => Poset a -> Poset b -> Bool +isOrderIso poseta posetb = (not . null) (orderIsos01 poseta posetb) + +-- Find all order isomorphisms between two posets +-- This algorithm is faster to find all isomorphisms, if there are many +-- (It may be that it is faster to find any, for large posets, but the break-even point seems to be quite big) +orderIsos posetA@(Poset (_,poa)) posetB@(Poset (_,pob)) + | map length heightPartA /= map length heightPartB = [] + | otherwise = dfs [] heightPartA heightPartB + where heightPartA = heightPartitionDAG (hasseDigraph posetA) + heightPartB = heightPartitionDAG (hasseDigraph posetB) + dfs xys [] [] = [xys] + dfs xys ([]:las) ([]:lbs) = dfs xys las lbs + dfs xys ((x:xs):las) (ys:lbs) = + concat [ dfs ((x,y):xys) (xs:las) (L.delete y ys : lbs) + | y <- ys, and [ (x `poa` x', x' `poa` x) == (y `pob` y', y' `pob` y) | (x',y') <- xys ] ] +-- A variant on this algorithm would use the Hasse digraph rather than the partial order in the test on the last line +-- This might be faster, depending how expensive the partial order comparison function is +-- In effect though, it would then be a DAG isomorphism function + +-- The order automorphisms of a poset +orderAuts1 poset = orderIsos poset poset +-- This returns all automorphisms +-- What we really want is to return generators of the permutation group + + + + + +pairs (x:xs) = map (x,) xs ++ pairs xs -- TupleSections +pairs [] = [] + +-- |A linear extension of a poset is a linear ordering of the elements which extends the partial order. +-- Equivalently, it is an ordering [x1..xn] of the underlying set, such that if xi <= xj then i <= j. +isLinext (Poset (set,po)) set' = all (\(x,y) -> not (y `po` x)) (pairs set') + + +-- |Linear extensions of a poset +linexts (Poset (set,po)) = linexts' [[]] set + where linexts' lss (r:rs) = + let lss' = [ lts ++ [r] ++ gts + | ls <- lss, + let ls' = takeWhile (not . (r `po`)) ls, + (lts,gts) <- zip (inits ls') (tails ls), + all (not . (`po` r)) gts ] + in linexts' lss' rs + linexts' lss [] = lss + + diff --git a/Math/QuantumAlgebra/TensorCategory.hs b/Math/QuantumAlgebra/TensorCategory.hs --- a/Math/QuantumAlgebra/TensorCategory.hs +++ b/Math/QuantumAlgebra/TensorCategory.hs @@ -14,7 +14,10 @@ id_ :: Ob c -> Ar c source, target :: Ar c -> Ob c (>>>) :: Ar c -> Ar c -> Ar c +-- Note that the class Category defined in Control.Category is about categories whose objects are Haskell types, +-- whereas we want the objects to be values of a single type. + -- Kassel p282 -- The following is actually definition of a strict tensor category class Category c => TensorCategory c where @@ -34,6 +37,7 @@ instance (TensorCategory c, Eq (Ar c), Show (Ar c)) => Num (Ar c) where (*) = tar -} + -- SYMMETRIC GROUPOID diff --git a/Math/Test/TAlgebras/TMatrix.hs b/Math/Test/TAlgebras/TMatrix.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TAlgebras/TMatrix.hs @@ -0,0 +1,32 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +{-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +module Math.Test.TAlgebras.TMatrix where + +import Test.QuickCheck + +import Math.Algebra.Field.Base + +import Math.Algebras.VectorSpace +import Math.Algebras.TensorProduct +import Math.Algebras.Matrix + +import Math.Test.TAlgebras.TVectorSpace +import Math.Test.TAlgebras.TStructures + +import Math.Algebras.Structures -- not really needed + +instance Arbitrary Mat2 where + arbitrary = elements [E2 1 1, E2 1 2, E2 2 1, E2 2 2] + +instance Arbitrary Mat2' where + arbitrary = elements [E2' 1 1, E2' 1 2, E2' 2 1, E2' 2 2] + + +prop_Algebra_Mat2 (k,x,y,z) = prop_Algebra (k,x,y,z) + where types = (k,x,y,z) :: (Q, Vect Q Mat2, Vect Q Mat2, Vect Q Mat2) + +prop_Coalgebra_Mat2' x = prop_Coalgebra x + where types = x :: Vect Q Mat2' diff --git a/Math/Test/TAlgebras/TQuaternions.hs b/Math/Test/TAlgebras/TQuaternions.hs --- a/Math/Test/TAlgebras/TQuaternions.hs +++ b/Math/Test/TAlgebras/TQuaternions.hs @@ -1,39 +1,72 @@ -- Copyright (c) 2010, David Amos. All rights reserved. {-# LANGUAGE FlexibleInstances, TypeSynonymInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} module Math.Test.TAlgebras.TQuaternions where import Test.QuickCheck +import Math.Algebra.Field.Base import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct import Math.Algebras.Quaternions +import Math.Test.TAlgebras.TVectorSpace import Math.Test.TAlgebras.TStructures +import Math.Algebras.Structures -- not really needed + instance Arbitrary HBasis where arbitrary = elements [One,I,J,K] +-- TVectorSpace defines an Arbitrary instance for Vect k b, given Arbitrary instances for k and b +{- instance Arbitrary (Quaternion Integer) where arbitrary = do ts <- arbitrary :: Gen [(HBasis, Integer)] return $ nf $ V ts - +-} prop_Algebra_Quaternion (k,x,y,z) = prop_Algebra (k,x,y,z) - where types = (k,x,y,z) :: (Integer, Quaternion Integer, Quaternion Integer, Quaternion Integer) + where types = (k,x,y,z) :: (Q, Quaternion Q, Quaternion Q, Quaternion Q) +-- (Integer, Quaternion Integer, Quaternion Integer, Quaternion Integer) +prop_Coalgebra_DualQuaternion x = prop_Coalgebra x + where types = x :: Vect Q (Dual HBasis) + +conjH = linear conjH' + where conjH' One = 1 + conjH' I = -i + conjH' J = -j + conjH' K = -k + +normH x = x * conjH x + + +-- The following property fails: conjugation is not an algebra morphism +-- It fails to commute with mult: conjH (i*j) /= conjH i * conjH j +nonprop_AlgebraMorphism_ConjH = prop_AlgebraMorphism conjH + +-- The following property also fails: norm is not an algebra morphism +-- It fails to commute with unit: conjH (unit (-1)) /= unit (-1) +nonprop_AlgebraMorphism_NormH = prop_AlgebraMorphism normH + + +{- prop_Coalgebra_Quaternion x = prop_Coalgebra x where types = x :: Quaternion Integer -- Fails - the algebra and coalgebra structures I've given are not compatible prop_Bialgebra_Quaternion (k,x,y) = prop_Bialgebra (k,x,y) where types = (k,x,y) :: (Integer, Quaternion Integer, Quaternion Integer) - +-} {- prop_FrobeniusRelation_Quaternion (x,y) = prop_FrobeniusRelation (x,y) where types = (x,y) :: (Quaternion Integer, Quaternion Integer) -- !! fails, because the counit we have given is not a Frobenius form -} +instance Algebra2 Integer (Quaternion Integer) where + unit2 k = unit k + mult2 xy = mult xy diff --git a/Math/Test/TAlgebras/TStructures.hs b/Math/Test/TAlgebras/TStructures.hs --- a/Math/Test/TAlgebras/TStructures.hs +++ b/Math/Test/TAlgebras/TStructures.hs @@ -1,6 +1,7 @@ -- Copyright (c) 2010, David Amos. All rights reserved. {-# LANGUAGE NoMonomorphismRestriction #-} +{-# LANGUAGE TypeFamilies, RankNTypes, MultiParamTypeClasses #-} module Math.Test.TAlgebras.TStructures where @@ -48,25 +49,31 @@ -- ALGEBRAS prop_Algebra (k,x,y,z) = - mult (x `te` mult (y `te` z)) == mult (mult (x `te` y) `te` z) && -- associativity - smultL k x == mult (unit k `te` x) && -- left unit - -- mult (k' `te` x) == (mult . (unit' `tf` id)) (k' `te` x) && -- left unit - smultR x k == mult (x `te` unit k) -- && -- right unit - -- mult (x `te` k') == (mult . (id `tf` unit')) (x `te` k') -- right unit - where k' = V [( (),k)] + (mult . (id `tf` mult)) (x `te` (y `te` z)) == + (mult . (mult `tf` id)) ((x `te` y) `te` z) && -- associativity + unitOutL (k' `te` x) == (mult . (unit' `tf` id)) (k' `te` x) && -- left unit + unitOutR (x `te` k') == (mult . (id `tf` unit')) (x `te` k') -- right unit + -- mult (x `te` mult (y `te` z)) == mult (mult (x `te` y) `te` z) && -- associativity + -- smultL k x == mult (unit k `te` x) && -- left unit + -- smultR x k == mult (x `te` unit k) -- && -- right unit + where k' = k *> return () -- additionally, unit and mult must be linear prop_Commutative (x,y) = let xy = x `te` y in (mult . twist) xy == mult xy +prop_Algebra_DSum (k,(a1,a2,a3),(b1,b2,b3)) = prop_Algebra (k, a1 `dsume` b1, a2 `dsume` b2, a3 `dsume` b3) +prop_Algebra_TProd (k,(a1,a2,a3),(b1,b2,b3)) = prop_Algebra (k, a1 `te` b1, a2 `te` b2, a3 `te` b3) + + -- COALGEBRAS prop_Coalgebra x = ((comult `tf` id) . comult) x == (assocL . (id `tf` comult) . comult) x && -- coassociativity - ((counit' `tf` id) . comult) x == V [((),1)] `te` x && -- left counit - ((id `tf` counit') . comult) x == x `te` V [((),1)] -- right counit + ((counit' `tf` id) . comult) x == unitInL x && -- left counit + ((id `tf` counit') . comult) x == unitInR x -- right counit -- additionally, counit and comult must be linear prop_Cocommutative x = @@ -75,20 +82,17 @@ -- MORPHISMS -prop_AlgebraMorphism f (k,l,x,y) = - prop_Linear f (k,l,x,y) && - -- (f . unit) k == unit k && +prop_AlgebraMorphism f (k,x,y) = + (f . unit) k == unit k && (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y) -- in this version we supply z of the intended return type of f, -- so that we can make sure we select the correct instance for f polymorphic in return type prop_AlgebraMorphism' f (k,l,x,y,z) = - prop_Linear f (k,l,x,y) && (f . unit) k + z == unit k + z && (f . mult) (x `te` y) == (mult . (f `tf` f)) (x `te` y) prop_CoalgebraMorphism f x = - -- prop_Linear f (k,l,x,y) && (counit . f) x == counit x && ( (f `tf` f) . comult) x == (comult . f) x @@ -107,7 +111,7 @@ prop_Bialgebra2 (k,xy) = (comult . unit') k' + xy == ((unit' `tf` unit') . iso) k' + xy where iso = fmap (\ () -> ((),()) ) -- the isomorphism k ~= k tensor k - k' = unit k :: Trivial Integer -- inject into the trivial algebra + k' = wrap k -- inject into the trivial algebra -- the +xy is just to force the other expression to be of the right type prop_Bialgebra3 (x,y) = @@ -149,6 +153,17 @@ prop_Module_Unit (k,m) = (action . (unit' `tf` id)) k' == -} + + +-- ALTERNATIVE DEFINITION OF ALGEBRA + +type TensorProd k u v = + (u ~ Vect k a, v ~ Vect k b) => Vect k (Tensor a b) + +class Algebra2 k a where + unit2 :: k -> a + mult2 :: TensorProd k a a -> a + -- FROBENIUS ALGEBRAS diff --git a/Math/Test/TAlgebras/TTensorProduct.hs b/Math/Test/TAlgebras/TTensorProduct.hs --- a/Math/Test/TAlgebras/TTensorProduct.hs +++ b/Math/Test/TAlgebras/TTensorProduct.hs @@ -9,7 +9,7 @@ import Math.Algebras.VectorSpace import Math.Algebras.TensorProduct import Math.Algebra.Field.Base -import Math.Test.TAlgebras.TVectorSpace hiding (i1, i2) +import Math.Test.TAlgebras.TVectorSpace import Prelude as P import Control.Category as C @@ -35,7 +35,25 @@ In the definition of `it': it = e1 `te` e2 :: TensorProd Q En En -} +-- check that tf is linear +prop_Linear_tf ((f,g),k,(a1,a2,b1,b2)) = prop_Linear (linfun f `tf` linfun g) (k, a1 `te` b1, a2 `te` b2) + where types = (f,g,k,a1,a2,b1,b2) :: (LinFun Q ABasis SBasis, LinFun Q BBasis TBasis, Q, + Vect Q ABasis, Vect Q ABasis, Vect Q BBasis, Vect Q BBasis) + +-- check that tensor product is a functor, as required +prop_TensorFunctor ((f1,f2,g1,g2),(a,b)) = + (P.id `tf` P.id) (a `te` b) == P.id (a `te` b) && + ((f' P.. f) `tf` (g' P.. g)) (a `te` b) == ((f' `tf` g') P.. (f `tf` g)) (a `te` b) + where f = linfun f1 + f' = linfun f2 + g = linfun g1 + g' = linfun g2 + types = (f1,f2,g1,g2,a,b) :: (LinFun Q ABasis ABasis, LinFun Q ABasis ABasis, + LinFun Q BBasis BBasis, LinFun Q BBasis BBasis, + Vect Q ABasis, Vect Q BBasis) + + -- Now test eg -- > quickCheck (\x -> (distrL . undistrL) x == id x) -- but need to make x be of interesting type (not just () ) @@ -65,16 +83,13 @@ instance Num k => Arrow (Linear k) where arr f = Linear (fmap f) -- requires nf call afterwards - first (Linear f) = Linear $ \(V ts) -> V $ - concat [let V us = x *> te (f $ return a) (return c) in us | ((a,c),x) <- ts] - second (Linear f) = Linear $ \(V ts) -> V $ - concat [let V us = x *> te (return c) (f $ return a) in us | ((c,a),x) <- ts] + first (Linear f) = Linear f *** Linear P.id + second (Linear f) = Linear P.id *** Linear f Linear f *** Linear g = Linear (f `tf2` g) where tf2 f g (V ts) = V $ concat [let V us = x *> te (f $ return a) (g $ return b) in us | ((a,b), x) <- ts] -- can't use tf, as it uses add, which assumes Ord instance - -- hence we should call nf afterwards - -- !! What about &&& + Linear f &&& Linear g = (Linear f *** Linear g) C.. Linear (\a -> a `te` a) {- -- The following are morally correct, but don't work because they require Ord instance diff --git a/Math/Test/TAlgebras/TVectorSpace.hs b/Math/Test/TAlgebras/TVectorSpace.hs --- a/Math/Test/TAlgebras/TVectorSpace.hs +++ b/Math/Test/TAlgebras/TVectorSpace.hs @@ -30,6 +30,11 @@ arbitrary = do n <- arbitrary :: Gen Int return (E n) +instance Arbitrary b => Arbitrary (Dual b) where + arbitrary = fmap Dual arbitrary +-- arbitrary = do b <- arbitrary :: Gen b -- ScopedTypeVariables +-- return (Dual b) + instance Arbitrary Q where arbitrary = do n <- arbitrary :: Gen Integer d <- arbitrary :: Gen Integer diff --git a/Math/Test/TCombinatorics/TDigraph.hs b/Math/Test/TCombinatorics/TDigraph.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCombinatorics/TDigraph.hs @@ -0,0 +1,93 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +module Math.Test.TCombinatorics.TDigraph where + +import Test.HUnit +import Control.Monad (when, unless) + +import Math.Combinatorics.Digraph +import Math.Combinatorics.Poset + + +testlistDigraph = TestList [ + testlistIsDagIsoPositive, + testlistIsDagIsoNegative, + testlistIsoRepDAGIsIso, + testlistIsoRepDAGPositive, + testlistIsoRepDAGNegative + ] + + +testcaseIsDagIsoPositive desc dag1 dag2 = TestCase $ assertBool desc $ isDagIso dag1 dag2 + +testlistIsDagIsoPositive = TestList [ + testcaseIsDagIsoPositive "D 30 ~= D 42" (hasseDigraph $ posetD 30) (hasseDigraph $ posetD 42), + testcaseIsDagIsoPositive "D 60 ~= D 90" (hasseDigraph $ posetD 30) (hasseDigraph $ posetD 42), + testcaseIsDagIsoPositive "D 30 ~= B 3" (hasseDigraph $ posetD 30) (hasseDigraph $ posetB 3), + testcaseIsDagIsoPositive "B 2 ~= 2 * 2" (hasseDigraph $ posetB 2) (hasseDigraph $ dprod (chainN 2) (chainN 2)) + ] + +testcaseIsDagIsoNegative desc dag1 dag2 = TestCase $ assertBool desc $ not (isDagIso dag1 dag2) + +testlistIsDagIsoNegative = TestList [ + testcaseIsDagIsoNegative "D 20 ~/= D 30" (hasseDigraph $ posetD 20) (hasseDigraph $ posetD 30), + testcaseIsDagIsoNegative "Subposets B4 - 1" (hasseDigraph $ subposet (posetB 4) (/= [1])) + (hasseDigraph $ subposet (posetB 4) (/= [1,2])), + testcaseIsDagIsoNegative "Subposets B4 - 2" (hasseDigraph $ subposet (posetB 4) (`notElem` [[1],[1,2]])) + (hasseDigraph $ subposet (posetB 4) (`notElem` [[1],[2,3]])) + ] + + +-- test that the isoRepDAG is isomorphic to the DAG +testcaseIsoRepDAGIsIso desc dag = TestCase $ assertBool desc $ isDagIso dag (isoRepDAG dag) + +testlistIsoRepDAGIsIso = TestList [ + testcaseIsoRepDAGIsIso "D 30" (hasseDigraph $ posetD 30), + testcaseIsoRepDAGIsIso "B 4 - [1,2]" (hasseDigraph $ subposet (posetB 4) (/= [1,2])) + ] + + +testcaseIsoRepDAGPositive desc dag1 dag2 = TestCase (assertEqual desc (isoRepDAG dag1) (isoRepDAG dag2)) + +testlistIsoRepDAGPositive = TestList [ + testcaseIsoRepDAGPositive "D 30 ~= D 42" (hasseDigraph $ posetD 30) (hasseDigraph $ posetD 42), + testcaseIsoRepDAGPositive "D 60 ~= D 90" (hasseDigraph $ posetD 30) (hasseDigraph $ posetD 42), + testcaseIsoRepDAGPositive "D 30 ~= B 3" (hasseDigraph $ posetD 30) (hasseDigraph $ posetB 3), + testcaseIsoRepDAGPositive "B 2 ~= 2 * 2" (hasseDigraph $ posetB 2) (hasseDigraph $ dprod (chainN 2) (chainN 2)) + ] + + +assertNotEqual desc val1 val2 = + when (val1 == val2) (assertFailure desc) + -- unless (val1 /= val2) (assertFailure desc) + +testcaseIsoRepDAGNegative desc dag1 dag2 = TestCase (assertNotEqual desc (isoRepDAG dag1) (isoRepDAG dag2)) + +testlistIsoRepDAGNegative = TestList [ + testcaseIsoRepDAGNegative "Subposets B4 - 1" (hasseDigraph $ subposet (posetB 4) (/= [1])) + (hasseDigraph $ subposet (posetB 4) (/= [1,2])), + testcaseIsoRepDAGNegative "Subposets B4 - 2" (hasseDigraph $ subposet (posetB 4) (`notElem` [[1],[1,2]])) + (hasseDigraph $ subposet (posetB 4) (`notElem` [[1],[2,3]])) + ] + + + +allDags n = [DG [1..n] es | es <- powerset (pairs [1..n])] + +-- > all (uncurry (==)) [(dag1 `isDagIso` dag2, isoRepDAG dag1 == isoRepDAG dag2) | (dag1,dag2) <- pairs (allDags 4)] + +{- +-- Following tests no longer valid, as isoRepDAG doesn't produce same representative as isoRepDAG1 +testcaseIsoRepDAG desc dag = TestCase (assertEqual desc (isoRepDAG1 dag) (isoRepDAG dag)) + +testlistIsoRepDAG = TestList [ + testcaseIsoRepDAG "posetB 3" (hasseDigraph (posetB 3)), + testcaseIsoRepDAG "posetP 3" (hasseDigraph (posetP 3)), + testcaseIsoRepDAG "dual (chainN 5)" (hasseDigraph (dual (chainN 5))), + testcaseIsoRepDAG "antiChainN 5" (hasseDigraph (antichainN 5)), + testcaseIsoRepDAG "posetD 60" (hasseDigraph (posetD 60)), + testcaseIsoRepDAG "dprod (posetB 2) (chainN 3)" (hasseDigraph (dprod (posetB 2) (chainN 3))), + testcaseIsoRepDAG "DG ['a'..'e'] [('a','d'),('a','e'),('b','c'),('b','d'),('d','e')]" + (DG ['a'..'e'] [('a','d'),('a','e'),('b','c'),('b','d'),('d','e')]) + ] +-} diff --git a/Math/Test/TCombinatorics/TIncidenceAlgebra.hs b/Math/Test/TCombinatorics/TIncidenceAlgebra.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCombinatorics/TIncidenceAlgebra.hs @@ -0,0 +1,38 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +module Math.Test.TCombinatorics.TIncidenceAlgebra where + +import Test.HUnit + +import Math.Algebra.Field.Base + +import Math.Combinatorics.Digraph +import Math.Combinatorics.Poset +import Math.Combinatorics.IncidenceAlgebra + +testlistIncidenceAlgebra = TestList [ + testlistMuReference, + testlistMuInverse + ] + + + +-- test that the calculated mu matches reference definition +testcaseMuReference desc poset muref = TestCase $ assertEqual desc muref (muIA poset) + +testlistMuReference = TestList [ + testcaseMuReference "chainN 3" (chainN 3) (muC 3), + testcaseMuReference "posetB 3" (posetB 3) (muB 3), + testcaseMuReference "posetL 3 f3" (posetL 3 f3) (muL 3 f3) + ] + + +-- test that muIA is multiplicative inverse of zetaIA +testcaseMuInverse desc poset = TestCase (assertEqual desc (unitIA poset) (muIA poset * zetaIA poset)) + +testlistMuInverse = TestList [ + testcaseMuInverse "chainN 3" (chainN 3), + testcaseMuInverse "antichainN 3" (antichainN 3), + testcaseMuInverse "posetB 3" (posetB 3), + testcaseMuInverse "posetP 3" (posetP 3) + ] diff --git a/Math/Test/TCombinatorics/TPoset.hs b/Math/Test/TCombinatorics/TPoset.hs new file mode 100644 --- /dev/null +++ b/Math/Test/TCombinatorics/TPoset.hs @@ -0,0 +1,42 @@ +-- Copyright (c) 2011, David Amos. All rights reserved. + +module Math.Test.TCombinatorics.TPoset where + +import Test.HUnit + +import Math.Combinatorics.Digraph +import Math.Combinatorics.Poset + +testlistPoset = TestList [ + testlistIsOrderIsoPositive, + testlistIsOrderIsoNegative, + testlistHasseDigraph + ] + + + +testcaseIsOrderIsoPositive desc p1 p2 = TestCase (assertBool desc (isOrderIso p1 p2)) + +testlistIsOrderIsoPositive = TestList [ + testcaseIsOrderIsoPositive "D 30 ~= D 42" (posetD 30) (posetD 42), + testcaseIsOrderIsoPositive "D 60 ~= D 90" (posetD 30) (posetD 42), + testcaseIsOrderIsoPositive "D 30 ~= B 3" (posetD 30) (posetB 3), + testcaseIsOrderIsoPositive "B 2 ~= 2 * 2" (posetB 2) (dprod (chainN 2) (chainN 2)) + ] + + +testcaseIsOrderIsoNegative desc p1 p2 = TestCase (assertBool desc (not (isOrderIso p1 p2))) + +testlistIsOrderIsoNegative = TestList [ + testcaseIsOrderIsoNegative "D 20 ~/= D 30" (posetD 20) (posetD 30) + ] + + +testcaseHasseDigraph desc poset = TestCase $ assertEqual desc poset (reachabilityPoset $ hasseDigraph poset) + +testlistHasseDigraph = TestList [ + testcaseHasseDigraph "chain 3" (chainN 3), + testcaseHasseDigraph "antichain 3" (antichainN 3), + testcaseHasseDigraph "posetB 3" (posetB 3), + testcaseHasseDigraph "posetB 4" (posetB 4) + ] diff --git a/Math/Test/TGraph.hs b/Math/Test/TGraph.hs --- a/Math/Test/TGraph.hs +++ b/Math/Test/TGraph.hs @@ -31,7 +31,7 @@ && all (uncurry (==)) graphPropsTestsInt graphPropsTestsBool = - [(isConnected nullGraph, True)] ++ + -- [(isConnected nullGraph, True)] ++ [(isConnected (c n), True) | n <- [3..8] ] ++ [(isConnected $ complement $ k n, False) | n <- [3..6] ]