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