HaskellForMaths 0.3.2 → 0.3.3

raw patch · 27 files changed

+1739/−127 lines, 27 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• HaskellForMaths.cabal +7/−2
• Math/Algebra/Field/Base.hs +2/−0
• Math/Algebra/Field/Extension.hs +9/−5
• Math/Algebra/NonCommutative/TensorAlgebra.hs +4/−1
• Math/Algebras/Commutative.hs +12/−13
• Math/Algebras/Matrix.hs +9/−10
• Math/Algebras/NonCommutative.hs +3/−2
• Math/Algebras/Quaternions.hs +52/−15
• Math/Algebras/Structures.hs +54/−18
• Math/Algebras/TensorAlgebra.hs +299/−12
• Math/Algebras/TensorProduct.hs +21/−2
• Math/Algebras/VectorSpace.hs +58/−8
• Math/Combinatorics/Digraph.hs +288/−0
• Math/Combinatorics/Graph.hs +32/−14
• Math/Combinatorics/GraphAuts.hs +13/−0
• Math/Combinatorics/IncidenceAlgebra.hs +292/−0
• Math/Combinatorics/Poset.hs +282/−0
• Math/QuantumAlgebra/TensorCategory.hs +4/−0
• Math/Test/TAlgebras/TMatrix.hs +32/−0
• Math/Test/TAlgebras/TQuaternions.hs +36/−3
• Math/Test/TAlgebras/TStructures.hs +29/−14
• Math/Test/TAlgebras/TTensorProduct.hs +22/−7
• Math/Test/TAlgebras/TVectorSpace.hs +5/−0
• Math/Test/TCombinatorics/TDigraph.hs +93/−0
• Math/Test/TCombinatorics/TIncidenceAlgebra.hs +38/−0
• Math/Test/TCombinatorics/TPoset.hs +42/−0
• Math/Test/TGraph.hs +1/−1

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

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

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

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

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

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

@@ -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)]

@@ -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)+-}

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

@@ -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 &#x2295; V &#x2295; V&#x2297;V &#x2295; V&#x2297;V&#x2297;V &#x2295; ...+-- 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&#x2297;v - v&#x2297;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&#x2297;v - v&#x2297;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&#x2297;u and sums of products u&#x2297;v + v&#x2297;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&#x2297;u and sums of products u&#x2297;v + v&#x2297;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'+-}++++

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

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

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

@@ -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<y
 data Graph a = G [a] [[a]] deriving (Eq,Ord,Show)
 
+instance Functor Graph where
+    -- |If f is not order-preserving, then you should call nf afterwards
+    fmap f (G vs es) = G (map f vs) (map (map f) es)
+-- could use DeriveFunctor to derive this Functor instance
+
+-- |Convert a graph to normal form. The input is assumed to be a valid graph apart from order
+nf :: Ord a => 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
 

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

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

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

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

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

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

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

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

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

@@ -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')])+    ] +-}

@@ -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)+    ] 

@@ -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)+    ]

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