FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ClosestPair2D.cc
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 
38 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
39 
40 const unsigned int twopow31 = 2147483648U;
41 
42 using namespace std;
43 
44 //----------------------------------------------------------------------
45 /// takes a point and sets a shuffle with the given shift...
46 void 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
64 bool 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 //----------------------------------------------------------------------
78 void 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 //----------------------------------------------------------------------=
181 void 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 //----------------------------------------------------------------------
195 inline 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 //----------------------------------------------------------------------
206 inline 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 //----------------------------------------------------------------------
217 void 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 //----------------------------------------------------------------------
233 void 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 //----------------------------------------------------------------------
299 void 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 //----------------------------------------------------------------------
349 unsigned 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 //----------------------------------------------------------------------
371 unsigned 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 //----------------------------------------------------------------------
403 void 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 //----------------------------------------------------------------------
434 void 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 
498 FASTJET_END_NAMESPACE
499 
bool floor_ln2_less(unsigned x, unsigned y)
returns true if floor(ln_base2(x)) < floor(ln_base2(y)), using Chan's neat trick...