Difference between revisions of "Graham Scan Implementation"

From HaskellWiki
Jump to: navigation, search
m (Testing: add graph of points)
m (fixup Pt definition)
 
(3 intermediate revisions by the same user not shown)
Line 2: Line 2:
   
 
<haskell>
 
<haskell>
import Data.Ord (comparing)
 
 
import Data.List (sortBy)
 
import Data.List (sortBy)
--Graham Scan exercise
 
  +
import Data.Function (on)
   
--Direction type
 
  +
data Direction = LEFT | RIGHT | STRAIGHT
data Direction = LeftTurn
 
  +
deriving (Show, Eq)
| RightTurn
 
| Straight
 
deriving (Show, Eq)
 
   
--Point type
 
  +
data Pt = Pt (Double, Double)
data Point = Point (Double, Double)
 
  +
deriving (Show, Eq, Ord)
deriving (Show)
 
   
--some points
 
  +
isTurned :: Pt -> Pt -> Pt -> Direction
p0 = Point (2.1,2.0)
 
  +
isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of
p1 = Point (4.2,2.0)
 
  +
EQ -> STRAIGHT
p2 = Point (0.5,2.5)
 
  +
LT -> RIGHT
p3 = Point (3.2,3.5)
 
  +
GT -> LEFT
p4 = Point (1.2,4.0)
 
  +
where sign = compare ((bx - ax) * (cy - ay)) ((cx - ax) * (by - ay))
p5 = Point (0.7,4.7)
 
p6 = Point (1.0,1.0)
 
p7 = Point (3.0,5.2)
 
p8 = Point (4.0,4.0)
 
p9 = Point (3.5,1.5)
 
pA = Point (0.5,1.0)
 
points = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,pA]
 
   
-- Actually, I'd leave it as EQ, GT, LT. Then, actually,
 
  +
gScan :: [Pt] -> [Pt]
-- if you wanted to sort points rotationally around a single point,
 
  +
gScan pts
-- sortBy (dir x) would actually work. --wasserman.louis@gmail.com
 
  +
| length pts >= 3 = scan [pt0] rests
--Get direction of single set of line segments
 
  +
| otherwise = pts
dir :: Point -> Point -> Point -> Direction
 
  +
where
dir (Point (ax, ay)) (Point (bx, by)) (Point (cx, cy)) = case sign of
 
  +
-- Find the most bottom-left point pt0
EQ -> Straight
 
  +
pt0 = foldr bottomLeft (Pt (1/0, 1/0)) pts where
GT -> LeftTurn
+
bottomLeft pa pb = case ord of
LT -> RightTurn
+
LT -> pa
where sign = compare ((bx - ax) * (cy - ay)) ((by - ay) * (cx - ax))
+
GT -> pb
  +
EQ -> pa
  +
where ord = (compare `on` (\ (Pt (x, y)) -> (y, x))) pa pb
   
--Get a list of Directions from a list of Points
 
  +
-- Sort other points based on angle
dirlist :: [Point] -> [Direction]
 
  +
rests = tail (sortBy (compare `on` compkey pt0) pts) where
dirlist (x:y:z:xs) = dir x y z : dirlist (y:z:xs)
 
  +
compkey (Pt (x0, y0)) (Pt (x, y)) = (atan2 (y - y0) (x - x0),
dirlist _ = []
 
  +
abs (x - x0))
   
--Compare Y axes
 
  +
-- Scan the points to find out convex
sortByY :: [Point] -> [Point]
 
  +
-- -- handle the case that all points are collinear
sortByY xs = sortBy lowestY xs
 
  +
scan [p0] (p1:ps)
where lowestY (Point(x1,y1)) (Point (x2,y2)) = compare (y1,x1) (y2,x2)
 
  +
| isTurned pz p0 p1 == STRAIGHT = [pz, p0]
--get COT of line defined by two points and the x-axis
 
  +
where pz = last ps
pointAngle :: Point -> Point -> Double
 
pointAngle (Point (x1, y1)) (Point (x2, y2)) = (x2 - x1) / (y2 - y1)
 
   
--compare based on point angle
 
  +
scan (x:xs) (y:z:rsts) = case isTurned x y z of
pointOrdering :: Point -> Point -> Ordering
 
  +
RIGHT -> scan xs (x:z:rsts)
pointOrdering a b = compare (pointAngle a b) 0.0
 
  +
STRAIGHT -> scan (x:xs) (z:rsts) -- skip collinear points
  +
LEFT -> scan (y:x:xs) (z:rsts)
   
--Sort by angle
 
  +
scan xs [z] = z : xs
sortByAngle :: [Point] -> [Point]
 
sortByAngle ps = bottomLeft : sortBy (compareAngles bottomLeft) (tail (sortedPs))
 
where sortedPs = sortByY ps
 
bottomLeft = head (sortedPs)
 
 
 
 
--Compare angles
 
compareAngles :: Point -> Point -> Point -> Ordering
 
compareAngles = comparing . pointAngle
 
 
--Graham Scan
 
gscan :: [Point] -> [Point]
 
gscan ps = scan (sortByAngle ps)
 
where scan (x:y:z:xs) = if dir x y z == LeftTurn
 
then scan (x:z:xs)
 
else x: scan (y:z:xs)
 
scan [x,y] = [x,y] -- there's no shame in a pattern match
 
