Fork me on GitHub

source: svn/trunk/external/fastjet/internal/ClosestPair2D.hh@ 1204

Last change on this file since 1204 was 859, checked in by Pavel Demin, 12 years ago

update fastjet to version 3.0.3

File size: 8.1 KB
Line 
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
39FASTJET_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
48class ClosestPair2D : public ClosestPair2DBase {
49public:
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
99private:
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
191class ClosestPair2D::Point {
192public:
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...
218inline 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...
226inline int ClosestPair2D::_ID(const Point * point) const {
227 return point - &(_points[0]);
228}
229
230
231//
232inline unsigned int ClosestPair2D::size() {
233 return _points.size() - _available_points.size();
234}
235
236
237
238FASTJET_END_NAMESPACE
239
240#endif // __FASTJET_CLOSESTPAIR2D__HH__
Note: See TracBrowser for help on using the repository browser.