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

{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, NoMonomorphismRestriction, ScopedTypeVariables, DeriveFunctor #-}

-- |A module defining the following Combinatorial Hopf Algebras, together with coalgebra or Hopf algebra morphisms between them:
--
-- * SSym, the Malvenuto-Reutnenauer Hopf algebra of permutations
--
-- * YSym, the (dual of the) Loday-Ronco Hopf algebra of binary trees
--
-- * QSym, the Hopf algebra of quasi-symmetric functions (having a basis indexed by compositions)
--
-- * Sh, the Shuffle Hopf algebra
module Math.Combinatorics.CombinatorialHopfAlgebra where

-- Sources:

-- Structure of the Malvenuto-Reutenauer Hopf algebra of permutations
-- Marcelo Aguiar and Frank Sottile
-- http://www.math.tamu.edu/~sottile/research/pdf/SSym.pdf

-- Structure of the Loday-Ronco Hopf algebra of trees
-- Marcelo Aguiar and Frank Sottile
-- http://www.math.tamu.edu/~sottile/research/pdf/Loday.pdf

-- Hopf Structures on the Multiplihedra
-- Stefan Forcey, Aaron Lauve and Frank Sottile
-- http://www.math.tamu.edu/~sottile/research/pdf/MSym.pdf


import Data.List as L
import Data.Maybe (fromJust)
import qualified Data.Set as S

import Math.Core.Field
import Math.Core.Utils

import Math.Algebras.VectorSpace hiding (E)
import Math.Algebras.TensorProduct
import Math.Algebras.Structures

import Math.Combinatorics.Poset

-- import Math.Algebra.Group.PermutationGroup
import Math.CommutativeAlgebra.Polynomial


-- SHUFFLE ALGEBRA
-- This is just the tensor algebra, but with shuffle product (and deconcatenation coproduct)

-- |A basis for the shuffle algebra. As a vector space, the shuffle algebra is identical to the tensor algebra.
-- However, we consider a different algebra structure, based on the shuffle product. Together with the
-- deconcatenation coproduct, this leads to a Hopf algebra structure.
newtype Shuffle a = Sh [a] deriving (Eq,Ord,Show)

-- |Construct a basis element of the shuffle algebra
sh :: [a] -> Vect Q (Shuffle a)
sh = return . Sh

shuffles (x:xs) (y:ys) = map (x:) (shuffles xs (y:ys)) ++ map (y:) (shuffles (x:xs) ys)
shuffles xs [] = [xs]
shuffles [] ys = [ys]

instance (Eq k, Num k, Ord a) => Algebra k (Shuffle a) where
    unit x = x *> return (Sh [])
    mult = linear mult'
        where mult' (Sh xs, Sh ys) = sumv [return (Sh zs) | zs <- shuffles xs ys]

deconcatenations xs = zip (inits xs) (tails xs)

instance (Eq k, Num k, Ord a) => Coalgebra k (Shuffle a) where
    counit = unwrap . linear counit' where counit' (Sh xs) = if null xs then 1 else 0
    comult = linear comult'
        where comult' (Sh xs) = sumv [return (Sh us, Sh vs) | (us, vs) <- deconcatenations xs]

instance (Eq k, Num k, Ord a) => Bialgebra k (Shuffle a) where {}

instance (Eq k, Num k, Ord a) => HopfAlgebra k (Shuffle a) where
    antipode = linear (\(Sh xs) -> (-1)^length xs *> return (Sh (reverse xs)))


-- SSYM: PERMUTATIONS
-- (This is permutations considered as combinatorial objects rather than as algebraic objects)

-- Permutations with shifted shuffle product
-- This is the Malvenuto-Reutenauer Hopf algebra of permutations, SSym.
-- It is neither commutative nor co-commutative

-- ssymF xs is the fundamental basis F_xs (Aguiar and Sottile)

-- |The fundamental basis for the Malvenuto-Reutenauer Hopf algebra of permutations, SSym.
newtype SSymF = SSymF [Int] deriving (Eq)

instance Ord SSymF where
    compare (SSymF xs) (SSymF ys) = compare (length xs, xs) (length ys, ys)

instance Show SSymF where
    show (SSymF xs) = "F " ++ show xs

-- |Construct a fundamental basis element in SSym.
-- The list of ints must be a permutation of [1..n], eg [1,2], [3,4,2,1].
ssymF :: [Int] -> Vect Q SSymF
ssymF xs | L.sort xs == [1..n] = return (SSymF xs)
         | otherwise = error "Not a permutation of [1..n]"
         where n = length xs

-- so this is a candidate mult. It is associative and SSymF [] is obviously a left and right identity
-- (need quickcheck properties to prove that)
shiftedConcat (SSymF xs) (SSymF ys) = let k = length xs in SSymF (xs ++ map (+k) ys)

prop_Associative f (x,y,z) = f x (f y z) == f (f x y) z

-- > quickCheck (prop_Associative shiftedConcat)
-- +++ OK, passed 100 tests.


instance (Eq k, Num k) => Algebra k SSymF where
    unit x = x *> return (SSymF [])
    mult = linear mult'
        where mult' (SSymF xs, SSymF ys) =
                  let k = length xs
                  in sumv [return (SSymF zs) | zs <- shuffles xs (map (+k) ys)]


-- standard permutation, also called flattening, eg [6,2,5] -> [3,1,2]
flatten xs = let mapping = zip (L.sort xs) [1..]
        in [y | x <- xs, let Just y = lookup x mapping] 

instance (Eq k, Num k) => Coalgebra k SSymF where
    counit = unwrap . linear counit' where counit' (SSymF xs) = if null xs then 1 else 0
    comult = linear comult'
        where comult' (SSymF xs) = sumv [return (SSymF (st us), SSymF (st vs)) | (us, vs) <- deconcatenations xs]
              st = flatten

instance (Eq k, Num k) => Bialgebra k SSymF where {}

instance (Eq k, Num k) => HopfAlgebra k SSymF where
    antipode = linear antipode'
        where antipode' (SSymF []) = return (SSymF [])
              antipode' x@(SSymF xs) = (negatev . mult . (id `tf` antipode) . removeTerm (SSymF [],x) . comult . return) x
        -- This expression for antipode is derived from mult . (id `tf` antipode) . comult == unit . counit
        -- It's possible because this is a graded, connected Hopf algebra. (connected means the counit is projection onto the grade 0 part)
-- It would be nicer to have an explicit expression for antipode.
{-
instance (Eq k, Num k) => HopfAlgebra k SSymF where
    antipode = linear antipode'
        where antipode' (SSymF v) = sumv [lambda v w *> return (SSymF w) | w <- L.permutations v]
              lambda v w = length [s | s <- powerset [1..n-1],  odd (length s), descentSet (w^-1 * v_s) `isSubset` s]
                         - length [s | s <- powerset [1..n-1],  even (length s), descentSet (w^-1 * v_s) `isSubset` s]
-}

-- |An alternative \"monomial\" basis for the Malvenuto-Reutenauer Hopf algebra of permutations, SSym.
-- This basis is related to the fundamental basis by Mobius inversion in the poset of permutations with the weak order.
newtype SSymM = SSymM [Int] deriving (Eq)

instance Ord SSymM where
    compare (SSymM xs) (SSymM ys) = compare (length xs, xs) (length ys, ys)

instance Show SSymM where
    show (SSymM xs) = "M " ++ show xs

-- |Construct a monomial basis element in SSym.
-- The list of ints must be a permutation of [1..n], eg [1,2], [3,4,2,1].
ssymM :: [Int] -> Vect Q SSymM
ssymM xs | L.sort xs == [1..n] = return (SSymM xs)
         | otherwise = error "Not a permutation of [1..n]"
         where n = length xs

inversions xs = let ixs = zip [1..] xs
                in [(i,j) | ((i,xi),(j,xj)) <- pairs ixs, xi > xj]

weakOrder xs ys = inversions xs `isSubsetAsc` inversions ys

mu (set,po) x y = mu' x y 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

-- |Convert an element of SSym represented in the monomial basis to the fundamental basis
toSSymF :: (Eq k, Num k) => Vect k SSymM -> Vect k SSymF
toSSymF = linear toSSymF'
    where toSSymF' (SSymM u) = sumv [mu (set,po) u v *> return (SSymF v) | v <- set, po u v]
              where set = L.permutations u
                    po = weakOrder

-- |Convert an element of SSym represented in the fundamental basis to the monomial basis
toSSymM :: (Eq k, Num k) => Vect k SSymF -> Vect k SSymM
toSSymM = linear toSSymM'
    where toSSymM' (SSymF u) = sumv [return (SSymM v) | v <- set, po u v]
              where set = L.permutations u
                    po = weakOrder

-- (p,q)-shuffles: permutations of [1..p+q] having at most one descent, at position p
-- denoted S^{(p,q)} in Aguiar&Sottile
-- (Grassmannian permutations?)
-- pqShuffles p q = [u++v | u <- combinationsOf p [1..n], let v = [1..n] `diffAsc` u] where n = p+q

-- The inverse of a (p,q)-shuffle.
-- The special form of (p,q)-shuffles makes an O(n) algorithm possible
-- pqInverse :: Int -> Int -> [Int] -> [Int]
{-
-- incorrect
pqInverse p q xs = pqInverse' [1..p] [p+1..p+q] xs
    where pqInverse' (l:ls) (r:rs) (x:xs) =
              if x <= p then l : pqInverse' ls (r:rs) xs else r : pqInverse' (l:ls) rs xs
          pqInverse' ls rs _ = ls ++ rs -- one of them is null
-}
-- pqInverseShuffles p q = shuffles [1..p] [p+1..p+q]


instance (Eq k, Num k) => Algebra k SSymM where
    unit x = x *> return (SSymM [])
    mult = toSSymM . mult . (toSSymF `tf` toSSymF)

{-
mult2 = linear mult'
    where mult' (SSymM u, SSymM v) = sumv [alpha u v w *> return (SSymM w) | w <- L.permutations [1..p+q] ]
                                     where p = length u; q = length v

alpha u v w = length [z | z <- pqInverseShuffles p q, let uv = shiftedConcat u v,
                          uv * z `weakOrder` w, u and v are maximal, ie no transposition of adjacents in either also works]
    where p = length u
          q = length v
-- so we need to define (*) for permutations in row form
-}

instance (Eq k, Num k) => Coalgebra k SSymM where
    counit = unwrap . linear counit' where counit' (SSymM xs) = if null xs then 1 else 0
    -- comult = (toSSymM `tf` toSSymM) . comult . toSSymF
    comult = linear comult'
        where comult' (SSymM xs) = sumv [return (SSymM (flatten ys), SSymM (flatten zs))
                                        | (ys,zs) <- deconcatenations xs,
                                          minimum (infinity:ys) > maximum (0:zs)] -- ie deconcatenations at a global descent
              infinity = maxBound :: Int

