FastJet 3.0.5
DnnPlane.cc
00001 //STARTHEADER
00002 // $Id: DnnPlane.cc 2577 2011-09-13 15:11:38Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00031 
00032 #include<set>
00033 #include<list>
00034 #include "fastjet/internal/DnnPlane.hh"
00035 using namespace std;
00036 
00037 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00038 
00039 
00040 /// Initialiser from a set of points on an Eta-Phi plane, where both
00041 /// eta and phi can have arbitrary ranges
00042 DnnPlane::DnnPlane(const vector<EtaPhi> & input_points, 
00043                    const bool & verbose ) {
00044 
00045   _verbose = verbose;
00046   int n = input_points.size();
00047   
00048   // construct Voronoi diagram in such a way as to get the vertex handles
00049   // and remember to set CGAL info with the index of the vertex
00050   SuperVertex sv;
00051   for (int i = 0; i < n; i++) {
00052     sv.vertex = 
00053        _TR.insert(Point(input_points[i].first, input_points[i].second));
00054 
00055     // we are not up to dealing with coincident vertices, so make 
00056     // sure the user knows!
00057     _CrashIfVertexPresent(sv.vertex, i);
00058     
00059     // we need to assicate an index to each vertex -- thus when we get
00060     // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
00061     // able to figure out which particle it corresponded to.
00062     sv.vertex->info() = i;
00063     _supervertex.push_back(sv);    
00064   }
00065 
00066   // label infinite vertex info with negative index 
00067   _TR.infinite_vertex()->info() = INFINITE_VERTEX;
00068 
00069   // set up the structure that holds nearest distances and neighbours
00070   for (int j = 0; j < n; j++) {_SetNearest(j);}
00071 
00072 }
00073 
00074 
00075 //----------------------------------------------------------------------
00076 /// Crashes if the given vertex handle already exists. Otherwise
00077 /// it does the bookkeeping for future such tests
00078 void DnnPlane::_CrashIfVertexPresent(
00079         const Vertex_handle & vertex, const int & its_index) {
00080   if (!_crash_on_coincidence) return;
00081 
00082   // vertices that do not have the same geometric position as any
00083   // other vertex so far added have info().val() == NEW_VERTEX -- this
00084   // is ensured by the InitialisedInt class, which forms the "info"
00085   // part of our
00086   // CGAL::Triangulation_vertex_base_with_info_2<InitialisedInt,K>
00087   //
00088   // If the vertex coincides with one that already exists, then
00089   // info().val() it's info().val() will have been updated (in
00090   // DNN:DNN) to be equal to a vertex "index".
00091   if (vertex->info().val() != NEW_VERTEX) {
00092     ostringstream err;
00093     err << "ERROR in DnnPlane::_CrashIfVertexPresent"
00094          <<endl << "Point "<<its_index<<" coincides with point "
00095          <<vertex->info().val() << endl;
00096     throw DnnError(err.str());
00097   } 
00098 }
00099 
00100 
00101 //----------------------------------------------------------------------
00102 /// remove the points labelled by the vector indices_to_remove, and
00103 /// add the points specified by the vector points_to_add
00104 /// (corresponding indices will be calculated automatically); the
00105 /// idea behind this routine is that the points to be added will
00106 /// somehow be close to the one or other of the points being removed
00107 /// and this can be used by the implementation to provide hints for
00108 /// inserting the new points in whatever structure it is using.  In a
00109 /// kt-algorithm the points being added will be a result of a
00110 /// combination of the points to be removed -- hence the proximity
00111 /// is (more or less) guaranteed.
00112 void DnnPlane::RemoveAndAddPoints(
00113                           const vector<int> & indices_to_remove,
00114                           const vector<EtaPhi> & points_to_add,
00115                           vector<int> & indices_added,
00116                           vector<int> & indices_of_updated_neighbours) {
00117 
00118 
00119   // build set of UNION of Voronoi neighbours of a pair of nearest
00120   // neighbours
00121   set<int> NeighbourUnion;
00122   // later on it will be convenient to have access to a set (rather
00123   // than vector) of indices being removed
00124   set<int> indices_removed;
00125 
00126   // for each of the indices to be removed add the voronoi neighbourhood to
00127   // the NeighbourUnion set.
00128   for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
00129     int index = indices_to_remove[ir];
00130     indices_removed.insert(index);
00131     if (_verbose) cout << " Starting  RemoveAndAddPoints" << endl; 
00132     if (_verbose) cout << " point " << index << endl;                    
00133     // have a circulators that will go round the Voronoi neighbours of
00134     // _supervertex[index1].vertex
00135     Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
00136     Vertex_circulator done = vc;
00137     do  {
00138       // if a neighbouring vertex not the infinite vertex, then add it
00139       // to our union of neighbouring vertices.
00140       if (_verbose) cout << "examining " << vc->info().val() << endl;
00141       if (vc->info().val() != INFINITE_VERTEX) {
00142         // NB: from it=1 onwards occasionally it might already have
00143         // been inserted -- but double insertion still leaves only one
00144         // copy in the set, so there's no problem
00145         NeighbourUnion.insert(vc->info().val());
00146         if (_verbose) cout << "inserted " << vc->info().val() << endl;
00147       } 
00148     } while (++vc != done);
00149   }
00150   
00151   if (_verbose) {
00152     set<int>::iterator it = NeighbourUnion.begin();
00153     cout << "Union of neighbours of combined points" << endl;
00154     for ( ; it != NeighbourUnion.end(); ++it ) {
00155       cout << *it << endl;
00156     }
00157   }
00158 
00159   // update set, triangulation and supervertex info
00160   for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
00161     int index = indices_to_remove[ir];
00162 
00163     // NeighbourUnion should not contain the points to be removed
00164     // (because later we will assume they still exist).
00165     NeighbourUnion.erase(indices_to_remove[ir]);
00166     
00167     // points to be removed should also be eliminated from the
00168     // triangulation and the supervertex structure should be updated
00169     // to reflect the fact that the points are no longer valid.
00170     _TR.remove(_supervertex[index].vertex);
00171     _supervertex[index].vertex = NULL;
00172   }
00173 
00174   // add new point: give a "hint" to the inserter that
00175   // the new point should be added close to old points -- the easiest way 
00176   // of getting this is to take a point from the NeighbourUnion AFTER we have
00177   // removed point1, point2, and to get one of its incident faces.
00178   // 
00179   // This hinting improves speed by c. 25% for 10^4 points because it
00180   // avoids the need for a costly (sqrt{N}) location step (at least
00181   // with a non-hierarchical triangulation -- with a hierarchical one,
00182   // this step could be done away with, though there will still be a
00183   // cost of O(ln N) to pay.
00184   // 
00185   // For some reason inserting the point before the two removals
00186   // slows things down by c. 25%. This importance of the order
00187   // is not understood.
00188   //
00189   // At some point it might be worth trying to select the "nearest"
00190   // of the various points in the neighbour union to avoid large 
00191   // steps in cases where we have 0..2pi periodicity and the first member
00192   // of the neighbour union happens to be on the wrong side.
00193   Face_handle face;
00194   if (indices_to_remove.size() > 0) {
00195     // face can only be found if there were points to remove in first place
00196     face = _TR.incident_faces(
00197                            _supervertex[*NeighbourUnion.begin()].vertex);}
00198   // make sure the output arrays are empty
00199   indices_added.clear();
00200   indices_of_updated_neighbours.clear();
00201   for (size_t ia = 0; ia < points_to_add.size(); ia++) {
00202     SuperVertex sv;
00203     _supervertex.push_back(sv);
00204     int index = _supervertex.size()-1;
00205     indices_added.push_back(index);
00206 
00207     if (indices_to_remove.size() > 0) {
00208       // be careful of using face (for location hinting) only when it exists
00209       _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
00210                                   points_to_add[ia].second),face);}
00211     else { 
00212       _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
00213                                                     points_to_add[ia].second));
00214     }
00215     // we are not up to dealing with coincident vertices, so make 
00216     // sure the user knows!
00217     _CrashIfVertexPresent(_supervertex[index].vertex, index);
00218     _supervertex[index].vertex->info() = index;
00219     
00220     // first find nearest neighbour of "newpoint" (shorthand for
00221     // _supervertex[index].vertex); while we're at it, for each of the
00222     // voronoi neighbours, "D", of newpoint, examine whether newpoint is
00223     // closer to "D" than D's current nearest neighbour -- when this
00224     // occurs, put D into indices_of_updated_neighbours.
00225     // 
00226     // manually put newpoint on indices_of_updated_neighbours
00227     indices_of_updated_neighbours.push_back(index);
00228     _SetAndUpdateNearest(index, indices_of_updated_neighbours);
00229   }
00230 
00231   // for Voronoi neighbours j of any of the removed points for which
00232   // one of those removed points was the nearest neighbour,
00233   // redetermine the nearest neighbour of j and add j onto the vector
00234   // of indices_of_updated_neighbours.
00235   set<int>::iterator it2 = NeighbourUnion.begin();
00236   for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
00237     int j = *it2;
00238     // the if avoids the vertex at infinity, which gets a negative index
00239     if( j != INFINITE_VERTEX ) {
00240       // this is where we check if the nearest neighbour of j was one
00241       // of the removed points
00242       if (indices_removed.count(_supervertex[j].NNindex)) {
00243         if (_verbose) cout << "j " << j << endl;
00244         _SetNearest(j);
00245         indices_of_updated_neighbours.push_back(j);
00246         if (_verbose) cout << "NN of " << j << " : " 
00247                           << _supervertex[j].NNindex
00248                           << ", dist = " << _supervertex[j].NNdistance <<endl;
00249       }
00250     }
00251   }
00252 
00253 }
00254 
00255 
00256 //----------------------------------------------------------------------
00257 /// Determines the index and distance of the nearest neighbour to 
00258 /// point j and puts the information into the _supervertex entry for j.
00259 void DnnPlane::_SetNearest (const int & j) {
00260   Vertex_handle current = _supervertex[j].vertex;
00261   Vertex_circulator vc = _TR.incident_vertices(current);
00262   Vertex_circulator done = vc;
00263   double dist;
00264   double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
00265   Vertex_handle nearest = _TR.infinite_vertex();
00266 
00267   // when there is only one finite point left in the triangulation, 
00268   // there are no triangles. Presumably this is why voronoi returns
00269   // NULL for the incident vertex circulator. Check if this is
00270   // happening before circulating over it... (Otherwise it crashes
00271   // when looking for neighbours of last point)
00272   if (vc != NULL) do { 
00273     if ( vc->info().val() != INFINITE_VERTEX) {
00274       // find distance between j and its Voronoi neighbour (vc)
00275       if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
00276       dist = _euclid_distance(current->point(), vc->point());
00277       // check if j is closer to vc than vc's currently registered
00278       // nearest neighbour (and update things if it is)
00279       if (dist < mindist) {
00280         mindist = dist; nearest = vc; 
00281         if (_verbose) cout << "nearer ";
00282       } 
00283       if (_verbose) cout << vc->point() << "; "<< dist << endl;
00284     }
00285   } while (++vc != done); // move on to next Voronoi neighbour
00286   // set j's supervertex info about nearest neighbour
00287   _supervertex[j].NNindex = nearest->info().val();
00288   _supervertex[j].NNdistance = mindist;
00289 }
00290 
00291 //----------------------------------------------------------------------
00292 /// Determines and stores the nearest neighbour of j, and where
00293 /// necessary update the nearest-neighbour info of Voronoi neighbours
00294 /// of j;
00295 ///
00296 /// For each voronoi neighbour D of j if the distance between j and D
00297 /// is less than D's own nearest neighbour, then update the
00298 /// nearest-neighbour info in D; push D's index onto 
00299 /// indices_of_updated_neighbours
00300 ///
00301 /// Note that j is NOT pushed onto indices_of_updated_neighbours --
00302 /// if you want it there, put it there yourself.
00303 ///
00304 /// NB: note that we have _SetAndUpdateNearest as a completely
00305 ///     separate routine from _SetNearest because we want to
00306 ///     use one single ciruclation over voronoi neighbours to find the
00307 ///     nearest neighbour and to update the voronoi neighbours if need
00308 ///     be.
00309 void DnnPlane::_SetAndUpdateNearest(
00310                           const int & j, 
00311                           vector<int> & indices_of_updated_neighbours) {
00312 
00313   Vertex_handle current = _supervertex[j].vertex;
00314   Vertex_circulator vc = _TR.incident_vertices(current);
00315   Vertex_circulator done = vc;
00316   double dist;
00317   double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
00318   Vertex_handle nearest = _TR.infinite_vertex();
00319 
00320   // when there is only one finite point left in the triangulation, 
00321   // there are no triangles. Presumably this is why voronoi returns
00322   // NULL for the incident vertex circulator. Check if this is
00323   // happening before circulating over it... (Otherwise it crashes
00324   // when looking for neighbours of last point)
00325   if (vc != NULL) do { 
00326     if (vc->info().val() != INFINITE_VERTEX) {
00327       if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
00328       // find distance between j and its Voronoi neighbour (vc)
00329       dist = _euclid_distance(current->point(), vc->point());
00330       // update the mindist if we are closer than anything found so far
00331       if (dist < mindist) {
00332         mindist = dist; nearest = vc; 
00333         if (_verbose) cout << "nearer ";
00334       } 
00335       // find index corresponding to vc for easy manipulation
00336       int vcindx = vc->info().val();
00337       if (_verbose) cout << vc->point() << "; "<< dist << endl;
00338       // check if j is closer to vc than vc's currently registered
00339       // nearest neighbour (and update things if it is)
00340       if (dist < _supervertex[vcindx].NNdistance) {
00341         if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
00342         _supervertex[vcindx].NNdistance = dist;
00343         _supervertex[vcindx].NNindex = j;
00344         indices_of_updated_neighbours.push_back(vcindx);
00345       }
00346     }
00347   } while (++vc != done); // move on to next Voronoi neighbour
00348   // set j's supervertex info about nearest neighbour
00349   _supervertex[j].NNindex = nearest->info().val();
00350   _supervertex[j].NNdistance = mindist;
00351 }
00352 
00353 FASTJET_END_NAMESPACE
00354 
00355 #endif //  DROP_CGAL
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends