-- Copyright 2009 Mikael Vejdemo Johansson <mik@stanford.edu>
-- Released under a BSD license

-- | The module 'OperadGB' carries the implementations of the Buchberger algorithm and most utility functions 
-- related to this.

module Math.Operad.OperadGB where

import Prelude hiding (mapM, sequence)
import Data.List (sort, sortBy, findIndex, nub, (\\), permutations)
import Data.Ord
import Data.Foldable (foldMap, Foldable)
import Control.Monad hiding (mapM)
import Data.Maybe

#if defined USE_MAPOPERAD
import Math.Operad.MapOperad
#elif defined USE_POLYBAG
import Math.Operad.PolyBag
#else
import Math.Operad.MapOperad
#endif

import Math.Operad.OrderedTree

--import Debug.Trace

-- * Fundamental data types and instances

-- | The number of internal vertices of a tree.
operationDegree :: (Ord a, Show a) => DecoratedTree a -> Int
operationDegree (DTLeaf _) = 0
operationDegree vertex = 1 + sum (map operationDegree (subTrees vertex))

-- | A list of operation degrees occurring in the terms of the operad element
operationDegrees :: (Ord a, Show a, TreeOrdering t, Num n) => OperadElement a n t -> [Int]
operationDegrees op = foldMonomials (\(k,_) lst -> lst ++ [(operationDegree . dt $ k)]) op

-- | The maximal operation degree of an operadic element
maxOperationDegree :: (Ord a, Show a, TreeOrdering t, Num n) => OperadElement a n t -> Int
maxOperationDegree = maximum . operationDegrees

-- | Check that an element of a free operad is homogenous
isHomogenous :: (Ord a, Show a, TreeOrdering t, Num n) => OperadElement a n t -> Bool
isHomogenous m = let 
    trees = getTrees m
--    equalityCheck :: OrderedTree a t -> OrderedTree a t -> Bool
    equalityCheck s t = arityDegree (dt s) == arityDegree (dt t) &&
                        operationDegree (dt s) == operationDegree (dt t)
  in and $ zipWith (equalityCheck) trees (tail trees)

-- * Free operad

-- ** Operadic compositions

-- | Composition in the shuffle operad
shuffleCompose :: (Ord a, Show a) => Int -> Shuffle -> DecoratedTree a -> DecoratedTree a -> DecoratedTree a
shuffleCompose i sh s t | not (isPermutation sh) = error "shuffleCompose: sh needs to be a permutation\n"
                        | (nLeaves s) + (nLeaves t) - 1 /= length sh = 
                            error $ "Permutation permutes the wrong number of things:" ++ show i ++ " " ++ show sh ++ " " ++ show s ++ " " ++ show t ++ "\n"
                        | not (isShuffleIPQ sh i (nLeaves t-1)) = 
                            error $ "Need a correct pointed shuffle permutation!\n" ++ 
                                  show sh ++ " is not in Sh" ++ show i ++ 
                                           "(" ++ show (nLeaves t-1) ++ "," ++ show (nLeaves s-i) ++ ")\n"
                        | otherwise = symmetricCompose i sh s t

-- | Composition in the non-symmetric operad. We compose s o_i t. 
nsCompose :: (Ord a, Show a) => Int -> DecoratedTree a -> DecoratedTree a -> DecoratedTree a
nsCompose i s t = if i-1 > nLeaves s then error "Composition point too large"
                  else let
                      pS = rePackLabels s
                      lookupList = zip (leafOrder s) (leafOrder pS)
                      idx = if isNothing newI then error "Index not in tree" else fromJust newI where newI = lookup i lookupList
                      trees = map DTLeaf [1..nLeaves s]
                      newTrees = take (idx-1) trees ++ [t] ++ drop idx trees
                 in
                   if length newTrees /= nLeaves s then error "Shouldn't happen" 
                   else 
                     nsComposeAll s newTrees