instance (Eq k, Num k) => Bialgebra k SSymM where {}

instance (Eq k, Num k) => HopfAlgebra k SSymM where
    antipode = toSSymM . antipode . toSSymF


-- YSYM: PLANAR BINARY TREES
-- These are really rooted planar binary trees.
-- It's because they're planar that we can distinguish left and right child branches.
-- (Non-planar would be if we considered trees where left and right children are swapped relative to one another as the same tree)
-- It is neither commutative nor co-commutative

-- |A type for (rooted) planar binary trees. The basis elements of the Loday-Ronco Hopf algebra are indexed by these.
--
-- Although the trees are labelled, we're really only interested in the shapes of the trees, and hence in the type PBT ().
-- The Algebra, Coalgebra and HopfAlgebra instances all ignore the labels.
-- However, it is convenient to allow labels, as they can be useful for seeing what is going on, and they also make it possible
-- to define various ways to create trees from lists of labels.
data PBT a = T (PBT a) a (PBT a) | E deriving (Eq, Show, Functor)

instance Ord a => Ord (PBT a) where
    compare u v = compare (shapeSignature u, prefix u) (shapeSignature v, prefix v)

-- |The fundamental basis for (the dual of) the Loday-Ronco Hopf algebra of binary trees, YSym.
newtype YSymF a = YSymF (PBT a) deriving (Eq, Ord, Functor)

