Rotating caliper in Haskell
I am trying to implement rotating calipers in Haskell from Wikipedia . The only difference with Wikipedia is , i am calculating square of maximum width of convex polygon rather than minimum width to test the implementation of rotating calipers. It seems that this implementation is not correct because i got 97 for last test case of TFOSS rather than 98. Could some one please tell me what is wrong with this implementation. In case of indentation problem , i have posted the code on ideone.
Thank Youimport Data.List
import Data.Array
import Data.Maybe
data Point a = P a a deriving ( Show , Ord , Eq )
data Vector a = V a a deriving ( Show , Ord , Eq )
data Turn = S | L | R deriving ( Show , Eq , Ord , Enum )
--start of convex hull
compPoint :: ( Num a , Ord a ) => Point a -> Point a -> Ordering
compPoint ( P x1 y1 ) ( P x2 y2 )
| compare x1 x2 == EQ = compare y1 y2
| otherwise = compare x1 x2
sortPoint :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortPoint xs = sortBy ( \ x y -> compPoint x y ) xs
findTurn :: ( Num a , Ord a , Eq a ) => Point a -> Point a -> Point a -> Turn
findTurn ( P x0 y0 ) ( P x1 y1 ) ( P x2 y2 )
| ( y1 - y0 ) * ( x2- x0 ) < ( y2 - y0 ) * ( x1 - x0 ) = L
| ( y1 - y0 ) * ( x2- x0 ) == ( y2 - y0 ) * ( x1 - x0 ) = S
| otherwise = R
hullComputation :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ]
hullComputation [x] ( z:ys ) = hullComputation [z,x] ys
hullComputation xs [] = xs
hullComputation ( y : x : xs ) ( z : ys )
| findTurn x y z == R = hullComputation ( x:xs ) ( z : ys )
| findTurn x y z == S = hullComputation ( x:xs ) ( z : ys )
| otherwise = hullComputation ( z : y : x : xs ) ys
convexHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
convexHull [] = []
convexHull [ p ] = [ p ]
convexHull [ p1 , p2 ] = [ p1 , p2 ]
convexHull xs = final where
txs = sortPoint xs
( x : y : ys ) = txs
lhull = hullComputation [y,x] ys
( x': y' : xs' ) = reverse txs
uhull = hullComputation [ y' , x' ] xs'
final = ( init lhull ) ++ ( init uhull )
--end of convex hull
--dot product for getting angle
angVectors :: ( Num a , Ord a , Floating a ) => Vector a -> Vector a -> a
angVectors ( V ax ay ) ( V bx by ) = theta where
dot = ax * bx + ay * by
a = sqrt $ ax ^ 2 + ay ^ 2
b = sqrt $ bx ^ 2 + by ^ 2
theta = acos $ dot / a / b
--start of rotating caliper part http://en.wikipedia.org/wiki/Rotating_calipers
--rotate the vector x y by angle t
rotVector :: ( Num a , Ord a , Floating a ) => Vector a -> a -> Vector a
rotVector ( V x y ) t = V ( x * cos t - y * sin t ) ( x * sin t + y * cos t )
--square of dist between two points
distPoints :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> a
distPoints ( P x1 y1 ) ( P x2 y2 ) = ( x1 - x2 ) ^ 2 + ( y1 - y2 ) ^ 2
--rotating caliipers
rotCal :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a -> Int -> Int -> Vector a -> Vector a -> a -> Int -> a
rotCal arr ang pa pb ca@( V ax ay ) cb@( V bx by ) dia n
| ang > pi = dia
| otherwise = rotCal arr ang' pa' pb' ca' cb' dia' n where
P x1 y1 = arr !! pa
P x2 y2 = arr !! ( mod ( pa + 1 ) n )
P x3 y3 = arr !! pb
P x4 y4 = arr !! ( mod ( pb + 1 ) n )
t1 = angVectors ca ( V ( x2 - x1 ) ( y2 - y1 ) )
t2 = angVectors cb ( V ( x4 - x3 ) ( y4 - y3 ) )
ca' = rotVector ca $ min t1 t2
cb' = rotVector cb $ min t1 t2
ang' = ang + min t1 t2
pa' = if t1 < t2 then mod ( pa + 1 ) n else pa
pb' = if t1 >= t2 then mod ( pb + 1 ) n else pb
dia' = max dia $ distPoints ( arr !! pa' ) ( arr !! pb' )
--dia' = max dia $ if t1 < t2 then distPoints ( arr !! pa' ) ( arr !! pb ) else distPoints ( arr !! pb' ) ( arr !! pa )
solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> String
solve [] = "0"
solve [ p ] = "0"
solve [ p1 , p2 ] = show $ distPoints p1 p2
solve [ p1 , p2 , p3 ] = show $ max ( distPoints p1 p2 ) $ max ( distPoints p2 p3 ) ( distPoints p3 p1 )
solve arr = show $ rotCal arr' 0 pa pb ( V 1 0 ) ( V (-1) 0 ) dia n where
arr' = convexHull arr
y1 = minimumBy ( \( P _ y1 ) ( P _ y2 ) -> compare y1 y2 ) arr'
y2 = maximumBy ( \( P _ y1 ) ( P _ y2 ) -> compare y1 y2 ) arr'
pa = fromJust . findIndex ( \ t -> t == y1 ) $ arr'
pb = fromJust . findIndex ( \ t -> t == y2 ) $ arr'
dia = distPoints ( arr' !! pa ) ( arr' !! pb )
n = length arr'
--end of rotating caliper
--spoj code for testing
final :: ( Num a , Ord a , Floating a ) => [ Point a ] -> String
final [] = "0"
final [ p ] = "0"
final [ 开发者_JAVA技巧p1 , p2 ] = show $ distPoints p1 p2
final [ p1 , p2 , p3 ] = show $ max ( distPoints p1 p2 ) $ max ( distPoints p2 p3 ) ( distPoints p3 p1 )
final arr = solve . convexHull $ arr
format :: ( Num a , Ord a , Floating a ) => [ Int ] -> [ [ Point a ]]
format [] = []
format (x:xs ) = t : format b where
( a , b ) = splitAt ( 2 * x ) xs
t = helpFormat a where
helpFormat [] = []
helpFormat ( x' : y' : xs' ) = ( P ( fromIntegral x' ) ( fromIntegral y' ) ) : helpFormat xs'
readD :: String -> Int
readD = read
main = interact $ unlines . map final . format . concat . ( map . map ) readD . map words . tail . lines
--end of spoj code
I am not going to figure out where the mistake is in your code.
I am going to tell you about some simple debugging techniques.
Load your code into ghci, run the code interactively, and check the results are as you expect.
$ ghci ghci> :load your-program.hs ghci> compPoint (P 0 0) (P 0 0) EQ ghci>
Try calling
compPoint
with different arguments until you are satisfied it is correct. Then move onto the next function.Use Test.QuickCheck.
This is essentially automating the previous method.
ghci> :load your-program.hs ghci> :m +Test.QuickCheck ghci Test.QuickCheck> let prop_equalPointsAreEqual x y = EQ == compPoint (P x y) (P x y) ghci Test.QuickCheck> quickCheck prop_equalPointsAreEqual
...and test more complicated properties until you are satisfied
compPoint
is correct. Then move onto the next function.Google for a QuickCheck tutorial.
If you prefer to print out intermediate values as a means of debugging, then use
trace
and/ortraceShow
from Debug.Trace.
N.B. My examples start by testing the lower level functions and working up, but you may prefer to start at the upper level and work down.
I don't know what's wrong with your code, but I made it a bit simpler.
import Data.List
import Data.Array
import Data.Maybe
import Data.Monoid
data Point a = P a a deriving (Show, Ord, Eq)
--data Vector a = V a a deriving (Show, Ord, Eq)
--data Turn = S | L | R deriving (Show, Eq, Ord, Enum)
-- L is LT, S is EQ, R is GT
-- The is really the same as just compare on Point
compPoint :: (Ord a) => Point a -> Point a -> Ordering
compPoint (P x1 y1) (P x2 y2) = compare x1 x2 `mappend` compare y1 y2
sortPoint :: (Ord a) => [Point a] -> [Point a]
sortPoint = sortBy compPoint
-- simpler sortPoint = sort
findTurn :: (Num a, Ord a) => Point a -> Point a -> Point a -> Ordering
findTurn (P x0 y0) (P x1 y1) (P x2 y2) =
compare ((y1 - y0) * (x2 - x0)) ((y2 - y0) * (x1 - x0))
hullComputation :: (Num a, Ord a) => [Point a] -> [Point a] -> [Point a]
hullComputation [x] (z:ys) = hullComputation [z,x] ys
hullComputation xs [] = xs
hullComputation (y : x : xs) (z : ys) =
case findTurn x y z of
GT -> hullComputation (x : xs) (z : ys)
EQ -> hullComputation (x : xs) (z : ys) -- same as above
LT -> hullComputation (z : y : x : xs) ys
精彩评论