Maximum Triangle Area Wrong Answer [ Haskell ]

mukesh_iiitm 2011-08-04 22:19:00 UTC #1

Hello all , I am trying to solve this problem [ spoj.pl/problems/MTRIAREA ] but i am getting wrong answer .I am using the rotating caliper technique explained at http://en.wikipedia.org/wiki/Rotating_calipers . Could some one please tell me what is wrong with this approach . It would be also great if some one accepted can post some test case for this problem .
Thank You
[bbone=haskell,263]
import Data.List
import Data.Array
import Data.Maybe
import Data.Function
import Text.Printf
import qualified Data.ByteString.Char8 as BS

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

findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
findMinx xs = sortBy ( \x y -> compPoint x y ) xs

compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a -> Ordering
compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( ( y1 - y0 ) * ( x2 - x0 ) ) ( ( y2 - y0) * ( x1 - x0 ) )

sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) 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

findHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ]
findHull [x] ( z : ys ) = findHull [ z , x ] ys --incase of second point on line from x to z
findHull xs [] = xs
findHull ( y : x : xs ) ( z : ys )
| findTurn x y z == R = findHull ( x : xs ) ( z:ys )
| findTurn x y z == S = findHull ( x : xs ) ( z:ys )
| otherwise = findHull ( z : y : x : xs ) ys

convexHull : Num a , Ord a ) => [ Point a ] -> [ Point a ]
convexHull xs = reverse . findHull [ y , x ] $ ys where
( x : y : ys ) = sortByangle . findMinx $ xs

--end of convex hull

--start of rotating caliper part http://en.wikipedia.org/wiki/Rotating_calipers
--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

--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 )

--area of triangle

computeArea :: ( Num a , Ord a , Floating a ) => Point a -> Point a -> Point a -> a
computeArea ( P x1 y1 ) ( P x2 y2 ) ( P x3 y3 ) = 0.5 * abs ( ( x1 * y2 + x2 * y3 + x3 * y1 ) - ( y1 * x2 + y2 * x3 + y3 * x1 ) )

--rotating caliipers

rotCal :: ( Num a , Ord a , Floating a ) => Array Int ( 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 ) area n
| ang > pi = area
| otherwise = rotCal arr ang' pa' pb' ca' cb' area' 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' , pb' , pre ) = if t1 < t2 then ( mod ( pa + 1 ) n , pb , pa ) else ( pa , mod ( pb + 1 ) n , pb )
area' = max area $ computeArea ( arr ! pre ) ( arr ! pa' ) ( arr ! pb' )

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
solve [] = 0
solve [ p ] = 0
solve [ p1 , p2 ] = 0
solve arr = rotCal arr' 0 pa pb ( V 1 0 ) ( V (-1) 0 ) area $ n where
y1 = minimumBy ( on compare fN ) arr
y2 = maximumBy ( on compare fN ) arr
pa = fromJust . findIndex ( == y1 ) $ arr
pb = fromJust . findIndex ( == y2 ) $ arr
n = length arr
arr' = listArray ( 0 , n ) arr
area = computeArea ( arr' ! pa ) ( arr' ! ( mod ( pa + 1 ) n ) ) ( arr' ! pb ) --compute first area
fN ( P x y ) = y

--end of rotating caliper

--spoj code for testing but time limit exceed

final :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
final [] = 0
final [ p ] = 0
final [ p1 , p2 ] = 0
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 :: BS.ByteString -> Int
readD = fst . fromJust . BS.readInt

main = BS.interact $ BS.unlines . map ( BS.pack . ( printf "%.2f" :: Double -> String ) . final ) . format . concat . map ( map readD . BS.words ) . init . BS.lines

--end of spoj code
[/bbone]

system 2011-08-06 17:17:35 UTC #2

Answer : 32214600.50

mukesh_iiitm 2011-08-06 20:56:26 UTC #3

Thanks vipul . Now i am suspecting that if this algorithm is able to solve this problem .Any one who is trying to solve this problem , rotating caliper algorithm assumes that one side of triangle will always overlap with convex polygon . May be ternary search could be helpful.

system 2011-08-06 21:22:32 UTC #4

Hint : Rotate calliper on vertices not on edges

mukesh_iiitm 2011-08-06 22:55:24 UTC #5

I did not get you. Initially i chose two caliper parallel to x axis , one passing from minimum y coordinate ( i + 0 * j ) [ say pa ] and second one passing through maximum y coordinate ( - i + 0 * j ) [ say pb ] . Now i computed the angle between vector ( pa , pa + 1 ) and ( i + 0 * j ) [ say patheta ] and angle between ( pb , pb + 1 ) with ( -i + 0 j ) [ say pbtheta ] . Now rotate the calipers by minimum of these two angles so technically , at every iterating one of caliper will overlap with side of polygon . Could you explain what do you mean by " Rotate calliper on vertices not on edges " .

system 2011-08-06 23:51:24 UTC #6

Required triangle's edges may not coincide with that of convex hull but the 3 points will coincide with its vertices.
1. Choose a,b,c as first three points of convex hull( initial area = area of triangle abc )
2. Keep selecting next point to c as new c till area increases.
3. Now, select next point to b as new b, if area increases go to step 2.
4. Repeat 2,3 for all a and keep a track of maximum area found so far.

mukesh_iiitm 2011-08-10 16:00:51 UTC #7

HI Vipul , i implemented the algorithm explained by you [ also same algorithm ideone.com/XoGhN ] .
Thank you
[bbone=haskell,267]
import Data.List
import Data.Array
import Data.Maybe
import qualified Data.ByteString.Char8 as BS
import Text.Printf

data Point a = P a a deriving ( Show , Ord , Eq )
data Turn = S | L | R deriving ( Show , Ord , Eq , Enum )

--start of graham scan

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

findMinx :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
findMinx xs = sortBy ( \x y -> compPoint x y ) xs

compAngle ::(Num a , Ord a ) => Point a -> Point a -> Point a -> Ordering
compAngle ( P x1 y1 ) ( P x2 y2 ) ( P x0 y0 ) = compare ( ( y1 - y0 ) * ( x2 - x0 ) ) ( ( y2 - y0) * ( x1 - x0 ) )

sortByangle :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ]
sortByangle (z:xs) = z : sortBy ( \x y -> compAngle x y z ) xs

findHull :: ( Num a , Ord a ) => [ Point a ] -> [ Point a ] -> [ Point a ]
findHull [x] ( z : ys ) = findHull [ z , x ] ys --incase of second point on line from x to z
findHull xs [] = xs
findHull ( y : x : xs ) ( z:ys )
| findTurn x y z == R = findHull ( x : xs ) ( z:ys )
| findTurn x y z == S = findHull ( x : xs ) ( z:ys )
| otherwise = findHull ( z : y : x : xs ) ys

convexHull : Num a , Ord a ) => [ Point a ] -> [ Point a ]
convexHull xs = reverse . findHull [y,x] $ ys where
(x:y:ys) = sortByangle.findMinx $ xs