instance Show a => Show (YSymF a) where
    show (YSymF t) = "F(" ++ show t ++ ")"

-- |Construct the element of YSym in the fundamental basis indexed by the given tree
ysymF :: PBT a -> Vect Q (YSymF a)
ysymF t = return (YSymF t)

{-
depth (T l x r) = 1 + max (depth l) (depth r)
depth E = 0
-}
nodecount (T l x r) = 1 + nodecount l + nodecount r
nodecount E = 0

-- in fact leafcount t = 1 + nodecount t (easiest to see with a picture)
leafcount (T l x r) = leafcount l + leafcount r
leafcount E = 1

prefix E = []
prefix (T l x r) = x : prefix l ++ prefix r

-- The shape signature uniquely identifies the shape of a tree.
-- Trees with distinct shapes have distinct signatures.
-- In addition, if sorting on shapeSignature, smaller trees sort before larger trees,
-- and leftward leaning trees sort before rightward leaning trees
shapeSignature t = shapeSignature' (nodeCountTree t)
    where shapeSignature' E = [0] -- not [], otherwise we can't distinguish T (T E () E) () E from T E () (T E () E)
          shapeSignature' (T l x r) = x : shapeSignature' r ++ shapeSignature' l

nodeCountTree E = E
nodeCountTree (T l _ r) = T l' n r'
    where l' = nodeCountTree l
          r' = nodeCountTree r
          n = 1 + (case l' of E -> 0; T _ lc _ -> lc) + (case r' of E -> 0; T _ rc _ -> rc)

