Brent Yorgey byorgey at seas.upenn.edu
Mon Aug 9 05:44:43 EDT 2010

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

> 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