-- | Composition in the symmetric operad
symmetricCompose :: (Ord a, Show a) => Int -> Shuffle -> DecoratedTree a -> DecoratedTree a -> DecoratedTree a
symmetricCompose i sh s t = if not (isPermutation sh) then error "symmetricCompose: sh needs to be a permutation\n"
                            else if (nLeaves s) + (nLeaves t) - 1 /= length sh then error "Permutation permutes the wrong number of things.\n"
                            else fmap ((sh!!) . (subtract 1)) $  nsCompose i s t

-- | Non-symmetric composition in the g(s;t1,...,tk) style. 
nsComposeAll :: (Ord a, Show a) => DecoratedTree a -> [DecoratedTree a] -> DecoratedTree a
nsComposeAll s trees = if nLeaves s /= length trees then error "NS: Need as many trees as leaves\n"
                        else if length trees == 0 then s
                        else let 
                            treesArities = map nLeaves trees
                            packedTrees = map rePackLabels trees
                            offSets = (0:) $ scanl1 (+) treesArities
                            newTrees = zipWith (\t n -> fmap (+n) t) packedTrees offSets
                          in
                            rePackLabels $ glueTrees $ fmap ((newTrees!!) . (subtract 1)) $ rePackLabels s
                        
-- | Verification for a shuffle used for the g(s;t1,..,tk) style composition in the shuffle operad.
checkShuffleAll :: Shuffle -> [Int] -> Bool
checkShuffleAll sh blockL = let
      checkOrders :: Shuffle -> [Int] -> Bool
      checkOrders [] _ = True
      checkOrders _ [] = True
      checkOrders restSh restBlock = (isSorted (take (head restBlock) restSh)) && 
                                     (length restSh <= head restBlock || 
                                      (head restSh) < (head (drop (head restBlock) restSh))) &&
                                     checkOrders (drop (head restBlock) restSh) (tail restBlock)
    in
      sum blockL == length sh &&
      checkOrders sh blockL
                            
-- | Sanity check for permutations. 
isPermutation :: Shuffle -> Bool
isPermutation sh = and ((zipWith (==) [1..]) (sort sh))

-- | Shuffle composition in the g(s;t1,...,tk) style. 
shuffleComposeAll :: (Ord a, Show a) => Shuffle -> DecoratedTree a -> [DecoratedTree a] -> DecoratedTree a
shuffleComposeAll sh s trees = if nLeaves s /= length trees then error "Shuffle: Need as many trees as leaves\n"
                               else if sum (map nLeaves trees) /= length sh then error "Permutation permutes the wrong number of things.\n"
                               else if not (isPermutation sh) then error "shuffleComposeAll: sh needs to be a permutation\n"
                               else if length trees == 0 then s
                                    else if not (checkShuffleAll sh (map nLeaves trees)) then error "Bad shuffle"
                                         else let
      newTree = nsComposeAll s trees
   in
            fmap ((sh!!) . (subtract 1)) newTree

-- | Symmetric composition in the g(s;t1,...,tk) style. 
symmetricComposeAll :: (Ord a, Show a) => Shuffle -> DecoratedTree a -> [DecoratedTree a] -> DecoratedTree a
symmetricComposeAll sh s trees = if nLeaves s /= length trees then error "Symm: Need as many trees as leaves\n"
                               else if sum (map nLeaves trees) /= length sh then error "Permutation permutes the wrong number of things.\n"
                               else if not (isPermutation sh) then error "sh needs to be a permutation"
                               else if length trees == 0 then s
                                         else let
      newTree = nsComposeAll s trees
   in
            fmap ((sh!!) . (subtract 1)) newTree




-- ** Divisibility among trees

-- | Data type to move the results of finding an embedding point for a subtree in a larger tree
-- around. The tree is guaranteed to have exactly one corolla tagged Nothing, the subtrees on top of
-- that corolla sorted by minimal covering leaf in the original setting, and the shuffle carried around
-- is guaranteed to restore the original leaf labels before the search.
type Embedding a = DecoratedTree (Maybe a)

