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
#ifdef TRACE
import Debug.Trace
import Math.Operad.PPrint
#endif
operationDegree :: (Ord a, Show a) => DecoratedTree a -> Int
operationDegree (DTLeaf _) = 0
operationDegree vertex = 1 + sum (map operationDegree (subTrees vertex))
operationDegrees :: (Ord a, Show a, TreeOrdering t, Num n) => OperadElement a n t -> [Int]
operationDegrees op = foldMonomials (\(k,_) lst -> lst ++ [(operationDegree . dt $ k)]) op
maxOperationDegree :: (Ord a, Show a, TreeOrdering t, Num n) => OperadElement a n t -> Int
maxOperationDegree element = if null op then 0 else maximum op where op = operationDegrees element
isHomogenous :: (Ord a, Show a, TreeOrdering t, Num n) => OperadElement a n t -> Bool
isHomogenous m = let
trees = getTrees m
equalityCheck s t = arityDegree (dt s) == arityDegree (dt t) &&
operationDegree (dt s) == operationDegree (dt t)
in and $ zipWith (equalityCheck) trees (tail trees)
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 t1)) =
error $ "Need a correct pointed shuffle permutation!\n" ++
show sh ++ " is not in Sh" ++ show i ++
"(" ++ show (nLeaves t1) ++ "," ++ show (nLeaves si) ++ ")\n"
| otherwise = symmetricCompose i sh s t
nsCompose :: (Ord a, Show a) => Int -> DecoratedTree a -> DecoratedTree a -> DecoratedTree a
nsCompose i s t = if i1 > 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 (idx1) trees ++ [t] ++ drop idx trees
in
if length newTrees /= nLeaves s then error "Shouldn't happen"
else
nsComposeAll s newTrees
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
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
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
isPermutation :: Shuffle -> Bool
isPermutation sh = and ((zipWith (==) [1..]) (sort sh))
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
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
type Embedding a = DecoratedTree (Maybe a)
divides :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> Bool
divides s t = not . null $ findAllEmbeddings s t
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 (i1) (map toJustTree $ subTrees t) ++ [tree] ++ drop i (map toJustTree $ subTrees t)))
in map glueTree ems
in rootFind ++ concatMap reGlue (zip [1..] subFinds)
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))
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
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]
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
planarTree :: (Ord a, Show a) => DecoratedTree a -> Bool
planarTree (DTLeaf _) = True
planarTree (DTVertex _ subs) = all planarTree subs && isSorted (map minimalLeaf subs)
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
let rootedLCMs = if (operationDegree s) > n || (operationDegree t) > n then [] else findRootedLCM s t
childLCMs = map (findSmallBoundedLCM (n1) 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 (i1) (subTrees t) ++ [leaf (minimalLeaf (subTrees t !! (i1)))] ++ drop i (subTrees t))
in concatMap (\emt -> accumulateTrees [(leaf i,emt)] [template]) ems
zippedChildLCMs = zip [1..] childLCMs
filter ((<=n) . operationDegree) rootedLCMs ++ (concatMap reGlue zippedChildLCMs)
findAllLCM :: (Ord a, Show a) => DecoratedTree a -> DecoratedTree a -> [DecoratedTree a]
findAllLCM s t = (findSmallBoundedLCM maxBound s t) ++ (findSmallBoundedLCM maxBound t s)
findAllBoundedLCM :: (Ord a, Show a) => Int -> DecoratedTree a -> DecoratedTree a -> [DecoratedTree a]
findAllBoundedLCM n s t = (findSmallBoundedLCM n s t) ++ (findSmallBoundedLCM n t s)
rePackLabels :: (Ord a, Show a, Ord b) => PreDecoratedTree a b -> DecoratedTree a
rePackLabels tree = fmap (fromJust . (flip lookup (zip (sort (foldMap (:[]) tree)) [1..]))) tree
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
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)
equivalentOrders :: [Int] -> [Int] -> Bool
equivalentOrders o1 o2 = if length o1 /= length o2 then False
else and $ zipWith (==) (zipWith compare o1 (tail o1)) (zipWith compare o2 (tail o2))
subTreeHasNothing :: (Ord a, Show a) => DecoratedTree (Maybe a) -> Bool
subTreeHasNothing (DTLeaf _) = False
subTreeHasNothing (DTVertex Nothing _) = True
subTreeHasNothing (DTVertex (Just _) sts) = any subTreeHasNothing sts
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
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)
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)
findAllSPolynomials :: (Ord a, Show a, TreeOrdering t, Fractional n) =>
[OperadElement a n t] -> [OperadElement a n t] -> [OperadElement a n t]
findAllSPolynomials = findInitialSPolynomials maxBound
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 $ findAllBoundedLCM n lmg1 lmg2
#ifdef TRACE
trace ("Found LCM: \n\t" ++ pp lmg1 ++ ", \n\t" ++ pp lmg2 ++ ": \n\t" ++ pp gamma ++ "\n") (return ())
#endif
mg1 <- findAllEmbeddings lmg1 gamma
mg2 <- findAllEmbeddings lmg2 gamma
return $ (applyReconstruction mg1 g1) (cf12 .*. (applyReconstruction mg2 g2))
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
stepOperadicBuchberger :: (Ord a, Show a, TreeOrdering t, Fractional n) =>
[OperadElement a n t] -> [OperadElement a n t] -> [OperadElement a n t]
stepOperadicBuchberger oldGb newGb = stepInitialOperadicBuchberger maxBound oldGb newGb
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 $
filter (not . isZero) $
do
spol <- findInitialSPolynomials maxD oldGb newGb
guard $ maxOperationDegree spol <= maxD
let red =
#ifdef TRACE
trace ("Reducing S-polynomial: " ++ pp spol ++ "\n") $
#endif
reduceCompletely spol (oldGb ++ newGb)
guard $ not . isZero $ red
return red
operadicBuchberger :: (Ord a, Show a, TreeOrdering t, Fractional n) => [OperadElement a n t] -> [OperadElement a n t]
operadicBuchberger gb = nub $ initialOperadicBuchberger maxBound gb
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 =
#ifdef TRACE
trace ("Stepping through\n") $
#endif
stepInitialOperadicBuchberger maxOD oldgb new
gbo = reduceBasis $ oldgb ++ new
gbc = (reduceBasis (gbn ++ gbo)) \\ gbo
in
operadicBuchbergerAcc gbo gbc
in nub $ operadicBuchbergerAcc [] gb
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
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 (k1)
let degC = nLeaves gen
degT = nLeaves tree
i <- [1..degT]
shuffle <- allShuffles i (degC 1) (degT i)
return $ shuffleCompose i shuffle tree gen
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 do
b <- basisElements' generators divisors (maxDegree1)
gen <- generators
let degC = nLeaves gen
degT = nLeaves b
i <- [1..degT]
shuffle <- allShuffles i (degC1) (degTi)
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
changeOrder :: (Ord a, Show a, TreeOrdering s, TreeOrdering t) => t -> OrderedTree a s -> OrderedTree a t
changeOrder o' (OT t _) = OT t o'