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

module Math.Algebra.Group.StringRewriting where

import Data.List as L
import Data.Maybe (catMaybes)

-- REWRITING

-- |Given a list of rewrite rules of the form (left,right), and a word,
-- rewrite it by repeatedly replacing any left substring in the word by the corresponding right
rewrite :: (Eq a) => [([a], [a])] -> [a] -> [a]
rewrite rules word = rewrite' rules word where
    rewrite' (r:rs) xs =
        case rewrite1 r xs of
        Nothing -> rewrite' rs xs
        Just ys -> rewrite' rules ys
    rewrite' [] xs = xs

rewrite1 (l,r) xs =
    case xs `splitSubstring` l of
    Nothing -> Nothing
    Just (a,b) -> Just (a++r++b)

-- given a string x and a substring b, find if possible (a,c) such that xs = abc
splitSubstring xs ys = splitSubstring' [] xs where
    splitSubstring' ls [] = Nothing
    splitSubstring' ls (r:rs) =
        if ys `L.isPrefixOf` (r:rs)
        then Just (reverse ls, drop (length ys) (r:rs))
        else splitSubstring' (r:ls) rs
-- there might be a more efficient way to do this


-- KNUTH-BENDIX

-- given two strings x,y, find if possible a,b,c with x=ab y=bc
findOverlap xs ys = findOverlap' [] xs ys where
    findOverlap' as [] cs = Nothing -- (reverse as, [], cs)
    findOverlap' as (b:bs) cs =
        if (b:bs) `L.isPrefixOf` cs
        then Just (reverse as, b:bs, drop (length (b:bs)) cs)
        else findOverlap' (b:as) bs cs
-- there might be a more efficient way to do this

-- note that findOverlap "abab" "abab" won't find the partial overlap ("ab","ab","ab")

