something something homology groups

STL.hs

module STL
where

import qualified Data.List as L

data K3 f = K3 f f f
  deriving (Show, Eq)

data Void

instance Eq Void where
  (==) = const . const True

instance Num f => Num (K3 f) where
  K3 x1 x2 x3 + K3 y1 y2 y3 = K3 (x1 + y1) (x2 + y2) (x3 + y3)
  K3 l1 l2 l3 * K3 x1 x2 x3 = K3 (l1 * x1) (l2 * x2) (l3 * x3)
  negate (K3 x1 x2 x3) = K3 (negate x1) (negate x2) (negate x3)
  abs (K3 x1 x2 x3) = K3 (abs x1) (abs x2) (abs x3)
  signum (K3 x1 x2 x3) = K3 (signum x1) (signum x2) (signum x3)
  fromInteger n = K3 (fromInteger n) (fromInteger n) (fromInteger n)

(.*) :: Num f => f -> K3 f -> K3 f
l .* v = K3 l l l * v

data Simplex2 v f = Simplex2 (v f) (v f) (v f)
  deriving Show

instance (Eq f, Eq (v f)) => Eq (Simplex2 v f) where
  Simplex2 p1 p2 p3 == Simplex2 q1 q2 q3
    = (p1, p2, p3) == (q1, q2, q3) ||
      (p2, p3, p1) == (q1, q2, q3) ||
      (p3, p1, p2) == (q1, q2, q3) ||
      (p1, p3, p2) == (q1, q2, q3) ||
      (p2, p1, p3) == (q1, q2, q3) ||
      (p3, p2, p1) == (q1, q2, q3)

data Simplex1 v f = Simplex1 (v f) (v f)
  deriving Show

instance (Eq f, Eq (v f)) => Eq (Simplex1 v f) where
  Simplex1 p1 p2 == Simplex1 q1 q2
    = (p1, p2) == (q1, q2) || (p2, p1) == (q1, q2)

data Simplex0 v f = Simplex0 (v f)
  deriving (Show, Eq)

class Simplex s where
  generalize :: s v f -> GenSimplex v f
  dim :: s v f -> Int
  dim s = let (GenSimplex vs) = generalize s in length vs - 1

instance Simplex Simplex2 where
  generalize (Simplex2 p1 p2 p3) = GenSimplex [p1, p2, p3]

instance Simplex Simplex1 where
  generalize (Simplex1 p1 p2) = GenSimplex [p1, p2]

instance Simplex Simplex0 where
  generalize (Simplex0 p1) = GenSimplex [p1]

data GenSimplex v f = GenSimplex [v f]
  deriving Show

instance Simplex GenSimplex where
  generalize = id

instance (Eq f, Eq (v f)) => Eq (GenSimplex v f) where
  GenSimplex left == GenSimplex right
    = right `elem` (L.permutations left)

new3Simplex2 :: (f, f, f) -> (f, f, f) -> (f, f, f) -> Simplex2 K3 f
new3Simplex2 (x1, x2, x3) (y1, y2, y3) (z1, z2, z3)
  = Simplex2 (K3 x1 x2 x3) (K3 y1 y2 y3) (K3 z1 z2 z3)

new3Simplex1 :: (f, f, f) -> (f, f, f) -> Simplex1 K3 f
new3Simplex1 (x1, x2, x3) (y1, y2, y3)
  = Simplex1 (K3 x1 x2 x3) (K3 y1 y2 y3)

boundary1 :: Simplex2 v f -> [Simplex1 v f]
boundary1 (Simplex2 p1 p2 p3) = [Simplex1 p1 p2, Simplex1 p2 p3, Simplex1 p3 p1]

boundary0 :: Simplex1 v f -> [Simplex0 v f]
boundary0 (Simplex1 p1 p2) = [Simplex0 p1, Simplex0 p2]

skeleton1 :: [Simplex2 v f] -> [Simplex1 v f]
skeleton1 (s:sx3) = (skeleton1 sx3) ++ (boundary1 s)
skeleton1 [] = []

skeleton0 :: [Simplex1 v f] -> [Simplex0 v f]
skeleton0 (s:sx2) = (skeleton0 sx2) ++ (boundary0 s)
skeleton0 [] = []

skeleton2to0 :: [Simplex2 v f] -> [Simplex0 v f]
skeleton2to0 = skeleton0 . skeleton1

class Mod t where
  (+~) :: Eq x => t x -> t x -> t x
  (*~) :: Eq x => Integer -> t x -> t x
  unit :: Eq x => x -> t x
  zero :: t x
  (-~) :: Eq x => t x -> t x -> t x
  x -~ y = x +~ ((-1) *~ y)

infixl 6 +~
infixl 7 *~

data FreeAbGrp t = FreeAbGrp [t] (t -> Integer)

unitF :: Eq x => x -> FreeAbGrp x
unitF = unit

instance Show t => Show (FreeAbGrp t) where
  show (FreeAbGrp (x:xs) f)
    = "(" ++ show (f x) ++ ") *~ unitF (" ++ show x ++ ")" ++ case xs of
      [] -> ""
      x  -> " +~ " ++ show (FreeAbGrp xs f)