leafCountTree E = E
leafCountTree (T l _ r) = T l' n r'
    where l' = leafCountTree l
          r' = leafCountTree r
          n = (case l' of E -> 1; T _ lc _ -> lc) + (case r' of E -> 1; T _ rc _ -> rc)

-- A tree that counts nodes in left and right subtrees
lrCountTree E = E
lrCountTree (T l _ r) = T l' (lc,rc) r'
    where l' = lrCountTree l
          r' = lrCountTree r
          lc = case l' of E -> 0; T _ (llc,lrc) _ -> 1 + llc + lrc
          rc = case r' of E -> 0; T _ (rlc,rrc) _ -> 1 + rlc + rrc

shape :: PBT a -> PBT ()
shape t = fmap (\_ -> ()) t

-- label the nodes of a tree in infix order while preserving its shape
numbered t = numbered' 1 t
    where numbered' _ E = E
          numbered' i (T l x r) = let k = nodecount l in T (numbered' i l) (i+k) (numbered' (i+k+1) r)
-- could also pair the numbers with the input labels


splits E = [(E,E)]
splits (T l x r) = [(u, T v x r) | (u,v) <- splits l] ++ [(T l x u, v) | (u,v) <- splits r]

instance (Eq k, Num k, Ord a) => Coalgebra k (YSymF a) where
    counit = unwrap . linear counit' where counit' (YSymF E) = 1; counit' (YSymF (T _ _ _)) = 0
    comult = linear comult'
        where comult' (YSymF t) = sumv [return (YSymF u, YSymF v) | (u,v) <- splits t]
              -- using sumv rather than sum to avoid requiring Show a
    -- so again this is a kind of deconcatenation coproduct

multisplits 1 t = [ [t] ]
multisplits 2 t = [ [u,v] | (u,v) <- splits t ]
multisplits n t = [ u:ws | (u,v) <- splits t, ws <- multisplits (n-1) v ]

graft [t] E = t
graft ts (T l x r) = let (ls,rs) = splitAt (leafcount l) ts
                     in T (graft ls l) x (graft rs r)

instance (Eq k, Num k, Ord a) => Algebra k (YSymF a) where
    unit x = x *> return (YSymF E)
    mult = linear mult'
        where mult' (YSymF t, YSymF u) = sumv [return (YSymF (graft ts u)) | ts <- multisplits (leafcount u) t]
              -- using sumv rather than sum to avoid requiring Show a

instance (Eq k, Num k, Ord a) => Bialgebra k (YSymF a) where {}

instance (Eq k, Num k, Ord a) => HopfAlgebra k (YSymF a) where
    antipode = linear antipode'
        where antipode' (YSymF E) = return (YSymF E)
              antipode' x = (negatev . mult . (id `tf` antipode) . removeTerm (YSymF E,x) . comult . return) x


-- |An alternative "monomial" basis for (the dual of) the Loday-Ronco Hopf algebra of binary trees, YSym.
newtype YSymM = YSymM (PBT ()) deriving (Eq, Ord)

instance Show YSymM where
    show (YSymM t) = "M(" ++ show t ++ ")"

-- |Construct the element of YSym in the monomial basis indexed by the given tree
ysymM :: PBT () -> Vect Q YSymM
ysymM t = return (YSymM t)


trees 0 = [E]
trees n = [T l () r | i <- [0..n-1], l <- trees (n-1-i), r <- trees i]

-- |The covering relation for the Tamari partial order on binary trees
covers E = []
covers (T t@(T u x v) y w) = [T t' y w | t' <- covers t]
                          ++ [T t y w' | w' <- covers w]
                          ++ [T u y (T v x w)]
                          -- Note that this preserves the descending property, and hence the bijection with permutations
                          -- If we were to swap x and y, we would preserve the binary search tree property instead (if our trees had it)
covers (T E x u) = [T E x u' | u' <- covers u]  

-- |The up-set of a binary tree in the Tamari partial order
tamariUpSet t = upSet' [] [t]
    where upSet' interior boundary =
              if null boundary
              then interior
              else let interior' = setUnionAsc interior boundary
                       boundary' = toSet $ concatMap covers boundary
                   in upSet' interior' boundary'

-- tamariOrder1 u v = v `elem` upSet u

tamariOrder u v = weakOrder (minPerm u) (minPerm v)
-- It should be possible to unpack this to be a statement purely about trees, but probably not worth

-- |Convert an element of YSym represented in the monomial basis to the fundamental basis
toYSymF :: (Eq k, Num k) => Vect k YSymM -> Vect k (YSymF ())
toYSymF = linear toYSymF'
    where toYSymF' (YSymM t) = sumv [mu (set,po) t s *> return (YSymF s) | s <- set]
              where po = tamariOrder
                    set = tamariUpSet t -- [s | s <- trees (nodecount t), t `tamariOrder` s]

-- |Convert an element of YSym represented in the fundamental basis to the monomial basis
toYSymM :: (Eq k, Num k) => Vect k (YSymF ()) -> Vect k YSymM
toYSymM = linear toYSymM'
    where toYSymM' (YSymF t) = sumv [return (YSymM s) | s <- tamariUpSet t]
                            -- sumv [return (YSymM s) | s <- trees (nodecount t), t `tamariOrder` s]


instance (Eq k, Num k) => Algebra k YSymM where
    unit x = x *> return (YSymM E)
    mult = toYSymM . mult . (toYSymF `tf` toYSymF)

instance (Eq k, Num k) => Coalgebra k YSymM where
    counit = unwrap . linear counit' where counit' (YSymM E) = 1; counit' (YSymM (T _ _ _)) = 0
    -- comult = (toYSymM `tf` toYSymM) . comult . toYSymF
    comult = linear comult'
        where comult' (YSymM t) = sumv [return (YSymM r, YSymM s) | (rs,ss) <- deconcatenations (underDecomposition t),
                                                                    let r = foldl under E rs, let s = foldl under E ss]