-- of this type!
 
scan _ = []
 
 
</haskell>
 
</haskell>
   
== Testing ==
 
The above implementation has an error, which can be found using
 
   
  +
== Test ==
 
<haskell>
 
<haskell>
 
import Test.QuickCheck
 
import Test.QuickCheck
Line 90: Line 56:
   
   
gscan' = map (\(Point x) -> x) . gscan . map Point
+
convex' = map (\(Pt x) -> x) . convex . map Pt
   
prop_gscan n = forAll (replicateM n arbitrary) (\xs -> sort (gscan' xs) == sort (gscan' (gscan' xs)))
+
prop_convex n = forAll (replicateM n arbitrary) (\xs -> sort (convex' xs) == sort (convex' (convex' xs)))
   
 
</haskell>
 
</haskell>
   
This leads to a counterexample for lists of 5 elements
 
  +
The results are:
   
 
<haskell>
 
<haskell>
> quickCheck (prop_gscan 3)
+
> quickCheck (prop_convex 0)
 
+++ OK, passed 100 tests.
 
+++ OK, passed 100 tests.
> quickCheck (prop_gscan 4)
+
> quickCheck (prop_convex 1)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 2)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 3)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 4)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 5)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 10)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 100)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 1000)
  +
+++ OK, passed 100 tests.
  +
> quickCheck (prop_convex 10000)
 
+++ OK, passed 100 tests.
 
+++ OK, passed 100 tests.
> quickCheck (prop_gscan 5)
 
*** Failed! Falsifiable (after 10 tests):
 
[(-2.4444173639942894,13.729627457461254),(65.72912810263666,5.955962930412828),(-34.288098030422404,80.6230134460068),(-6.446932942713564,-11.632835144720378),(2.861905401095031,1.1159493836896193)]
 
 
> let f = [(-2.4444173639942894,13.729627457461254),(65.72912810263666,5.955962930412828),(-34.288098030422404,80.6230134460068),(-6.446932942713564,-11.632835144720378),(2.861905401095031,1.1159493836896193)]
 
> gscan' f
 
[(-6.446932942713564,-11.632835144720378),(-34.288098030422404,80.6230134460068),(-2.4444173639942894,13.729627457461254),(65.72912810263666,5.955962930412828)]
 
> gscan' (gscan' f)
 
[(-6.446932942713564,-11.632835144720378),(-34.288098030422404,80.6230134460068),(65.72912810263666,5.955962930412828)]
 
> gscan' (gscan' (gscan' f))
 
[(-6.446932942713564,-11.632835144720378),(-34.288098030422404,80.6230134460068),(65.72912810263666,5.955962930412828)]
 
 
</haskell>
 
</haskell>
 
The versions of <tt>f</tt> with differing numbers of gscan applied are:
 
[[File:Graham_scan_implementation_failure.png]]
 

Latest revision as of 14:36, 5 January 2015

Descriptions of this problem can be found in Real World Haskell, Chapter 3

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

data Direction = LEFT | RIGHT | STRAIGHT
               deriving (Show, Eq)

data Pt = Pt (Double, Double)
        deriving (Show, Eq, Ord)

isTurned :: Pt -> Pt -> Pt -> Direction
isTurned (Pt (ax, ay)) (Pt (bx, by)) (Pt (cx, cy)) = case sign of
    EQ -> STRAIGHT
    LT -> RIGHT
    GT -> LEFT
    where sign = compare ((bx - ax) * (cy - ay)) ((cx - ax) * (by - ay))

gScan :: [Pt] -> [Pt]
gScan pts 
    | length pts >= 3 = scan [pt0] rests
    | otherwise       = pts
    where 
        -- Find the most bottom-left point pt0
        pt0 = foldr bottomLeft (Pt (1/0, 1/0)) pts where
            bottomLeft pa pb = case ord of
                               LT -> pa
                               GT -> pb
                               EQ -> pa
                       where ord = (compare `on` (\ (Pt (x, y)) -> (y, x))) pa pb

        -- Sort other points based on angle
        rests = tail (sortBy (compare `on` compkey pt0) pts) where
            compkey (Pt (x0, y0)) (Pt (x, y)) = (atan2 (y - y0) (x - x0),
                                       abs (x - x0))

        -- Scan the points to find out convex
        -- -- handle the case that all points are collinear
        scan [p0] (p1:ps)
            | isTurned pz p0 p1 == STRAIGHT = [pz, p0]
            where pz = last ps

        scan (x:xs) (y:z:rsts) = case isTurned x y z of
            RIGHT    -> scan xs (x:z:rsts)
            STRAIGHT -> scan (x:xs) (z:rsts) -- skip collinear points
            LEFT     -> scan (y:x:xs) (z:rsts)

        scan xs [z] = z : xs


Test

import Test.QuickCheck
import Control.Monad
import Data.List


convex' = map (\(Pt x) -> x) . convex . map Pt

prop_convex n = forAll (replicateM n arbitrary) (\xs -> sort (convex' xs) == sort (convex' (convex' xs)))

The results are:

> quickCheck (prop_convex 0)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 1)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 2)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 3)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 4)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 5)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 10)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 100)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 1000)
+++ OK, passed 100 tests.
> quickCheck (prop_convex 10000)
+++ OK, passed 100 tests.