instance Eq t => Eq (FreeAbGrp t) where
  x == y = ((xs L.\\ ys) `L.union` (ys L.\\ xs)) == []
    where
      xs = toTupleList x
      ys = toTupleList y

instance Mod FreeAbGrp where
  x +~ FreeAbGrp [] _ = x
  FreeAbGrp [] _ +~ y = y
  FreeAbGrp (r:rs) rf +~ FreeAbGrp (s:ss) sf
    = FreeAbGrp (L.nub $ s:r:rs) f +~ FreeAbGrp ss sf
    where
      f t = if t == s then rf t + sf t else rf t
  unit x = FreeAbGrp [x] (\s -> if s == x then 1 else 0)
  n *~ FreeAbGrp rs f = FreeAbGrp rs ((*n) . f)
  zero = FreeAbGrp [] (const 0)

toTupleList :: FreeAbGrp t -> [(Integer, t)]
toTupleList (FreeAbGrp (x:xs) f) = (f x, x):toTupleList (FreeAbGrp xs f)
toTupleList _ = []

testObj :: [Simplex2 K3 Float]
testObj = L.nub [
    -- green
    new3Simplex2 (0, 10, 0) (15, 5, 0) (0, 0, 0)
  , new3Simplex2 (35, 15, 0) (40, 0, 0) (30, 0, 0)
  , new3Simplex2 (40, 40, 0) (40, 30, 0) (25, 35, 0)
  , new3Simplex2 (10, 40, 0) (0, 40, 0) (5, 25, 0)
    -- blue
  , new3Simplex2 (15, 5, 0) (30, 0, 0) (0, 0, 0)
  , new3Simplex2 (40, 30, 0) (40, 0, 0) (35, 15, 0)
  , new3Simplex2 (40, 40, 0) (25, 35, 0) (10, 40, 0)
  , new3Simplex2 (0, 40, 0) (5, 25, 0) (0, 10, 0)
    -- fuchsia
  , new3Simplex2 (30, 0, 0) (35, 15, 0) (15, 5, 0)
  , new3Simplex2 (35, 15, 0) (40, 30, 0) (25, 35, 0)
  , new3Simplex2 (25, 35, 0) (10, 40, 0) (5, 25, 0)
  , new3Simplex2 (5, 25, 0) (15, 5, 0) (0, 10, 0)
    -- aqua
  , new3Simplex2 (25, 35, 0) (35, 15, 0) (5, 25, 0)
    -- yellow
  , new3Simplex2 (5, 25, 0) (35, 15, 0) (15, 5, 0)
  ]

triangle :: [Simplex1 K3 Float]
triangle = [
    Simplex1 (K3 0 0 0) (K3 2 0 0)
  , Simplex1 (K3 0 0 0) (K3 1 1 0)
  , Simplex1 (K3 1 1 0) (K3 2 0 0)
  ]

fsum :: Eq t => [FreeAbGrp t] -> FreeAbGrp t
fsum = foldl (+~) zero

d2 :: (Num f, Eq f, Eq (v f)) => FreeAbGrp (Simplex2 v f) -> FreeAbGrp (Simplex1 v f)
d2 (FreeAbGrp [] _) = zero
d2 (FreeAbGrp (s:ss) f)
  = d2 (FreeAbGrp ss f) +~ f s *~ (fsum . (unit <$>) . boundary1) s

d1 :: (Num f, Eq f, Eq (v f)) => FreeAbGrp (Simplex1 v f) -> FreeAbGrp (Simplex0 v f)
d1 (FreeAbGrp [] _) = zero
d1 (FreeAbGrp (s:ss) f)
  = d1 (FreeAbGrp ss f) +~ f s *~ (fsum . (unit <$>) . boundary0) s

d0 :: FreeAbGrp (Simplex0 v f) -> FreeAbGrp Void
d0 _ = FreeAbGrp [] (const 0)

data Matrix t = Matrix [[t]] deriving Eq

instance Show t => Show (Matrix t) where
  show (Matrix lines)
    = "[" ++ (foldl (++) "" . L.intersperse "\n ") (show <$> L.transpose lines) ++ "]"

instance Functor Matrix where
  fmap f (Matrix xs) = Matrix $ fmap (fmap f) xs

matrix :: (Eq t, Eq u) => [t] -> [u] -> (FreeAbGrp t -> FreeAbGrp u) -> Matrix Integer
matrix [] _ _ = Matrix []
matrix (sb:sbasis) dbasis hom
  = let FreeAbGrp _ f = hom (unit sb)
        Matrix rem = matrix sbasis dbasis hom
    in Matrix $ [f db | db <- dbasis] : rem

homologyMat :: (Eq t, Eq u) => [t] -> [u] -> (FreeAbGrp t -> FreeAbGrp u) -> Matrix Integer
homologyMat sbasis dbasis hom = (`mod` 2) <$> matrix sbasis dbasis hom