fastjet 2.4.5
DnnPlane.cc
Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: DnnPlane.cc 310 2006-08-24 10:52:50Z cacciari $
00003 //
00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam
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, write to the Free Software
00026 //  Foundation, Inc.:
00027 //      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00028 //----------------------------------------------------------------------
00029 //ENDHEADER
00030 
00031 
00032 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00033 
00034 #include<set>
00035 #include<list>
00036 #include "fastjet/internal/DnnPlane.hh"
00037 using namespace std;
00038 
00039 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00040 
00041 
00044 DnnPlane::DnnPlane(const vector<EtaPhi> & input_points, 
00045                    const bool & verbose ) {
00046 
00047   _verbose = verbose;
00048   int n = input_points.size();
00049   
00050   // construct Voronoi diagram in such a way as to get the vertex handles
00051   // and remember to set CGAL info with the index of the vertex
00052   SuperVertex sv;
00053   for (int i = 0; i < n; i++) {
00054     sv.vertex = 
00055        _TR.insert(Point(input_points[i].first, input_points[i].second));
00056 
00057     // we are not up to dealing with coincident vertices, so make 
00058     // sure the user knows!
00059     _CrashIfVertexPresent(sv.vertex, i);
00060     
00061     // we need to assicate an index to each vertex -- thus when we get
00062     // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
00063     // able to figure out which particle it corresponded to.
00064     sv.vertex->info() = i;
00065     _supervertex.push_back(sv);    
00066   }
00067 
00068   // label infinite vertex info with negative index 
00069   _TR.infinite_vertex()->info() = INFINITE_VERTEX;
00070 
00071   // set up the structure that holds nearest distances and neighbours
00072   for (int j = 0; j < n; j++) {_SetNearest(j);}
00073 
00074 }
00075 
00076 
00077 //----------------------------------------------------------------------
00080 void DnnPlane::_CrashIfVertexPresent(
00081         const Vertex_handle & vertex, const int & its_index) {
00082   if (!_crash_on_coincidence) return;
00083 
00084   // vertices that do not have the same geometric position as any
00085   // other vertex so far added have info().val() == NEW_VERTEX -- this
00086   // is ensured by the InitialisedInt class, which forms the "info"
00087   // part of our
00088   // CGAL::Triangulation_vertex_base_with_info_2<InitialisedInt,K>
00089   //
00090   // If the vertex coincides with one that already exists, then
00091   // info().val() it's info().val() will have been updated (in
00092   // DNN:DNN) to be equal to a vertex "index".
00093   if (vertex->info().val() != NEW_VERTEX) {
00094     ostringstream err;
00095     err << "ERROR in DnnPlane::_CrashIfVertexPresent"
00096          <<endl << "Point "<<its_index<<" coincides with point "
00097          <<vertex->info().val() << endl;
00098     throw DnnError(err.str());
00099   } 
00100 }
00101 
00102 
00103 //----------------------------------------------------------------------
00114 void DnnPlane::RemoveAndAddPoints(
00115                           const vector<int> & indices_to_remove,
00116                           const vector<EtaPhi> & points_to_add,
00117                           vector<int> & indices_added,
00118                           vector<int> & indices_of_updated_neighbours) {
00119 
00120 
00121   // build set of UNION of Voronoi neighbours of a pair of nearest
00122   // neighbours
00123   set<int> NeighbourUnion;
00124   // later on it will be convenient to have access to a set (rather
00125   // than vector) of indices being removed
00126   set<int> indices_removed;
00127 
00128   // for each of the indices to be removed add the voronoi neighbourhood to
00129   // the NeighbourUnion set.
00130   for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
00131     int index = indices_to_remove[ir];
00132     indices_removed.insert(index);
00133     if (_verbose) cout << " Starting  RemoveAndAddPoints" << endl; 
00134     if (_verbose) cout << " point " << index << endl;                    
00135     // have a circulators that will go round the Voronoi neighbours of
00136     // _supervertex[index1].vertex
00137     Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
00138     Vertex_circulator done = vc;
00139     do  {
00140       // if a neighbouring vertex not the infinite vertex, then add it
00141       // to our union of neighbouring vertices.
00142       if (_verbose) cout << "examining " << vc->info().val() << endl;
00143       if (vc->info().val() != INFINITE_VERTEX) {
00144         // NB: from it=1 onwards occasionally it might already have
00145         // been inserted -- but double insertion still leaves only one
00146         // copy in the set, so there's no problem
00147         NeighbourUnion.insert(vc->info().val());
00148         if (_verbose) cout << "inserted " << vc->info().val() << endl;
00149       } 
00150     } while (++vc != done);
00151   }
00152   
00153   if (_verbose) {
00154     set<int>::iterator it = NeighbourUnion.begin();
00155     cout << "Union of neighbours of combined points" << endl;
00156     for ( ; it != NeighbourUnion.end(); ++it ) {
00157       cout << *it << endl;
00158     }
00159   }
00160 
00161   // update set, triangulation and supervertex info
00162   for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
00163     int index = indices_to_remove[ir];
00164 
00165     // NeighbourUnion should not contain the points to be removed
00166     // (because later we will assume they still exist).
00167     NeighbourUnion.erase(indices_to_remove[ir]);
00168     
00169     // points to be removed should also be eliminated from the
00170     // triangulation and the supervertex structure should be updated
00171     // to reflect the fact that the points are no longer valid.
00172     _TR.remove(_supervertex[index].vertex);
00173     _supervertex[index].vertex = NULL;
00174   }
00175 
00176   // add new point: give a "hint" to the inserter that
00177   // the new point should be added close to old points -- the easiest way 
00178   // of getting this is to take a point from the NeighbourUnion AFTER we have
00179   // removed point1, point2, and to get one of its incident faces.
00180   // 
00181   // This hinting improves speed by c. 25% for 10^4 points because it
00182   // avoids the need for a costly (sqrt{N}) location step (at least
00183   // with a non-hierarchical triangulation -- with a hierarchical one,
00184   // this step could be done away with, though there will still be a
00185   // cost of O(ln N) to pay.
00186   // 
00187   // For some reason inserting the point before the two removals
00188   // slows things down by c. 25%. This importance of the order
00189   // is not understood.
00190   //
00191   // At some point it might be worth trying to select the "nearest"
00192   // of the various points in the neighbour union to avoid large 
00193   // steps in cases where we have 0..2pi periodicity and the first member
00194   // of the neighbour union happens to be on the wrong side.
00195   Face_handle face;
00196   if (indices_to_remove.size() > 0) {
00197     // face can only be found if there were points to remove in first place
00198     face = _TR.incident_faces(
00199                            _supervertex[*NeighbourUnion.begin()].vertex);}
00200   // make sure the output arrays are empty
00201   indices_added.clear();
00202   indices_of_updated_neighbours.clear();
00203   for (size_t ia = 0; ia < points_to_add.size(); ia++) {
00204     SuperVertex sv;
00205     _supervertex.push_back(sv);
00206     int index = _supervertex.size()-1;
00207     indices_added.push_back(index);
00208 
00209     if (indices_to_remove.size() > 0) {
00210       // be careful of using face (for location hinting) only when it exists
00211       _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
00212                                   points_to_add[ia].second),face);}
00213     else { 
00214       _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, 
00215                                                     points_to_add[ia].second));
00216     }
00217     // we are not up to dealing with coincident vertices, so make 
00218     // sure the user knows!
00219     _CrashIfVertexPresent(_supervertex[index].vertex, index);
00220     _supervertex[index].vertex->info() = index;
00221     
00222     // first find nearest neighbour of "newpoint" (shorthand for
00223     // _supervertex[index].vertex); while we're at it, for each of the
00224     // voronoi neighbours, "D", of newpoint, examine whether newpoint is
00225     // closer to "D" than D's current nearest neighbour -- when this
00226     // occurs, put D into indices_of_updated_neighbours.
00227     // 
00228     // manually put newpoint on indices_of_updated_neighbours
00229     indices_of_updated_neighbours.push_back(index);
00230     _SetAndUpdateNearest(index, indices_of_updated_neighbours);
00231   }
00232 
00233   // for Voronoi neighbours j of any of the removed points for which
00234   // one of those removed points was the nearest neighbour,
00235   // redetermine the nearest neighbour of j and add j onto the vector
00236   // of indices_of_updated_neighbours.
00237   set<int>::iterator it2 = NeighbourUnion.begin();
00238   for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
00239     int j = *it2;
00240     // the if avoids the vertex at infinity, which gets a negative index
00241     if( j != INFINITE_VERTEX ) {
00242       // this is where we check if the nearest neighbour of j was one
00243       // of the removed points
00244       if (indices_removed.count(_supervertex[j].NNindex)) {
00245         if (_verbose) cout << "j " << j << endl;
00246         _SetNearest(j);
00247         indices_of_updated_neighbours.push_back(j);
00248         if (_verbose) cout << "NN of " << j << " : " 
00249                           << _supervertex[j].NNindex
00250                           << ", dist = " << _supervertex[j].NNdistance <<endl;
00251       }
00252     }
00253   }
00254 
00255 }
00256 
00257 
00258 //----------------------------------------------------------------------
00261 void DnnPlane::_SetNearest (const int & j) {
00262   Vertex_handle current = _supervertex[j].vertex;
00263   Vertex_circulator vc = _TR.incident_vertices(current);
00264   Vertex_circulator done = vc;
00265   double dist;
00266   double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
00267   Vertex_handle nearest = _TR.infinite_vertex();
00268 
00269   // when there is only one finite point left in the triangulation, 
00270   // there are no triangles. Presumably this is why voronoi returns
00271   // NULL for the incident vertex circulator. Check if this is
00272   // happening before circulating over it... (Otherwise it crashes
00273   // when looking for neighbours of last point)
00274   if (vc != NULL) do { 
00275     if ( vc->info().val() != INFINITE_VERTEX) {
00276       // find distance between j and its Voronoi neighbour (vc)
00277       if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
00278       dist = _euclid_distance(current->point(), vc->point());
00279       // check if j is closer to vc than vc's currently registered
00280       // nearest neighbour (and update things if it is)
00281       if (dist < mindist) {
00282         mindist = dist; nearest = vc; 
00283         if (_verbose) cout << "nearer ";
00284       } 
00285       if (_verbose) cout << vc->point() << "; "<< dist << endl;
00286     }
00287   } while (++vc != done); // move on to next Voronoi neighbour
00288   // set j's supervertex info about nearest neighbour
00289   _supervertex[j].NNindex = nearest->info().val();
00290   _supervertex[j].NNdistance = mindist;
00291 }
00292 
00293 //----------------------------------------------------------------------
00311 void DnnPlane::_SetAndUpdateNearest(
00312                           const int & j, 
00313                           vector<int> & indices_of_updated_neighbours) {
00314 
00315   Vertex_handle current = _supervertex[j].vertex;
00316   Vertex_circulator vc = _TR.incident_vertices(current);
00317   Vertex_circulator done = vc;
00318   double dist;
00319   double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
00320   Vertex_handle nearest = _TR.infinite_vertex();
00321 
00322   // when there is only one finite point left in the triangulation, 
00323   // there are no triangles. Presumably this is why voronoi returns
00324   // NULL for the incident vertex circulator. Check if this is
00325   // happening before circulating over it... (Otherwise it crashes
00326   // when looking for neighbours of last point)
00327   if (vc != NULL) do { 
00328     if (vc->info().val() != INFINITE_VERTEX) {
00329       if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
00330       // find distance between j and its Voronoi neighbour (vc)
00331       dist = _euclid_distance(current->point(), vc->point());
00332       // update the mindist if we are closer than anything found so far
00333       if (dist < mindist) {
00334         mindist = dist; nearest = vc; 
00335         if (_verbose) cout << "nearer ";
00336       } 
00337       // find index corresponding to vc for easy manipulation
00338       int vcindx = vc->info().val();
00339       if (_verbose) cout << vc->point() << "; "<< dist << endl;
00340       // check if j is closer to vc than vc's currently registered
00341       // nearest neighbour (and update things if it is)
00342       if (dist < _supervertex[vcindx].NNdistance) {
00343         if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
00344         _supervertex[vcindx].NNdistance = dist;
00345         _supervertex[vcindx].NNindex = j;
00346         indices_of_updated_neighbours.push_back(vcindx);
00347       }
00348     }
00349   } while (++vc != done); // move on to next Voronoi neighbour
00350   // set j's supervertex info about nearest neighbour
00351   _supervertex[j].NNindex = nearest->info().val();
00352   _supervertex[j].NNdistance = mindist;
00353 }
00354 
00355 FASTJET_END_NAMESPACE
00356 
00357 #endif //  DROP_CGAL
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines