Fwd: [Haskell-beginners] Optimization Problem

Mon Aug 9 21:06:46 EDT 2010

I believe the improvement was in avoiding the call to 'sqrt', rather
than substituting (^2) by a mult.
Although that last substitution is always welcomed in C, GHC could
easily make it for you.

El lun, 09-08-2010 a las 16:23 +0530, 200901104 at daiict.ac.in escribió:
> You were right. I have changed the code to 
> 
> distance' :: Point →  Point →  Double
> distance' (x1, y1) (x2, y2) = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2)
> 
> inDisc :: Disc →  Point →  Bool
> inDisc (Disc r p1 _) p2 = (distance' p1 p2) ≤ r*r
> 
> There was not any performance improvement by changing ** to ^.
> Why is * faster that **2 or ^2?
> 
> Note: The change reduced run time to 1.24 sec.
> 
> ----- Original Message -----
> From: "Brent Yorgey" <byorgey at seas.upenn.edu>
> To: beginners at haskell.org
> Sent: Monday, August 9, 2010 3:14:43 PM
> Subject: Re: [Haskell-beginners] Optimization Problem
> 
> 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
> 
> _______________________________________________
> Beginners mailing list
> Beginners at haskell.org
> http://www.haskell.org/mailman/listinfo/beginners
> _______________________________________________
> Beginners mailing list
> Beginners at haskell.org
> http://www.haskell.org/mailman/listinfo/beginners