gistfile1.hs

module Math.Homology where

import Data.List  ((\\), subsequences)
import Data.Maybe (fromMaybe)
import Data.Tuple (swap)

newtype Simplex a = Simplex { getVertices :: [a] } deriving (Show, Eq)

simplexDim :: Simplex a -> Int
simplexDim = (subtract 1) . length . getVertices

facets :: Simplex a -> [Simplex a]
facets = map Simplex . deletions . getVertices
         where deletions [] = []
               deletions (x:xs) = xs : map (x:) (deletions xs)

newtype Chain r b = Chain { getCoefficients :: [(r, b)] } deriving (Show)

boundary :: (Num r) => Simplex a -> Chain r (Simplex a)
boundary = Chain . zip (cycle [1, -1]) . facets

newtype Complex a = Complex { getSimplices :: [Simplex a] } deriving (Show)

complexDim :: Complex a -> Int
complexDim = maximum . map simplexDim . getSimplices

-- All simplices in a complex that have a specified dimension
simplicesWithDim :: Complex a -> Int -> [Simplex a]
simplicesWithDim c n = filter (\s -> simplexDim s == n) (getSimplices c)

-- Add a simplex into a complex, recursively adding all faces of the simplex
addSimplex :: (Eq a) => Simplex a -> Complex a -> Complex a
addSimplex (Simplex []) c = c
addSimplex s c | s `elem` getSimplices c = c
               | otherwise               = foldr addSimplex newComplex (facets s)
                 where newComplex        = Complex (s : getSimplices c)

type Matrix r = [[r]]

matrix :: Int -> Int -> ((Int, Int) -> r) -> Matrix r
matrix n m f = map genRow [1..n]
               where
                 genRow i = map (\j -> f (i, j)) [1..m]

-- Given to bases and the image of each element for the first base, produce a Matrix
matrixFromAction :: (Eq a, Num r) => [a] -> [a] -> (a -> Chain r a) -> Matrix r
matrixFromAction n m f = matrix (length n) (length m) gen
                         where gen coord = fromMaybe 0 (coefficient coord)
                               coefficient (i, j) = lookup (n!!(i-1)) (assoc (m!!(j-1)))
                               assoc s = map swap (getCoefficients $ f s)

-- Copied from online --

removeZeroRowsAndColumns :: (Eq r, Num r) => Matrix r -> Matrix r
removeZeroRowsAndColumns = removeZeroRows . removeZeroColumns
  where
    removeZeroColumns x = map (drop n) x where
      n = foldl1 min (map countZeros x)
      countZeros = length . fst . break (0 /=)
    removeZeroRows = filter (or . map (0 /=))

rank :: (Eq r, Num r) => Matrix r -> Int
rank m | cleaned  == [] = 0
       | otherwise      = 1 + rank (map (cancelZero pivot) (zeros++rest))
  where cleaned = removeZeroRowsAndColumns m
        (zeros, (pivot:rest)) = break ((0 /=) . head) cleaned
        cancelZero (a:as) (b:bs) = zipWith (\x y -> a*y-b*x) as bs

------------------------

-- Find the boundary map at a specified dimension
d :: (Eq a) => Complex a -> Int -> Matrix Rational
d c n = matrixFromAction cn cn' boundary
        where cn  = simplicesWithDim c (n-1)
              cn' = simplicesWithDim c n

-- Calculate nth betti number of a complex
betti :: (Eq a) => Complex a -> Int -> Int
betti c n = kernelDim - imageDim
  where count = length $ simplicesWithDim c n
        kernelDim = count - (rank $ d c n)
        imageDim = rank $ d c (n+1)

-- Calculate all betti numbers of a complex
bettis :: (Eq a) => Complex a -> [Int]
bettis c = map (betti c) [0..(complexDim c)]

complexFromLists :: [[a]] -> Complex a
complexFromLists = Complex . map Simplex

ball :: Int -> Complex Int
ball n = complexFromLists $ subsequences [0..n] \\ [[]]

sphere :: Int -> Complex Int
sphere n = complexFromLists $ subsequences [0..(n+1)] \\ [[0..(n+1)], []]

torus :: Complex Int
torus = complexFromLists
        [[1],[2],[3],[4],[5],[6],[7],[8],[9],
         [1,2],[2,3],[1,3],
         [5,9],[8,9],[5,8],
         [4,6],[6,7],[4,7],

         [1,5],[4,5],[1,4],
         [2,9],[6,9],[2,6],
         [3,8],[7,8],[3,7],

         [1,9],[3,9],[1,8],
         [4,9],[6,8],[5,7],
         [1,6],[2,7],[1,7],

         [1,5,9],[1,2,9],[2,3,9],[3,8,9],[1,3,8],[1,5,8],
         [4,5,9],[4,6,9],[6,8,9],[6,7,8],[5,7,8],[4,5,7],
         [1,4,6],[1,2,6],[2,6,7],[2,3,7],[1,3,7],[1,4,7]]