Fork me on GitHub

Changeset 49234af in git for external/fastjet/DnnPlane.cc


Ignore:
Timestamp:
Dec 9, 2014, 1:27:13 PM (10 years ago)
Author:
Michele <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
37deb3b, 9e991f8
Parents:
f6b6ee7 (diff), e7e90df (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'TestFastJet310b1'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/fastjet/DnnPlane.cc

    rf6b6ee7 r49234af  
    1 //STARTHEADER
    2 // $Id$
     1//FJSTARTHEADER
     2// $Id: DnnPlane.cc 3442 2014-07-24 07:20:49Z salam $
    33//
    4 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
     4// Copyright (c) 2005-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
    55//
    66//----------------------------------------------------------------------
    … …  
    1313//
    1414//  The algorithms that underlie FastJet have required considerable
    15 //  development and are described in hep-ph/0512210. If you use
     15//  development. They are described in the original FastJet paper,
     16//  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
    1617//  FastJet as part of work towards a scientific publication, please
    17 //  include a citation to the FastJet paper.
     18//  quote the version you use and include a citation to the manual and
     19//  optionally also to hep-ph/0512210.
    1820//
    1921//  FastJet is distributed in the hope that it will be useful,
    … …  
    2527//  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
    2628//----------------------------------------------------------------------
    27 //ENDHEADER
     29//FJENDHEADER
    2830
    2931
    … …  
    3335#include<list>
    3436#include "fastjet/internal/DnnPlane.hh"
     37
    3538using namespace std;
    3639
    3740FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
     41
     42const double DnnPlane::DISTANCE_FOR_CGAL_CHECKS=1.0e-12; 
    3843
    3944
    … …  
    5358       _TR.insert(Point(input_points[i].first, input_points[i].second));
    5459
    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;
     60    // check if we are dealing with coincident vertices
     61    int coinciding_index = _CheckIfVertexPresent(sv.vertex, i);
     62    if (coinciding_index == i){
     63      // we need to associate an index to each vertex -- thus when we get
     64      // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
     65      // able to figure out which particle it corresponded to.
     66      sv.vertex->info() = sv.coincidence = i;
     67    } else {
     68      //cout << "  coincident with " << coinciding_index << endl;
     69      // the new vertex points to the already existing one and we
     70      // record the coincidence
     71      //
     72      // Note that we must not only set the coincidence of the
     73      // currently-added particle, the one it coincides with also
     74      // needs be updated (taking into account that it might already
     75      // coincide with another one)
     76      //
     77      // An example may help. Say coinciding_index = i1 and we're adding i2==i.
     78      // Then _sv[i2].coincidence = i1; _sv[i1].coincidence = i2. In both
     79      // cases sv.vertex->info() == i1;
     80      //
     81      // Later on we add i3; we find out that its coinciding index is i1;
     82      // so we set _sv[i3].coincidence = i2 and sv[i1].coincidence = i3.
     83      //
     84      // This gives us the structure
     85      //  _supervertex[i1].coincidence == in
     86      //  _supervertex[i2].coincidence == i1
     87      //  ...
     88      //  _supervertex[in].coincidence == in-1
     89      //
     90      sv.coincidence = _supervertex[coinciding_index].coincidence; // handles cases with previous coincidences
     91      _supervertex[coinciding_index].coincidence = i;
     92    }
     93     
    6394    _supervertex.push_back(sv);   
    6495  }
    … …  
    76107/// Crashes if the given vertex handle already exists. Otherwise
    77108/// it does the bookkeeping for future such tests
    78 void DnnPlane::_CrashIfVertexPresent(
    79         const Vertex_handle & vertex, const int & its_index) {
    80   if (!_crash_on_coincidence) return;
    81 
     109int DnnPlane::_CheckIfVertexPresent(
     110        const Vertex_handle & vertex, const int its_index) {
    82111  // vertices that do not have the same geometric position as any
    83112  // other vertex so far added have info().val() == NEW_VERTEX -- this
    … …  
    90119  // DNN:DNN) to be equal to a vertex "index".
    91120  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   }
     121    if (_crash_on_coincidence){
     122      ostringstream err;
     123      err << "Error: DnnPlane::_CheckIfVertexPresent"
     124          << "Point "<<its_index<<" coincides with point "
     125          <<vertex->info().val() << endl;
     126      throw DnnError(err.str());
     127    }
     128    return vertex->info().val();
     129  }
     130
     131  return its_index;
    98132}
    99133
    … …  
    116150                          vector<int> & indices_of_updated_neighbours) {
    117151
     152  if (_verbose) cout << "Starting  DnnPlane::RemoveAndAddPoints" << endl;
    118153
    119154  // build set of UNION of Voronoi neighbours of a pair of nearest
    … …  
    124159  set<int> indices_removed;
    125160
    126   // for each of the indices to be removed add the voronoi neighbourhood to
    127   // the NeighbourUnion set.
     161  // for each of the indices to be removed add the voronoi
     162  // neighbourhood to the NeighbourUnion set as well as the coinciding
     163  // points that had the current point as coincidence before.
    128164  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
    129165    int index = indices_to_remove[ir];
    130166    indices_removed.insert(index);
    131     if (_verbose) cout << " Starting  RemoveAndAddPoints" << endl;
    132     if (_verbose) cout << " point " << index << endl;                   
     167    if (_verbose) cout << "  scheduling point " << index << " for removal" << endl;
     168
     169    if (_supervertex[index].coincidence != index){
     170      // we have a coincidence
     171      //
     172      // The only one of the coincident points that has to be
     173      // inserted in the neighbourhood list (and thus updated) is the
     174      // one that has 'index' as coincidence.
     175      int new_index = _supervertex[index].coincidence;
     176      while (_supervertex[new_index].coincidence != index)
     177        new_index = _supervertex[new_index].coincidence;
     178      if (_verbose) cout << "  inserted coinciding " << new_index << " to neighbours union" << endl;
     179      NeighbourUnion.insert(new_index);
     180
     181      // if this is the point among the coiciding ones that holds the
     182      // CGAL vertex, then also insert the CGAL neighbours, otherwise
     183      // just skip that step.
     184      if (index != _supervertex[index].vertex->info().val()) continue;
     185    }
     186
    133187    // have a circulators that will go round the Voronoi neighbours of
    134188    // _supervertex[index1].vertex
    135189    Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
    136190    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);
     191    if (vc != NULL){ // a safety check in case there is no Voronoi
     192                     // neighbour (which may happen e.g. if we just
     193                     // have a bunch of coincident points)
     194      do  {
     195        // if a neighbouring vertex is not the infinite vertex, then add it
     196        // to our union of neighbouring vertices.
     197        if (_verbose) cout << "examining " << vc->info().val() << endl;
     198        if (vc->info().val() != INFINITE_VERTEX) {
     199          // NB: from it=1 onwards occasionally it might already have
     200          // been inserted -- but double insertion still leaves only one
     201          // copy in the set, so there's no problem
     202          NeighbourUnion.insert(vc->info().val());
     203          if (_verbose) cout << "  inserted " << vc->info().val() << " to neighbours union" << endl;
     204        }
     205      } while (++vc != done);
     206    }
    149207  }
    150208 
    … …  
    160218  for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
    161219    int index = indices_to_remove[ir];
     220    if (_verbose) cout << "  removing " << index << endl;
    162221
    163222    // NeighbourUnion should not contain the points to be removed
    164223    // (because later we will assume they still exist).
    165224    NeighbourUnion.erase(indices_to_remove[ir]);
    166    
     225
     226    // first deal with  coincidences
     227    if (_supervertex[index].coincidence != index){
     228      int new_index = _supervertex[index].coincidence;
     229
     230      // if this is the point among the coiciding ones that "owns" the
     231      // CGAL vertex we need to re-label the CGAL vertex so that it
     232      // points to the coincident particle and set the current one to
     233      // NULL
     234      //
     235      // This can be done only on the first point as they all share
     236      // the same value
     237      //
     238      // Note that this has to be done before the following step since
     239      // it will alter the coincidence information
     240      if (index == _supervertex[index].vertex->info().val())
     241        _supervertex[new_index].vertex->info() = new_index;
     242
     243      // we need to browse the coincidences until we end the loop, at
     244      // which point we reset the coincidence of the point that has
     245      // the current one as a coincidence
     246      while (_supervertex[new_index].coincidence != index)
     247        new_index = _supervertex[new_index].coincidence;
     248      _supervertex[new_index].coincidence = _supervertex[index].coincidence;
     249
     250      // remove the coincidence on the point being removed and mark it
     251      // as removed
     252      _supervertex[index].coincidence = index;
     253      _supervertex[index].vertex = NULL;
     254
     255      continue;
     256    }
     257
    167258    // points to be removed should also be eliminated from the
    168259    // triangulation and the supervertex structure should be updated
    … …  
    192283  // of the neighbour union happens to be on the wrong side.
    193284  Face_handle face;
    194   if (indices_to_remove.size() > 0) {
     285  //if (indices_to_remove.size() > 0) { // GS: use NeighbourUnion instead
     286                                        //     (safe also in case of coincidences)
     287  if (NeighbourUnion.size() > 0) {
    195288    // face can only be found if there were points to remove in first place
    196289    face = _TR.incident_faces(
    … …  
    204297    int index = _supervertex.size()-1;
    205298    indices_added.push_back(index);
    206 
    207     if (indices_to_remove.size() > 0) {
     299    if (_verbose) cout << "  adding " << index << endl;
     300
     301    //if (indices_to_remove.size() > 0) {
     302    if (NeighbourUnion.size() > 0) {
    208303      // be careful of using face (for location hinting) only when it exists
    209304      _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
    … …  
    213308                                                    points_to_add[ia].second));
    214309    }
    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;
     310
     311    // check if this leads to a coincidence
     312    int coinciding_index = _CheckIfVertexPresent(_supervertex[index].vertex, index);
     313    if (coinciding_index == index){
     314      // we need to associate an index to each vertex -- thus when we get
     315      // a vertex (e.g. as a nearest neighbour) from CGAL, we will be
     316      // able to figure out which particle it corresponded to.
     317      _supervertex[index].vertex->info() = _supervertex[index].coincidence = index;
     318    } else {
     319      if (_verbose) cout << "  coinciding with vertex " << coinciding_index << endl;
     320      // the new vertex points to an already existing one and we
     321      // record the coincidence
     322      //
     323      // we also update the NN of the coinciding particle (to avoid
     324      // having to loop over the list of coinciding neighbours later)
     325      // This is done first as it allows us to check if this is a new
     326      // coincidence or a coincidence added to a particle that was
     327      // previously "alone"
     328      _supervertex[coinciding_index].NNindex = index;
     329      _supervertex[coinciding_index].NNdistance = 0.0;
     330      indices_of_updated_neighbours.push_back(coinciding_index);
     331
     332      // Note that we must not only set the coincidence of the
     333      // currently-added particle, the one it coincides with also
     334      // needs be updated (taking into account that it might already
     335      // coincide with another one)
     336      _supervertex[index].coincidence = _supervertex[coinciding_index].coincidence; // handles cases with previous coincidences
     337      _supervertex[coinciding_index].coincidence = index;
     338
     339    }
    219340   
    220341    // first find nearest neighbour of "newpoint" (shorthand for
    … …  
    227348    indices_of_updated_neighbours.push_back(index);
    228349    _SetAndUpdateNearest(index, indices_of_updated_neighbours);
     350
     351    //cout << "Added: " << setprecision(20) << " ("
     352    //     << points_to_add[ia].first << "," << points_to_add[ia].second
     353    //     << ") with index " << index << endl;
    229354  }
    230355
    … …  
    251376  }
    252377
     378  if (_verbose) cout << "Leaving  DnnPlane::RemoveAndAddPoints" << endl;
    253379}
    254 
    255380
    256381//----------------------------------------------------------------------
    257382/// Determines the index and distance of the nearest neighbour to
    258383/// point j and puts the information into the _supervertex entry for j.
    259 void DnnPlane::_SetNearest (const int & j) {
     384void DnnPlane::_SetNearest (const int j) {
     385  // first deal with the cases where we have a coincidence
     386  if (_supervertex[j].coincidence != j){
     387    _supervertex[j].NNindex = _supervertex[j].coincidence;
     388    _supervertex[j].NNdistance = 0.0;
     389    return;
     390  }
     391
     392  // The code below entirely uses CGAL distance comparisons to compute
     393  // the nearest neighbour. It has the mais drawback to induice a
     394  // 10-20% time penalty so we switched to our own comparison (which
     395  // only turns to CGAL for dangerous situations)
     396  //
     397  //  Vertex_handle current = _supervertex[j].vertex;
     398  //  Vertex_circulator vc = _TR.incident_vertices(current);
     399  //  Vertex_circulator done = vc;
     400  //  Vertex_handle nearest = _TR.infinite_vertex();
     401  //  double mindist = HUGE_DOUBLE;
     402  //
     403  //   // when there is only one finite point left in the triangulation,
     404  //   // there are no triangles. Presumably this is why voronoi returns
     405  //   // NULL for the incident vertex circulator. Check if this is
     406  //   // happening before circulating over it... (Otherwise it crashes
     407  //   // when looking for neighbours of last point)
     408  //   if (vc != NULL){
     409  //     // initialise the nearest vertex handle to the first incident
     410  //     // vertex that is not INFINITE_VERTEX
     411  //     while (vc->info().val() == INFINITE_VERTEX){
     412  //       vc++;
     413  //       if (vc==done) break; // if vc==done, then INFINITE_VERTEX is the
     414  //                       // only element in the neighbourhood
     415  //     }
     416  //
     417  //     // if there is just the infinite vertex, we have vc->info().val()
     418  //     // == INFINITE_VERTEX and nothing has to be done
     419  //     // otherwise, use the current vc as an initialisation
     420  //     if (vc->info().val() != INFINITE_VERTEX){
     421  //       nearest = vc; // initialisation to the first non-infinite vertex
     422  //
     423  //       // and loop over the following ones
     424  //       while (++vc != done){
     425  //    // we should not compare with the infinite vertex
     426  //    if (vc->info().val() == INFINITE_VERTEX) continue;
     427  //
     428  //    if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
     429  //    // use CGAL's distance comparison to check if 'vc' is closer to
     430  //    // 'current' than the nearest so far (we include the == case for
     431  //    // safety though it should not matter in this precise case)
     432  //    if (CGAL::compare_distance_to_point(current->point(), vc->point(), nearest->point())!=CGAL::LARGER){
     433  //      nearest = vc;
     434  //      if (_verbose) cout << "nearer";
     435  //    }
     436  //       }
     437  //
     438  //       // now compute the distance
     439  //       //
     440  //       // Note that since we're always using CGAL to compare distances
     441  //       // (and never the distance computed using _euclid_distance) we
     442  //       // should not worry about rounding errors in mindist
     443  //       mindist = _euclid_distance(current->point(), nearest->point());
     444  //     }
     445  //   }
     446  //
     447  //  // set j's supervertex info about nearest neighbour
     448  //  _supervertex[j].NNindex = nearest->info().val();
     449  //  _supervertex[j].NNdistance = mindist;
     450
    260451  Vertex_handle current = _supervertex[j].vertex;
    261452  Vertex_circulator vc = _TR.incident_vertices(current);
    … …  
    264455  double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double?
    265456  Vertex_handle nearest = _TR.infinite_vertex();
    266 
     457 
    267458  // when there is only one finite point left in the triangulation,
    268459  // there are no triangles. Presumably this is why voronoi returns
    … …  
    274465      // find distance between j and its Voronoi neighbour (vc)
    275466      if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
    276       dist = _euclid_distance(current->point(), vc->point());
     467
    277468      // check if j is closer to vc than vc's currently registered
    278469      // nearest neighbour (and update things if it is)
    279       if (dist < mindist) {
    280         mindist = dist; nearest = vc;
    281         if (_verbose) cout << "nearer ";
     470      if (_is_closer_to(current->point(), vc->point(), nearest, dist, mindist)){
     471        nearest = vc;
     472        if (_verbose) cout << "nearer ";
    282473      }
    283474      if (_verbose) cout << vc->point() << "; "<< dist << endl;
    284475    }
    285476  } while (++vc != done); // move on to next Voronoi neighbour
     477 
    286478  // set j's supervertex info about nearest neighbour
    287479  _supervertex[j].NNindex = nearest->info().val();
    … …  
    291483//----------------------------------------------------------------------
    292484/// Determines and stores the nearest neighbour of j, and where
    293 /// necessary update the nearest-neighbour info of Voronoi neighbours
     485/// necessary updates the nearest-neighbour info of Voronoi neighbours
    294486/// of j;
    295487///
    … …  
    304496/// NB: note that we have _SetAndUpdateNearest as a completely
    305497///     separate routine from _SetNearest because we want to
    306 ///     use one single ciruclation over voronoi neighbours to find the
     498///     use one single circulation over voronoi neighbours to find the
    307499///     nearest neighbour and to update the voronoi neighbours if need
    308500///     be.
    309501void DnnPlane::_SetAndUpdateNearest(
    310                           const int & j,
     502                          const int j,
    311503                          vector<int> & indices_of_updated_neighbours) {
     504  //cout << "SetAndUpdateNearest for point " << j << endl;
     505  // first deal with coincidences
     506  if (_supervertex[j].coincidence != j){
     507    _supervertex[j].NNindex = _supervertex[j].coincidence;
     508    _supervertex[j].NNdistance = 0.0;
     509    //cout << "  set to coinciding point " << _supervertex[j].coincidence << endl;
     510    return;
     511  }
    312512
    313513  Vertex_handle current = _supervertex[j].vertex;
    … …  
    326526    if (vc->info().val() != INFINITE_VERTEX) {
    327527      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());
     528
    330529      // 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 ";
     530      if (_is_closer_to(current->point(), vc->point(), nearest, dist, mindist)){
     531        nearest = vc;
     532        if (_verbose) cout << "nearer ";
    334533      }
     534
    335535      // find index corresponding to vc for easy manipulation
    336536      int vcindx = vc->info().val();
    337537      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) {
     538
     539      if (_is_closer_to_with_hint(vc->point(), current->point(),
     540                                  _supervertex[_supervertex[vcindx].NNindex].vertex,
     541                                  dist, _supervertex[vcindx].NNdistance)){
    341542        if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
    342         _supervertex[vcindx].NNdistance = dist;
    343543        _supervertex[vcindx].NNindex = j;
    344544        indices_of_updated_neighbours.push_back(vcindx);
    345545      }
     546
     547      // original code without the use of CGAL distance in potentially
     548      // dangerous cases
     549      //
     550      // // check if j is closer to vc than vc's currently registered
     551      // // nearest neighbour (and update things if it is)
     552      // //
     553      // // GS: originally, the distance test below was a strict <. It
     554      // //     has to be <= because if the two distances are ==, it is
     555      // //     possible that the old NN is no longer connected to vc in
     556      // //     the triangulation, and we are sure that the newly
     557      // //     inserted point (j) is (since we loop over j's
     558      // //     neighbouring points in the triangulation).
     559      // if (dist <= _supervertex[vcindx].NNdistance) {
     560      //        if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
     561      //        _supervertex[vcindx].NNdistance = dist;
     562      //        _supervertex[vcindx].NNindex = j;
     563      //        indices_of_updated_neighbours.push_back(vcindx);
     564      // }
    346565    }
    347566  } while (++vc != done); // move on to next Voronoi neighbour
    348567  // set j's supervertex info about nearest neighbour
     568  //cout << "  set to point " << nearest->info().val() << endl;
    349569  _supervertex[j].NNindex = nearest->info().val();
    350570  _supervertex[j].NNdistance = mindist;
Note: See TracChangeset for help on using the changeset viewer.