00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/internal/ClosestPair2D.hh"
00032
00033 #include<limits>
00034 #include<iostream>
00035 #include<iomanip>
00036
00037 FASTJET_BEGIN_NAMESPACE
00038
00039 const unsigned int huge_unsigned = 4294967295U;
00040 const unsigned int twopow31 = 2147483648U;
00041
00042 using namespace std;
00043
00044
00046 void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle,
00047 unsigned int shift) {
00048
00049 Coord2D renorm_point = (point.coord - _left_corner)/_range;
00050
00051
00052 assert(renorm_point.x >=0);
00053 assert(renorm_point.x <=1);
00054 assert(renorm_point.y >=0);
00055 assert(renorm_point.y <=1);
00056
00057 shuffle.x = static_cast<unsigned int>(twopow31 * renorm_point.x) + shift;
00058 shuffle.y = static_cast<unsigned int>(twopow31 * renorm_point.y) + shift;
00059 shuffle.point = &point;
00060 }
00061
00062
00064 bool ClosestPair2D::Shuffle::operator<(const Shuffle & q) const {
00065
00066 if (floor_ln2_less(x ^ q.x, y ^ q.y)) {
00067
00068 return (y < q.y);
00069 } else {
00070
00071 return (x < q.x);
00072 }
00073 }
00074
00075
00076
00077
00078 void ClosestPair2D::_initialize(const std::vector<Coord2D> & positions,
00079 const Coord2D & left_corner,
00080 const Coord2D & right_corner,
00081 unsigned int max_size) {
00082
00083 unsigned int n_positions = positions.size();
00084 assert(max_size >= n_positions);
00085
00086
00087
00088
00089 _points.resize(max_size);
00090
00091
00092 for (unsigned int i = n_positions; i < max_size; i++) {
00093 _available_points.push(&(_points[i]));
00094 }
00095
00096 _left_corner = left_corner;
00097 _range = max((right_corner.x - left_corner.x),
00098 (right_corner.y - left_corner.y));
00099
00100
00101
00102 vector<Shuffle> shuffles(n_positions);
00103 for (unsigned int i = 0; i < n_positions; i++) {
00104
00105 _points[i].coord = positions[i];
00106 _points[i].neighbour_dist2 = numeric_limits<double>::max();
00107 _points[i].review_flag = 0;
00108
00109
00110 _point2shuffle(_points[i], shuffles[i], 0);
00111 }
00112
00113
00114 for (unsigned ishift = 0; ishift < _nshift; ishift++) {
00115
00116
00117 _shifts[ishift] = static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift);
00118 if (ishift == 0) {_rel_shifts[ishift] = 0;}
00119 else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];}
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 _cp_search_range = 30;
00131 _points_under_review.reserve(_nshift * _cp_search_range);
00132
00133
00134 for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
00135
00136
00137 if (ishift > 0) {
00138 unsigned rel_shift = _rel_shifts[ishift];
00139 for (unsigned int i = 0; i < shuffles.size(); i++) {
00140 shuffles[i] += rel_shift; }
00141 }
00142
00143
00144 sort(shuffles.begin(), shuffles.end());
00145
00146
00147 _trees[ishift] = auto_ptr<Tree>(new Tree(shuffles, max_size));
00148
00149
00150 circulator circ = _trees[ishift]->somewhere(), start=circ;
00151
00152 unsigned int CP_range = min(_cp_search_range, n_positions-1);
00153 do {
00154 Point * this_point = circ->point;
00155
00156 this_point->circ[ishift] = circ;
00157
00158 circulator other = circ;
00159 for (unsigned i=0; i < CP_range; i++) {
00160 ++other;
00161 double dist2 = this_point->distance2(*other->point);
00162 if (dist2 < this_point->neighbour_dist2) {
00163 this_point->neighbour_dist2 = dist2;
00164 this_point->neighbour = other->point;
00165 }
00166 }
00167 } while (++circ != start);
00168
00169 }
00170
00171
00172 vector<double> mindists2(n_positions);
00173 for (unsigned int i = 0; i < n_positions; i++) {
00174 mindists2[i] = _points[i].neighbour_dist2;}
00175
00176 _heap = auto_ptr<MinHeap>(new MinHeap(mindists2, max_size));
00177 }
00178
00179
00180
00181 void ClosestPair2D::closest_pair(unsigned int & ID1, unsigned int & ID2,
00182 double & distance2) const {
00183 ID1 = _heap->minloc();
00184 ID2 = _ID(_points[ID1].neighbour);
00185 distance2 = _points[ID1].neighbour_dist2;
00186 if (ID1 > ID2) swap(ID1,ID2);
00187 }
00188
00189
00190
00191 inline void ClosestPair2D::_add_label(Point * point, unsigned int review_flag) {
00192
00193
00194
00195 if (point->review_flag == 0) _points_under_review.push_back(point);
00196
00197
00198 point->review_flag |= review_flag;
00199 }
00200
00201
00202 inline void ClosestPair2D::_set_label(Point * point, unsigned int review_flag) {
00203
00204
00205
00206 if (point->review_flag == 0) _points_under_review.push_back(point);
00207
00208
00209 point->review_flag = review_flag;
00210 }
00211
00212
00213 void ClosestPair2D::remove(unsigned int ID) {
00214
00215
00216
00217 Point * point_to_remove = & (_points[ID]);
00218
00219
00220 _remove_from_search_tree(point_to_remove);
00221
00222
00223
00224 _deal_with_points_to_review();
00225 }
00226
00227
00228
00229 void ClosestPair2D::_remove_from_search_tree(Point * point_to_remove) {
00230
00231
00232
00233
00234 _available_points.push(point_to_remove);
00235
00236
00237 _set_label(point_to_remove, _remove_heap_entry);
00238
00239
00240
00241
00242
00243 unsigned int CP_range = min(_cp_search_range, size()-1);
00244
00245
00246 for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
00247
00248
00249 circulator removed_circ = point_to_remove->circ[ishift];
00250 circulator right_end = removed_circ.next();
00251
00252 _trees[ishift]->remove(removed_circ);
00253
00254
00255 circulator left_end = right_end, orig_right_end = right_end;
00256 for (unsigned int i = 0; i < CP_range; i++) {left_end--;}
00257
00258 if (size()-1 < _cp_search_range) {
00259
00260
00261
00262
00263 left_end--; right_end--;
00264 }
00265
00266
00267
00268
00269 do {
00270 Point * left_point = left_end->point;
00271
00272
00273
00274 if (left_point->neighbour == point_to_remove) {
00275
00276 _add_label(left_point, _review_neighbour);
00277 } else {
00278
00279 double dist2 = left_point->distance2(*right_end->point);
00280 if (dist2 < left_point->neighbour_dist2) {
00281 left_point->neighbour = right_end->point;
00282 left_point->neighbour_dist2 = dist2;
00283
00284 _add_label(left_point, _review_heap_entry);
00285 }
00286 }
00287 ++right_end;
00288 } while (++left_end != orig_right_end);
00289 }
00290
00291 }
00292
00293
00294
00295 void ClosestPair2D::_deal_with_points_to_review() {
00296
00297
00298
00299 unsigned int CP_range = min(_cp_search_range, size()-1);
00300
00301
00302
00303 while(_points_under_review.size() > 0) {
00304
00305 Point * this_point = _points_under_review.back();
00306
00307 _points_under_review.pop_back();
00308
00309 if (this_point->review_flag & _remove_heap_entry) {
00310
00311 assert(!(this_point->review_flag ^ _remove_heap_entry));
00312 _heap->remove(_ID(this_point));
00313 }
00314
00315 else {
00316 if (this_point->review_flag & _review_neighbour) {
00317 this_point->neighbour_dist2 = numeric_limits<double>::max();
00318
00319 for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
00320 circulator other = this_point->circ[ishift];
00321
00322 for (unsigned i=0; i < CP_range; i++) {
00323 ++other;
00324 double dist2 = this_point->distance2(*other->point);
00325 if (dist2 < this_point->neighbour_dist2) {
00326 this_point->neighbour_dist2 = dist2;
00327 this_point->neighbour = other->point;
00328 }
00329 }
00330 }
00331 }
00332
00333
00334 _heap->update(_ID(this_point), this_point->neighbour_dist2);
00335 }
00336
00337
00338 this_point->review_flag = 0;
00339
00340 }
00341
00342 }
00343
00344
00345 unsigned int ClosestPair2D::insert(const Coord2D & new_coord) {
00346
00347
00348 assert(_available_points.size() > 0);
00349 Point * new_point = _available_points.top();
00350 _available_points.pop();
00351
00352
00353 new_point->coord = new_coord;
00354
00355
00356 _insert_into_search_tree(new_point);
00357
00358
00359
00360 _deal_with_points_to_review();
00361
00362
00363 return _ID(new_point);
00364 }
00365
00366
00367 unsigned int ClosestPair2D::replace(unsigned int ID1, unsigned int ID2,
00368 const Coord2D & position) {
00369
00370
00371 Point * point_to_remove = & (_points[ID1]);
00372 _remove_from_search_tree(point_to_remove);
00373
00374 point_to_remove = & (_points[ID2]);
00375 _remove_from_search_tree(point_to_remove);
00376
00377
00378
00379 Point * new_point = _available_points.top();
00380 _available_points.pop();
00381
00382
00383 new_point->coord = position;
00384
00385
00386 _insert_into_search_tree(new_point);
00387
00388
00389
00390 _deal_with_points_to_review();
00391
00392
00393 return _ID(new_point);
00394
00395 }
00396
00397
00398
00399 void ClosestPair2D::replace_many(
00400 const std::vector<unsigned int> & IDs_to_remove,
00401 const std::vector<Coord2D> & new_positions,
00402 std::vector<unsigned int> & new_IDs) {
00403
00404
00405 for (unsigned int i = 0; i < IDs_to_remove.size(); i++) {
00406 _remove_from_search_tree(& (_points[IDs_to_remove[i]]));
00407 }
00408
00409
00410 new_IDs.resize(0);
00411 for (unsigned int i = 0; i < new_positions.size(); i++) {
00412 Point * new_point = _available_points.top();
00413 _available_points.pop();
00414
00415 new_point->coord = new_positions[i];
00416
00417 _insert_into_search_tree(new_point);
00418
00419 new_IDs.push_back(_ID(new_point));
00420 }
00421
00422
00423
00424 _deal_with_points_to_review();
00425
00426 }
00427
00428
00429
00430 void ClosestPair2D::_insert_into_search_tree(Point * new_point) {
00431
00432
00433 _set_label(new_point, _review_heap_entry);
00434
00435
00436 new_point->neighbour_dist2 = numeric_limits<double>::max();
00437
00438
00439 unsigned int CP_range = min(_cp_search_range, size()-1);
00440
00441 for (unsigned ishift = 0; ishift < _nshift; ishift++) {
00442
00443 Shuffle new_shuffle;
00444 _point2shuffle(*new_point, new_shuffle, _shifts[ishift]);
00445
00446
00447 circulator new_circ = _trees[ishift]->insert(new_shuffle);
00448 new_point->circ[ishift] = new_circ;
00449
00450
00451
00452 circulator right_edge = new_circ; right_edge++;
00453 circulator left_edge = new_circ;
00454 for (unsigned int i = 0; i < CP_range; i++) {left_edge--;}
00455
00456
00457 do {
00458 Point * left_point = left_edge->point;
00459 Point * right_point = right_edge->point;
00460
00461
00462
00463 double new_dist2 = left_point->distance2(*new_point);
00464 if (new_dist2 < left_point->neighbour_dist2) {
00465 left_point->neighbour_dist2 = new_dist2;
00466 left_point->neighbour = new_point;
00467 _add_label(left_point, _review_heap_entry);
00468 }
00469
00470
00471
00472 new_dist2 = new_point->distance2(*right_point);
00473 if (new_dist2 < new_point->neighbour_dist2) {
00474 new_point->neighbour_dist2 = new_dist2;
00475 new_point->neighbour = right_point;
00476 }
00477
00478
00479
00480
00481
00482
00483
00484 if (left_point->neighbour == right_point) {
00485 _add_label(left_point, _review_neighbour);
00486 }
00487
00488
00489 right_edge++;
00490 } while (++left_edge != new_circ);
00491 }
00492 }
00493
00494 FASTJET_END_NAMESPACE
00495