module Math.Topology.CubeCmplx.DirCubeCmplx (
T, Vertex, vertex, coords, vertexUnsafe, vertexVectorUnsafe, vertexToList,
vertexPtWise, vAdd, vSub, vSubUnsafe, vMin, vMax, vGT, vLT, vDim,
VertSpan, vertSpan, vsFst, vsSnd, vsUnsafe, vsVert, vsFstList, vsSndList,
vsCoords, vsCoordsUnsafe, vsDim, vsIsCell, vsFatten, vsCornerPairs,
vsCornerVerts, vsBdry,
CubeCell, minVert, maxVert, cell, cellUnsafe, cellDim, cellVertsUnsafe,
cellVerts, spanTopCells, vertToCell, inSpan, vInSpan, inBdry, spanBdryCells,
nCubes, nCubeVerts, nCubeCells, nCubeProperCells, nCubeBdry, nCubeKSkels,
verts, subCells, properSubCells, bdry, kSkel, isSubCell, isPropSubCell,
opFaceUnsafe,
genToNonGen, nonGenToGen,
CubeCmplx, cells, cmplxEmpty, cmplxNull, cmplxSize, cmplxApply, vsCmplx,
cmplxDelCell, cmplxDelCells, cmplxAddCells, cmplxUnions, cmplxFilter,
cmplxHullUnsafe, cmplxFilterSpan, cmplxFilterSpans, cellNhd,
swissFlag, sqPairFwd, sqPairBack, torus3d, genusTwo3d,
lazyProd
) where
import Data.Int (Int8)
import Data.Maybe (fromJust)
import Data.List (transpose, groupBy, sortBy)
import Data.Ord (comparing)
import Data.Function (on)
import Control.Monad (liftM, guard)
import Data.Hashable (Hashable, hashWithSalt, hash)
import Control.DeepSeq (deepseq, NFData(..))
import qualified Data.HashSet as S
(HashSet, fromList, filter, toList, union, empty, unions,
difference, map, size, singleton, null, foldr, delete, member)
import qualified Data.HashMap.Strict as M
(HashMap, empty, null, insertWith, fromListWith, filter,
fromList, lookup, toList)
import qualified Data.Vector.Unboxed as V
((//), sum, toList, all, fromList, Vector, zipWith, length, map, update,
(!), replicate, elemIndex, elemIndices, (++), zip, enumFromN, drop, take,
singleton, Unbox, accumulate)
import Data.Bits ((.&.), (.|.), xor)
import Control.Parallel.Strategies
(rdeepseq, parBuffer, withStrategy, parList, dot, evalTuple3, r0)
import Control.Arrow ((***))
import Test.QuickCheck (Arbitrary, arbitrary, suchThat, choose, vectorOf,
resize)
lazyProd :: [[a]] -> [[a]]
lazyProd [] = [[]]
lazyProd [x] = map (:[]) x
lazyProd (x1:x2:xs) = concat . concat $
[[[y1:y2:yn | y1<-x1] | y2<-x2] | yn <- (lazyProd xs)]
bitFill :: V.Vector T -> V.Vector T -> V.Vector T -> V.Vector T
bitFill v m f = V.accumulate (+) v $ V.zip maskIndices f
where maskIndices = V.elemIndices 1 m
type T = Int8
data Vertex = Vertex { coords :: V.Vector T, _hash :: Int } deriving (Eq, Ord)
instance Show Vertex where show v = show (V.toList $ coords v)
instance Arbitrary Vertex
where arbitrary = do l <- choose (1,5)
ts <- vectorOf l (arbitrary `suchThat`
(\t -> t >= 0 && t <= 63))
return (vertexUnsafe ts)
instance Hashable Vertex where hashWithSalt s v = s + (_hash v)
instance NFData Vertex where rnf v = (rnf $ coords v) `seq`
(rnf $ _hash v) `seq` ()
vertex :: [T] -> Maybe Vertex
vertex ts | null ts = Nothing
| any (< 0) ts = Nothing
| otherwise = Just . vertexUnsafe $ ts
vertexUnsafe :: [T] -> Vertex
vertexUnsafe ts = Vertex { coords = V.fromList ts, _hash = hash ts }
vertexVectorUnsafe :: V.Vector T -> Vertex
vertexVectorUnsafe v = Vertex { coords = v, _hash = hash $ V.toList v }
vertexToList :: Vertex -> [T]
vertexToList = V.toList . coords
vertexPtWise :: (T -> T -> T) -> Vertex -> Vertex -> Vertex
vertexPtWise f v1 v2 = vertexVectorUnsafe $
V.zipWith (\x y -> if (f x y) < 0 then 0 else f x y)
(coords v1) (coords v2)
vAdd :: Vertex -> Vertex -> Vertex
vAdd = vertexPtWise (+)
vSub :: Vertex -> Vertex -> Vertex
vSub = vertexPtWise ()
vSubUnsafe :: Vertex -> Vertex -> Vertex
vSubUnsafe v1 v2 = vertexVectorUnsafe $ V.zipWith () (coords v1) (coords v2)
vMin :: Vertex -> Vertex -> Vertex
vMin = vertexPtWise (min)
vMax :: Vertex -> Vertex -> Vertex
vMax = vertexPtWise (max)
vLT :: Vertex -> Vertex -> Bool
vLT v1 v2 = V.all (==True) $ V.zipWith (<=) (coords v1) (coords v2)
vGT :: Vertex -> Vertex -> Bool
vGT = flip vLT
vDim :: Vertex -> Int
vDim = V.length . coords
data VertSpan = VertSpan { vsFst :: !Vertex, vsSnd :: !Vertex }
deriving (Show, Eq, Ord)
instance NFData VertSpan
instance Arbitrary VertSpan
where arbitrary = do v1 <- arbitrary
v2 <- (resize 6 arbitrary) `suchThat`
((== vDim v1).vDim)
return $ vsUnsafe v1 (v1 `vAdd` v2)
vertSpan :: Vertex -> Vertex -> Maybe VertSpan
vertSpan v1 v2
| (v1 `vLT` v2) && (vDim v1 == vDim v2) = Just $ VertSpan v1 v2
| otherwise = Nothing
vsUnsafe :: Vertex -> Vertex -> VertSpan
vsUnsafe = VertSpan
vsVert :: Vertex -> VertSpan
vsVert v = vsUnsafe v v
vsFstList :: VertSpan -> [T]
vsFstList = vertexToList . vsFst
vsSndList :: VertSpan -> [T]
vsSndList = vertexToList . vsSnd
vsCoords :: [T] -> [T] -> Maybe VertSpan
vsCoords t1 t2 = do v1 <- vertex t1; v2 <- vertex t2; vertSpan v1 v2
vsCoordsUnsafe :: [T] -> [T] -> VertSpan
vsCoordsUnsafe t1 t2 = vsUnsafe (vertexUnsafe t1) (vertexUnsafe t2)
vsDim :: VertSpan -> Int
vsDim vs = V.sum $ V.zipWith (\x y -> if x /= y then 1 else 0)
(coords $ vsFst vs) (coords $ vsSnd vs)
vsIsCell :: VertSpan -> Bool
vsIsCell vs = V.all (flip elem [0,1]) . coords $
(vsSnd vs) `vSubUnsafe` (vsFst vs)
vsFatten :: VertSpan -> VertSpan
vsFatten vs = vsUnsafe ((vsFst vs) `vSub` d) ((vsSnd vs) `vAdd` d)
where c = head $ spanTopCells vs
d = (maxVert c) `vSub` (minVert c)
vsCornerPairs :: VertSpan -> S.HashSet (CubeCell, Vertex)
vsCornerPairs vs
| vsDim vs == 0 = S.singleton $ (cellUnsafe (vsFstList vs) (vsSndList vs),
vertexUnsafe (vsFstList vs))
| otherwise = S.fromList $ zip cells corners
where coordSpans = transpose [vsFstList vs, vsSndList vs]
coordRans = map (\cs -> enumFromTo (head cs) (last cs)) coordSpans
coordRans' = map (\cs -> enumFromThenTo (last cs) (last cs1)
(head cs)) coordSpans
possCoords = zipWith (\l1 l2 -> [l1, reverse l2])
(map (take 2) coordRans) (map (take 2) coordRans')
cells = map (\[x,y] -> cellUnsafe x y) . map transpose $
lazyProd possCoords
corners = map vertexUnsafe . lazyProd $
map (\[x,y] -> [head x, last y]) possCoords
vsCornerVerts :: VertSpan -> S.HashSet Vertex
vsCornerVerts = S.map snd . vsCornerPairs
newtype BitVector = BitVector { bitVect :: V.Vector T } deriving (Show)
instance Arbitrary BitVector where
arbitrary = do l <- choose (1,7)
bs <- vectorOf l (choose (0,1))
return . BitVector $ V.fromList bs
data CubeCell = CubeCell { _minVert :: !Vertex, _maxVert :: !Vertex } deriving (Eq)
instance NFData CubeCell
instance Hashable CubeCell
where hashWithSalt s c = hashWithSalt s (_minVert c, _maxVert c)
instance Ord CubeCell where
c1 <= c2 = (minVert c1, maxVert c1) <= (minVert c2, maxVert c2)
instance Show CubeCell
where show c = "(" ++ show (cellDim c) ++ ","
++ show (minVert c) ++ "," ++ show (maxVert c) ++ ")"
instance Arbitrary CubeCell
where arbitrary = do v1 <- arbitrary
v2 <- (liftM (vertexVectorUnsafe . bitVect) $ arbitrary)
`suchThat` ((== vDim v1) . vDim)
return $ cellVertsUnsafe v1 (v1 `vAdd` v2)
minVert :: CubeCell -> Vertex
minVert c = _minVert c
maxVert :: CubeCell -> Vertex
maxVert c = _maxVert c
cellVertsUnsafe :: Vertex -> Vertex -> CubeCell
cellVertsUnsafe v1 v2 = CubeCell v1 v2
cellUnsafe :: [T] -> [T] -> CubeCell
cellUnsafe t1 t2 = cellVertsUnsafe (vertexUnsafe t1) (vertexUnsafe t2)
cellVerts :: Vertex -> Vertex -> Maybe CubeCell
cellVerts v1 v2 = do vs <- vertSpan v1 v2
guard (vsIsCell vs)
return $ cellVertsUnsafe v1 v2
cell :: [T] -> [T] -> Maybe CubeCell
cell t1 t2 = do v1 <- vertex t1; v2 <- vertex t2; cellVerts v1 v2
cellDim :: CubeCell -> Int
cellDim c = fromEnum . V.sum . coords $ maxVert c `vSubUnsafe` minVert c
spanTopCells :: VertSpan -> [CubeCell]
spanTopCells = map pairUp . vertSpans
where pairUp [a,b] = cellUnsafe a b
vertSpans vs = map transpose . lazyProd .
map (pairs . uncurry enumFromTo) $
zip (vsFstList vs) (vsSndList vs)
pairs [] = []
pairs [x] = [[x,x]]
pairs xs = zipWith (\a b -> [a,b]) xs (tail xs)
vertToCell :: Vertex -> CubeCell
vertToCell v = cellVertsUnsafe v v
inSpan :: CubeCell -> VertSpan -> Bool
inSpan c vs = (vsFst vs `vLT` minVert c) && (maxVert c `vLT` vsSnd vs)
vInSpan :: Vertex -> VertSpan -> Bool
vInSpan v vs = (vertToCell v) `inSpan` vs
data VertType = Min | Max | Neither deriving (Show,Eq)
inBdry :: CubeCell -> VertSpan -> Bool
inBdry c vs = any (==True) $
zipWith (\a b -> a == b && a /= Neither)
(vertBdryCmpts vs $ minVert c)
(vertBdryCmpts vs $ maxVert c)
where vertBdryCmpts vs v = zipWith3 cmp (vsFstList vs) (vsSndList vs)
(vertexToList v)
cmp min max i | i == min = Min
| i == max = Max
| otherwise = Neither
vsBdry :: VertSpan -> [VertSpan]
vsBdry vs = map (uncurry vsUnsafe) (fstSnd fst ++ fstSnd snd)
where ranges = V.zip (coords $ vsFst vs) (coords $ vsSnd vs)
modVec f i = V.take i ranges
V.++ (V.singleton . (\t -> (t,t)) . f $ ranges V.! i)
V.++ V.drop (i+1) ranges
fstSnd f = zip (map (vertexVectorUnsafe .
V.map fst . modVec f) [0..V.length ranges1])
(map (vertexVectorUnsafe .
V.map snd . modVec f) [0..V.length ranges1])
spanBdryCells :: VertSpan -> [[CubeCell]]
spanBdryCells = map spanTopCells . vsBdry
nCubes :: [CubeCell]
nCubes = map gen [0..]
where gen n = cellUnsafe (replicate n 0) (replicate n 1)
nCubeVerts :: Int -> [CubeCell]
nCubeVerts n | n < 0 = []
| otherwise = nCubesVerts !! n
nCubesVerts = map nCubeVerts' [0..]
where nCubeVerts' 0 = map (vertToCell . vertexUnsafe) [[0]]
nCubeVerts' n = map (vertToCell . vertexUnsafe) . lazyProd $
replicate n [0,1]
nCubeCells :: Int -> [CubeCell]
nCubeCells n | n < 0 = []
| otherwise = nCubesCells !! n
nCubesCells = map nCubeCells' [0..]
where nCubeCells' n = [cellVertsUnsafe v1 v2 |
v1 <- map minVert $ nCubeVerts n,
v2 <- map minVert $ nCubeVerts n, v1 `vLT` v2]
nCubeProperCells :: Int -> [CubeCell]
nCubeProperCells n = filter ((/= n) . cellDim) . nCubeCells $ n
nCubeBdry :: Int -> [CubeCell]
nCubeBdry n | n < 0 = []
| otherwise = nCubesBdry !! n
nCubesBdry = map nCubeBdry' [0..]
where nCubeBdry' n = concat . spanBdryCells $ vsCoordsUnsafe
(replicate n 0) (replicate n 1)
nCubeKSkels :: Int -> Int -> [CubeCell]
nCubeKSkels n k | n < 0 || k < 0 = []
| k > n = [nCubes !! n]
| otherwise = nCubesKSkels !! n !! k
nCubesKSkels = map nCubeKSkels' [0..]
where nCubeKSkels' = groupBy ((==) `on` cellDim) .
sortBy (comparing cellDim) . nCubeCells
genToNonGen :: CubeCell -> CubeCell -> CubeCell
genToNonGen c g = cellVertsUnsafe l u
where bitMask = coords $ maxVert c `vSubUnsafe` minVert c
minc = coords $ minVert c
l = vertexVectorUnsafe $ bitFill minc bitMask (coords $ minVert g)
u = vertexVectorUnsafe $ bitFill minc bitMask (coords $ maxVert g)
nonGenToGen :: CubeCell -> CubeCell -> CubeCell
nonGenToGen c s = cellUnsafe (zipWith (V.!) (repeat $ locMin) indices)
(zipWith (V.!) (repeat $ locMax) indices)
where locMin = coords $ minVert s `vSubUnsafe` minVert c
locMax = coords $ maxVert s `vSubUnsafe` minVert c
bitMask = coords $ maxVert c `vSubUnsafe` minVert c
indices = V.toList . V.elemIndices 1 $ bitMask
lookupSubCells :: [[CubeCell]] -> CubeCell -> [CubeCell]
lookupSubCells l c = map (genToNonGen c) $ l !! cellDim c
verts :: CubeCell -> [Vertex]
verts c = map minVert $ lookupSubCells nCubesVerts c
subCells :: CubeCell -> [CubeCell]
subCells = lookupSubCells nCubesCells
properSubCells :: CubeCell -> [CubeCell]
properSubCells = lookupSubCells (map nCubeProperCells [0..])
bdry :: CubeCell -> [CubeCell]
bdry = lookupSubCells nCubesBdry
kSkel :: Int -> CubeCell -> [CubeCell]
kSkel k c | k < 0 = []
| otherwise = map (genToNonGen c) gs
where gs = nCubeKSkels (cellDim c) k
isSubCell :: CubeCell -> CubeCell -> Bool
isSubCell s c = inSpan s $ vsUnsafe (minVert c) (maxVert c)
isPropSubCell :: CubeCell -> CubeCell -> Bool
isPropSubCell s c = (isSubCell s c) && (cellDim c /= cellDim s)
genOpFaces :: [M.HashMap CubeCell CubeCell]
genOpFaces = map opFaces [0..]
where differ v1 v2 = V.zipWith xor (V.zipWith (.&.) v1 v2)
(V.zipWith (.|.) v1 v2)
invert v1 v2 = V.map (xor 1) $ differ v1 v2
index v1 v2 = fromJust $ V.elemIndex 1 $ invert v1 v2
newVal v1 v2 = (index v1 v2, 1 v1 V.! index v1 v2)
newVerts v1 v2 = map (vertexUnsafe . V.toList .
flip (V.//) [newVal v1 v2]) [v1, v2]
opVerts c = newVerts (coords $ minVert c) (coords $ maxVert c)
opFace c = cellVertsUnsafe (head $ opVerts c) (last $ opVerts c)
opFaces n = M.fromList . zip (nCubesBdry !! n) $
map (opFace) (nCubesBdry !! n)
opFaceUnsafe :: CubeCell -> CubeCell -> CubeCell
opFaceUnsafe c f = let g = fromJust $ M.lookup f' (genOpFaces !! (cellDim c))
in genToNonGen c g
where f' = nonGenToGen c f
newtype CubeCmplx = CubeCmplx { cells :: S.HashSet CubeCell } deriving (Show,Eq)
instance NFData CubeCmplx where rnf cx = rnf (cells cx)
instance Arbitrary CubeCmplx
where arbitrary = do vs <- arbitrary `suchThat` ((<= 3).vsDim)
let cx = vsCmplx vs
let cs = zip (cycle [1..100]) $ S.toList (cells cx)
return . CubeCmplx . S.fromList . map snd .
filter ((>=10) . fst) $ cs
cmplxEmpty :: CubeCmplx
cmplxEmpty = CubeCmplx { cells = S.empty }
cmplxNull :: CubeCmplx -> Bool
cmplxNull cx = S.null $ cells cx
cmplxSize :: CubeCmplx -> Int
cmplxSize cx = S.size $ cells cx
cmplxApply :: CubeCmplx -> (CubeCell -> S.HashSet CubeCell) -> CubeCmplx
cmplxApply cx f = CubeCmplx . S.unions . map f . S.toList $ cells cx
vsCmplx :: VertSpan -> CubeCmplx
vsCmplx vs = CubeCmplx { cells = S.fromList $ spanTopCells vs }
cmplxDelCell :: CubeCmplx -> CubeCell -> CubeCmplx
cmplxDelCell cx c = CubeCmplx { cells = S.delete c (cells cx) }
cmplxDelCells :: CubeCmplx -> S.HashSet CubeCell -> CubeCmplx
cmplxDelCells cx cs = CubeCmplx { cells = S.difference (cells cx) cs }
cmplxAddCells :: CubeCmplx -> S.HashSet CubeCell -> CubeCmplx
cmplxAddCells cx cs = CubeCmplx { cells = S.union cs (cells cx) }
cmplxUnions :: [CubeCmplx] -> CubeCmplx
cmplxUnions = CubeCmplx . S.unions . map cells
cmplxFilter :: (CubeCell -> Bool) -> CubeCmplx -> CubeCmplx
cmplxFilter f cx = CubeCmplx . S.filter f $ cells cx
cmplxHullUnsafe :: CubeCmplx -> VertSpan
cmplxHullUnsafe cx = vsUnsafe minv maxv
where (f,s) = unzip . map (\c -> (minVert c, maxVert c)) . S.toList $ cells cx
minv = foldr vMin (vertexUnsafe $ replicate (vDim $ head f)
(maxBound :: T)) f
maxv = foldr vMax (vertexUnsafe $ replicate (vDim $ head f)
(minBound :: T)) s
cmplxFilterSpan :: CubeCmplx -> VertSpan -> CubeCmplx
cmplxFilterSpan cx vs = cmplxFilter (flip inSpan vs) cx
cmplxFilterSpans :: CubeCmplx -> [VertSpan] -> [(CubeCmplx, VertSpan)]
cmplxFilterSpans cx vss = withStrategy (parBuffer 100 rdeepseq) $
zip (map (cmplxFilterSpan cx) vss) vss
cellNhd :: CubeCmplx -> CubeCell -> CubeCmplx
cellNhd cx c = cmplxFilterSpan cx $ vsUnsafe minv maxv
where minv = (minVert c) `vSub`
(vertexVectorUnsafe $ V.replicate (vDim (minVert c)) 1)
maxv = (maxVert c) `vAdd`
(vertexVectorUnsafe $ V.replicate (vDim (minVert c)) 1)
swissFlag :: (CubeCmplx, [VertSpan])
swissFlag = (cx, [vsVert $ vertexUnsafe [1,1], vsVert $ vertexUnsafe [6,6]])
where cx = cmplxDelCells (vsCmplx $ vsCoordsUnsafe [1,1] [6,6]) $
S.fromList $ [cellUnsafe [2,3] [3,4], cellUnsafe [3,2] [4,3],
cellUnsafe [3,3] [4,4], cellUnsafe [4,3] [5,4],
cellUnsafe [3,4] [4,5]]
sqPairFwd :: (CubeCmplx, [VertSpan])
sqPairFwd = (cx, [vsVert $ vertexUnsafe [1,1], vsVert $ vertexUnsafe [6,6]])
where cx = cmplxDelCells (vsCmplx $ vsCoordsUnsafe [1,1] [6,6]) $
S.fromList $ [cellUnsafe [2,2] [3,3], cellUnsafe [4,4] [5,5]]
sqPairBack :: (CubeCmplx, [VertSpan])
sqPairBack = (cx, [vsVert $ vertexUnsafe [1,1], vsVert $ vertexUnsafe [6,6]])
where cx = cmplxDelCells (vsCmplx $ vsCoordsUnsafe [1,1] [6,6]) $
S.fromList $ [cellUnsafe [2,4] [3,5], cellUnsafe [4,2] [5,3]]
torus3d :: (CubeCmplx, [VertSpan])
torus3d = (cx, [vsVert $ vertexUnsafe [1,1,1], vsVert $ vertexUnsafe [4,4,2]])
where cx = cmplxDelCells (vsCmplx $ vsCoordsUnsafe [1,1,1] [4,4,2]) $
S.fromList $ [cellUnsafe [2,2,1] [3,3,2]]
genusTwo3d :: (CubeCmplx, [VertSpan])
genusTwo3d = (cx, [vsVert $ vertexUnsafe [1,1,1], vsVert $ vertexUnsafe [4,6,2]])
where cx = cmplxDelCells (vsCmplx $ vsCoordsUnsafe [1,1,1] [4,6,2]) $
S.fromList $ [cellUnsafe [2,2,1] [3,3,2],
cellUnsafe [2,4,1] [3,5,2]]