-- {-# LINE 1 "delaunay.hsc" #-}
{-# LANGUAGE ForeignFunctionInterface #-}
module Geometry.Delaunay.CDelaunay
  ( 
    cTessellationToTessellation
  , c_tessellation 
  )
  where
import           Control.Monad              ( (<$!>) )
import qualified Data.HashMap.Strict.InsOrd as H
import           Data.IntMap.Strict         ( fromAscList, (!) )
import qualified Data.IntSet                as IS
import           Data.List                  ( findIndex )
import           Data.Maybe                 ( fromJust, isJust )
import           Data.Tuple.Extra           ( both, fst3, snd3, thd3, (&&&) )
import           Geometry.Delaunay.Types    ( Tessellation(..),
                                              Tile(..),
                                              TileFacet(..),
                                              Simplex(..),
                                              Site(..) )
import           Foreign  ( Ptr,
                            Storable(pokeByteOff, poke, peek, alignment, sizeOf, peekByteOff),
                            peekArray )
import           Foreign.C.Types            ( CInt, CDouble(..), CUInt(..) )
import           Geometry.Qhull.Types       ( Family(Family, None), IndexPair(Pair) )

data CSite = CSite {
    __id             :: CUInt
  , __neighsites     :: Ptr CUInt
  , __nneighsites    :: CUInt
  , __neighridgesids :: Ptr CUInt
  , __nneighridges   :: CUInt
  , __neightiles     :: Ptr CUInt
  , __nneightiles    :: CUInt
}

instance Storable CSite where
    sizeOf    __ = (56)
-- {-# LINE 29 "delaunay.hsc" #-}
    alignment __ = 8
-- {-# LINE 30 "delaunay.hsc" #-}
    peek ptr = do
      id'              <- (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
-- {-# LINE 32 "delaunay.hsc" #-}
      neighsites'      <- (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
-- {-# LINE 33 "delaunay.hsc" #-}
      nneighsites'     <- (\hsc_ptr -> peekByteOff hsc_ptr 16) ptr
-- {-# LINE 34 "delaunay.hsc" #-}
      neighridgesids'  <- (\hsc_ptr -> peekByteOff hsc_ptr 24) ptr
-- {-# LINE 35 "delaunay.hsc" #-}
      nneighridges'    <- (\hsc_ptr -> peekByteOff hsc_ptr 32) ptr
-- {-# LINE 36 "delaunay.hsc" #-}
      neightiles'      <- (\hsc_ptr -> peekByteOff hsc_ptr 40) ptr
-- {-# LINE 37 "delaunay.hsc" #-}
      nneightiles'     <- (\hsc_ptr -> peekByteOff hsc_ptr 48) ptr
-- {-# LINE 38 "delaunay.hsc" #-}
      return CSite { __id             = id'
                   , __neighsites     = neighsites'
                   , __nneighsites    = nneighsites'
                   , __neighridgesids = neighridgesids'
                   , __nneighridges   = nneighridges'
                   , __neightiles     = neightiles'
                   , __nneightiles    = nneightiles'
                  }
    poke ptr (CSite r1 r2 r3 r4 r5 r6 r7)
      = do
          (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1
-- {-# LINE 49 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2
-- {-# LINE 50 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 16) ptr r3
-- {-# LINE 51 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 24) ptr r4
-- {-# LINE 52 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 32) ptr r5
-- {-# LINE 53 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 40) ptr r6
-- {-# LINE 54 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 48) ptr r7
-- {-# LINE 55 "delaunay.hsc" #-}

cSiteToSite :: [[Double]] -> CSite -> IO (Int, Site, [(Int,Int)])
cSiteToSite sites csite = do
  let id'          = fromIntegral $ __id csite
      nneighsites  = fromIntegral $ __nneighsites csite
      nneighridges = fromIntegral $ __nneighridges csite
      nneightiles  = fromIntegral $ __nneightiles csite
      point        = sites !! id'
  neighsites <- (<$!>) (map fromIntegral)
                       (peekArray nneighsites (__neighsites csite))
  neighridges <- (<$!>) (map fromIntegral)
                        (peekArray nneighridges (__neighridgesids csite))
  neightiles <- (<$!>) (map fromIntegral)
                       (peekArray nneightiles (__neightiles csite))
  return ( id'
         , Site {
                  _point          = point
                , _neighsitesIds  = IS.fromAscList neighsites
                , _neighfacetsIds = IS.fromAscList neighridges
                , _neightilesIds  = IS.fromAscList neightiles
                }
         , map (\j -> (id', j)) (filterAscList id' neighsites) )
  where
    filterAscList :: Int -> [Int] -> [Int]
    filterAscList n list =
      let i = findIndex (> n) list in
      if isJust i
        then drop (fromJust i) list
        else []

data CSimplex = CSimplex {
    __sitesids :: Ptr CUInt
  , __center   :: Ptr CDouble
  , __radius   :: CDouble
  , __volume   :: CDouble
}

instance Storable CSimplex where
    sizeOf    __ = (32)
-- {-# LINE 85 "delaunay.hsc" #-}
    alignment __ = 8
-- {-# LINE 86 "delaunay.hsc" #-}
    peek ptr = do
      sitesids'    <- (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
-- {-# LINE 88 "delaunay.hsc" #-}
      center'      <- (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
-- {-# LINE 89 "delaunay.hsc" #-}
      radius'      <- (\hsc_ptr -> peekByteOff hsc_ptr 16) ptr
-- {-# LINE 90 "delaunay.hsc" #-}
      volume'      <- (\hsc_ptr -> peekByteOff hsc_ptr 24) ptr
-- {-# LINE 91 "delaunay.hsc" #-}
      return CSimplex { __sitesids    = sitesids'
                      , __center      = center'
                      , __radius      = radius'
                      , __volume      = volume'
                    }
    poke ptr (CSimplex r1 r2 r3 r4)
      = do
          (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1
-- {-# LINE 99 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2
-- {-# LINE 100 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 16) ptr r3
-- {-# LINE 101 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 24) ptr r4
-- {-# LINE 102 "delaunay.hsc" #-}

cSimplexToSimplex :: [[Double]] -> Int -> CSimplex -> IO Simplex
cSimplexToSimplex sites simplexdim csimplex = do
  let radius      = cdbl2dbl $ __radius csimplex
      volume      = cdbl2dbl $ __volume csimplex
      dim         = length (head sites)
  sitesids <- (<$!>) (map fromIntegral)
                     (peekArray simplexdim (__sitesids csimplex))
  let points = fromAscList
               (zip sitesids (map ((!!) sites) sitesids))
  center <- (<$!>) (map cdbl2dbl) (peekArray dim (__center csimplex))
  return Simplex { _vertices'       = points
                 , _circumcenter = center
                 , _circumradius = radius
                 , _volume'       = volume }
  where
    cdbl2dbl :: CDouble -> Double
    cdbl2dbl x = if isNaN x then 0/0 else realToFrac x

data CSubTile = CSubTile {
    __id'        :: CUInt
  , __subsimplex :: CSimplex
  , __ridgeOf1   :: CUInt
  , __ridgeOf2   :: CInt
  , __normal     :: Ptr CDouble
  , __offset     :: CDouble
}

instance Storable CSubTile where
    sizeOf    __ = (72)
-- {-# LINE 132 "delaunay.hsc" #-}
    alignment __ = 8
-- {-# LINE 133 "delaunay.hsc" #-}
    peek ptr = do
      id'       <- (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
-- {-# LINE 135 "delaunay.hsc" #-}
      simplex'  <- (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
-- {-# LINE 136 "delaunay.hsc" #-}
      ridgeOf1' <- (\hsc_ptr -> peekByteOff hsc_ptr 40) ptr
-- {-# LINE 137 "delaunay.hsc" #-}
      ridgeOf2' <- (\hsc_ptr -> peekByteOff hsc_ptr 44) ptr
-- {-# LINE 138 "delaunay.hsc" #-}
      normal'   <- (\hsc_ptr -> peekByteOff hsc_ptr 48) ptr
-- {-# LINE 139 "delaunay.hsc" #-}
      offset'   <- (\hsc_ptr -> peekByteOff hsc_ptr 56) ptr
-- {-# LINE 140 "delaunay.hsc" #-}
      return CSubTile { __id'        = id'
                      , __subsimplex = simplex'
                      , __ridgeOf1   = ridgeOf1'
                      , __ridgeOf2   = ridgeOf2'
                      , __normal     = normal'
                      , __offset     = offset' }
    poke ptr (CSubTile r1 r2 r3 r4 r5 r6)
      = do
          (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1
-- {-# LINE 149 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2
-- {-# LINE 150 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 40) ptr r3
-- {-# LINE 151 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 44) ptr r4
-- {-# LINE 152 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 48) ptr r5
-- {-# LINE 153 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 56) ptr r6
-- {-# LINE 154 "delaunay.hsc" #-}

cSubTiletoTileFacet :: [[Double]] -> CSubTile -> IO (Int, TileFacet)
cSubTiletoTileFacet points csubtile = do
  let dim        = length (head points)
      ridgeOf1   = fromIntegral $ __ridgeOf1 csubtile
      ridgeOf2   = fromIntegral $ __ridgeOf2 csubtile
      ridgeOf    = if ridgeOf2 == -1 then [ridgeOf1] else [ridgeOf1, ridgeOf2]
      id'        = fromIntegral $ __id' csubtile
      subsimplex = __subsimplex csubtile
      offset     = realToFrac $ __offset csubtile
  simplex <- cSimplexToSimplex points dim subsimplex
  normal <- (<$!>) (map realToFrac) (peekArray dim (__normal csubtile))
  return (id', TileFacet { _subsimplex = simplex
                         , _facetOf    = IS.fromAscList ridgeOf
                         , _normal'     = normal
                         , _offset'     = offset })

data CTile = CTile {
    __id''        :: CUInt
  , __simplex     :: CSimplex
  , __neighbors   :: Ptr CUInt
  , __nneighbors  :: CUInt
  , __ridgesids   :: Ptr CUInt
  , __nridges     :: CUInt
  , __family      :: CInt
  , __orientation :: CInt
}

instance Storable CTile where
    sizeOf    __ = (80)
-- {-# LINE 184 "delaunay.hsc" #-}
    alignment __ = 8
-- {-# LINE 185 "delaunay.hsc" #-}
    peek ptr = do
      id'         <- (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
-- {-# LINE 187 "delaunay.hsc" #-}
      simplex'    <- (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
-- {-# LINE 188 "delaunay.hsc" #-}
      neighbors'  <- (\hsc_ptr -> peekByteOff hsc_ptr 40) ptr
-- {-# LINE 189 "delaunay.hsc" #-}
      nneighbors' <- (\hsc_ptr -> peekByteOff hsc_ptr 48) ptr
-- {-# LINE 190 "delaunay.hsc" #-}
      ridgesids'  <- (\hsc_ptr -> peekByteOff hsc_ptr 56) ptr
-- {-# LINE 191 "delaunay.hsc" #-}
      nridges'    <- (\hsc_ptr -> peekByteOff hsc_ptr 64) ptr
-- {-# LINE 192 "delaunay.hsc" #-}
      family'     <- (\hsc_ptr -> peekByteOff hsc_ptr 68) ptr
-- {-# LINE 193 "delaunay.hsc" #-}
      orient      <- (\hsc_ptr -> peekByteOff hsc_ptr 72) ptr
-- {-# LINE 194 "delaunay.hsc" #-}
      return CTile { __id''        = id'
                   , __simplex     = simplex'
                   , __neighbors   = neighbors'
                   , __nneighbors  = nneighbors'
                   , __ridgesids   = ridgesids'
                   , __nridges     = nridges'
                   , __family      = family'
                   , __orientation = orient
                  }
    poke ptr (CTile r1 r2 r3 r4 r5 r6 r7 r8)
      = do
          (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1
-- {-# LINE 206 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2
-- {-# LINE 207 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 40) ptr r3
-- {-# LINE 208 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 48) ptr r4
-- {-# LINE 209 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 56) ptr r5
-- {-# LINE 210 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 64) ptr r6
-- {-# LINE 211 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 68) ptr r7
-- {-# LINE 212 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 72) ptr r8
-- {-# LINE 213 "delaunay.hsc" #-}

cTileToTile :: [[Double]] -> CTile -> IO (Int, Tile)
cTileToTile points ctile = do
  let id'        = fromIntegral $ __id'' ctile
      csimplex   = __simplex ctile
      nneighbors = fromIntegral $ __nneighbors ctile
      nridges    = fromIntegral $ __nridges ctile
      family     = __family ctile
      orient     = __orientation ctile
      dim        = length (head points)
  simplex <- cSimplexToSimplex points (dim+1) csimplex
  neighbors <- (<$!>) (map fromIntegral)
                      (peekArray nneighbors (__neighbors ctile))
  ridgesids <- (<$!>) (map fromIntegral)
                      (peekArray nridges (__ridgesids ctile))
  return (id', Tile {  _simplex      = simplex
                     , _neighborsIds = IS.fromAscList neighbors
                     , _facetsIds    = IS.fromAscList ridgesids
                     , _family'       = if family == -1
                                        then None
                                        else Family (fromIntegral family)
                     , _toporiented  = orient == 1 })

data CTessellation = CTessellation {
    __sites     :: Ptr CSite
  , __tiles     :: Ptr CTile
  , __ntiles    :: CUInt
  , __subtiles  :: Ptr CSubTile
  , __nsubtiles :: CUInt
}

instance Storable CTessellation where
    sizeOf    __ = (40)
-- {-# LINE 246 "delaunay.hsc" #-}
    alignment __ = 8
-- {-# LINE 247 "delaunay.hsc" #-}
    peek ptr = do
      sites'     <- (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
-- {-# LINE 249 "delaunay.hsc" #-}
      tiles'     <- (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
-- {-# LINE 250 "delaunay.hsc" #-}
      ntiles'    <- (\hsc_ptr -> peekByteOff hsc_ptr 16) ptr
-- {-# LINE 251 "delaunay.hsc" #-}
      subtiles'  <- (\hsc_ptr -> peekByteOff hsc_ptr 24) ptr
-- {-# LINE 252 "delaunay.hsc" #-}
      nsubtiles' <- (\hsc_ptr -> peekByteOff hsc_ptr 32) ptr
-- {-# LINE 253 "delaunay.hsc" #-}
      return CTessellation {
                     __sites     = sites'
                   , __tiles     = tiles'
                   , __ntiles    = ntiles'
                   , __subtiles  = subtiles'
                   , __nsubtiles = nsubtiles'
                  }
    poke ptr (CTessellation r1 r2 r3 r4 r5)
      = do
          (\hsc_ptr -> pokeByteOff hsc_ptr 0) ptr r1
-- {-# LINE 263 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 8) ptr r2
-- {-# LINE 264 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 16) ptr r3
-- {-# LINE 265 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 24) ptr r4
-- {-# LINE 266 "delaunay.hsc" #-}
          (\hsc_ptr -> pokeByteOff hsc_ptr 32) ptr r5
-- {-# LINE 267 "delaunay.hsc" #-}

foreign import ccall unsafe "tessellation" c_tessellation
  :: Ptr CDouble -- sites
  -> CUInt       -- dim
  -> CUInt       -- nsites
  -> CUInt       -- 0/1, point at infinity
  -> CUInt       -- 0/1, include degenerate
  -> CDouble     -- volume threshold
  -> Ptr CUInt   -- exitcode
  -> IO (Ptr CTessellation)

cTessellationToTessellation :: [[Double]] -> CTessellation -> IO Tessellation
cTessellationToTessellation vertices ctess = do
  let ntiles    = fromIntegral $ __ntiles ctess
      nsubtiles = fromIntegral $ __nsubtiles ctess
      nsites    = length vertices
  sites''    <- peekArray nsites (__sites ctess)
  tiles''    <- peekArray ntiles (__tiles ctess)
  subtiles'' <- peekArray nsubtiles (__subtiles ctess)
  sites'     <- mapM (cSiteToSite vertices) sites''
  let sites = fromAscList (map (fst3 &&& snd3) sites')
      edgesIndices = concatMap thd3 sites'
      edges = map (toPair &&& both (_point . ((!) sites))) edgesIndices
  tiles'     <- mapM (cTileToTile vertices) tiles''
  subtiles'  <- mapM (cSubTiletoTileFacet vertices) subtiles''
  return Tessellation
         { _sites      = sites
         , _tiles      = fromAscList tiles'
         , _tilefacets = fromAscList subtiles'
         , _edges'     = H.fromList edges }
  where
    toPair (i,j) = Pair i j