Difference between revisions of "Graham Scan Implementation"
From HaskellWiki
(add example where this code doesn't work) 
m (fixup Pt definition) 

(4 intermediate revisions by 2 users 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 bottomleft point pt0 

−  EQ > Straight 

+  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 

−  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 xaxis 

+  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:  
−  +  convex' = map (\(Pt x) > x) . convex . map Pt 

−  +  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 ( 
+  > quickCheck (prop_convex 0) 
+++ OK, passed 100 tests. 
+++ OK, passed 100 tests. 

−  > quickCheck ( 
+  > 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> 
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 bottomleft 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.