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

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

#elif defined USE_POLYBAG
#else
#endif

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

-- | 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, bounded in total operation degree. findSmallBoundedLCM :: (Ord a, Show a) => Int -> DecoratedTree a -> DecoratedTree a -> [DecoratedTree a] findSmallBoundedLCM _ (DTLeaf _) _ = [] findSmallBoundedLCM _ _ (DTLeaf _) = [] findSmallBoundedLCM 0 _ _ = [] findSmallBoundedLCM n 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 = if (operationDegree s) > n || (operationDegree t) > n then [] else findRootedLCM s t
childLCMs = map (findSmallBoundedLCM (n-1) 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 filter ((<=n) . operationDegree) 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 = (findSmallBoundedLCM maxBound s t) ++ (findSmallBoundedLCM maxBound t s) -- | Finds all small common multiples of s and t, bounded in total operation degree. findAllBoundedLCM :: (Ord a, Show a) => Int -> DecoratedTree a -> DecoratedTree a -> [DecoratedTree a] findAllBoundedLCM n s t = (findSmallBoundedLCM n s t) ++ (findSmallBoundedLCM n 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) =>
findAllSPolynomials = findInitialSPolynomials maxBound

-- | Finds all S polynomials for which the operationdegree stays bounded.
findInitialSPolynomials :: (Ord a, Show a, TreeOrdering t, Fractional n) =>
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$ findAllBoundedLCM n lmg1 lmg2
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) =>

-- | 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) =>
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 = nub $initialOperadicBuchberger maxBound 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 let gbn = stepInitialOperadicBuchberger maxOD oldgb new gbo = nub$ oldgb ++ new
gbc = gbn \\ gbo
in
in nub $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'