Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/src/ClosestPair2D.cc@ 38

Last change on this file since 38 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

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