source: HiSusy/trunk/Delphes-3.0.0/external/fastjet/DnnPlane.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 14.5 KB
Line 
1//STARTHEADER
2// $Id: DnnPlane.cc 859 2012-11-28 01:49:23Z pavel $
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
30#ifndef DROP_CGAL // in case we do not have the code for CGAL
31
32#include<set>
33#include<list>
34#include "fastjet/internal/DnnPlane.hh"
35using namespace std;
36
37FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
38
39
40/// Initialiser from a set of points on an Eta-Phi plane, where both
41/// eta and phi can have arbitrary ranges
42DnnPlane::DnnPlane(const vector<EtaPhi> & input_points, 
43                   const bool & verbose ) {
44
45  _verbose = verbose;
46  int n = input_points.size();
47 
48  // construct Voronoi diagram in such a way as to get the vertex handles
49  // and remember to set CGAL info with the index of the vertex
50  SuperVertex sv;
51  for (int i = 0; i < n; i++) {
52    sv.vertex = 
53       _TR.insert(Point(input_points[i].first, input_points[i].second));
54
55    // we are not up to dealing with coincident vertices, so make
56    // sure the user knows!
57    _CrashIfVertexPresent(sv.vertex, i);
58   
59    // we need to assicate an index to each vertex -- thus when we get
60    // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
61    // able to figure out which particle it corresponded to.
62    sv.vertex->info() = i;
63    _supervertex.push_back(sv);   
64  }
65
66  // label infinite vertex info with negative index
67  _TR.infinite_vertex()->info() = INFINITE_VERTEX;
68
69  // set up the structure that holds nearest distances and neighbours
70  for (int j = 0; j < n; j++) {_SetNearest(j);}
71
72}
73
74
75//----------------------------------------------------------------------
76/// Crashes if the given vertex handle already exists. Otherwise
77/// it does the bookkeeping for future such tests
78void DnnPlane::_CrashIfVertexPresent(
79        const Vertex_handle & vertex, const int & its_index) {
80  if (!_crash_on_coincidence) return;
81
82  // vertices that do not have the same geometric position as any
83  // other vertex so far added have info().val() == NEW_VERTEX -- this
84  // is ensured by the InitialisedInt class, which forms the "info"
85  // part of our
86  // CGAL::Triangulation_vertex_base_with_info_2<InitialisedInt,K>
87  //
88  // If the vertex coincides with one that already exists, then
89  // info().val() it's info().val() will have been updated (in
90  // DNN:DNN) to be equal to a vertex "index".
91  if (vertex->info().val() != NEW_VERTEX) {
92    ostringstream err;
93    err << "ERROR in DnnPlane::_CrashIfVertexPresent"
94         <<endl << "Point "<<its_index<<" coincides with point "
95         <<vertex->info().val() << endl;
96    throw DnnError(err.str());
97  } 
98}
99
100
101//----------------------------------------------------------------------
102/// remove the points labelled by the vector indices_to_remove, and
103/// add the points specified by the vector points_to_add
104/// (corresponding indices will be calculated automatically); the
105/// idea behind this routine is that the points to be added will
106/// somehow be close to the one or other of the points being removed
107/// and this can be used by the implementation to provide hints for
108/// inserting the new points in whatever structure it is using.  In a
109/// kt-algorithm the points being added will be a result of a
110/// combination of the points to be removed -- hence the proximity
111/// is (more or less) guaranteed.
112void DnnPlane::RemoveAndAddPoints(
113                          const vector<int> & indices_to_remove,
114                          const vector<EtaPhi> & points_to_add,
115                          vector<int> & indices_added,
116                          vector<int> & indices_of_updated_neighbours) {
117
118
119  // build set of UNION of Voronoi neighbours of a pair of nearest
120  // neighbours
121  set<int> NeighbourUnion;
122  // later on it will be convenient to have access to a set (rather
123  // than vector) of indices being removed
124  set<int> indices_removed;
125
126  // for each of the indices to be removed add the voronoi neighbourhood to
127  // the NeighbourUnion set.
128  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
129    int index = indices_to_remove[ir];
130    indices_removed.insert(index);
131    if (_verbose) cout << " Starting  RemoveAndAddPoints" << endl; 
132    if (_verbose) cout << " point " << index << endl;                   
133    // have a circulators that will go round the Voronoi neighbours of
134    // _supervertex[index1].vertex
135    Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
136    Vertex_circulator done = vc;
137    do  {
138      // if a neighbouring vertex not the infinite vertex, then add it
139      // to our union of neighbouring vertices.
140      if (_verbose) cout << "examining " << vc->info().val() << endl;
141      if (vc->info().val() != INFINITE_VERTEX) {
142        // NB: from it=1 onwards occasionally it might already have
143        // been inserted -- but double insertion still leaves only one
144        // copy in the set, so there's no problem
145        NeighbourUnion.insert(vc->info().val());
146        if (_verbose) cout << "inserted " << vc->info().val() << endl;
147      } 
148    } while (++vc != done);
149  }
150 
151  if (_verbose) {
152    set<int>::iterator it = NeighbourUnion.begin();
153    cout << "Union of neighbours of combined points" << endl;
154    for ( ; it != NeighbourUnion.end(); ++it ) {
155      cout << *it << endl;
156    }
157  }
158
159  // update set, triangulation and supervertex info
160  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
161    int index = indices_to_remove[ir];
162
163    // NeighbourUnion should not contain the points to be removed
164    // (because later we will assume they still exist).
165    NeighbourUnion.erase(indices_to_remove[ir]);
166   
167    // points to be removed should also be eliminated from the
168    // triangulation and the supervertex structure should be updated
169    // to reflect the fact that the points are no longer valid.
170    _TR.remove(_supervertex[index].vertex);
171    _supervertex[index].vertex = NULL;
172  }
173
174  // add new point: give a "hint" to the inserter that
175  // the new point should be added close to old points -- the easiest way
176  // of getting this is to take a point from the NeighbourUnion AFTER we have
177  // removed point1, point2, and to get one of its incident faces.
178  //
179  // This hinting improves speed by c. 25% for 10^4 points because it
180  // avoids the need for a costly (sqrt{N}) location step (at least
181  // with a non-hierarchical triangulation -- with a hierarchical one,
182  // this step could be done away with, though there will still be a
183  // cost of O(ln N) to pay.
184  //
185  // For some reason inserting the point before the two removals
186  // slows things down by c. 25%. This importance of the order
187  // is not understood.
188  //
189  // At some point it might be worth trying to select the "nearest"
190  // of the various points in the neighbour union to avoid large
191  // steps in cases where we have 0..2pi periodicity and the first member
192  // of the neighbour union happens to be on the wrong side.
193  Face_handle face;
194  if (indices_to_remove.size() > 0) {
195    // face can only be found if there were points to remove in first place
196    face = _TR.incident_faces(
197                           _supervertex[*NeighbourUnion.begin()].vertex);}
198  // make sure the output arrays are empty
199  indices_added.clear();
200  indices_of_updated_neighbours.clear();
201  for (size_t ia = 0; ia < points_to_add.size(); ia++) {
202    SuperVertex sv;
203    _supervertex.push_back(sv);
204    int index = _supervertex.size()-1;
205    indices_added.push_back(index);
206
207    if (indices_to_remove.size() > 0) {
208      // be careful of using face (for location hinting) only when it exists
209      _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
210                                  points_to_add[ia].second),face);}
211    else { 
212      _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
213                                                    points_to_add[ia].second));
214    }
215    // we are not up to dealing with coincident vertices, so make
216    // sure the user knows!
217    _CrashIfVertexPresent(_supervertex[index].vertex, index);
218    _supervertex[index].vertex->info() = index;
219   
220    // first find nearest neighbour of "newpoint" (shorthand for
221    // _supervertex[index].vertex); while we're at it, for each of the
222    // voronoi neighbours, "D", of newpoint, examine whether newpoint is
223    // closer to "D" than D's current nearest neighbour -- when this
224    // occurs, put D into indices_of_updated_neighbours.
225    //
226    // manually put newpoint on indices_of_updated_neighbours
227    indices_of_updated_neighbours.push_back(index);
228    _SetAndUpdateNearest(index, indices_of_updated_neighbours);
229  }
230
231  // for Voronoi neighbours j of any of the removed points for which
232  // one of those removed points was the nearest neighbour,
233  // redetermine the nearest neighbour of j and add j onto the vector
234  // of indices_of_updated_neighbours.
235  set<int>::iterator it2 = NeighbourUnion.begin();
236  for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
237    int j = *it2;
238    // the if avoids the vertex at infinity, which gets a negative index
239    if( j != INFINITE_VERTEX ) {
240      // this is where we check if the nearest neighbour of j was one
241      // of the removed points
242      if (indices_removed.count(_supervertex[j].NNindex)) {
243        if (_verbose) cout << "j " << j << endl;
244        _SetNearest(j);
245        indices_of_updated_neighbours.push_back(j);
246        if (_verbose) cout << "NN of " << j << " : " 
247                          << _supervertex[j].NNindex
248                          << ", dist = " << _supervertex[j].NNdistance <<endl;
249      }
250    }
251  }
252
253}
254
255
256//----------------------------------------------------------------------
257/// Determines the index and distance of the nearest neighbour to
258/// point j and puts the information into the _supervertex entry for j.
259void DnnPlane::_SetNearest (const int & j) {
260  Vertex_handle current = _supervertex[j].vertex;
261  Vertex_circulator vc = _TR.incident_vertices(current);
262  Vertex_circulator done = vc;
263  double dist;
264  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
265  Vertex_handle nearest = _TR.infinite_vertex();
266
267  // when there is only one finite point left in the triangulation,
268  // there are no triangles. Presumably this is why voronoi returns
269  // NULL for the incident vertex circulator. Check if this is
270  // happening before circulating over it... (Otherwise it crashes
271  // when looking for neighbours of last point)
272  if (vc != NULL) do { 
273    if ( vc->info().val() != INFINITE_VERTEX) {
274      // find distance between j and its Voronoi neighbour (vc)
275      if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
276      dist = _euclid_distance(current->point(), vc->point());
277      // check if j is closer to vc than vc's currently registered
278      // nearest neighbour (and update things if it is)
279      if (dist < mindist) {
280        mindist = dist; nearest = vc; 
281        if (_verbose) cout << "nearer ";
282      } 
283      if (_verbose) cout << vc->point() << "; "<< dist << endl;
284    }
285  } while (++vc != done); // move on to next Voronoi neighbour
286  // set j's supervertex info about nearest neighbour
287  _supervertex[j].NNindex = nearest->info().val();
288  _supervertex[j].NNdistance = mindist;
289}
290
291//----------------------------------------------------------------------
292/// Determines and stores the nearest neighbour of j, and where
293/// necessary update the nearest-neighbour info of Voronoi neighbours
294/// of j;
295///
296/// For each voronoi neighbour D of j if the distance between j and D
297/// is less than D's own nearest neighbour, then update the
298/// nearest-neighbour info in D; push D's index onto
299/// indices_of_updated_neighbours
300///
301/// Note that j is NOT pushed onto indices_of_updated_neighbours --
302/// if you want it there, put it there yourself.
303///
304/// NB: note that we have _SetAndUpdateNearest as a completely
305///     separate routine from _SetNearest because we want to
306///     use one single ciruclation over voronoi neighbours to find the
307///     nearest neighbour and to update the voronoi neighbours if need
308///     be.
309void DnnPlane::_SetAndUpdateNearest(
310                          const int & j, 
311                          vector<int> & indices_of_updated_neighbours) {
312
313  Vertex_handle current = _supervertex[j].vertex;
314  Vertex_circulator vc = _TR.incident_vertices(current);
315  Vertex_circulator done = vc;
316  double dist;
317  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
318  Vertex_handle nearest = _TR.infinite_vertex();
319
320  // when there is only one finite point left in the triangulation,
321  // there are no triangles. Presumably this is why voronoi returns
322  // NULL for the incident vertex circulator. Check if this is
323  // happening before circulating over it... (Otherwise it crashes
324  // when looking for neighbours of last point)
325  if (vc != NULL) do { 
326    if (vc->info().val() != INFINITE_VERTEX) {
327      if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
328      // find distance between j and its Voronoi neighbour (vc)
329      dist = _euclid_distance(current->point(), vc->point());
330      // update the mindist if we are closer than anything found so far
331      if (dist < mindist) {
332        mindist = dist; nearest = vc; 
333        if (_verbose) cout << "nearer ";
334      } 
335      // find index corresponding to vc for easy manipulation
336      int vcindx = vc->info().val();
337      if (_verbose) cout << vc->point() << "; "<< dist << endl;
338      // check if j is closer to vc than vc's currently registered
339      // nearest neighbour (and update things if it is)
340      if (dist < _supervertex[vcindx].NNdistance) {
341        if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
342        _supervertex[vcindx].NNdistance = dist;
343        _supervertex[vcindx].NNindex = j;
344        indices_of_updated_neighbours.push_back(vcindx);
345      }
346    }
347  } while (++vc != done); // move on to next Voronoi neighbour
348  // set j's supervertex info about nearest neighbour
349  _supervertex[j].NNindex = nearest->info().val();
350  _supervertex[j].NNdistance = mindist;
351}
352
353FASTJET_END_NAMESPACE
354
355#endif //  DROP_CGAL
Note: See TracBrowser for help on using the repository browser.