29 #include "fastjet/internal/ClosestPair2D.hh" 36 FASTJET_BEGIN_NAMESPACE
38 const unsigned int huge_unsigned = 4294967295U;
39 const unsigned int twopow31 = 2147483648U;
45 void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle,
48 Coord2D renorm_point = (point.coord - _left_corner)/_range;
51 assert(renorm_point.x >=0);
52 assert(renorm_point.x <=1);
53 assert(renorm_point.y >=0);
54 assert(renorm_point.y <=1);
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;
63 bool ClosestPair2D::Shuffle::operator<(
const Shuffle & q)
const {
77 void ClosestPair2D::_initialize(
const std::vector<Coord2D> & positions,
78 const Coord2D & left_corner,
79 const Coord2D & right_corner,
80 unsigned int max_size) {
82 unsigned int n_positions = positions.size();
83 assert(max_size >= n_positions);
88 _points.resize(max_size);
91 for (
unsigned int i = n_positions; i < max_size; i++) {
92 _available_points.push(&(_points[i]));
95 _left_corner = left_corner;
96 _range = max((right_corner.x - left_corner.x),
97 (right_corner.y - left_corner.y));
101 vector<Shuffle> shuffles(n_positions);
102 for (
unsigned int i = 0; i < n_positions; i++) {
104 _points[i].coord = positions[i];
105 _points[i].neighbour_dist2 = numeric_limits<double>::max();
106 _points[i].review_flag = 0;
109 _point2shuffle(_points[i], shuffles[i], 0);
113 for (
unsigned ishift = 0; ishift < _nshift; ishift++) {
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];}
129 _cp_search_range = 30;
130 _points_under_review.reserve(_nshift * _cp_search_range);
133 for (
unsigned int ishift = 0; ishift < _nshift; ishift++) {
137 unsigned rel_shift = _rel_shifts[ishift];
138 for (
unsigned int i = 0; i < shuffles.size(); i++) {
139 shuffles[i] += rel_shift; }
143 sort(shuffles.begin(), shuffles.end());
146 _trees[ishift] = auto_ptr<Tree>(
new Tree(shuffles, max_size));
149 circulator circ = _trees[ishift]->somewhere(), start=circ;
151 unsigned int CP_range = min(_cp_search_range, n_positions-1);
153 Point * this_point = circ->point;
155 this_point->circ[ishift] = circ;
157 circulator other = circ;
158 for (
unsigned i=0; i < CP_range; i++) {
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;
166 }
while (++circ != start);
171 vector<double> mindists2(n_positions);
172 for (
unsigned int i = 0; i < n_positions; i++) {
173 mindists2[i] = _points[i].neighbour_dist2;}
175 _heap = auto_ptr<MinHeap>(
new MinHeap(mindists2, max_size));
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) std::swap(ID1,ID2);
190 inline void ClosestPair2D::_add_label(Point * point,
unsigned int review_flag) {
194 if (point->review_flag == 0) _points_under_review.push_back(point);
197 point->review_flag |= review_flag;
201 inline void ClosestPair2D::_set_label(Point * point,
unsigned int review_flag) {
205 if (point->review_flag == 0) _points_under_review.push_back(point);
208 point->review_flag = review_flag;
212 void ClosestPair2D::remove(
unsigned int ID) {
216 Point * point_to_remove = & (_points[ID]);
219 _remove_from_search_tree(point_to_remove);
223 _deal_with_points_to_review();
228 void ClosestPair2D::_remove_from_search_tree(Point * point_to_remove) {
233 _available_points.push(point_to_remove);
236 _set_label(point_to_remove, _remove_heap_entry);
242 unsigned int CP_range = min(_cp_search_range, size()-1);
245 for (
unsigned int ishift = 0; ishift < _nshift; ishift++) {
248 circulator removed_circ = point_to_remove->circ[ishift];
249 circulator right_end = removed_circ.next();
251 _trees[ishift]->remove(removed_circ);
254 circulator left_end = right_end, orig_right_end = right_end;
255 for (
unsigned int i = 0; i < CP_range; i++) {left_end--;}
257 if (size()-1 < _cp_search_range) {
262 left_end--; right_end--;
269 Point * left_point = left_end->point;
273 if (left_point->neighbour == point_to_remove) {
275 _add_label(left_point, _review_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;
283 _add_label(left_point, _review_heap_entry);
287 }
while (++left_end != orig_right_end);
294 void ClosestPair2D::_deal_with_points_to_review() {
298 unsigned int CP_range = min(_cp_search_range, size()-1);
302 while(_points_under_review.size() > 0) {
304 Point * this_point = _points_under_review.back();
306 _points_under_review.pop_back();
308 if (this_point->review_flag & _remove_heap_entry) {
310 assert(!(this_point->review_flag ^ _remove_heap_entry));
311 _heap->remove(_ID(this_point));
315 if (this_point->review_flag & _review_neighbour) {
316 this_point->neighbour_dist2 = numeric_limits<double>::max();
318 for (
unsigned int ishift = 0; ishift < _nshift; ishift++) {
319 circulator other = this_point->circ[ishift];
321 for (
unsigned i=0; i < CP_range; i++) {
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;
333 _heap->update(_ID(this_point), this_point->neighbour_dist2);
337 this_point->review_flag = 0;
344 unsigned int ClosestPair2D::insert(
const Coord2D & new_coord) {
347 assert(_available_points.size() > 0);
348 Point * new_point = _available_points.top();
349 _available_points.pop();
352 new_point->coord = new_coord;
355 _insert_into_search_tree(new_point);
359 _deal_with_points_to_review();
362 return _ID(new_point);
366 unsigned int ClosestPair2D::replace(
unsigned int ID1,
unsigned int ID2,
367 const Coord2D & position) {
370 Point * point_to_remove = & (_points[ID1]);
371 _remove_from_search_tree(point_to_remove);
373 point_to_remove = & (_points[ID2]);
374 _remove_from_search_tree(point_to_remove);
378 Point * new_point = _available_points.top();
379 _available_points.pop();
382 new_point->coord = position;
385 _insert_into_search_tree(new_point);
389 _deal_with_points_to_review();
392 return _ID(new_point);
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) {
404 for (
unsigned int i = 0; i < IDs_to_remove.size(); i++) {
405 _remove_from_search_tree(& (_points[IDs_to_remove[i]]));
410 for (
unsigned int i = 0; i < new_positions.size(); i++) {
411 Point * new_point = _available_points.top();
412 _available_points.pop();
414 new_point->coord = new_positions[i];
416 _insert_into_search_tree(new_point);
418 new_IDs.push_back(_ID(new_point));
423 _deal_with_points_to_review();
429 void ClosestPair2D::_insert_into_search_tree(Point * new_point) {
432 _set_label(new_point, _review_heap_entry);
435 new_point->neighbour_dist2 = numeric_limits<double>::max();
438 unsigned int CP_range = min(_cp_search_range, size()-1);
440 for (
unsigned ishift = 0; ishift < _nshift; ishift++) {
443 _point2shuffle(*new_point, new_shuffle, _shifts[ishift]);
446 circulator new_circ = _trees[ishift]->insert(new_shuffle);
447 new_point->circ[ishift] = new_circ;
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--;}
457 Point * left_point = left_edge->point;
458 Point * right_point = right_edge->point;
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);
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;
483 if (left_point->neighbour == right_point) {
484 _add_label(left_point, _review_neighbour);
489 }
while (++left_edge != new_circ);
493 FASTJET_END_NAMESPACE
bool floor_ln2_less(unsigned x, unsigned y)
returns true if floor(ln_base2(x)) < floor(ln_base2(y)), using Chan's neat trick...