4
\$\begingroup\$

Here is an implementation of Graham Scan algorithm for computing the convex hull of a finite set of points:

data Point = Point {
  xcoord :: Double,
  ycoord :: Double
} deriving (Show)

data Direction = ToRight | Straight | ToLeft  deriving (Show, Eq)

distance :: Point -> Point -> Double
distance a b = sqrt (((xcoord a) - (xcoord b))^2 + ((ycoord a) - (ycoord b))^2)

-- The direction function can be simplified after some algebraic analysis:
direction :: Point -> Point -> Point -> Direction
direction a b c
  | s < 0 = ToRight
  | s == 0 = Straight
  | s > 0 = ToLeft
  where
    s = ((xcoord c) - (xcoord a)) * ((ycoord a) - (ycoord b)) + ((ycoord c) - (ycoord a)) * ((xcoord b) - (xcoord a))

directions_ :: [Point] -> [Direction]
directions_ [] = []
directions_ [_] = []
directions_ [_, _] = []
directions_ (a:b:c:xs) = (direction a b c) : (directions_ (b:c:xs))

directions :: [Point] -> [Direction]
directions xs = ToLeft : ToLeft : (directions_ xs)


lowestYPoint :: [Point] -> Point
lowestYPoint [x] = x
lowestYPoint (x:xs)
  | (ycoord x) <= (ycoord lowestRest) = x
  | otherwise = lowestRest
  where
    lowestRest = lowestYPoint xs


cosineAB :: Point -> Point -> Double
cosineAB a b = ((xcoord b) - (xcoord a)) / (distance a b)

sortByAngle :: Point -> [Point] -> [Point]
sortByAngle x xs = sortBy comparePoints xs where
  comparePoints :: Point -> Point -> Ordering
  comparePoints a b
    | cA > cB = LT
    | cA < cB = GT
    | otherwise = case () of
      _ | xA < xB  -> LT
        | xA > xB  -> GT
        | xA == xB -> case () of
          _ | yA < yB   -> LT
            | yA > yB   -> GT
            | yA == yB  -> EQ
    where
      yA = (ycoord a)
      yB = (ycoord b)
      xA = (xcoord a)
      xB = (xcoord b)
      cA = (cosineAB x a)
      cB = (cosineAB x b)

grahamScanDirections :: [Direction] -> [Bool]
grahamScanDirections [] = []
grahamScanDirections [_] = [True]
grahamScanDirections (ToLeft:ToRight:xs) = False : (grahamScanDirections (ToLeft:xs))
grahamScanDirections (ToLeft:ToLeft:xs)  = True  : (grahamScanDirections (ToLeft:xs))
grahamScanDirections (ToLeft:Straight:xs) = False : (grahamScanDirections (ToLeft:xs))
grahamScanDirections _ = error "Impossible scenario"

filterByList :: [Bool] -> [a] -> [a]
filterByList [] _ = []
filterByList _ [] = []
filterByList (x:xs) (y:ys)
  | x = y : (filterByList xs ys)
  | otherwise = (filterByList xs ys)

grahamScan :: [Point] -> [Point]
grahamScan ps = filterByList bool_list ps1 where
  -- sort ps by their angles (or cosine values) formed with the lowest point
  ps1 = sortByAngle (lowestYPoint ps) ps
  bool_list = grahamScanDirections $ directions ps1

Test:

*Main> let ps = [(Point 1 0), (Point 1 2), (Point 0 3), (Point 2 3), (Point 2 2), (Point 3 2), (Point 2 1)]
*Main> grahamScan ps
[Point {xcoord = 1.0, ycoord = 0.0},Point {xcoord = 3.0, ycoord = 2.0},Point {xcoord = 2.0, ycoord = 3.0},Point {xcoord = 0.0, ycoord = 3.0}]

Any suggestion about correctness/performance or any other possible improvement is appreciated.

\$\endgroup\$
2
\$\begingroup\$

I almost always prefer case analysis to records and guards.‡ E.g.,

direction :: Point -> Point -> Point -> Direction
direction (Point ax ay) (Point bx by) (Point cx cy) =
    let
        s = (cx - ax) * (ay - by) + (cy - ay) * (bx - ax)
    in
        case s `compare` 0 of
            LT -> ToRight
            EQ -> Straight
            GT -> ToLeft

