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