Fwd: [Haskell-beginners] Optimization Problem

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"
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@haskell.org http://www.haskell.org/mailman/listinfo/beginners
_______________________________________________ Beginners mailing list Beginners@haskell.org http://www.haskell.org/mailman/listinfo/beginners

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@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"
To: beginners@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@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@haskell.org http://www.haskell.org/mailman/listinfo/beginners
_______________________________________________ Beginners mailing list Beginners@haskell.org http://www.haskell.org/mailman/listinfo/beginners _______________________________________________ Beginners mailing list Beginners@haskell.org http://www.haskell.org/mailman/listinfo/beginners

On Tuesday 10 August 2010 03:06:46, MAN wrote:
I believe the improvement was in avoiding the call to 'sqrt', rather than substituting (^2) by a mult.
That's the big one. sqrt is pretty expensive compared to a multiplication.
Although that last substitution is always welcomed in C, GHC could easily make it for you.
It could, but it doesn't (at least not in core). Most of the time it doesn't make a significant difference whether you use (x ^ 2) or (x * x), but sometimes it can make a big difference (we had an example earlier this year, but I can't seem to find it now).

On Tuesday 10 August 2010 04:02:32, Daniel Fischer wrote:
It could, but it doesn't (at least not in core). Most of the time it doesn't make a significant difference whether you use (x ^ 2) or (x * x), but sometimes it can make a big difference (we had an example earlier this year, but I can't seem to find it now).
Finally found it: http://www.haskell.org/pipermail/haskell-cafe/2010-May/077862.html shows an example where the use of (x ^ 2) instead of (x * x) prevented sharing of a common subexpression, thus creating an O(2^n) algorithm from what rightfully was an O(n) algorithm.
participants (3)
-
200901104@daiict.ac.in
-
Daniel Fischer
-
MAN