开发者

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 You

import 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.

  1. 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.

  2. 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.

  3. If you prefer to print out intermediate values as a means of debugging, then use trace and/or traceShow 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 
0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