instance (Eq k, Num k) => Bialgebra k YSymM where {}

instance (Eq k, Num k) => HopfAlgebra k YSymM where
    antipode = toYSymM . antipode . toYSymF 


-- QSYM: QUASI-SYMMETRIC FUNCTIONS
-- The following is the Hopf algebra QSym of quasi-symmetric functions
-- using the monomial basis (indexed by compositions)

-- compositions in ascending order
-- might be better to use bfs to get length order
-- |List the compositions of an integer n. For example, the compositions of 4 are [[1,1,1,1],[1,1,2],[1,2,1],[1,3],[2,1,1],[2,2],[3,1],[4]]
compositions :: Int -> [[Int]]
compositions 0 = [[]]
compositions n = [i:is | i <- [1..n], is <- compositions (n-i)]

-- can retrieve subsets of [1..n-1] from compositions n as follows
-- > map (tail . scanl (+) 0) (map init $ compositions 4)
-- [[],[3],[2],[2,3],[1],[1,3],[1,2],[1,2,3]]

-- quasi shuffles of two compositions
quasiShuffles :: [Int] -> [Int] -> [[Int]]
quasiShuffles (x:xs) (y:ys) = map (x:) (quasiShuffles xs (y:ys)) ++
                              map (y:) (quasiShuffles (x:xs) ys) ++
                              map ((x+y):) (quasiShuffles xs ys)
