{--
Author: Damek Davis
Website: http://www.math.ucla.edu/~damek/code.html
License: GNU GPL
--}

module HomologyZ2 (Homology, Complex, homology, example1, example2) where 
import Data.List
import qualified Data.Vector as IV  (unsafeIndex, unsafeFreeze, unsafeThaw, fromList, cons, head, tail, reverse, unsafeUpd, length, Vector, toList)
import qualified Data.Vector.Mutable as MIV (unsafeWrite, unsafeRead, clear, set, unsafeNew)
import Control.Monad.ST
import Control.Monad
import Data.STRef
import qualified Data.Set as Set

-- |A complex is a vector of triples.  Each index represents a generator, $x$, and the triple refers $(dx, d^{-1}x, alive?)$, where $dx$ is the list of generators that $x$ maps to, $d^{-1}x$ is the list of generators $y$ such that $x \in dy$, and alive? indicates whether x is still an element of the complex (alive? = true by default).
--
-- >>> example1
-- IV.fromList [([], [3], True), ([], [3], True), ([], [], True), ([0, 1, 2], [], True), ([], [3], True)]

type Complex = IV.Vector ([Int], [Int], Bool)

-- |A list of representatives of homology classes.

type Homology = [Int] 


-- |Compute the homology of a complex by the following method: Represent the complex by a directed graph, where the vertex set is the set of generators, and there is an edge from x to y if $y \in dx$. For each $x$ such that $dx \neq \emptyset$, pick $y \in dx$.  For each $z \in d^{-1}y$, and $w \in dx$, if there is an edge from $z$ to $w$, delete it, otherwise create an edge from $z$ to $w$. Finally, delete $x$, $y$, and all edges into and out of $x$ and $y$. Continue iterating this process until there are no edges left, then read off the homology (the list of elements that are still alive in the complex).
--
-- >>> homology example1
-- [1,2,4]
-- >>> homology example2
-- [2,3,4,5]

homology :: Complex -> Homology
homology cr = let v = runST $ do
                         vec <- IV.unsafeThaw cr
                         forM_ [0..IV.length cr - 1] $ \x -> do
                              (dx, d'1x, xAlive) <- MIV.unsafeRead vec x
                              when ((not $ null dx) && (xAlive)) $ do
                                      (dy, d'1y, yAlive) <- MIV.unsafeRead vec (head dx)
                                      forM_ d'1y $ \z -> do
                                         (dz, d'1z, zAlive) <- MIV.unsafeRead vec z
                                         when (zAlive && (z/=x)) $ MIV.unsafeWrite vec z (symmDiffList dz dx, d'1z, zAlive)
                                      forM_ dx $ \w -> do
                                         (dw, d'1w, wAlive) <- MIV.unsafeRead vec w
                                         when wAlive $ MIV.unsafeWrite vec w (dw, symmDiffList d'1y d'1w, wAlive)
                                      forM_ d'1x $ \z -> do
                                         (dz, d'1z, zAlive) <- MIV.unsafeRead vec z
                                         when zAlive $ MIV.unsafeWrite vec z (delete x dz, d'1z, zAlive) 
                                      forM_ dy $ \w -> do 
                                         (dw, d'1w, wAlive) <- MIV.unsafeRead vec w
                                         when wAlive $ MIV.unsafeWrite vec w (dw, delete (head dx) d'1w, wAlive) 
                                      MIV.unsafeWrite vec x ([], [], False)
                                      MIV.unsafeWrite vec (head dx) ([], [], False) 
                         IV.unsafeFreeze vec
               in [x | x <- [0..IV.length v - 1], triple3 (IV.unsafeIndex v x) ]


triple1 (x, y, z) = x
triple2 (x, y, z) = y
triple3 (x, y, z) = z

-- |Takes the symmetric difference of two lists.
--
-- >>> symmDiffList [1, 2, 3, 4] [2, 4, 5, 6, 7] 
-- [1, 3, 5, 6, 7]

symmDiffList :: (Eq a) => [a] -> [a] -> [a]
symmDiffList x y = (x `union` y) \\ (x `intersect` y)

-- |An example complex.
--
-- >>> example1
-- IV.fromList [([], [3], True), ([], [3], True), ([], [], True), ([0, 1, 2], [], True), ([], [3], True)]

example1 :: Complex
example1 = IV.fromList [([], [3], True), ([], [3], True), ([], [], True), ([0, 1, 2], [], True), ([], [3], True)]

-- |An example complex.
--
-- >>> example2
-- IV.fromList [([1,4],[],True), ([],[5,2,0],True), ([4,1],[],True), ([],[],True), ([],[5,2,0],True), ([1,4],[],True)]

example2 :: Complex
example2 = IV.fromList [([1,4],[],True),
           ([],[5,2,0],True),
           ([4,1],[],True),
           ([],[],True),
           ([],[5,2,0],True),
           ([1,4],[],True)]