[Haskell-beginners] Optimization Problem

[Brent Yorgey]

I don't know all that much about profiling/optimization.  But one
simple thing occurs to me -- does it help to change inDisc to

  inDisc (Disc r p1 _) p2 = (distsqr p1 p2) <= r*r

  distsqr :: Point -> Point -> Double
  distsqr (x1, y1) (x2, y2) = (x1 - x2)^2 + (y1 - y2)^2

Then you avoid using sqrt.  Also there's no need to use ** since you
are taking a positive integral power.  (^2) will optimize to just a
single multiplication, whereas I have no idea what (**2) does, it's
probably much slower.

You could even store the square of the radius in Disc, then you also
avoid having to compute r*r every time.

-Brent

On Mon, Aug 09, 2010 at 01:16:54PM +0530, 200901104 at daiict.ac.in wrote:
> Optimization Problem 
> 
> I am trying to make an algorithm which solves Smallest circle problem( http ://en. wikipedia .org/ wiki /Smallest_circle_problem). 
> My code takes ~3.2 sec to solve ~500 data sets of 2000 points apiece. I need it to be done in less than 1 sec. 
> After profiling it seems that ~1.76 (55%) sec goes to the inDisc function and ~1 (32%) sec goes into gcnts . So I am trying 
> to optimize these parts of the code. 
> 
> -- is the point p2 is inside disc? 
> inDisc :: Disc -> Point -> Bool 
> inDisc (Disc r p1 _) p2 = (distance p1 p2) <= r 
> 
> -- gives the double values in one line 
> gcnts :: IO [Double] 
> gcnts = do 
> line <- getLine 
> return (map read (words line)) 
> 
> And data and point are 
> 
> -- (x, y) 
> type Point = (Double, Double) 
> 
> - radius, center, know points in the disc 
> data Disc = Disc Double Point [Point] 
> deriving (Show) 
> 
> The relevant profiler part: 
> COST CENTRE MODULE %time % alloc 
> 
> inDisc Main 55.0 0.0 
> gcnts Main 32.5 58.1 
> 
> individual inherited 
> COST CENTRE MODULE no. entries %time % alloc %time % alloc 
> inDisc Main 247 1237577 35.6 0.0 35.6 0.0 
> inDisc Main 244 527216 11.9 0.0 11.9 0.0 
> inDisc Main 241 214076 7.5 0.0 7.5 0.0 
> gcnts Main 235 8043 31.3 58.1 31.3 58.1 
> gcnts Main 233 10 1.3 0.0 1.3 0.0 
> gcnts Main 231 1 0.0 0.0 0.0 0.0 
> 
> I have attached my code and profiler output. 
> 
> BSRK Aditya . 

> import Data.List
> 
> type Point = (Double, Double)
> data Disc = Disc Double Point [Point]
>             deriving (Show)
> 
> getRad :: Disc -> Double
> getRad (Disc r _ _) = r
> 
> getPoints :: Disc -> [Point]
> getPoints (Disc _ _ ps) = ps
> 
> getCen :: Disc -> Point
> getCen (Disc _ p _) = p
> 
> midpoint :: Point -> Point -> Point
> midpoint (x1, y1) (x2, y2) = ((x1+x2)/2, (y1+y2)/2)
> 
> distance :: Point -> Point -> Double
> distance (x1, y1) (x2, y2) = sqrt((x1 - x2)**2 + (y1 - y2)**2)
> 
> inDisc :: Disc -> Point -> Bool
> inDisc (Disc r p1 _) p2 = (distance p1 p2) <= r
> 
> appendPoint :: Disc -> Point -> Disc
> appendPoint (Disc r p ps) p2 = (Disc r p (p2:ps))
> 
> circumCenter :: Point -> Point -> Point -> Point
> circumCenter (ax, ay) (bx, by) (cx, cy)
>     =  (((ay**2+ax**2)*(by-cy)+(by**2+bx**2)*(cy-ay)+(cy**2+cx**2)*(ay-by))/d,
>        ((ay**2+ax**2)*(cx-bx)+(by**2+bx**2)*(ax-cx)+(cy**2+cx**2)*(bx-ax))/d)
>        where d = 2*(ax*(by-cy)+bx*(cy-ay)+cx*(ay-by))
> 
> mround x = (fromIntegral (round (100*x)))/100
> -- the real code
> minDisc :: [Point] -> Disc
> minDisc (p1:p2:[]) = Disc ((distance p1 p2)/2) (midpoint p1 p2) (p1:p2:[])
> minDisc (p1:p2:ps) = foldl' helper (minDisc (p2:p1:[])) ps
>    where helper d p
>            | (inDisc d p) = (appendPoint d p)
>            | otherwise = (minDiscwp (getPoints d) p)
> minDisc _ = Disc 0 (0, 0) []
> 
> minDiscwp :: [Point] -> Point -> Disc
> minDiscwp ps q = foldl' (helper q) (minDisc ((head ps):q:[])) (tail ps)
>     where helper q d p
>            | (inDisc d p) = (appendPoint d p)
>            | otherwise = (minDiscwp2 (init (getPoints d)) p q)
> 
> minDiscwp2 :: [Point] -> Point -> Point -> Disc
> minDiscwp2 ps q1 q2 = foldl' (helper q1 q2) (minDisc (q1:q2:[])) ps
>     where helper q1 q2 d p
>            | (inDisc d p) = (appendPoint d p)
>            | otherwise = Disc (distance cr p) cr (p:(getPoints d))
>                where cr = circumCenter q1 q2 p
> --IO
> 
> solver :: Double -> [(Double, Double, Double, Double)] -> Double
> solver t lst =  mround(getRad (minDisc (map (h1 t) lst)))
>     where h1 t (ix, iy, vx, vy) = (ix+vx*t, iy+vy*t)
> 
> gcnts :: IO [Double]
> gcnts = do
>           line <- getLine
>           return (map read (words line))
> 
> ecase :: Double -> IO ()
> ecase cnt
>   | cnt == 0 = do return ()
>   | otherwise = do (n:t:[]) <- gcnts
>                    pts <- gvecs n
>                    execute 1 (t + 1) pts
>                    ecase (cnt - 1)
> 
> gvecs :: Double -> IO [(Double, Double, Double, Double)]
> gvecs n
>   | n == 0 = do return ([])
>   | otherwise = do  (ix:iy:vx:vy:[]) <- gcnts
>                     nex <- gvecs (n-1)
>                     return ((ix, iy, vx, vy):nex)
> 
> execute :: Double -> Double -> [(Double, Double, Double, Double)] -> IO ()
> execute ct tt lsts
>    | ct == tt = do return ()
>    | otherwise = do putStrLn (show (solver ct lsts))
>                     execute (ct + 1) tt lsts
> 
> main :: IO ()
> main = do
>           (t:[]) <- gcnts
>           ecase t


> _______________________________________________
> Beginners mailing list
> Beginners at haskell.org
> http://www.haskell.org/mailman/listinfo/beginners