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