More delaunay optimisations
[packages/dph.git] / examples / quickhull / QuickHull.hs
1 -- This program is heavily GC bound unless we use a very large heap (eg, -H500M)
2
3 import GHC.Conc (par, pseq)
4
5 import System.IO
6 import Control.Exception (evaluate)
7 import System.Environment
8 import System.CPUTime
9 import System.Time
10 import System.Random
11
12 import qualified Data.Array.Parallel.Unlifted as U
13
14
15 -- Time
16 --
17
18
19 -- Random points generation
20 --
21
22 -- IMPORTANT: We use the same seed with the same random generator in all
23 -- quickhull codes. The asymptotic work complexity of quickhull
24 -- is between O (N) and O (N^2) depending on the input.
25 -- To compare benchmark results, they always need to use the same
26 -- input.
27
28 generatePoints :: Int -> [Point]
29 generatePoints n
30 = let rg = mkStdGen 42742 -- always use the same seed
31 in toPoints (take (2*n) (randomRs (-100, 100) rg))
32 where
33 toPoints [] = []
34 toPoints (x:y:pts) = Point x y : toPoints pts
35
36 loadPoints :: String -> IO [Point]
37 loadPoints file
38 = do
39 h <- openBinaryFile file ReadMode
40 upts <- U.hGet h
41 hClose h
42 convert (U.fsts upts) (U.snds upts)
43
44 convert :: U.Array Double -> U.Array Double -> IO [Point]
45 convert xs ys
46 = do
47 let pts = zipWith Point (U.toList xs) (U.toList ys)
48 evaluate $ nf pts
49 return pts
50
51
52 -- Benchmark
53 --
54
55 data Point = Point !Double !Double
56 data Line = Line Point Point
57
58 instance Show Point where
59 show (Point x y) = show (x, y)
60
61 nf (Point x y:xs) = x `seq` y `seq` nf xs
62 nf [] = ()
63
64 upper :: (a -> a -> Bool) -> [(a, b)] -> b
65 upper above = snd . foldl1 pick
66 where
67 pick left@(kl, _) right@(kr, _) | kl `above` kr = left
68 | otherwise = right
69
70 distance :: Point -> Line -> Double
71 distance (Point xo yo) (Line (Point x1 y1) (Point x2 y2))
72 = (x1-xo) * (y2 - yo) - (y1 - yo) * (x2 - xo)
73
74 hsplit :: [Point] -> Line -> [Point]
75 hsplit points line@(Line p1 p2)
76 | length packed < 2 = p1:packed
77 | otherwise = hsplit packed (Line p1 pm) ++ hsplit packed (Line pm p2)
78 where
79 cross = [ (distance p line, p) | p <- points ]
80 packed = [ p | (p, (c, _)) <- zip points cross, c > 0.0 ]
81
82 pm = upper (>) cross
83
84 quickHull :: [Point] -> [Point]
85 quickHull [] = []
86 quickHull points
87 = hsplit points (Line minx maxx) ++ hsplit points (Line maxx minx)
88 where
89 xs = [ (x, p) | p@(Point x y) <- points ]
90 minx = upper (<) xs
91 maxx = upper (>) xs
92
93
94 -- Parallel version
95
96 hsplitPar :: [Point] -> Line -> [Point]
97 hsplitPar points line@(Line p1 p2)
98 | length packed < 2 = p1:packed
99 | otherwise = let left = hsplitPar packed (Line p1 pm)
100 right = hsplitPar packed (Line pm p2)
101 in
102 right `par`
103 (left ++ right)
104 where
105 cross = [ (distance p line, p) | p <- points ]
106 packed = [ p | (p, (c, _)) <- zip points cross, c > 0.0 ]
107
108 pm = upper (>) cross
109
110 quickHullPar :: [Point] -> [Point]
111 quickHullPar [] = []
112 quickHullPar points
113 = let left = hsplitPar points (Line minx maxx)
114 right = hsplitPar points (Line maxx minx)
115 in
116 right `par`
117 (left ++ right)
118 where
119 xs = [ (x, p) | p@(Point x y) <- points ]
120 minx = upper (<) xs
121 maxx = upper (>) xs
122
123 -- OBSERVATION: If we use nf on 'right' in 'quickHullPar' and 'hsplitPar' (and maybe even
124 -- 'nf right `par` nf left `pseq` ...') the parallel GC takes a big hit and makes everything much
125 -- slower. (Keep in mind that even in the good case, this program spends 2/3 of its running time
126 -- in the GC.)
127
128 -- main
129 --
130
131 main :: IO ()
132 main
133 = do
134 [mode, args1, args2] <- getArgs
135 let runs = read args1
136 -- n = read args2
137 --
138 -- let pts = generatePoints n
139 -- eval pts `seq` return ()
140 pts <- loadPoints args2
141 let {-# NOINLINE oneRun #-} -- important to execute multiple runs
142 oneRun pts = do
143 t1 <- getTime
144 let res = case mode of
145 "seq" -> quickHull pts
146 "par" -> quickHullPar pts
147 _ -> error "mode must be 'seq' or 'par'"
148 evaluate $ nf res
149 t2 <- getTime
150 return (length res, fromTime (t2 `minus` t1))
151 results <- sequence (replicate runs (oneRun pts))
152
153 let (lens, times) = unzip results
154 (walls, cpus) = unzip times
155 putStrLn $ "Result length = " ++ show (head lens) ++ ": " ++
156 showWallCPU (minimum walls) (minimum cpus) ++ " " ++
157 showWallCPU (sum walls `div` toInteger runs)
158 (sum cpus `div` toInteger runs) ++ " " ++
159 showWallCPU (maximum walls) (maximum cpus)
160 where
161 showWallCPU wall cpu = show wall ++"/" ++ show cpu