quasiShuffles xs [] = [xs]
quasiShuffles [] ys = [ys]


-- |A type for the monomial basis for the quasi-symmetric functions, indexed by compositions.
newtype QSymM = QSymM [Int] deriving (Eq)

instance Ord QSymM where
    compare (QSymM xs) (QSymM ys) = compare (length xs, xs) (length ys, ys)

instance Show QSymM where
    show (QSymM xs) = "M " ++ show xs

-- |Construct the element of QSym in the monomial basis indexed by the given composition
qsymM :: [Int] -> Vect Q QSymM
qsymM = return . QSymM

instance (Eq k, Num k) => Algebra k QSymM where
    unit x = x *> return (QSymM [])
    mult = linear mult'
        where mult' (QSymM alpha, QSymM beta) = sum [return (QSymM gamma) | gamma <- quasiShuffles alpha beta]

instance (Eq k, Num k) => Coalgebra k QSymM where
    counit = unwrap . linear counit' where counit' (QSymM alpha) = if null alpha then 1 else 0
    comult = linear comult' where
        comult' (QSymM gamma) = sum [return (QSymM alpha, QSymM beta) | (alpha,beta) <- deconcatenations gamma]

instance (Eq k, Num k) => Bialgebra k QSymM where {}

instance (Eq k, Num k) => HopfAlgebra k QSymM where
    antipode = linear antipode' where
        antipode' (QSymM alpha) = (-1)^length alpha * sum [return (QSymM (reverse beta)) | beta <- coarsenings alpha]

coarsenings (x1:x2:xs) = coarsenings ((x1+x2):xs) ++ map (x1:) (coarsenings (x2:xs))
coarsenings xs = [xs] -- for xs a singleton or null

refinements (x:xs) = [y++ys | y <- compositions x, ys <- refinements xs]
refinements [] = [[]]


newtype QSymF = QSymF [Int] deriving (Eq)

instance Ord QSymF where
    compare (QSymF xs) (QSymF ys) = compare (length xs, xs) (length ys, ys)

instance Show QSymF where
    show (QSymF xs) = "F " ++ show xs

-- |Construct the element of QSym in the fundamental basis indexed by the given composition
qsymF :: [Int] -> Vect Q QSymF
qsymF = return . QSymF

-- |Convert an element of QSym represented in the monomial basis to the fundamental basis
toQSymF :: (Eq k, Num k) => Vect k QSymM -> Vect k QSymF
toQSymF = linear toQSymF'
    where toQSymF' (QSymM alpha) = sumv [(-1) ^ (length beta - length alpha) *> return (QSymF beta) | beta <- refinements alpha]

