Fork me on GitHub

source: git/external/fastjet/ClosestPair2D.cc@ 8cbe720

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8cbe720 was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

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