-- | Returns True if there is a subtree of @t@ isomorphic to s, respecting leaf orders. 
divides :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> Bool
divides s t = not . null $ findAllEmbeddings s t

-- | Finds all ways to embed s into t respecting leaf orders.
findAllEmbeddings :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> [Embedding a]
findAllEmbeddings _ (DTLeaf _) = []
findAllEmbeddings s t = let 
    rootFind = maybeToList $ findRootedEmbedding s t
    subFinds = map (findAllEmbeddings s) (subTrees t)
    reGlue (i, ems) = let
        glueTree tree = 
            (DTVertex 
             (Just $ vertexType t) 
             (take (i-1) (map toJustTree $ subTrees t) ++ [tree] ++ drop i (map toJustTree $ subTrees t)))
      in map glueTree ems
  in rootFind ++ concatMap reGlue (zip [1..] subFinds)
    
-- | Finds all ways to embed s into t, respecting leaf orders and mapping the root of s to the root of t.
findRootedEmbedding :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> Maybe (Embedding a)
findRootedEmbedding (DTLeaf _) t = Just (DTVertex Nothing [toJustTree t])
findRootedEmbedding (DTVertex _ _) (DTLeaf _) = Nothing
findRootedEmbedding s t = do
  guard $ vertexArity s == vertexArity t
  guard $ vertexType s == vertexType t
  guard $ equivalentOrders (map minimalLeaf (subTrees s)) (map minimalLeaf (subTrees t))
  let mTreeFinds = zipWith findRootedEmbedding (subTrees s) (subTrees t)
  guard $ all isJust mTreeFinds
  let treeFinds = map fromJust mTreeFinds
  guard $ all (isNothing . vertexType) treeFinds
  guard $ equivalentOrders (leafOrder s) (concatMap (map minimalLeaf . subTrees) treeFinds)
  return $ DTVertex Nothing (sortBy (comparing minimalLeaf) (concatMap subTrees treeFinds))

-- | Finds a large common divisor of two trees, such that it embeds into both trees, mapping its root 
-- to the roots of the trees, respectively. 
findRootedDecoratedGCD :: (Ord a, Show a) => 
         DecoratedTree a -> DecoratedTree a -> Maybe (PreDecoratedTree a (DecoratedTree a,DecoratedTree a))
findRootedDecoratedGCD (DTLeaf k) t = Just $ DTLeaf (DTLeaf k, t)
findRootedDecoratedGCD s (DTLeaf k) = Just $ DTLeaf (s, DTLeaf k)
findRootedDecoratedGCD s t = do
  guard $ vertexArity s == vertexArity t
  guard $ vertexType s == vertexType t
  let mrdGCDs = zipWith findRootedDecoratedGCD (subTrees s) (subTrees t)
  guard $ all isJust mrdGCDs
  let rdGCDs = map fromJust mrdGCDs
  return $ DTVertex (vertexType s) rdGCDs

-- | Finds all small common multiples of trees s and t, under the assumption that the common multiples shares
-- root with both trees.
findRootedLCM :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> [DecoratedTree a]
findRootedLCM s t = filter (\tree -> divides s tree && divides t tree) $
    if operationDegree s < operationDegree t then findRootedLCM t s 
    else 
        do
  let mrdGCD = findRootedDecoratedGCD s t
  guard $ isJust mrdGCD 
  let rdGCD = fromJust mrdGCD
      leafDecorations = foldMap (:[]) rdGCD
      rebuildRecipe = reverse . sortBy (comparing fst) $ filter (isLeaf . fst) leafDecorations
  accumulateTrees rebuildRecipe [s]

-- | Internal utility function. Reassembles a tree according to a "building recipe", and gives the orbit
-- of the resulting tree under the symmetric group action back.
accumulateTrees :: (Ord a, Show a) => 
                   [(DecoratedTree a, DecoratedTree a)] -> [DecoratedTree a] -> [DecoratedTree a]