-- |Convert an element of QSym represented in the fundamental basis to the monomial basis
toQSymM :: (Eq k, Num k) => Vect k QSymF -> Vect k QSymM
toQSymM = linear toQSymM'
    where toQSymM' (QSymF alpha) = sumv [return (QSymM beta) | beta <- refinements alpha] -- ie beta <- up-set of alpha

instance (Eq k, Num k) => Algebra k QSymF where
    unit x = x *> return (QSymF [])
    mult = toQSymF . mult . (toQSymM `tf` toQSymM)

instance (Eq k, Num k) => Coalgebra k QSymF where
    counit = unwrap . linear counit' where counit' (QSymF xs) = if null xs then 1 else 0
    comult = (toQSymF `tf` toQSymF) . comult . toQSymM

instance (Eq k, Num k) => Bialgebra k QSymF where {}

instance (Eq k, Num k) => HopfAlgebra k QSymF where
    antipode = toQSymF . antipode . toQSymM


-- QUASI-SYMMETRIC POLYNOMIALS

-- the above induces Hopf algebra structure on quasi-symmetric functions via
-- m_alpha -> sum [product (zipWith (^) (map x_ is) alpha | is <- combinationsOf k [] ] where k = length alpha

xvars n = [glexvar ("x" ++ show i) | i <- [1..n] ]

-- compare with Reynolds operator
-- so a basis for quasi-symmetric functions over xvars n consists of [quasiSymM xs is | m <- [0..], is <- compositions m]
quasiSymM xs is = sum [product (zipWith (^) xs' is) | xs' <- combinationsOf r xs]
    where r = length is


-- MAPS BETWEEN (POSETS AND) HOPF ALGEBRAS

-- A descending tree is one in which a child is always less than a parent.
descendingTree [] = E
descendingTree [x] = T E x E
descendingTree xs = T l x r
    where x = maximum xs
          (ls,_:rs) = L.break (== x) xs
          l = descendingTree ls
          r = descendingTree rs
-- This is a bijection from permutations to "ordered trees".
-- It is order-preserving on trees with the same nodecount.
-- We can recover the permutation by reading the node labels in infix order.
-- This is the map called lambda in Loday.pdf


-- |A Hopf algebra morphism from SSymF to YSymF
descendingTreeMap :: (Eq k, Num k) => Vect k SSymF -> Vect k (YSymF ())
descendingTreeMap = nf . fmap (YSymF . shape .  descendingTree')
    where descendingTree' (SSymF xs) = descendingTree xs
-- This is the map called Lambda in Loday.pdf, or tau in MSym.pdf
-- It is an algebra morphism.

-- One of the ideas in the MSym paper is to look at the intermediate result (fmap descendingTree' x),
-- which is an "ordered tree", and consider the map as factored through this

-- The map is surjective but not injective. The fibers tau^-1(t) are intervals in the weak order on permutations

-- "inverse" for descendingTree
-- These are the maps called gamma in Loday.pdf
minPerm t = minPerm' (lrCountTree t)
    where minPerm' E = []
          minPerm' (T l (lc,rc) r) = minPerm' l ++ [lc+rc+1] ++ map (+lc) (minPerm' r)

maxPerm t = maxPerm' (lrCountTree t)
    where maxPerm' E = []
          maxPerm' (T l (lc,rc) r) = map (+rc) (maxPerm' l) ++ [lc+rc+1] ++ maxPerm' r


-- The composition of [1..n] obtained by treating each left-facing leaf as a cut
-- Specifically, we visit the nodes in infix order, cutting after a node if it does not have an E as its right child
-- This is the map called L in Loday.pdf
leftLeafComposition E = []
leftLeafComposition t = cuts $ tail $ leftLeafs t
    where leftLeafs (T l x E) = leftLeafs l ++ [False]
          leftLeafs (T l x r) = leftLeafs l ++ leftLeafs r
          leftLeafs E = [True]
          cuts bs = case break id bs of
                    (ls,r:rs) -> (length ls + 1) : cuts rs
                    (ls,[]) -> [length ls]

leftLeafComposition' (YSymF t) = QSymF (leftLeafComposition t)

-- |A Hopf algebra morphism from YSymF to QSymF
leftLeafCompositionMap :: (Eq k, Num k) => Vect k (YSymF a) -> Vect k QSymF
leftLeafCompositionMap = nf . fmap leftLeafComposition'


-- The descent set of a permutation is [i | x_i > x_i+1], where we start the indexing from 1
descents [] = []
descents xs = map (+1) $ L.elemIndices True $ zipWith (>) xs (tail xs)

-- The composition of [1..n] obtained by treating each descent as a cut
descentComposition [] = []
descentComposition xs = dc $ zipWith (>) xs (tail xs) ++ [False]
    where dc bs = case break id bs of
                  (ls,r:rs) -> (length ls + 1) : dc rs
                  (ls,[]) -> [length ls]

-- |A Hopf algebra morphism from SSymF to QSymF
descentMap :: (Eq k, Num k) => Vect k SSymF -> Vect k QSymF
descentMap = nf . fmap (\(SSymF xs) -> QSymF (descentComposition xs))
-- descentMap == leftLeafCompositionMap . descendingTreeMap

underComposition (QSymF ps) = foldr under (SSymF []) [SSymF [1..p] | p <- ps]
    where under (SSymF xs) (SSymF ys) = let q = length ys
                                            zs = map (+q) xs ++ ys -- so it has a global descent at the split
                                        in SSymF zs
-- This is a poset morphism (indeed, it forms a Galois connection with descentComposition)
-- but it does not extend to a Hopf algebra morphism.
-- (It does extend to a coalgebra morphism.)
-- (It is picking the maximum permutation having a given descent composition,
-- so there's an element of arbitrariness to it.)
-- This is the map called Z (Zeta?) in Loday.pdf

{-
-- This is O(n^2), whereas an O(n) implementation should be possible
-- Also, we would really like the associated composition (obtained by treating each global descent as a cut)?
globalDescents xs = globalDescents' 0 [] xs
    where globalDescents' i ls (r:rs) = (if minimum (infinity:ls) > maximum (0:r:rs) then [i] else [])
                                     ++ globalDescents' (i+1) (r:ls) rs
          globalDescents' n _ [] = [n]
          infinity = maxBound :: Int
-- The idea is that this leads to a map from SSymM to QSymM

globalDescentComposition [] = []
globalDescentComposition (x:xs) = globalDescents' 1 x xs
    where globalDescents' i minl (r:rs) = if minl > maximum (r:rs)
                                          then i : globalDescents' 1 r rs
                                          else globalDescents' (i+1) r rs
          globalDescents' i _ [] = [i]

globalDescentMap :: (Eq k, Num k) => Vect k SSymM -> Vect k QSymM
globalDescentMap = nf . fmap (\(SSymM xs) -> QSymM (globalDescentComposition xs))
-}

-- A multiplication operation on trees
-- (Connected with their being cofree)
-- (intended to be used as infix)
under E t = t
under (T l x r) t = T l x (under r t)

isUnderIrreducible (T l x E) = True
isUnderIrreducible _ = False

underDecomposition (T l x r) = T l x E : underDecomposition r
underDecomposition E = []


-- GHC7.4.1 doesn't like the following type signature - a bug.
-- ysymmToSh :: (Eq k, Num k) => Vect k (YSymM) => Vect k (Shuffle (PBT ()))
ysymmToSh = fmap ysymmToSh'
    where ysymmToSh' (YSymM t) = Sh (underDecomposition t)
-- This is a coalgebra morphism (but not an algebra morphism)
-- It shows that YSym is co-free
{-
-- This one not working yet - perhaps it needs an nf, or to go via S/YSymF, or ...
ssymmToSh = nf . fmap ssymmToSh'
    where ssymmToSh' (SSymM xs) = (Sh . underDecomposition . shape . descendingTree) xs
-}