dph-examples: Simplify QuickHull eventlog tracing.
[packages/dph.git] / dph-examples / examples / spectral / QuickHull / handvec / Handvec.hs
1 {-# LANGUAGE MagicHash, BangPatterns, ParallelArrays, TypeOperators #-}
2
3 module Handvec ( hsplit_v, Point, Line )
4 where
5
6 import Points2D.Types
7
8 import Data.Array.Parallel as PA
9 import Data.Array.Parallel.Base as B
10 import Data.Array.Parallel.PArray ()
11 import Data.Array.Parallel.PArray.PData
12 import Data.Array.Parallel.PArray.PDataInstances
13 import Data.Array.Parallel.PArray.ScalarInstances
14 import Data.Array.Parallel.Lifted.Closure
15 import Data.Array.Parallel.Prelude as P'
16 import Data.Array.Parallel.Prelude.Double as D
17 import Data.Array.Parallel.Unlifted as U
18 import GHC.Exts
19
20 import Debug.Trace as Debug
21
22 import Prelude as P
23
24
25 hsplit_v :: PArray Point :-> Line :-> PArray Point
26 hsplit_v = closure2 hsplit_s hsplit_l
27 {-# NOINLINE hsplit_v #-}
28
29 -- | Scalar version of hsplit (generalised to lifted hsplit)
30 -- 1. construct a nested array out of 'points' array (single segment)
31 -- 2. construct a singleton array containing 'line'
32 --
33 -- UNUSED.
34 hsplit_s :: PArray Point -> Line -> PArray Point
35 hsplit_s ps@(PArray n# pts) line
36 = let segd = U.singletonSegd (I# n#)
37 pts_s = PArray 1# (PNested segd pts)
38 lines = PArray 1# (toPData line)
39 hull@(PArray _ (PNested segd' hull_pts))
40 = hsplit_l pts_s lines
41 !(I# n'#) = U.lengthSegd segd'
42 in PArray n'# hull_pts
43 where toPData ((x1,y1), (x2,y2)) = P_2 (P_2 (r x1) (r y1))
44 (P_2 (r x2) (r y2))
45 r x = PDouble (U.replicate 1 x)
46 {-# NOINLINE hsplit_s #-}
47
48
49 -- | The heart of QuickHull.
50 --
51 -- Wrapper for lifted hsplit with the type expected by the vectoriser
52 -- (which we still go through to avoid hand vectorising 'HandvecWrp.quickHull'.
53 hsplit_l :: PArray (PArray Point) -> PArray Line -> PArray (PArray Point)
54 hsplit_l (PArray _ points_s) (PArray nlines# lines)
55 = let result@(PNested segd _) = hsplit_l' nlines# points_s lines
56 !(I# nsegs#) = U.lengthSegd segd -- should be 2
57 in PArray nsegs# result
58 {-# INLINE hsplit_l #-}
59
60
61 -- | Actually the heart of QuickHull.
62 hsplit_l' :: Int# -> PData (PArray Point) -> PData Line -> PData (PArray Point)
63 hsplit_l' nlines# (PNested segd _) _
64 | Debug.traceEvent ("GANG[1/1] Issuing hsplit_l' " ++
65 "with " ++ show (I# nlines#) ++ " lines " ++
66 "and " ++ show (U.elementsSegd segd) ++ " points, " ++
67 "distributed as: " ++ show (U.toList $ U.lengthsSegd segd)) False = undefined
68 hsplit_l' 0# _ _ = PNested U.emptySegd (P_2 (PDouble U.empty) (PDouble U.empty))
69 hsplit_l' nlines# points_s -- E.g.: nlines#: 6, npts >= 13
70 lines@(P_2 (P_2 (PDouble line_x1s) (PDouble line_y1s))
71 (P_2 (PDouble line_x2s) (PDouble line_y2s)))
72 = let -- Find distance-like measures between each point and its respective line.
73 cross_s = calc_cross_s points_s lines
74
75 -- Throw out all the points which are below the line (i.e. distance <= 0).
76 above_s@(PNested aboveSegd (P_2 (PDouble xos) (PDouble yos)))
77 = calc_above_s points_s cross_s
78
79 -- For some of the lines there are no more points left (segment length 0).
80 -- The recursion for these lines has finished and we'll treat them separately.
81 (aboveSegdLens,_,_) = unpackSegd aboveSegd -- Segd ( [4, 5, 0, 2, 2, 0]
82 -- , [0, 4, 9, 9, 11, 13]
83 -- , 13 )
84
85 -- Empty segments selector (empty -> 1; otherwise -> 0)
86 selE = U.tagsToSel2 -- Sel2 ( [0, 0, 1, 0, 0, 1]
87 $ U.map B.fromBool -- , [0, 1, 0, 2, 3, 1]
88 $ U.zipWith (P.==) aboveSegdLens -- , 4, 2 )
89 (U.replicate (I# nlines#) 0)
90 tagsE = U.tagsSel2 selE -- [0, 0, 1, 0, 0, 1]
91
92 --- *** First take care of the lines for which the recursion hasn't finished *** ---
93 --- *** *** ---
94
95 -- Create segd for the points above those lines, which still have points above them.
96 -- This is equivalent to removing all empty segments from 'aboveSegd'.
97 segdNE = U.lengthsToSegd -- Segd ( [4, 5, 2, 2]
98 $ U.packByTag aboveSegdLens tagsE 0 -- , [0, 4, 9, 11]
99 -- , 13 )
100 (lensNE, idsNE, _) = unpackSegd segdNE
101
102 tagsE_r = U.replicate_s aboveSegd tagsE -- [0, 0, 0, 0, -- 4
103 -- 0, 0, 0, 0, 0, -- 5
104 -- -- 0
105 -- 0, 0, -- 2
106 -- 0, 0] -- 2
107 -- -- 0
108
109 -- WATCH OUT: Redundant packByTag. Basicaly remove all empty segments from 'aboveSegd',
110 -- without touching the data arrays
111 aboveNE@(PNested _ (P_2 (PDouble xos') (PDouble yos')))
112 = PNested segdNE (P_2 (PDouble $ U.packByTag xos tagsE_r 0)
113 (PDouble $ U.packByTag yos tagsE_r 0))
114
115 --- *** Prepare segd to call hsplit anew *** ---
116 --- *** *** ---
117
118 -- Number of non-empty segments
119 !nNE@(I# nNE#) = U.elementsSel2_0 selE -- 4
120
121 -- New segd with one element per non-empty segment of 'aboveSegd' (farthest point)
122 segdNEx1 = U.mkSegd (U.replicate nNE 1) -- Segd ( [1, 1, 1, 1]
123 (U.enumFromStepLen 0 1 nNE) -- , [0, 1, 3, 4]
124 nNE -- , 4 )
125 -- New segd with two elements per non-empty segment of 'aboveSegd'.
126 -- It's used indirectly in multiple places to denote the fact that we get two lines
127 -- for each line for which the recursion hasn't finished. They share the point farthest
128 -- from the old line.
129 segdNEx2 = segdNEx1 `U.plusSegd` segdNEx1 -- Segd ( [2, 2, 2, 2]
130 -- , [0, 2, 4, 8]
131 (_,_,eltsNEx2) = unpackSegd segdNEx2 -- , 8 )
132
133
134 -- Replicate each non-empty segment twice
135 lens_dbl = U.replicate_s segdNEx2 lensNE -- Segd ( [4, 4, 5, 5, 2, 2, 2, 2]
136 segd_dbl = U.lengthsToSegd lens_dbl -- , [0, 4, 8,13,18,20,22,24]
137 (_,_,elts_dbl)= unpackSegd segd_dbl -- , 26 )
138
139 -- Find indices to take from the points array to make it adhere to the doubled segd_dbl
140 -- (vsegs to the rescue?)
141 ids_to_take = U.enumFromStepLenEach elts_dbl -- 26 (length hint)
142 (U.replicate_s segdNEx2 idsNE) -- [0, 0, 4, 4, 9, 9,11,11]
143 (U.replicate eltsNEx2 1) -- [1, 1, 1, 1, 1, 1, 1, 1]
144 lens_dbl -- [4, 4, 5, 5, 2, 2, 2, 2]
145 -- ids_to_take: [0,1,2,3, 0,1,2,3, 4,5,6,7,8, 4,5,6,7,8, 9,10, 9,10, 11,12, 11,12]
146
147 -- Actually perform doubling of all segments.
148 -- Each line now has two identical segments for the points above it.
149 above_dbl = PNested segd_dbl (P_2 (PDouble $ U.bpermute xos' ids_to_take)
150 (PDouble $ U.bpermute yos' ids_to_take))
151
152 -- Take a moment to find points farthest from each line
153 fars@(P_2 (PDouble far_xs) (PDouble far_ys))
154 = calc_farthest_s points_s cross_s tagsE
155
156 -- Find lines to and from the farthest points (i.e. [:(p1, pm), (pm, p2):] in the original code).
157 -- Remember we are only dealing with lines for which there are still points above them.
158 -- The use of segments here is a convenient way to interleave the old line starts, with the
159 -- farthest points (also used as line starts in the next iteration).
160 -- segdNEx2 is not used directly and is merely a hint to segmented append.
161 line_1s' = P_2 (PDouble $ U.append_s segdNEx2
162 segdNEx1 (U.packByTag line_x1s tagsE 0)
163 segdNEx1 far_xs)
164
165 (PDouble $ U.append_s segdNEx2
166 segdNEx1 (U.packByTag line_y1s tagsE 0)
167 segdNEx1 far_ys)
168
169 -- Line ends are interleaved in a similar manner.
170 line_2s' = P_2 (PDouble $ U.append_s segdNEx2
171 segdNEx1 far_xs
172 segdNEx1 (U.packByTag line_x2s tagsE 0))
173
174 (PDouble $ U.append_s segdNEx2
175 segdNEx1 far_ys
176 segdNEx1 (U.packByTag line_y2s tagsE 0))
177 lines' = P_2 line_1s' line_2s' -- 8 lines in total
178
179 -- Finally make recursive call to compute convex hull from (when the recursion hasn't finished)
180 !(I# eltsNEx2#) = eltsNEx2
181 hullNE@(PNested hullNE_segd (P_2 (PDouble hullNE_xs) (PDouble hullNE_ys)))
182 = hsplit_l' eltsNEx2# above_dbl lines'
183
184 -- hsplit returns *only* line starts. In the example above we had 4 lines with points above them
185 (hullNE_lens,_,_) = unpackSegd hullNE_segd -- Segd: [ 2, 1, 2, 2, 1, 1, 2, 1 ]
186 -- ^ ^
187 -- Lines with finished recursion each contribute
188 -- one pt to hull. They go where ^ points.
189
190 --- *** Now deal with the empty segnments, i.e. lines which have no points above them *** ---
191 --- *** and will thus form part of the convex hull *** ---
192 !nE@(I# nE#)= U.elementsSel2_1 selE -- 2 lines have no more points above them
193
194 -- Prepare the final segd to hold both the lines that have no above points as well as the remainder
195 -- of the convex hull returned by the recursive call.
196 result_segd = U.lengthsToSegd -- Segd: [3, 4, 1, 2, 3, 1]
197 $ U.combine2 (U.tagsSel2 selE) -- [0, 0, 1, 0, 0, 1] (actually tagsE)
198 (U.repSel2 selE) -- Distribution across PEs
199 (U.sum_s segdNEx2 hullNE_lens) -- [3, 4, 2, 3]
200 (U.replicate nE 1) -- [1, 1]
201
202 result_sel = U.tagsToSel2 -- [0,0,0, 0,0,0,0, 1, 0,0, 0,0,0, 1]
203 $ (U.replicate_s result_segd tagsE)
204 result_tags = U.tagsSel2 result_sel -- As above
205 result_repsel = U.repSel2 result_sel -- Distribution across PEs
206
207 -- Combine data returned from 'hsplit' and the starts of the lines we left alone.
208 result_xs = U.combine2 result_tags
209 result_repsel
210 hullNE_xs
211 (U.packByTag line_x1s tagsE 1)
212 result_ys = U.combine2 result_tags
213 result_repsel
214 hullNE_ys
215 (U.packByTag line_y1s tagsE 1)
216
217 in (PNested result_segd (P_2 (PDouble result_xs) (PDouble result_ys)))
218 {-# INLINE hsplit_l' #-}
219
220 unpackSegd :: Segd -> (Array Int, Array Int, Int)
221 unpackSegd segd = ( U.lengthsSegd segd
222 , U.indicesSegd segd
223 , U.elementsSegd segd )
224 {-# INLINE unpackSegd #-}
225
226
227 -- | Find (relative) distances of points from line.
228 --
229 -- Each point can be above (positive distance) or below (negative distance)
230 -- a line as looking from the center of the convex hull.
231 --
232 -- Corresponds to 'cross' in the original program
233 -- > cross = [: distance p line | p <- points :]
234 calc_cross_s :: PData (PArray Point) -> PData Line -> PData (PArray Double)
235 calc_cross_s (PNested segd (P_2 (PDouble xos) (PDouble yos)))
236 (P_2 (P_2 (PDouble x1s) (PDouble y1s))
237 (P_2 (PDouble x2s) (PDouble y2s)))
238 = let crosses = U.zipWith6 distance xos
239 yos
240 (U.replicate_s segd x1s)
241 (U.replicate_s segd y1s)
242 (U.replicate_s segd x2s)
243 (U.replicate_s segd y2s)
244 in PNested segd (PDouble crosses)
245 {-# INLINE calc_cross_s #-}
246
247
248 -- | Calculate cross product between vectors formed between a point and
249 -- each of the two line ends.
250 --
251 -- I.e. |x1-xo| |x2-xo|
252 -- | | x | | = (x1-xo)(y2-yo)-(y1-yo)(x2-xo)
253 -- |y1-yo| |y2-yo|
254 --
255 -- Not changed from the original source version thanks to vectavoid
256 -- (except tuples are dissolved).
257 distance :: Double -> Double -- Point
258 -> Double -> Double -- Line start
259 -> Double -> Double -- Line end
260 -> Double -- Distance
261 distance xo yo x1 y1 x2 y2
262 = (x1 P.- xo) P.* (y2 P.- yo) P.- (y1 P.- yo) P.* (x2 P.- xo)
263 {-# INLINE distance #-}
264
265
266 -- | Find points above the lines given distance-like measures
267 --
268 -- Corresponds to 'packed' in the original program
269 -- > packed = [: p | (p,c) <- zipP points cross, c D.> 0.0 :]
270 calc_above_s :: PData (PArray Point) -> PData (PArray Double) -> PData (PArray Point)
271 calc_above_s (PNested segd (P_2 (PDouble xos) (PDouble yos))) -- points_s
272 (PNested _ (PDouble distances)) -- cross_s
273 = let length = U.elementsSegd segd
274
275 -- Compute selector for positive elements ((>0) -> 1; otherwise -> 0)
276 selAbove = U.tagsToSel2
277 $ U.map B.fromBool
278 $ U.zipWith (P.>) distances (U.replicate length 0.0)
279
280 -- Compute segd for just the positive elements
281 segdAbove = U.lengthsToSegd (U.count_s segd (tagsSel2 selAbove) 1)
282
283 -- Get the actual points corresponding to positive elements
284 tagsAbove = U.tagsSel2 selAbove
285 xosAbove = U.packByTag xos tagsAbove 1
286 yosAbove = U.packByTag yos tagsAbove 1
287
288 in (PNested segdAbove (P_2 (PDouble xosAbove) (PDouble yosAbove)))
289 {-# INLINE calc_above_s #-}
290
291
292 -- | For each line find a point fartherst from that line.
293 --
294 -- Each segment is a collection of points above a certain line.
295 -- The array of Doubles gives (relative) distances of points from the line.
296 --
297 -- Corresponds to 'pm' in the original program:
298 -- > pm = points !: maxIndexP cross
299 calc_farthest_s :: PData (PArray Point) -> PData (PArray Double) -> Array Tag -> PData Point
300 calc_farthest_s (PNested ptsSegd (P_2 (PDouble x1s) (PDouble y1s))) -- points_s
301 (PNested distSegd (PDouble distances)) -- cross_s
302 tags
303 = let -- Index the distances array, and find the indices corresponding to the
304 -- largest distance in each segment
305 indexed = U.zip (U.indices_s distSegd) distances
306 max_ids = U.fsts (U.fold1_s maxSnd distSegd indexed)
307
308 -- Find indices to take from the points array by offsetting from segment starts
309 ids = U.zipWith (P.+) (U.indicesSegd ptsSegd)
310 max_ids
311 max_xs = U.bpermute x1s ids
312 max_ys = U.bpermute y1s ids
313 -- however we are only interested in the ones which are from the non-empty segments
314 -- (TODO revise comment)
315 in P_2 (PDouble $ U.packByTag max_xs tags 0)
316 (PDouble $ U.packByTag max_ys tags 0)
317 {-# INLINE calc_farthest_s #-}
318
319
320 -- | Find pair with the bigest second element.
321 -- from dph-lifted-copy:D.A.P.Prelude.Int
322 maxSnd :: P.Ord b => (a, b) -> (a, b) -> (a, b)
323 maxSnd (i,x) (j,y) | x P.>= y = (i,x)
324 | P.otherwise = (j,y)
325 {-# INLINE maxSnd #-}