FastJet  3.1.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ClosestPair2D.hh
1 //FJSTARTHEADER
2 // $Id: ClosestPair2D.hh 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 #ifndef __FASTJET_CLOSESTPAIR2D__HH__
32 #define __FASTJET_CLOSESTPAIR2D__HH__
33 
34 #include<vector>
35 #include<stack>
36 #include<iostream>
37 #include "fastjet/internal/ClosestPair2DBase.hh"
38 #include "fastjet/internal/SearchTree.hh"
39 #include "fastjet/internal/MinHeap.hh"
40 
41 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
42 
43 //----------------------------------------------------------------------
44 /// \if internal_doc
45 /// @ingroup internal
46 /// \class ClosestPair2D
47 /// concrete implementation for finding closest pairs in 2D -- will
48 /// use Chan's (hopefully efficient) shuffle based structures
49 /// \endif
50 class ClosestPair2D : public ClosestPair2DBase {
51 public:
52  /// constructor from a vector of 2D positions -- number of objects
53  /// after insertion and deletion must never exceed positions.size();
54  /// objects are given IDs that correspond to their index in the vector
55  /// of positions
56  ClosestPair2D(const std::vector<Coord2D> & positions,
57  const Coord2D & left_corner, const Coord2D & right_corner) {
58  _initialize(positions, left_corner, right_corner, positions.size());
59  };
60 
61  /// constructor which allows structure to grow beyond positions.size(), up
62  /// to max_size
63  ClosestPair2D(const std::vector<Coord2D> & positions,
64  const Coord2D & left_corner, const Coord2D & right_corner,
65  const unsigned int max_size) {
66  _initialize(positions, left_corner, right_corner, max_size);
67  };
68 
69  /// provides the IDs of the closest pair as well as the distance between
70  /// them
71  void closest_pair(unsigned int & ID1, unsigned int & ID2,
72  double & distance2) const;
73 
74  /// removes the entry labelled by ID from the object;
75  void remove(unsigned int ID);
76 
77  /// inserts the position into the closest pair structure and returns the
78  /// ID that has been allocated for the object.
79  unsigned int insert(const Coord2D &);
80 
81  /// removes ID1 and ID2 and inserts position, returning the ID
82  /// corresponding to position...
83  virtual unsigned int replace(unsigned int ID1, unsigned int ID2,
84  const Coord2D & position);
85 
86  /// replaces IDs_to_remove with points at the new_positions
87  /// indicating the IDs allocated to the new points in new_IDs
88  virtual void replace_many(const std::vector<unsigned int> & IDs_to_remove,
89  const std::vector<Coord2D> & new_positions,
90  std::vector<unsigned int> & new_IDs);
91 
92  // mostly for checking how things are working...
93  inline void print_tree_depths(std::ostream & outdev) const {
94  outdev << _trees[0]->max_depth() << " "
95  << _trees[1]->max_depth() << " "
96  << _trees[2]->max_depth() << "\n";
97  };
98 
99  unsigned int size();
100 
101 private:
102 
103  void _initialize(const std::vector<Coord2D> & positions,
104  const Coord2D & left_corner, const Coord2D & right_corner,
105  const unsigned int max_size);
106 
107  static const unsigned int _nshift = 3;
108 
109  class Point; // will be defined below
110 
111  /// since sets of three objects will crop up repeatedly, useful
112  /// to have a triplet class?
113  template<class T> class triplet {
114  public:
115  inline const T & operator[](unsigned int i) const {return _contents[i];};
116  inline T & operator[](unsigned int i) {return _contents[i];};
117  private:
118  T _contents[_nshift];
119  };
120 
121 
122  /// class that will take care of ordering of shuffles for us
123  class Shuffle {
124  public:
125  unsigned int x, y;
126  Point * point;
127  bool operator<(const Shuffle &) const;
128  void operator+=(unsigned int shift) {x += shift; y+= shift;};
129  };
130 
131  typedef SearchTree<Shuffle> Tree;
132  typedef Tree::circulator circulator;
133  typedef Tree::const_circulator const_circulator;
134 
135 
136  triplet<std::auto_ptr<Tree> > _trees;
137  std::auto_ptr<MinHeap> _heap;
138  std::vector<Point> _points;
139  std::stack<Point *> _available_points;
140 
141  /// points that are "under review" in some way
142  std::vector<Point *> _points_under_review;
143 
144  // different statuses for review
145  static const unsigned int _remove_heap_entry = 1;
146  static const unsigned int _review_heap_entry = 2;
147  static const unsigned int _review_neighbour = 4;
148 
149  /// add a label to a point as to the nature of review needed
150  /// (includes adding it to list of points needing review) [doesn't
151  /// affect other labels already set for the point]
152  void _add_label(Point * point, unsigned int review_flag);
153 
154  /// sets the label for the point to be exclusively this
155  /// review flag (and adds it to list of points needing review
156  /// if not already there)
157  void _set_label(Point * point, unsigned int review_flag);
158 
159  /// for all entries of the _points_under_review[] vector, carry out
160  /// the actions indicated by its review flag; the points are
161  /// then removed from _points_under_review[] and their flags
162  /// set to zero
163  void _deal_with_points_to_review();
164 
165  /// carry out the search-tree related operations of point removal
166  void _remove_from_search_tree(Point * point_to_remove);
167 
168  /// carry out the search-tree related operations of point insertion
169  void _insert_into_search_tree(Point * new_point);
170 
171  /// takes a point and creates a shuffle with the given shift
172  void _point2shuffle(Point & , Shuffle & , unsigned int shift);
173 
174  /// pieces needed for converting coordinates to integer
175  Coord2D _left_corner;
176  double _range;
177 
178  int _ID(const Point *) const;
179 
180  triplet<unsigned int> _shifts; // absolute shifts
181  triplet<unsigned int> _rel_shifts; // shifts relative to previous shift
182 
183  unsigned int _cp_search_range;
184 };
185 
186 
187 //----------------------------------------------------------------------
188 /// \if internal_doc
189 /// @ingroup internal
190 /// \class ClosestPair2D::Point
191 /// class for representing all info needed about a point
192 /// \endif
193 class ClosestPair2D::Point {
194 public:
195  /// the point's coordinates
196  Coord2D coord;
197  /// a pointer to its closest neighbour in our structure
198  Point * neighbour;
199  /// the corresponding squared distance
200  double neighbour_dist2;
201  /// circulators for each of the shifts of the shuffles
202  triplet<circulator> circ;
203 
204  /// indicates that something special is currently happening to this point
205  unsigned int review_flag;
206 
207  /// returns the distance between two of these objects
208  double distance2(const Point & other) const {
209  return coord.distance2(other.coord);
210  };
211 
212  /// creates a shuffle for us with a given shift
213  //void set_shuffle(Shuffle & shuffle);
214 };
215 
216 
217 //----------------------------------------------------------------------
218 /// returns true if floor(ln_base2(x)) < floor(ln_base2(y)), using
219 /// Chan's neat trick...
220 inline bool floor_ln2_less(unsigned x, unsigned y) {
221  if (x>y) return false;
222  return (x < (x^y)); // beware of operator precedence...
223 }
224 
225 
226 //----------------------------------------------------------------------
227 /// returns the ID for the specified point...
228 inline int ClosestPair2D::_ID(const Point * point) const {
229  return point - &(_points[0]);
230 }
231 
232 
233 //
234 inline unsigned int ClosestPair2D::size() {
235  return _points.size() - _available_points.size();
236 }
237 
238 
239 
240 FASTJET_END_NAMESPACE
241 
242 #endif // __FASTJET_CLOSESTPAIR2D__HH__
bool floor_ln2_less(unsigned x, unsigned y)
returns true if floor(ln_base2(x)) < floor(ln_base2(y)), using Chan's neat trick...
bool operator<(SharedPtr< T > const &t, SharedPtr< U > const &u)
comparison: orgering
Definition: SharedPtr.hh:378