accumulateTrees [] partialTrees = partialTrees
accumulateTrees ((aLeaf,tree):rs) partialTrees = 
    if not $ isLeaf aLeaf then error "Should have a leaf" else
        let
            newTrees = do
              partialTree <- partialTrees
              let idx = minimalLeaf aLeaf
                  newTree = rePackLabels tree
                  packedPartialTree = rePackLabels partialTree
                  lookupList = zip (leafOrder partialTree) (leafOrder packedPartialTree)
                  i = fromJust $ lookup idx lookupList
              return $ nsCompose i packedPartialTree newTree
          in
            do
              t <- accumulateTrees rs newTrees
              p <- permutations [1..nLeaves t]
              let returnTree = relabelLeaves t p
              guard $ planarTree returnTree
              return $ returnTree

-- | Checks a tree for planarity.
planarTree :: (Ord a, Show a) => DecoratedTree a -> Bool
planarTree (DTLeaf _) = True
planarTree (DTVertex _ subs) = all planarTree subs && isSorted (map minimalLeaf subs)

-- | Finds all small common multiples of s and t such that t glues into s from above.
findSmallLCM :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> [DecoratedTree a]
findSmallLCM (DTLeaf _) _ = []
findSmallLCM _ (DTLeaf _) = []
findSmallLCM s t = nub $ filter (divides s) $ filter (isJust . findRootedEmbedding t) $ do
  -- find rLCMs of s and t.
  -- find LCMs of all subtrees of s with t
  -- for those, reglue the rest of t
  let rootedLCMs = findRootedLCM s t
      childLCMs = map (findSmallLCM s) (subTrees t)
      reGlue (i,ems) = if i > length (subTrees t) then error "Too high composition point, findSmallLCM:reGlue" else let
                    template = rePackLabels $  
                               DTVertex 
                               (vertexType t) 
                               (take (i-1) (subTrees t) ++ [leaf (minimalLeaf (subTrees t !! (i-1)))] ++ drop i (subTrees t))
                  in concatMap (\emt -> accumulateTrees [(leaf i,emt)] [template]) ems
      zippedChildLCMs = zip [1..] childLCMs
  rootedLCMs ++ (concatMap reGlue zippedChildLCMs)

-- | Finds all small common multiples of s and t.
findAllLCM :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> [DecoratedTree a]
findAllLCM s t = (findSmallLCM s t) ++ (findSmallLCM t s)

-- | Relabels a tree in the right order, but with entries from [1..]
rePackLabels :: (Ord a, Show a, Ord b) => PreDecoratedTree a b -> DecoratedTree a
rePackLabels tree = fmap (fromJust . (flip lookup (zip (sort (foldMap (:[]) tree)) [1..]))) tree

-- | Removes vertex type encapsulations.
fromJustTree :: (Ord a, Show a) => DecoratedTree (Maybe a) -> Maybe (DecoratedTree a)
fromJustTree (DTLeaf k) = Just (DTLeaf k)
fromJustTree (DTVertex Nothing _) = Nothing
fromJustTree (DTVertex (Just l) sts) = let
    newsts = map fromJustTree sts
  in
    if all isJust newsts then Just $ DTVertex l (map fromJust newsts)
    else Nothing

-- | Adds vertex type encapsulations. 
toJustTree :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree (Maybe a)
toJustTree (DTLeaf k) = DTLeaf k
toJustTree (DTVertex a sts) = DTVertex (Just a) (map toJustTree sts)

-- replace the following function by one that zips lists, sorts them once, and then unsplits them,
-- verifying both lists to be sorted afterwards?

-- | Verifies that two integer sequences correspond to the same total ordering of the entries.
equivalentOrders :: [Int] -> [Int] -> Bool
equivalentOrders ss ts = let
    sLookup :: [(Int,Int)]
    sLookup = zip (sort ss) [1..]
    tLookup :: [(Int,Int)]
    tLookup = zip (sort ts) [1..]
    sOrder = map (fromJust . flip lookup sLookup) ss
    tOrder = map (fromJust . flip lookup tLookup) ts
  in
    sOrder == tOrder