You have some redundant cases in the definition of directions_. My use of the as-pattern below might be a little too much flourish, but technically it saves you a few list construction operations in putting b and c back on the front of your version of xs.

directions_ :: [Point] -> [Direction]
directions_ (a:xs@(b:c:_)) = direction a b c : directions_ xs
directions_              _ = []

You could condense your code using some higher order functions from base. E.g.,

import Data.Function (on)
import Data.List (minimumBy)

lowestYPoint :: [Point] -> Point
lowestYPoint = minimumBy (compare `on` yCoord) -- Or `comparing` from `Data.Ord`

sortByAngle is a little abusive of case statements and guards, plus you end up performing a lot of comparisons when you should really just go by the ordering compare gives you back.

-- data Point ... deriving Ord
-- The default `Ord` instance works for your purposes here
--   but you may want to explicitly `instance Ord Point ...`

sortByAngle :: Point -> [Point] -> [Point]
sortByAngle p ps = sortBy comparePoints ps
    where
        comparePoints :: Point -> Point -> Ordering
        comparePoints a b =
            let
                cosPA = cosineAB p a
                cosPB = cosineAB p b
            in
                case cosPA `compare` cosPB of
                    LT -> GT
                    GT -> LT
                    EQ -> a `compare` b

filterByList is particularly smelly, pattern match those booleans!

filterByList :: [Bool] -> [a] -> [a]
filterbyList (True :xs) (y:ys) = y : filterByList xs ys
filterByList (False:xs) (_:ys) =     filterbyList xs ys
filterByList -- ...

I'd not write my own version of filter though.

filterByList :: [Bool] -> [a] -> [a]
filterByList ps as = map snd . filter fst $ zip ps as

‡ Why I prefer pattern matching to guards

  1. Precision
  2. Efficiency
  3. Consistency

Pattern matching repeats what is important about the value we're inspecting. The constructors and values we care about are directly on display, and the irrelevant portions are shown to be unimportant by naming them _. Consider the difference between these two fragments.

data Shape = Circle | Square | Rectangle Int Int
isSquare :: Shape -> Bool

func1 :: Shape -> a
func1 s | isSquare s = ...
        | otherwise  = ...

func2 :: Shape -> a
func2 Square                      = ...
func2 _                           = ...

Now, are func1 and func2 equivalent? That depends on whether isSquare checks to see if you've been passed a Rectangle with equal length sides.

If you compile with -fwarn-incomplete-patterns GHC will tell you about partial pattern matches in your functions as well, but if you forget to include a guard you're on your own. E.g., my version of direction will not compile if you don't include the EQ case, your version will fail at runtime if you leave out s == 0.

In the example of direction, the guard version will perform three comparisons on s. The default definition of Ord defines all of the comparators (<, <=, &c) in terms of compare, so in the worst case you'll have to call compare three times to figure out which guard matches. In the case version compare is called only once and its result gets reused, potentially saving a lot of time if the implementation of compare were costly.

And lastly, having a set style and preference leads to consistency, which is always a boon to readability. I did it, therefore I do it.

\$\endgroup\$
  • \$\begingroup\$ Note: though the last filterByList function is more elegant, it is unfortunately slower than the ugly one. \$\endgroup\$ – zigazou Mar 4 '15 at 7:27
  • \$\begingroup\$ Do you have a benchmark for that? Because it seems to be the other way around. Here is a gist compiling both versions with -O2 -ddump-simpl, you can see the Core for the elegant version turns into a tidy little foldr. There's a little teeny benchmark in there too, the elegant version runs >5 times faster. \$\endgroup\$ – bisserlis Mar 4 '15 at 19:40
  • \$\begingroup\$ Here are my outputs pastebin.com/5yxFDHJD As you can see the cost centre of filterByList increases when using the cleaner version. \$\endgroup\$ – zigazou Mar 5 '15 at 8:04
  • \$\begingroup\$ Dunno what to tell you, I'm using GHC 7.6.3 so it must be a compiler regression. I don't really understand your profiling output either because both reports reference graham.opt (lines 4 and 29) so it looks like you ran the same test twice. \$\endgroup\$ – bisserlis Mar 5 '15 at 13:57
  • \$\begingroup\$ You can download my tests here ouep.eu/pic/83125.zip just run make prof in the directory to compile and run the profiling. Concerning the profilings both saying graham.opt, it's only because I hand modified the same file to test the two functions. I am interested to see if it's a GHC regression from 7.6.3 to 7.8.3. Thanks :-) \$\endgroup\$ – zigazou Mar 6 '15 at 6:01
