module Math.Topology.CubeCmplx.CornerReduce (
cmplxCornersNaive, cmplxSpanIntCorners, cmplxCorners, cmplxCornersInt,
cmplxReduce', cmplxReduce
) where
import Data.List (transpose, groupBy, sortBy, partition)
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 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 Control.Parallel.Strategies
(rdeepseq, parBuffer, withStrategy, parList, dot, evalTuple3, r0)
import Control.Arrow ((***))
import Math.Topology.CubeCmplx.DirCubeCmplx
cmplxCornersNaive :: CubeCmplx -> S.HashSet (Vertex, CubeCell)
cmplxCornersNaive cx = S.fromList . map (\(v,(c,_)) -> (v,c)) . M.toList .
M.filter ((==1) . snd) $
M.fromListWith (\(c,n1) (_,n2) -> (c,n1+n2)) vcs
where cs = S.toList $ cells cx
vcs = concatMap (\c -> zip (verts c) (zip (repeat c) (repeat 1))) cs
cmplxCovSpans :: CubeCmplx -> [VertSpan]
cmplxCovSpans cx
| cmplxNull cx = []
| otherwise = vsCov3 $ cmplxCovHullUnsafe cx
where vsCov3 vs = map (\e -> vsUnsafe (vertexUnsafe $ map fst e)
(vertexUnsafe $ map snd e)) $ rProds vs
ranges vs = zip (vsFstList vs) (vsSndList vs)
f (n,m) = [(i,i+3) | i <- [n..m], i+3 <= m]
rProds vs = lazyProd $ map f (ranges vs)
cmplxSpanIntCorners :: CubeCmplx -> VertSpan -> S.HashSet (Vertex, CubeCell)
cmplxSpanIntCorners cx vs = S.filter (\(v,_) -> not $ (vertToCell v) `inBdry` vs)
(cmplxCornersNaive cx)
cmplxCorners :: CubeCmplx -> S.HashSet (Vertex, CubeCell)
cmplxCorners cx = S.unions . withStrategy (parBuffer 100 rdeepseq) .
map (uncurry cmplxSpanIntCorners) $
cmplxFilterSpans cx (cmplxCovSpans cx)
cmplxCornersInt :: CubeCmplx -> VertSpan -> S.HashSet (Vertex, CubeCell)
cmplxCornersInt cx vs
= S.filter (\(v,_) -> (v `vInSpan` vs) &&
(not $ (vertToCell v) `inBdry` vs)) (cmplxCorners cx)
cmplxCovHullUnsafe :: CubeCmplx -> VertSpan
cmplxCovHullUnsafe cx = vsUnsafe f' s'
where (f,s) = (vsFst $ cmplxHullUnsafe cx, vsSnd $ cmplxHullUnsafe cx)
d = s `vSubUnsafe` f
f' = f `vSub` (vertexUnsafe $ replicate (vDim f) 1)
s' = vertexPtWise (\s1 d1 -> if d1 < 2 then s1+(2d1) else s1+1) s d
genLx :: Int -> M.HashMap Vertex [(Vertex, CubeCell)]
genLx n = gensLx !! n
gensLx = map genLx' [0..] where
genLx' n = M.fromList . zip vs . zipWith zip (map xs vs) .
map (uncurry (zipWith cellVertsUnsafe)) $ zip (map ls vs) (map us vs)
where z = V.replicate n (0 :: T)
o = V.replicate n (1 :: T)
vs = map minVert $ nCubeVerts n
nx = V.toList . V.zip (V.enumFromN 0 n) . V.map (\e -> 1e) . coords
ls = map (vertexVectorUnsafe . V.update z . V.singleton) . nx
us = map (vertexVectorUnsafe . V.update o . V.singleton) . nx
xs = \v -> map (vertexVectorUnsafe . V.update (coords v)
. V.singleton) (nx v)
cellLx :: CubeCell -> Vertex -> S.HashSet (Vertex, CubeCell)
cellLx c x = let cs = case M.lookup gx $ genLx (cellDim c) of
Nothing -> S.empty
Just cs -> S.fromList cs in
S.map (\(v1,c1) -> (minVert $ genToNonGen c (vertToCell v1),
genToNonGen c c1)) cs
where gx = minVert $ nonGenToGen c (vertToCell x)
cmplxCornerUpdate :: CubeCmplx -> (Vertex, CubeCell)
-> (S.HashSet (Vertex, CubeCell),
S.HashSet CubeCell,
S.HashSet CubeCell)
cmplxCornerUpdate cx (x,c) = (cvs, newTops, notTops)
where delNhd = cmplxDelCell (cellNhd cx c) c
delNhdcs = S.toList $ cells delNhd
(notTops,newTops) = S.fromList *** S.fromList $
partition (\c -> any (==True) $
map (isSubCell c) delNhdcs) $
S.toList . S.map snd $ cellLx c x
cvs = S.filter (\(v1,_) -> (vertToCell v1) `isSubCell` c) .
cmplxCornersNaive $ cmplxAddCells delNhd newTops
uniqueCorners :: S.HashSet (Vertex, CubeCell) -> S.HashSet (Vertex, CubeCell)
uniqueCorners vcs = S.fromList . map (\(a,b) -> (b,a)) . M.toList $
S.foldr step M.empty vcs
where step (v,c) m = M.insertWith (\v1 v2 -> v1) c v m
cmplxCornersDelSerial ::
[VertSpan]
-> (CubeCmplx, S.HashSet (Vertex, CubeCell))
-> (CubeCmplx, S.HashSet (Vertex, CubeCell))
cmplxCornersDelSerial vs (cx,xs) = S.foldr step (cx,xs') corners
where vFilt = S.filter (\(v,_) -> not . any (==True) $
map (vInSpan v) vs)
xs' = vFilt xs
corners = uniqueCorners xs'
step (v,c) (cx',cs) = (if cellDim c == 0 then cx'
else cmplxAddCells (cmplxDelCell cx' c) $
(newTops cx' (v,c)),
vFilt $ if cellDim c == 0 then cs
else S.union (delVerts cs c)
(newCorn cx' (v,c)))
update cx' (v,c) = cmplxCornerUpdate cx' (v,c)
newTops cx' (v,c) = (\(_,s,_) -> s) $ update cx' (v,c)
newCorn cx' (v,c) = (\(f,_,_) -> f) $ update cx' (v,c)
delVerts cs c = S.difference cs . S.fromList $
zip (verts c) (repeat c)
cmplxReduceSerial :: CubeCmplx -> S.HashSet (Vertex, CubeCell) -> [VertSpan]
-> CubeCmplx
cmplxReduceSerial cx xs vs
| xs == xs' = cx'
| otherwise = cmplxReduceSerial cx' xs' vs
where (cx',xs') = cmplxCornersDelSerial vs (cx,xs)
cmplxCornerUpdates ::
CubeCmplx -> [(Vertex, CubeCell)]
-> Maybe (S.HashSet (Vertex, CubeCell),
S.HashSet CubeCell,
S.HashSet CubeCell)
cmplxCornerUpdates cx xss = if paraPoss
then Just (S.unions xs, S.unions cs, S.unions os)
else Nothing
where updates = map (cmplxCornerUpdate cx) xss
(xs,cs,os) = unzip3 updates
paraPoss = M.null . M.filter (>= 2) . M.fromListWith (+) $
zip (concatMap (verts . snd) xss) (repeat 1)
cmplxCornersDelPar :: Int
-> [VertSpan]
-> (CubeCmplx, S.HashSet (Vertex, CubeCell))
-> (CubeCmplx, S.HashSet (Vertex, CubeCell))
cmplxCornersDelPar frac vs (cx,xs)
= case cmplxCornerUpdates cx (S.toList corners) of
Just (ncs,nts,_) -> (cmplxAddCells (cmplxDelCells cx
(S.filter ((/=0).cellDim) .
S.map snd $ corners)) nts,
vFilt $ S.filter (\(v,c) -> (cellDim c == 0) ||
not (v `S.member` delVerts)) xs'
`S.union` ncs)
Nothing -> cmplxCornersDelSerial vs (cx,xs)
where vFilt = S.filter (\(v,_) -> not . any (==True) $
map (vInSpan v) vs)
xs' = vFilt xs
corners = S.fromList . cStrat . S.toList $ uniqueCorners xs'
delVerts = S.fromList . concatMap (verts.snd) $ S.toList corners
cStrat = map snd . filter ((==1).fst) . zip (cycle [1..frac])
cmplxReduceIter' :: Int
-> CubeCmplx
-> S.HashSet (Vertex, CubeCell)
-> [VertSpan]
-> [(CubeCmplx, S.HashSet (Vertex, CubeCell), [VertSpan])]
cmplxReduceIter' frac cx xs vs = iterate go (cx,xs,vs) where
go (cx',xs',vs') = (ncx,nxs,vs)
where (ncx,nxs) = cmplxCornersDelPar frac vs (cx',xs')
cmplxReduceIter :: Int
-> CubeCmplx
-> S.HashSet (Vertex, CubeCell)
-> [VertSpan]
-> CubeCmplx
cmplxReduceIter frac cx xs vs
| cmplxNull cx = cx
| otherwise = (\(ncx,_,_) -> ncx) . fst . head .
dropWhile (\((_,xs,_),(_,ys,_)) -> xs /= ys) $
zip rcx (tail rcx)
where rcx = cmplxReduceIter' frac cx xs vs
disjointCov :: VertSpan -> [VertSpan]
disjointCov vs = map (buildvs . transpose) . lazyProd .
map f $ zip (vsFstList vs) (vsSndList vs)
where f (i,j) | j <= i+2 = [[i,j]]
| otherwise = [[i,i+(ji) `div` 2], [(i+(ji) `div` 2)+1, j]]
buildvs [t1,t2] = vsCoordsUnsafe t1 t2
rSubProb :: CubeCmplx -> VertSpan
-> (CubeCmplx, S.HashSet (Vertex, CubeCell), [VertSpan])
rSubProb cx vs = (fatCmplx, corners, vsBdry vs)
where subCmplx = cmplxFilterSpan cx vs
corners = cmplxCornersInt subCmplx vs
c = head $ spanTopCells vs
d = (maxVert c) `vSub` (minVert c)
fatCmplx = cmplxFilterSpan cx $ vsUnsafe ((vsFst vs) `vSub` d)
((vsSnd vs) `vAdd` d)
rSubProbs :: CubeCmplx -> [(CubeCmplx, S.HashSet (Vertex, CubeCell), [VertSpan])]
rSubProbs cx
| cmplxNull cx = []
| otherwise = map (rSubProb cx) . disjointCov . vsFatten .
cmplxHullUnsafe $ cx
cmplxReduce' :: CubeCmplx -> [VertSpan] -> [CubeCmplx]
cmplxReduce' cx vs = iterate go cx where
go cx' = serStep . parStep $ cx'
where f (cx',xs',vs') = cmplxReduceIter 10 cx' xs' (vs ++ vs')
parStep = cmplxUnions . withStrategy (parList rdeepseq) .
map f . rSubProbs
serStep cx = (\(x,_,_) -> x) . head . drop 1 $
cmplxReduceIter' 10 cx (cmplxCorners cx) vs
cmplxReduce :: CubeCmplx -> [VertSpan] -> CubeCmplx
cmplxReduce cx vs
| cmplxNull cx = cx
| otherwise = fst . head . dropWhile (\(c,d) -> c /= d) $
zip rcx (tail rcx)
where rcx = cmplxReduce' cx vs