-- Knuth-Bendix algorithm
-- http://en.wikipedia.org/wiki/Knuth-Bendix_algorithm
-- Given a set of rules (assumed already reduced with respect to each other)
-- return a confluent rewrite system
knuthBendix1 rules = knuthBendix' rules pairs where
    pairs = [(lri,lrj) | lri <- rules, lrj <- rules, lri /= lrj]
    knuthBendix' rules [] = rules -- should reduce in some way
    knuthBendix' rules ( ((li,ri),(lj,rj)) : ps) =
        case findOverlap li lj of
        Nothing -> knuthBendix' rules ps
        Just (a,b,c) -> case ordpair (rewrite rules (ri++c)) (rewrite rules (a++rj)) of
                        Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
                        Just rule' -> let rules' = reduce rule' rules
                                          ps' = ps ++ [(rule',rule) | rule <- rules'] ++ [(rule,rule') | rule <- rules']
                                      in knuthBendix' (rule':rules') ps'
                        -- the new rule comes from seeing that
                        -- a ++ b ++ c == l1 ++ c -> r1 ++ c (by rule 1)
                        -- a ++ b ++ c == a ++ l2 -> a ++ r2 (by rule 2)
    reduce rule@(l,r) rules = filter (\(l',r') -> not (L.isInfixOf l l')) rules
        -- [rule' | rule'@(l',r') <- rules, not (l `L.isInfixOf` l')]

ordpair x y =
    case shortlex x y of
    LT -> Just (y,x)
    EQ -> Nothing
    GT -> Just (x,y)

shortlex x y = compare (length x, x) (length y, y)

-- for groups, where "letters" will take the form Either a a, we will want a different order, because we will want x^-1 -> x^3 to be the right way round


-- An optimisation - keep the rules ordered smallest first, and process the pairs smallest first
-- Appears to be significantly faster on average
knuthBendix2 rules = map snd $ knuthBendix' rules' pairs where
    rules' = L.sort $ map sizedRule rules
    pairs = L.sort [sizedPair sri srj | sri <- rules', srj <- rules', sri /= srj]
    knuthBendix' rules [] = rules
    knuthBendix' rules ( (s,(li,ri),(lj,rj)) : ps) =
        case findOverlap li lj of
        Nothing -> knuthBendix' rules ps
        Just (a,b,c) -> case ordpair (rewrite (map snd rules) (ri++c)) (rewrite (map snd rules) (a++rj)) of
                        Nothing -> knuthBendix' rules ps -- they both reduce to the same thing
                        Just rule' -> let rules' = reduce (snd rule') rules
                                          -- ps' = L.sort $ ps ++ [sizedPair rule' rule | rule <- rules'] ++ [sizedPair rule rule' | rule <- rules']
                                          ps' = merge ps $ merge [sizedPair rule' rule | rule <- rules'] [sizedPair rule rule' | rule <- rules']
                                     in knuthBendix' (L.insert rule' rules') ps'
    reduce rule@(l,r) rules = filter (\(s',(l',r')) -> not (L.isInfixOf l l')) rules
    -- reduce rule@(l,r) rules = [rule' | rule'@(s',(l',r')) <- rules, not (l `L.isInfixOf` l')]
    ordpair x y =
        let lx = length x; ly = length y in
            case compare (lx,x) (ly,y) of
            LT -> Just (ly,(y,x)); EQ -> Nothing; GT -> Just (lx,(x,y))
    sizedRule (rule@(l,r)) = (length l, rule)
    sizedPair (s1,r1) (s2,r2) = (s1+s2,r1,r2)

-- merge two ordered lists
merge (x:xs) (y:ys) =
    case compare x y of
    LT -> x : merge xs (y:ys)
    GT -> y : merge (x:xs) ys
    EQ -> error "" -- shouldn't happen in our case
merge xs ys = xs++ys

-- Another optimisation - at the stage where we remove some rules, we remove corresponding pairs too
-- Seems to perform about 25% faster on large problems (eg Coxeter groups A4-12, B4-12)
knuthBendix3 rules = knuthBendix' rules' pairs (length rules' + 1) where
    rules' = L.sort $ zipWith (\i (l,r) -> (length l,i,(l,r)) ) [1..] rules
    pairs = L.sort [sizedPair ri rj | ri <- rules', rj <- rules', ri /= rj]
    knuthBendix' rules [] k = map (\(s,i,r) -> r) rules
    knuthBendix' rules ( (s,(i,j),((li,ri),(lj,rj))) : ps) k =
        case findOverlap li lj of
        Nothing -> knuthBendix' rules ps k
        Just (a,b,c) -> case ordpair k (rewrite (map third rules) (ri++c)) (rewrite (map third rules) (a++rj)) of
                        Nothing -> knuthBendix' rules ps k -- they both reduce to the same thing
                        Just rule'@(_,_,(l,r)) ->
                            let (outrules,inrules) = L.partition (\(s',i',(l',r')) -> L.isInfixOf l l') rules
                                removedIndices = map second outrules
                                ps' = [p | p@(s,(i,j),(ri,rj)) <- ps, i `notElem` removedIndices, j `notElem` removedIndices]
                                ps'' = merge ps' $ merge [sizedPair rule' rule | rule <- inrules] [sizedPair rule rule' | rule <- inrules]
                            in knuthBendix' (L.insert rule' inrules) ps'' (k+1)
    ordpair k x y =
        let lx = length x; ly = length y in
            case compare (lx,x) (ly,y) of
            LT -> Just (ly,k,(y,x)); EQ -> Nothing; GT -> Just (lx,k,(x,y))
    second (s,i,r) = i
    third (s,i,r) = r
    sizedPair (si,i,ri) (sj,j,rj) = (si+sj,(i,j),(ri,rj))


-- |Implementation of the Knuth-Bendix algorithm. Given a list of relations, return a confluent rewrite system.
-- The algorithm is not guaranteed to terminate.
knuthBendix :: (Ord a) => [([a], [a])] -> [([a], [a])]
knuthBendix relations = knuthBendix3 (reduce [] rules) where
    rules = catMaybes [ordpair x y | (x,y) <- relations]
    reduce ls (r:rs) = reduce (r: reduce' r ls) (reduce' r rs)
    reduce ls [] = ls
    reduce' r rules = catMaybes [ordpair (rewrite [r] lhs) (rewrite [r] rhs) | (lhs,rhs) <- rules]


-- |Given generators and a confluent rewrite system, return (normal forms of) all elements
nfs :: (Ord a) => ([a], [([a], [a])]) -> [[a]]
nfs (gs,rs) = nfs' [[]] where
    nfs' [] = [] -- we have run out of words - this monoid is finite
    nfs' ws = let ws' = [g:w | g <- gs, w <- ws, not (any (`L.isPrefixOf` (g:w)) (map fst rs))]
              in ws ++ nfs' ws'

-- |Given generators and relations, return (normal forms of) all elements
elts :: (Ord a) => ([a], [([a], [a])]) -> [[a]]
elts (gs,rs) = nfs (gs, knuthBendix rs)


-- PRESENTATIONS FOR SOME STANDARD GROUPS
-- Would like to add a few more to this list

newtype SGen = S Int deriving (Eq,Ord)

instance Show SGen where
    show (S i) = "s" ++ show i

s_ i = S i
s1 = s_ 1
s2 = s_ 2
s3 = s_ 3

-- D L Johnson, Presentations of Groups, p62

-- symmetric group, generated by adjacent transpositions
_S n = (gs, r ++ s ++ t) where
    gs = map s_ [1..n-1]
    r = [([s_ i, s_ i],[]) | i <- [1..n-1]]
    s = [(concat $ replicate 3 [s_ i, s_ (i+1)],[]) | i <- [1..n-2]]
    t = [([s_ i, s_ j, s_ i, s_ j],[]) | i <- [1..n-1], j <- [i+2..n-1]]

-- braid presentation for Sn
_S' n = (gs, r ++ s ++ t) where
    gs = map s_ [1..n-1]
    r = [([s_ i, s_ i], []) | i <- [1..n-1]]
    s = [([s_ (i+1), s_ i, s_ (i+1)], [s_ i, s_ (i+1), s_ i] ) | i <- [1..n-2]]
    t = [([s_ i, s_ j, s_ i, s_ j], []) | i <- [1..n-1], j <- [i+2..n-1]]

-- http://en.wikipedia.org/wiki/Triangle_group
-- triangle groups - Johnson p127ff
tri l m n = ("abc", [("aa",""),("bb",""),("cc",""),("ab" ^ l,""),("bc" ^ n,""),("ca" ^ m,"" )])
    where xs ^ i = concat $ replicate i xs

-- von Dyck groups - Johnson p121ff
-- The subgroup of index 2 in the triangle group consisting of elts that preserve the orientation of the triangle
_D l m n = ("xy", [("x" ^ l,""), ("y" ^ m,""), ("xy" ^ n,"")])
    where xs ^ i = concat $ replicate i xs

-- Degenerate cases: n == 1 => cyclic group
-- l,2,2, l>=2 -> n-gon bipyramid - dihedral group
-- Spherical case: 1/l+1/m+1/n > 1
-- 3,3,2 -> tetrahedron; 4,3,2 -> octahedron; 5,3,2 -> icosahedron
-- Euclidean case: 1/l+1/m+1/n == 1
-- 3,3,3  4,4,2  6,3,2 
-- Hyperbolic case: 1/l+1/m+1/n < 1