2
\$\begingroup\$

Here are some improvements you may apply.

You can combine pattern matching and records. It allows you to avoid unnecessary numerous calls to xcoord and ycoord. It also makes your code easier to read.

You can force strictness on xcoord and ycoord fields.

These two improvements help with performance.

You can also do the following:

  • simplify pattern matching py reordering the patterns (see filterByList)
  • drop lots of parenthesis since space has bigger precedence than operators
  • simplify comparePoints with keeping only one level of guards
  • use more standard library functions (see lowestYPoint)

Most of these were given to me by hlint ;-)

import Data.List
import Data.Function (on)
import Control.DeepSeq

data Point = Point {
    xcoord :: !Double,
    ycoord :: !Double
} deriving Show

data Direction = ToRight | Straight | ToLeft  deriving (Show, Eq)

distance :: Point -> Point -> Double
distance (Point xA yA) (Point xB yB) = sqrt ((xA - xB)^2 + (yA - yB)^2)

direction :: Point -> Point -> Point -> Direction
direction (Point xA yA) (Point xB yB) (Point xC yC)
    | s <  0    = ToRight
    | s == 0    = Straight
    | otherwise = ToLeft
    where
        s = (xC - xA) * (yA - yB) + (yC - yA) * (xB - xA)

directions_ :: [Point] -> [Direction]
directions_ (a:b:c:xs) = direction a b c : directions_ (b:c:xs)
directions_ _ = []

directions :: [Point] -> [Direction]
directions xs = ToLeft : ToLeft : directions_ xs

lowestYPoint :: [Point] -> Point
lowestYPoint = minimumBy (compare `on` ycoord)

cosineAB :: Point -> Point -> Double
cosineAB a@(Point xA _) b@(Point xB _) = (xB - xA) / distance a b

comparePoints :: Point -> Point -> Point -> Ordering
comparePoints x a@(Point xA yA) b@(Point xB yB)
    | cA > cB   = LT
    | cA < cB   = GT
    | xA < xB   = LT
    | xA > xB   = GT
    | yA < yB   = LT
    | yA > yB   = GT
    | otherwise = EQ
    where
        cA = cosineAB x a
        cB = cosineAB x b

sortByAngle :: Point -> [Point] -> [Point]
sortByAngle x = sortBy (comparePoints x)

grahamScanDirections :: [Direction] -> [Bool]
grahamScanDirections [] = []
grahamScanDirections [_] = [True]
grahamScanDirections (ToLeft:ToRight:xs) = False : grahamScanDirections (ToLeft:xs)
grahamScanDirections (ToLeft:ToLeft:xs)  = True  : grahamScanDirections (ToLeft:xs)
grahamScanDirections (ToLeft:Straight:xs) = False : grahamScanDirections (ToLeft:xs)
grahamScanDirections _ = error "Impossible scenario"

filterByList :: [Bool] -> [a] -> [a]
filterByList (True :xs) (y:ys) = y : filterByList xs ys
filterByList (False:xs) (_:ys) = filterByList xs ys
filterByList _ _               = []

grahamScan :: [Point] -> [Point]
grahamScan ps = filterByList bool_list ps1
    where ps1       = sortByAngle (lowestYPoint ps) ps
          bool_list = grahamScanDirections $ directions ps1

pss = replicate 1000000
                [ Point 1 0
                , Point 1 2
                , Point 0 3
                , Point 2 3
                , Point 2 2
                , Point 3 2
                , Point 2 1
                ]

instance NFData Point where
    rnf a = a `seq` () 

main = (map grahamScan pss) `deepseq` putStrLn "Done"
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.