1 | //STARTHEADER |
---|
2 | // $Id: ClosestPair2D.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
3 | // |
---|
4 | // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez |
---|
5 | // |
---|
6 | //---------------------------------------------------------------------- |
---|
7 | // This file is part of FastJet. |
---|
8 | // |
---|
9 | // FastJet is free software; you can redistribute it and/or modify |
---|
10 | // it under the terms of the GNU General Public License as published by |
---|
11 | // the Free Software Foundation; either version 2 of the License, or |
---|
12 | // (at your option) any later version. |
---|
13 | // |
---|
14 | // The algorithms that underlie FastJet have required considerable |
---|
15 | // development and are described in hep-ph/0512210. If you use |
---|
16 | // FastJet as part of work towards a scientific publication, please |
---|
17 | // include a citation to the FastJet paper. |
---|
18 | // |
---|
19 | // FastJet is distributed in the hope that it will be useful, |
---|
20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
22 | // GNU General Public License for more details. |
---|
23 | // |
---|
24 | // You should have received a copy of the GNU General Public License |
---|
25 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>. |
---|
26 | //---------------------------------------------------------------------- |
---|
27 | //ENDHEADER |
---|
28 | |
---|
29 | #include "fastjet/internal/ClosestPair2D.hh" |
---|
30 | |
---|
31 | #include<limits> |
---|
32 | #include<iostream> |
---|
33 | #include<iomanip> |
---|
34 | #include<algorithm> |
---|
35 | |
---|
36 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
37 | |
---|
38 | const unsigned int huge_unsigned = 4294967295U; |
---|
39 | const unsigned int twopow31 = 2147483648U; |
---|
40 | |
---|
41 | using namespace std; |
---|
42 | |
---|
43 | //---------------------------------------------------------------------- |
---|
44 | /// takes a point and sets a shuffle with the given shift... |
---|
45 | void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle, |
---|
46 | unsigned int shift) { |
---|
47 | |
---|
48 | Coord2D renorm_point = (point.coord - _left_corner)/_range; |
---|
49 | // make sure the point is sensible |
---|
50 | //cerr << point.coord.x <<" "<<point.coord.y<<endl; |
---|
51 | assert(renorm_point.x >=0); |
---|
52 | assert(renorm_point.x <=1); |
---|
53 | assert(renorm_point.y >=0); |
---|
54 | assert(renorm_point.y <=1); |
---|
55 | |
---|
56 | shuffle.x = static_cast<unsigned int>(twopow31 * renorm_point.x) + shift; |
---|
57 | shuffle.y = static_cast<unsigned int>(twopow31 * renorm_point.y) + shift; |
---|
58 | shuffle.point = &point; |
---|
59 | } |
---|
60 | |
---|
61 | //---------------------------------------------------------------------- |
---|
62 | /// compares this shuffle with the other one |
---|
63 | bool ClosestPair2D::Shuffle::operator<(const Shuffle & q) const { |
---|
64 | |
---|
65 | if (floor_ln2_less(x ^ q.x, y ^ q.y)) { |
---|
66 | // i = 2 in Chan's algorithm |
---|
67 | return (y < q.y); |
---|
68 | } else { |
---|
69 | // i = 1 in Chan's algorithm |
---|
70 | return (x < q.x); |
---|
71 | } |
---|
72 | } |
---|
73 | |
---|
74 | |
---|
75 | |
---|
76 | //---------------------------------------------------------------------- |
---|
77 | void ClosestPair2D::_initialize(const std::vector<Coord2D> & positions, |
---|
78 | const Coord2D & left_corner, |
---|
79 | const Coord2D & right_corner, |
---|
80 | unsigned int max_size) { |
---|
81 | |
---|
82 | unsigned int n_positions = positions.size(); |
---|
83 | assert(max_size >= n_positions); |
---|
84 | |
---|
85 | //_points(positions.size()) |
---|
86 | |
---|
87 | // allow the points array to grow to the following size |
---|
88 | _points.resize(max_size); |
---|
89 | // currently unused points are immediately made available on the |
---|
90 | // stack |
---|
91 | for (unsigned int i = n_positions; i < max_size; i++) { |
---|
92 | _available_points.push(&(_points[i])); |
---|
93 | } |
---|
94 | |
---|
95 | _left_corner = left_corner; |
---|
96 | _range = max((right_corner.x - left_corner.x), |
---|
97 | (right_corner.y - left_corner.y)); |
---|
98 | |
---|
99 | // initialise the coordinates for the points and create the zero-shifted |
---|
100 | // shuffle array |
---|
101 | vector<Shuffle> shuffles(n_positions); |
---|
102 | for (unsigned int i = 0; i < n_positions; i++) { |
---|
103 | // set up the points |
---|
104 | _points[i].coord = positions[i]; |
---|
105 | _points[i].neighbour_dist2 = numeric_limits<double>::max(); |
---|
106 | _points[i].review_flag = 0; |
---|
107 | |
---|
108 | // create shuffle with 0 shift. |
---|
109 | _point2shuffle(_points[i], shuffles[i], 0); |
---|
110 | } |
---|
111 | |
---|
112 | // establish what our shifts will be |
---|
113 | for (unsigned ishift = 0; ishift < _nshift; ishift++) { |
---|
114 | // make sure we use double-precision for calculating the shifts |
---|
115 | // since otherwise we will hit integer overflow. |
---|
116 | _shifts[ishift] = static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift); |
---|
117 | if (ishift == 0) {_rel_shifts[ishift] = 0;} |
---|
118 | else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];} |
---|
119 | } |
---|
120 | //_shifts[0] = 0; |
---|
121 | //_shifts[1] = static_cast<unsigned int>((twopow31*1.0)/3.0); |
---|
122 | //_shifts[2] = static_cast<unsigned int>((twopow31*2.0)/3.0); |
---|
123 | //_rel_shifts[0] = 0; |
---|
124 | //_rel_shifts[1] = _shifts[1]; |
---|
125 | //_rel_shifts[2] = _shifts[2]-_shifts[1]; |
---|
126 | |
---|
127 | // and how we will search... |
---|
128 | //_cp_search_range = 49; |
---|
129 | _cp_search_range = 30; |
---|
130 | _points_under_review.reserve(_nshift * _cp_search_range); |
---|
131 | |
---|
132 | // now initialise the three trees |
---|
133 | for (unsigned int ishift = 0; ishift < _nshift; ishift++) { |
---|
134 | |
---|
135 | // shift the shuffles if need be. |
---|
136 | if (ishift > 0) { |
---|
137 | unsigned rel_shift = _rel_shifts[ishift]; |
---|
138 | for (unsigned int i = 0; i < shuffles.size(); i++) { |
---|
139 | shuffles[i] += rel_shift; } |
---|
140 | } |
---|
141 | |
---|
142 | // sort the shuffles |
---|
143 | sort(shuffles.begin(), shuffles.end()); |
---|
144 | |
---|
145 | // and create the search tree |
---|
146 | _trees[ishift] = auto_ptr<Tree>(new Tree(shuffles, max_size)); |
---|
147 | |
---|
148 | // now we look for the closest-pair candidates on this tree |
---|
149 | circulator circ = _trees[ishift]->somewhere(), start=circ; |
---|
150 | // the actual range in which we search |
---|
151 | unsigned int CP_range = min(_cp_search_range, n_positions-1); |
---|
152 | do { |
---|
153 | Point * this_point = circ->point; |
---|
154 | //cout << _ID(this_point) << " "; |
---|
155 | this_point->circ[ishift] = circ; |
---|
156 | // run over all points within _cp_search_range of this_point on tree |
---|
157 | circulator other = circ; |
---|
158 | for (unsigned i=0; i < CP_range; i++) { |
---|
159 | ++other; |
---|
160 | double dist2 = this_point->distance2(*other->point); |
---|
161 | if (dist2 < this_point->neighbour_dist2) { |
---|
162 | this_point->neighbour_dist2 = dist2; |
---|
163 | this_point->neighbour = other->point; |
---|
164 | } |
---|
165 | } |
---|
166 | } while (++circ != start); |
---|
167 | //cout << endl<<endl; |
---|
168 | } |
---|
169 | |
---|
170 | // now initialise the heap object... |
---|
171 | vector<double> mindists2(n_positions); |
---|
172 | for (unsigned int i = 0; i < n_positions; i++) { |
---|
173 | mindists2[i] = _points[i].neighbour_dist2;} |
---|
174 | |
---|
175 | _heap = auto_ptr<MinHeap>(new MinHeap(mindists2, max_size)); |
---|
176 | } |
---|
177 | |
---|
178 | |
---|
179 | //----------------------------------------------------------------------= |
---|
180 | void ClosestPair2D::closest_pair(unsigned int & ID1, unsigned int & ID2, |
---|
181 | double & distance2) const { |
---|
182 | ID1 = _heap->minloc(); |
---|
183 | ID2 = _ID(_points[ID1].neighbour); |
---|
184 | distance2 = _points[ID1].neighbour_dist2; |
---|
185 | if (ID1 > ID2) swap(ID1,ID2); |
---|
186 | } |
---|
187 | |
---|
188 | |
---|
189 | //---------------------------------------------------------------------- |
---|
190 | inline void ClosestPair2D::_add_label(Point * point, unsigned int review_flag) { |
---|
191 | |
---|
192 | // if it's not already under review, then put it on the list of |
---|
193 | // points needing review |
---|
194 | if (point->review_flag == 0) _points_under_review.push_back(point); |
---|
195 | |
---|
196 | // OR the point's current flag with the requested review flag |
---|
197 | point->review_flag |= review_flag; |
---|
198 | } |
---|
199 | |
---|
200 | //---------------------------------------------------------------------- |
---|
201 | inline void ClosestPair2D::_set_label(Point * point, unsigned int review_flag) { |
---|
202 | |
---|
203 | // if it's not already under review, then put it on the list of |
---|
204 | // points needing review |
---|
205 | if (point->review_flag == 0) _points_under_review.push_back(point); |
---|
206 | |
---|
207 | // SET the point's current flag to the requested review flag |
---|
208 | point->review_flag = review_flag; |
---|
209 | } |
---|
210 | |
---|
211 | //---------------------------------------------------------------------- |
---|
212 | void ClosestPair2D::remove(unsigned int ID) { |
---|
213 | |
---|
214 | //cout << "While removing " << ID <<endl; |
---|
215 | |
---|
216 | Point * point_to_remove = & (_points[ID]); |
---|
217 | |
---|
218 | // remove this point from the search tree |
---|
219 | _remove_from_search_tree(point_to_remove); |
---|
220 | |
---|
221 | // the above statement labels certain points as needing "review" -- |
---|
222 | // deal with them... |
---|
223 | _deal_with_points_to_review(); |
---|
224 | } |
---|
225 | |
---|
226 | |
---|
227 | //---------------------------------------------------------------------- |
---|
228 | void ClosestPair2D::_remove_from_search_tree(Point * point_to_remove) { |
---|
229 | |
---|
230 | // add this point to the list of "available" points (this is |
---|
231 | // relevant also from the point of view of determining the range |
---|
232 | // over which we circulate). |
---|
233 | _available_points.push(point_to_remove); |
---|
234 | |
---|
235 | // label it so that it goes from the heap |
---|
236 | _set_label(point_to_remove, _remove_heap_entry); |
---|
237 | |
---|
238 | // establish the range over which we search (a) for points that have |
---|
239 | // acquired a new neighbour and (b) for points which had ID as their |
---|
240 | // neighbour; |
---|
241 | |
---|
242 | unsigned int CP_range = min(_cp_search_range, size()-1); |
---|
243 | |
---|
244 | // then, for each shift |
---|
245 | for (unsigned int ishift = 0; ishift < _nshift; ishift++) { |
---|
246 | //cout << " ishift = " << ishift <<endl; |
---|
247 | // get the circulator for the point we'll remove and its successor |
---|
248 | circulator removed_circ = point_to_remove->circ[ishift]; |
---|
249 | circulator right_end = removed_circ.next(); |
---|
250 | // remove the point |
---|
251 | _trees[ishift]->remove(removed_circ); |
---|
252 | |
---|
253 | // next find the point CP_range points to the left |
---|
254 | circulator left_end = right_end, orig_right_end = right_end; |
---|
255 | for (unsigned int i = 0; i < CP_range; i++) {left_end--;} |
---|
256 | |
---|
257 | if (size()-1 < _cp_search_range) { |
---|
258 | // we have a smaller range now than before -- but when seeing who |
---|
259 | // could have had ID as a neighbour, we still need to use the old |
---|
260 | // range for seeing how far back we search (but new separation between |
---|
261 | // points). [cf CCN28-42] |
---|
262 | left_end--; right_end--; |
---|
263 | } |
---|
264 | |
---|
265 | // and then for each left-end point: establish if the removed |
---|
266 | // point was its neighbour [in which case a new neighbour must be |
---|
267 | // found], otherwise see if the right-end point is a closer neighbour |
---|
268 | do { |
---|
269 | Point * left_point = left_end->point; |
---|
270 | |
---|
271 | //cout << " visited " << setw(3)<<_ID(left_point)<<" (its neighbour was "<< setw(3)<< _ID(left_point->neighbour)<<")"<<endl; |
---|
272 | |
---|
273 | if (left_point->neighbour == point_to_remove) { |
---|
274 | // we'll deal with it later... |
---|
275 | _add_label(left_point, _review_neighbour); |
---|
276 | } else { |
---|
277 | // check to see if right point has become its closest neighbour |
---|
278 | double dist2 = left_point->distance2(*right_end->point); |
---|
279 | if (dist2 < left_point->neighbour_dist2) { |
---|
280 | left_point->neighbour = right_end->point; |
---|
281 | left_point->neighbour_dist2 = dist2; |
---|
282 | // NB: (LESSER) REVIEW NEEDED HERE TOO... |
---|
283 | _add_label(left_point, _review_heap_entry); |
---|
284 | } |
---|
285 | } |
---|
286 | ++right_end; |
---|
287 | } while (++left_end != orig_right_end); |
---|
288 | } // ishift... |
---|
289 | |
---|
290 | } |
---|
291 | |
---|
292 | |
---|
293 | //---------------------------------------------------------------------- |
---|
294 | void ClosestPair2D::_deal_with_points_to_review() { |
---|
295 | |
---|
296 | // the range in which we carry out searches for new neighbours on |
---|
297 | // the search tree |
---|
298 | unsigned int CP_range = min(_cp_search_range, size()-1); |
---|
299 | |
---|
300 | // now deal with the points that are "under review" in some way |
---|
301 | // (have lost their neighbour, or need their heap entry updating) |
---|
302 | while(_points_under_review.size() > 0) { |
---|
303 | // get the point to be considered |
---|
304 | Point * this_point = _points_under_review.back(); |
---|
305 | // remove it from the list |
---|
306 | _points_under_review.pop_back(); |
---|
307 | |
---|
308 | if (this_point->review_flag & _remove_heap_entry) { |
---|
309 | // make sure no other flags are on (it wouldn't be consistent?) |
---|
310 | assert(!(this_point->review_flag ^ _remove_heap_entry)); |
---|
311 | _heap->remove(_ID(this_point)); |
---|
312 | } |
---|
313 | // check to see if the _review_neighbour flag is on |
---|
314 | else { |
---|
315 | if (this_point->review_flag & _review_neighbour) { |
---|
316 | this_point->neighbour_dist2 = numeric_limits<double>::max(); |
---|
317 | // among all three shifts |
---|
318 | for (unsigned int ishift = 0; ishift < _nshift; ishift++) { |
---|
319 | circulator other = this_point->circ[ishift]; |
---|
320 | // among points within CP_range |
---|
321 | for (unsigned i=0; i < CP_range; i++) { |
---|
322 | ++other; |
---|
323 | double dist2 = this_point->distance2(*other->point); |
---|
324 | if (dist2 < this_point->neighbour_dist2) { |
---|
325 | this_point->neighbour_dist2 = dist2; |
---|
326 | this_point->neighbour = other->point; |
---|
327 | } |
---|
328 | } |
---|
329 | } |
---|
330 | } |
---|
331 | |
---|
332 | // for any non-zero review flag we'll have to update the heap |
---|
333 | _heap->update(_ID(this_point), this_point->neighbour_dist2); |
---|
334 | } |
---|
335 | |
---|
336 | // "delabel" the point |
---|
337 | this_point->review_flag = 0; |
---|
338 | |
---|
339 | } |
---|
340 | |
---|
341 | } |
---|
342 | |
---|
343 | //---------------------------------------------------------------------- |
---|
344 | unsigned int ClosestPair2D::insert(const Coord2D & new_coord) { |
---|
345 | |
---|
346 | // get hold of a point |
---|
347 | assert(_available_points.size() > 0); |
---|
348 | Point * new_point = _available_points.top(); |
---|
349 | _available_points.pop(); |
---|
350 | |
---|
351 | // set the point's coordinate |
---|
352 | new_point->coord = new_coord; |
---|
353 | |
---|
354 | // now find it's neighbour in the search tree |
---|
355 | _insert_into_search_tree(new_point); |
---|
356 | |
---|
357 | // sort out other points that may have been affected by this, |
---|
358 | // and/or for which the heap needs to be updated |
---|
359 | _deal_with_points_to_review(); |
---|
360 | |
---|
361 | // |
---|
362 | return _ID(new_point); |
---|
363 | } |
---|
364 | |
---|
365 | //---------------------------------------------------------------------- |
---|
366 | unsigned int ClosestPair2D::replace(unsigned int ID1, unsigned int ID2, |
---|
367 | const Coord2D & position) { |
---|
368 | |
---|
369 | // deletion from tree... |
---|
370 | Point * point_to_remove = & (_points[ID1]); |
---|
371 | _remove_from_search_tree(point_to_remove); |
---|
372 | |
---|
373 | point_to_remove = & (_points[ID2]); |
---|
374 | _remove_from_search_tree(point_to_remove); |
---|
375 | |
---|
376 | // insertion into tree |
---|
377 | // get hold of a point |
---|
378 | Point * new_point = _available_points.top(); |
---|
379 | _available_points.pop(); |
---|
380 | |
---|
381 | // set the point's coordinate |
---|
382 | new_point->coord = position; |
---|
383 | |
---|
384 | // now find it's neighbour in the search tree |
---|
385 | _insert_into_search_tree(new_point); |
---|
386 | |
---|
387 | // the above statement labels certain points as needing "review" -- |
---|
388 | // deal with them... |
---|
389 | _deal_with_points_to_review(); |
---|
390 | |
---|
391 | // |
---|
392 | return _ID(new_point); |
---|
393 | |
---|
394 | } |
---|
395 | |
---|
396 | |
---|
397 | //---------------------------------------------------------------------- |
---|
398 | void ClosestPair2D::replace_many( |
---|
399 | const std::vector<unsigned int> & IDs_to_remove, |
---|
400 | const std::vector<Coord2D> & new_positions, |
---|
401 | std::vector<unsigned int> & new_IDs) { |
---|
402 | |
---|
403 | // deletion from tree... |
---|
404 | for (unsigned int i = 0; i < IDs_to_remove.size(); i++) { |
---|
405 | _remove_from_search_tree(& (_points[IDs_to_remove[i]])); |
---|
406 | } |
---|
407 | |
---|
408 | // insertion into tree |
---|
409 | new_IDs.resize(0); |
---|
410 | for (unsigned int i = 0; i < new_positions.size(); i++) { |
---|
411 | Point * new_point = _available_points.top(); |
---|
412 | _available_points.pop(); |
---|
413 | // set the point's coordinate |
---|
414 | new_point->coord = new_positions[i]; |
---|
415 | // now find it's neighbour in the search tree |
---|
416 | _insert_into_search_tree(new_point); |
---|
417 | // record the ID |
---|
418 | new_IDs.push_back(_ID(new_point)); |
---|
419 | } |
---|
420 | |
---|
421 | // the above statement labels certain points as needing "review" -- |
---|
422 | // deal with them... |
---|
423 | _deal_with_points_to_review(); |
---|
424 | |
---|
425 | } |
---|
426 | |
---|
427 | |
---|
428 | //---------------------------------------------------------------------- |
---|
429 | void ClosestPair2D::_insert_into_search_tree(Point * new_point) { |
---|
430 | |
---|
431 | // this point will have to have it's heap entry reviewed... |
---|
432 | _set_label(new_point, _review_heap_entry); |
---|
433 | |
---|
434 | // set the current distance to "infinity" |
---|
435 | new_point->neighbour_dist2 = numeric_limits<double>::max(); |
---|
436 | |
---|
437 | // establish how far we will be searching; |
---|
438 | unsigned int CP_range = min(_cp_search_range, size()-1); |
---|
439 | |
---|
440 | for (unsigned ishift = 0; ishift < _nshift; ishift++) { |
---|
441 | // create the shuffle |
---|
442 | Shuffle new_shuffle; |
---|
443 | _point2shuffle(*new_point, new_shuffle, _shifts[ishift]); |
---|
444 | |
---|
445 | // insert it into the tree |
---|
446 | circulator new_circ = _trees[ishift]->insert(new_shuffle); |
---|
447 | new_point->circ[ishift] = new_circ; |
---|
448 | |
---|
449 | // now get hold of the right and left edges of the region we will be |
---|
450 | // looking at (cf CCN28-43) |
---|
451 | circulator right_edge = new_circ; right_edge++; |
---|
452 | circulator left_edge = new_circ; |
---|
453 | for (unsigned int i = 0; i < CP_range; i++) {left_edge--;} |
---|
454 | |
---|
455 | // now |
---|
456 | do { |
---|
457 | Point * left_point = left_edge->point; |
---|
458 | Point * right_point = right_edge->point; |
---|
459 | |
---|
460 | // see if the new point is closer to the left-edge than the latter's |
---|
461 | // current neighbour |
---|
462 | double new_dist2 = left_point->distance2(*new_point); |
---|
463 | if (new_dist2 < left_point->neighbour_dist2) { |
---|
464 | left_point->neighbour_dist2 = new_dist2; |
---|
465 | left_point->neighbour = new_point; |
---|
466 | _add_label(left_point, _review_heap_entry); |
---|
467 | } |
---|
468 | |
---|
469 | // see if the right-point is closer to the new point than it's current |
---|
470 | // neighbour |
---|
471 | new_dist2 = new_point->distance2(*right_point); |
---|
472 | if (new_dist2 < new_point->neighbour_dist2) { |
---|
473 | new_point->neighbour_dist2 = new_dist2; |
---|
474 | new_point->neighbour = right_point; |
---|
475 | } |
---|
476 | |
---|
477 | // if the right-edge point was the left-edge's neighbour, then |
---|
478 | // then it's just gone off-radar and the left-point will need to |
---|
479 | // have its neighbour recalculated [actually, this is overdoing |
---|
480 | // it a little, since right point may be an less "distant" |
---|
481 | // (circulator distance) in one of the other shifts -- but not |
---|
482 | // sure how to deal with this...] |
---|
483 | if (left_point->neighbour == right_point) { |
---|
484 | _add_label(left_point, _review_neighbour); |
---|
485 | } |
---|
486 | |
---|
487 | // shift the left and right edges until left edge hits new_circ |
---|
488 | right_edge++; |
---|
489 | } while (++left_edge != new_circ); |
---|
490 | } |
---|
491 | } |
---|
492 | |
---|
493 | FASTJET_END_NAMESPACE |
---|
494 | |
---|