QuickHullIO.hs

{-# LANGUAGE BangPatterns, PatternGuards, RankNTypes #-}

module QuickHullIO
	(quickHull)
where
import Data.Function
import Control.Monad
import Control.Exception
import Control.Concurrent
import Control.Concurrent.MVar
import Control.Monad.ST
import GHC.Conc
import Data.IORef
import Data.List
import Data.Ord
import Data.Vector.Unboxed			(Vector)
import qualified Data.Vector.Unboxed		as V
import qualified Data.Vector.Unboxed.Mutable	as MV
import qualified Data.Vector.Generic		as G
import Debug.Trace

type Point	= (Double, Double)
type Line	= (Point, Point)


-- | Compute the convex hull of a vector of points.
quickHull :: Vector Point -> IO (Vector Point)
quickHull !points
  | V.length points == 0	
  = return points

  | otherwise
  = do	-- Find the left and right-most points.
	let (minx, maxx) 	= minmax points

	-- Hull points get written to the vector in this IORef.
	hullRef	<- newIORef V.empty

	-- Fork off computations to handle half of the points each.
	-- For uniformly distributed points this first iteration takes most of the time.
	parIO 	[ hsplit hullRef points minx maxx
		, hsplit hullRef points maxx minx]

	-- Grab the finished hull points.
	hull	<- readIORef hullRef

	-- We've got the hull points, but they can appear in arbitrary order.
	-- Do a rubbish via-lists merge phase so that they appear clockwise around the edge.
	-- This isn't too expensive if there aren't many points on the hull.
	let (above, below) 
		= V.unstablePartition 
			(\p -> distance minx maxx p > 0)
			hull
	
	let aboveSorted	= V.fromList $ sortBy (comparing fst) $ V.toList above
	let belowSorted	= V.fromList $ sortBy (comparing fst) $ V.toList below
	let hull' = aboveSorted V.++ V.reverse belowSorted

	return hull'
	

hsplit :: IORef (Vector Point) -> Vector Point -> Point -> Point -> IO ()
{-# INLINE hsplit #-}
hsplit hullRef !points !p1@(!p1X, !p1Y) !p2@(!p2X, !p2Y)
	-- we've found one.
	| V.length packed == 0
	= addHullPoint hullRef p1
	
	-- do the two new segments in parallel.
	| V.length packed > 1000
	= parIO
		[ hsplit hullRef packed p1 pm
		, hsplit hullRef packed pm p2 ]
		
	| otherwise
	= do	hsplit hullRef packed p1 pm
		hsplit hullRef packed pm p2

	where	(packed, pm)	= parPackPoints points p1X p1Y p2X p2Y
	

-- | Copy points from the input vector that are on the left of the line into a
--	new buffer. While we're doing this, determine the point that is furthest
--	from the line.
--
--	If we have a big enough vector then split it in two and do both halves
--	in parallel. Doing this requires a copy afterwards to join the two
--	results back together. It's a trade off between decreased FP load and 
--	increased memory traffic. 
--
parPackPoints 
	:: Vector Point 
	-> Double -> Double
	-> Double -> Double
	-> ( Vector Point
	   , Point)
	
{-# INLINE parPackPoints #-}
parPackPoints !points !p1X !p1Y !p2X !p2Y
 |   numCapabilities == 1
  || V.length points < 1000
 = packPoints p1X p1Y p2X p2Y points

 | otherwise
 = let	
	numSegments	= numCapabilities

	-- Total number of points to process.
	lenPoints	= V.length points

	-- How many points to process in each segment.
	lenSeg		= lenPoints `div` numSegments

 	-- If the total number of points doesn't divide evenly into segments
	-- then there may be an odd number. Make sure to get the rest into the last segment.
	splitPacked count ixStart 
	    | count == 0	= []

	    | count == 1	
	    = let points'		= V.unsafeSlice ixStart (lenPoints - ixStart) points
	      	  result@(packed', _)	= packPoints p1X p1Y p2X p2Y points'
	      in  packed' `pseq` (result : [])

	    | otherwise	
	    = let points'		= V.unsafeSlice ixStart lenSeg points
	          result@(packed', _)	= packPoints p1X p1Y p2X p2Y points'
		  rest			= splitPacked (count - 1) (ixStart + lenSeg)
	      in  packed' `par` rest `par` (result : rest)

	results	= splitPacked numSegments 0
	vResult	= concatVectors $ map fst results
	pMax	= selectFurthest p1X p1Y p2X p2Y results
	
   in	(vResult, pMax)


selectFurthest 
 	:: Double -> Double 
	-> Double -> Double
	-> [(Vector Point, Point)] 
	-> Point
	
selectFurthest !p1X !p1Y !p2X !p2Y ps
 = go (0, 0) 0 ps

 where	go pMax !distMax []	
	 = pMax

	go pMax !distMax ((packed, pm):rest)
	 | V.length packed == 0
	 = go pMax distMax rest
	
	 | otherwise
	 , dist		<-  distance (p1X, p1Y) (p2X, p2Y) pm 
  	 = if dist > distMax
		then go pm   dist    rest
		else go pMax distMax rest


packPoints 
	:: Double -> Double 		-- First point on dividing line.
	-> Double -> Double 		-- Second point on dividing line.
	-> Vector Point 		-- Source points.
	-> ( Vector Point		-- Packed vector containing only points on the left of the line.
	   , Point)			-- The point on the left that was furthest from the line.

{-# INLINE packPoints #-}
packPoints !p1X !p1Y !p2X !p2Y !points 
 = let
	result	
	 = G.create 
 	 $ do	packed		 <- MV.new (V.length points + 1)
		(pMax, ixPacked) <- fill points packed p1X p1Y p2X p2Y 0 0

		-- We stash the maximum point on the end of the vector to get
		-- it through the create call.
		MV.unsafeWrite packed ixPacked pMax
		return $ MV.unsafeSlice 0 (ixPacked + 1) packed
	
   in	( V.unsafeSlice 0 (V.length result - 1) result
	, result `V.unsafeIndex` (V.length result - 1))
			

fill 	:: forall s
	.  Vector Point 		-- Source points.
	-> MV.MVector s Point 		-- Vector to write packed points into.
	-> Double -> Double 		-- First point on dividing line.
	-> Double -> Double		-- Second poitn on dividing line.
	-> Int 				-- Index into source points to start reading from.
	-> Int				-- Index into packed points to start writing to.
	-> ST s 
		( Point			-- Furthest point from the line that was found.
		, Int)			-- The number of packed points written.

{-# INLINE fill #-}
fill !points !packed !p1X !p1Y !p2X !p2Y !ixPoints' !ixPacked'
 = go (0, 0) 0 ixPoints' ixPacked'
 where go pMax !distMax !ixPoints !ixPacked
	| ixPoints >= V.length points	
	= do	return (pMax, ixPacked)
		
	| p	<- points `V.unsafeIndex` ixPoints
	, d	<- distance (p1X, p1Y) (p2X, p2Y) p
	, d > 0
	= do	MV.unsafeWrite packed ixPacked p
		if d > distMax
		 then	go p    d       (ixPoints + 1) (ixPacked + 1)
		 else	go pMax distMax (ixPoints + 1) (ixPacked + 1)
			
	| otherwise
	= go pMax distMax (ixPoints + 1) ixPacked


minmax :: Vector Point -> (Point, Point)
{-# INLINE minmax #-}
minmax !vec
 = go first first 0
 where	first	= vec V.! 0

	go pMin@(!minX, !minY) pMax@(!maxX, !maxY) !ix
	  | ix >= V.length vec	= (pMin, pMax)

	  | (x, y)	<- vec `V.unsafeIndex` ix
	  = if       x < minX then go (x, y) pMax   (ix + 1)
	    else if  x > maxX then go pMin   (x, y) (ix + 1)
	    else go pMin pMax (ix + 1)
	

distance :: Point -> Point -> Point -> Double
{-# INLINE distance #-}
distance (x1, y1) (x2, y2) (xo, yo)
  = (x1-xo) * (y2 - yo) - (y1 - yo) * (x2 - xo)


addHullPoint :: IORef (Vector Point) -> Point -> IO ()
addHullPoint hullRef p
 = atomicModifyIORef hullRef
 $ \hull -> (V.singleton p V.++ hull, ())



-- Can't find an equivalent for this in Control.Concurrent.
parIO :: [IO ()] -> IO ()
parIO stuff
 = do	mVars	<- replicateM (length stuff) newEmptyMVar
	zipWithM_ (\c v -> forkIO $ c `finally` putMVar v ()) stuff mVars
	mapM_ readMVar mVars
	

-- We really want a function in the vector library for this.
concatVectors :: [Vector Point] -> Vector Point
{-# NOINLINE concatVectors #-}
concatVectors vectors
 = G.create
 $ do	let len	= sum $ map V.length vectors
	vOut	<- MV.new len
	go vectors vOut 0
	return vOut

 where	{-# INLINE go #-}
	go [] _ _
	 = return ()

	go (vSrc:vsSrc) vDest !ixStart	
	 = do	let lenSrc	= V.length vSrc
		let vDestSlice	= MV.unsafeSlice ixStart lenSrc vDest
		V.copy vDestSlice vSrc
		go vsSrc vDest (ixStart + lenSrc)