--end of graham scan

caltmp :: ( Num a , Ord a , Floating a ) => Int -> Int -> Int -> Array Int ( Point a ) -> a
caltmp a b c arr = area where
P x1 y1 = arr ! a
P x2 y2 = arr ! b
P x3 y3 = arr ! c
area = 0.5 * ( abs $ ( x1 * y2 + x2 * y3 + x3 * y1 ) - ( y1 * x2 + y2 * x3 + y3 * x1 ) )

calArea :: ( Num a , Ord a , Floating a ) => Int -> Int -> Int -> Int -> Int -> a -> Array Int ( Point a ) -> a
calArea a b c k n area arr
| k >= n = area
| area' >= area = calArea a b ( mod ( c + 1 ) n ) k n area' arr --area a b ( c + 1 ) is greater than area a b c
| area'' >= area = calArea a ( mod ( b + 1 ) n ) c k n area'' arr --check if area a ( b + 1 ) c is greater area a b c
| otherwise = calArea a' b c ( k + 1 ) n area''' arr --repeat with next a
where
area' = caltmp a b ( mod ( c + 1 ) n ) arr
area'' = caltmp a ( mod ( b + 1 ) n ) c arr
a' = mod ( a + 1 ) n
area''' = max area $ caltmp a' b c arr

solve :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
solve [] = 0
solve [ p ] = 0
solve [ p1 , p2 ] = 0
solve arr =
calArea 0 1 2 0 n area arr' where
--maximum $ map ( \x -> calArea x ( mod ( x + 1 ) n ) ( mod ( x + 2 ) n ) 0 n ( caltmp x ( mod ( x + 1 ) n ) ( mod ( x + 2 ) n ) arr' ) arr' ) [0 .. pred n ] where
n = length arr
arr' = listArray ( 0 , n ) arr
area = caltmp 0 1 2 arr'

final :: ( Num a , Ord a , Floating a ) => [ Point a ] -> a
final [] = 0
final [ p ] = 0
final [ p1 , p2 ] = 0
final arr = solve . convexHull $ arr

readD :: BS.ByteString -> Int
readD = fst . fromJust . BS.readInt

main = BS.interact $ BS.unlines . map ( BS.pack . ( printf "%.2f" :: Double -> String ) . final ) . format . concat . map ( map readD . BS.words ) . init . BS.lines
[/bbone]