-- | Returns True if any of the vertices in the given tree has been tagged.
subTreeHasNothing :: (Ord a, Show a) => DecoratedTree (Maybe a) -> Bool
subTreeHasNothing (DTLeaf _) = False
subTreeHasNothing (DTVertex Nothing _) = True
subTreeHasNothing (DTVertex (Just _) sts) = any subTreeHasNothing sts

-- | The function that mimics resubstitution of a new tree into the hole left by finding embedding,
-- called m_\alpha,\beta in Dotsenko-Khoroshkin. This version only attempts to resubstitute the tree
-- at the root, bailing out if not possible.
reconstructNode :: (Ord a, Show a) => DecoratedTree a -> Embedding a -> Maybe (DecoratedTree a)
reconstructNode sub super = if isJust (vertexType super) then Nothing
                            else if (nLeaves sub) /= (vertexArity super) then Nothing
                            else let 
                                newSubTrees = map fromJustTree (subTrees super)
                           in
                             if any isNothing newSubTrees then Nothing
                                else let
                                    newTrees = map fromJust newSubTrees
                                    leafs = concatMap leafOrder newTrees
                                    newTree = nsComposeAll sub newTrees
                               in
                                 Just $ fmap ((leafs!!) . (subtract 1)) newTree

-- | The function that mimics resubstitution of a new tree into the hole left by finding embedding,
-- called m_\alpha,\beta in Dotsenko-Khoroshkin. This version recurses down in the tree in order
-- to find exactly one hole, and substitute the tree sub into it.
reconstructTree :: (Ord a, Show a) => DecoratedTree a -> Embedding a -> Maybe (DecoratedTree a)
reconstructTree sub super = if isLeaf super then Nothing
                            else if isNothing (vertexType super) then reconstructNode sub super
                            else if (1/=) . sum . map fromEnum $ map subTreeHasNothing (subTrees super) then Nothing
                            else
                                      let
                                          fromMaybeBy f s t = if isNothing t then f s else t
                                          subtreesSuper = subTrees super
                                          reconstructions = map (reconstructTree sub) subtreesSuper
                                          results = zipWith (fromMaybeBy fromJustTree) subtreesSuper reconstructions
                                      in
                                        if any isNothing results then Nothing
                                           else Just $ DTVertex (fromJust $ vertexType super) (map fromJust results)


-- ** Groebner basis methods

-- | Applies the reconstruction map represented by em to all trees in the operad element op. Any operad element that 
-- fails the reconstruction (by having the wrong total arity, for instance) will be silently dropped. We recommend
-- you apply this function only to homogenous operad elements, but will not make that check.
applyReconstruction :: (Ord a, Show a, TreeOrdering t, Num n) => Embedding a -> OperadElement a n t -> OperadElement a n t
applyReconstruction em m = let
      reconstructor (ordT, n) = do
        newDT <- reconstructTree (dt ordT) em
        return $ (ot newDT, n)
    in oe $ mapMaybe reconstructor (toList m)

-- | Finds all S polynomials for a given list of operad elements. 
findAllSPolynomials :: (Ord a, Show a, TreeOrdering t, Fractional n) => 
                       [OperadElement a n t] -> [OperadElement a n t] -> [OperadElement a n t]
findAllSPolynomials oldGb newGb = nub . map (\o -> (1/leadingCoefficient o) .*. o) . filter (not . isZero) $ do
  g1 <- oldGb ++ newGb
  g2 <- newGb
  let lmg1 = leadingMonomial g1
      lmg2 = leadingMonomial g2
      cf12 = (leadingCoefficient g1) / (leadingCoefficient g2)
  gamma <- nub $ findAllLCM lmg1 lmg2
  mg1 <- findAllEmbeddings lmg1 gamma
  mg2 <- findAllEmbeddings lmg2 gamma
  return $ (applyReconstruction mg1 g1) - (cf12 .*. (applyReconstruction mg2 g2))

-- | Finds all S polynomials for which the operationdegree stays bounded.
findInitialSPolynomials :: (Ord a, Show a, TreeOrdering t, Fractional n) =>
                           Int -> [OperadElement a n t] -> [OperadElement a n t] -> [OperadElement a n t]
findInitialSPolynomials n oldGb newGb = nub . map (\o -> (1/leadingCoefficient o) .*. o) . filter (not . isZero) $ do
  g1 <- oldGb ++ newGb
  g2 <- newGb
  let lmg1 = leadingMonomial g1
      lmg2 = leadingMonomial g2
      cf12 = (leadingCoefficient g1) / (leadingCoefficient g2)
  gamma <- nub $ findAllLCM lmg1 lmg2
  guard $ (operationDegree gamma) <= n
  mg1 <- findAllEmbeddings lmg1 gamma
  mg2 <- findAllEmbeddings lmg2 gamma
  return $ (applyReconstruction mg1 g1) - (cf12 .*. (applyReconstruction mg2 g2))

-- | Reduce g with respect to f and the embedding em: lt f -> lt g.
reduceOE :: (Ord a, Show a, TreeOrdering t, Fractional n) => Embedding a -> OperadElement a n t -> OperadElement a n t -> OperadElement a n t
reduceOE em f g = if not (divides (leadingMonomial f) (leadingMonomial g)) 
                  then g 
                  else let
                      cgf = (leadingCoefficient g) / (leadingCoefficient f)
                      ret = g - (cgf .*. (applyReconstruction em f))
                    in
                      if isZero ret then ret else (1/leadingCoefficient ret) .*. ret

reduceCompletely :: (Ord a, Show a, TreeOrdering t, Fractional n) => OperadElement a n t -> [OperadElement a n t] -> OperadElement a n t
reduceCompletely op [] = op
reduceCompletely op gb = if isZero op 
                         then op 
                         else let
                             divisorIdx = findIndex (flip divides (leadingMonomial op)) (map leadingMonomial gb)
                           in
                             if isNothing divisorIdx then op
                             else 
                               let 
                                 g1 = gb!!(fromJust divisorIdx) 
                                 em = head $ findAllEmbeddings (leadingMonomial g1) (leadingMonomial op)
                                 o1 = reduceOE em g1 op
                                in 
                                 reduceCompletely o1 gb

-- | Perform one iteration of the Buchberger algorithm: generate all S-polynomials. Reduce all S-polynomials.
-- Return anything that survived the reduction.
stepOperadicBuchberger :: (Ord a, Show a, TreeOrdering t, Fractional n) => 
                          [OperadElement a n t] -> [OperadElement a n t] -> [OperadElement a n t]
stepOperadicBuchberger oldGb newGb = nub $ do
  spol <- findAllSPolynomials oldGb newGb
  let red = reduceCompletely spol (oldGb ++ newGb)
  guard $ not . isZero $ red
  return red

-- | Perform one iteration of the Buchberger algorithm: generate all S-polynomials. Reduce all S-polynomials.
-- Return anything that survived the reduction. Keep the occurring operation degrees bounded. 
stepInitialOperadicBuchberger :: (Ord a, Show a, TreeOrdering t, Fractional n) => 
                          Int -> [OperadElement a n t] -> [OperadElement a n t] -> [OperadElement a n t]
stepInitialOperadicBuchberger maxD oldGb newGb = nub $ do
  spol <- findInitialSPolynomials maxD oldGb newGb
  guard $ maxOperationDegree spol <= maxD
  let red = reduceCompletely spol (oldGb ++ newGb)
  guard $ not . isZero $ red
  return red

-- | Perform the entire Buchberger algorithm for a given list of generators. Iteratively run the single iteration
-- from 'stepOperadicBuchberger' until no new elements are generated.
--
-- DO NOTE: This is entirely possible to get stuck in an infinite loop. It is not difficult to write down generators
-- such that the resulting Groebner basis is infinite. No checking is performed to catch this kind of condition.
operadicBuchberger :: (Ord a, Show a, TreeOrdering t, Fractional n) => [OperadElement a n t] -> [OperadElement a n t]
operadicBuchberger gb = let
    operadicBuchbergerAcc oldgb [] = oldgb
    operadicBuchbergerAcc oldgb new = operadicBuchbergerAcc 
                                      (newGB) 
                                      ((nub $ stepOperadicBuchberger oldgb new) \\ newGB) 
                                          where newGB = gb ++ new
  in operadicBuchbergerAcc [] gb

-- | Perform the entire Buchberger algorithm for a given list of generators. Iteratively run the single iteration
-- from 'stepOperadicBuchberger' until no new elements are generated. While doing this, maintain an upper bound
-- on the operation degree of any elements occurring.
--
initialOperadicBuchberger :: (Ord a, Show a, TreeOrdering t, Fractional n) =>
                             Int -> [OperadElement a n t] -> [OperadElement a n t]
initialOperadicBuchberger maxOD gb = let
    operadicBuchbergerAcc oldgb [] = oldgb
    operadicBuchbergerAcc oldgb new = if minimum (map maxOperationDegree new) > maxOD then oldgb 
                                   else
                                       operadicBuchbergerAcc 
                                       (newGB) 
                                       ((nub $ stepInitialOperadicBuchberger maxOD gb new) \\ newGB) 
                                           where newGB = gb ++ new
  in operadicBuchbergerAcc [] gb 

-- | Reduces a list of elements with respect to all other elements occurring in that same list.
reduceBasis :: (Ord a, Show a, TreeOrdering t, Fractional n) => [OperadElement a n t] -> [OperadElement a n t]
reduceBasis gb = map (\g -> reduceCompletely g (gb \\ [g])) gb

-- ** Low degree bases

-- | All trees composed from the given generators of operation degree n.
allTrees :: (Ord a, Show a) => 
            [DecoratedTree a] -> Int -> [DecoratedTree a]
allTrees _ 0 = []
allTrees gens 1 = gens
allTrees gens k = nub $ do
  gen <- gens
  tree <- allTrees gens (k-1)
  let degC = nLeaves gen
      degT = nLeaves tree
  i <- [1..degT]
  shuffle <- allShuffles i (degC - 1) (degT - i)
  return $ shuffleCompose i shuffle tree gen

-- | Generate basis trees for a given Groebner basis up through degree 'maxDegree'. 'divisors' is expected
-- to contain the leading monomials in the Groebner basis.
basisElements :: (Ord a, Show a) => 
                 [DecoratedTree a] -> [DecoratedTree a] -> Int -> [DecoratedTree a]
basisElements generators divisors maxDegree = nub $
    if maxDegree <= 0 then [] else if maxDegree == 1 then generators
--    else if null divisors then allTrees generators maxDegree
         else do
  b <- basisElements' generators divisors (maxDegree-1)
  gen <- generators
  let degC = nLeaves gen
      degT = nLeaves b
  i <- [1..degT]
  shuffle <- allShuffles i (degC-1) (degT-i)
  let newB = shuffleCompose i shuffle b gen
  guard $ not $ any (flip divides newB) divisors
  return newB

basisElements' :: (Ord a, Show a) => 
                 [DecoratedTree a] -> [DecoratedTree a] -> Int -> [DecoratedTree a]
basisElements' generators divisors maxDegree = if null divisors then allTrees generators maxDegree
                                              else do
  b <- allTrees generators maxDegree
  guard $ not $ any (flip divides b) divisors
  return b 

-- | Change the monomial order used for a specific tree. Use this in conjunction with mapMonomials
-- in order to change monomial order for an entire operad element. 
changeOrder :: (Ord a, Show a, TreeOrdering s, TreeOrdering t) => t -> OrderedTree a s -> OrderedTree a t
changeOrder o' (OT t _) = OT t o'