FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MinHeap.cc
1 //STARTHEADER
2 // $Id: MinHeap.cc 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 #include "fastjet/internal/MinHeap.hh"
30 #include<iostream>
31 #include<cmath>
32 #include<limits>
33 
34 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
35 
36 using namespace std;
37 
38 //----------------------------------------------------------------------
39 /// construct the MinHeap; structure will be as follows:
40 /// . _heap[0].minloc points to globally smallest entry
41 /// _heap[1].minloc points to smallest entry in one half of heap
42 /// _heap[2].minloc points to smallest entry in other half of heap
43 ///
44 /// . for _heap[i], the "parent" is to be found at (i-1)/2
45 void MinHeap::_initialise(const std::vector<double> & values){
46 
47  // fill the high-range of the heap with the largest possible value
48  // (minloc of each entry is itself)
49  for (unsigned i = values.size(); i < _heap.size(); i++) {
50  _heap[i].value = std::numeric_limits<double>::max();
51  _heap[i].minloc = &(_heap[i]);
52  }
53 
54  // fill the rest of the heap with the actual values
55  // (minloc of each entry is itself)
56  for (unsigned i = 0; i < values.size(); i++) {
57  _heap[i].value = values[i];
58  _heap[i].minloc = &(_heap[i]);
59  }
60 
61  // now adjust the minlocs so that everything is OK...
62  for (unsigned i = _heap.size()-1; i > 0; i--) {
63  ValueLoc * parent = &(_heap[(i-1)/2]);
64  ValueLoc * here = &(_heap[i]);
65  if (here->minloc->value < parent->minloc->value) {
66  parent->minloc = here->minloc;
67  }
68  }
69  //cout << minloc() << " "<<sqrt(minval())<<endl;
70  //cout << sqrt(_heap[47].value)<<endl;
71  //cout << sqrt(_heap[48].value)<<endl;
72  //cout << sqrt(_heap[25].value)<<endl;
73 }
74 
75 
76 //----------------------------------------------------------------------
77 void MinHeap::update(unsigned int loc, double new_value) {
78 
79 
80  assert(loc < _heap.size());
81  ValueLoc * start = &(_heap[loc]);
82 
83  // if the minloc is somewhere below us and our value is no smaller
84  // than the previous value, we can't possibly change the minloc
85  if (start->minloc != start && !(new_value < start->minloc->value)) {
86  start->value = new_value;
87  //std::cout << " had easy exit\n";
88  return;
89  }
90 
91  // update the value and put a temporary location
92  start->value = new_value;
93  start->minloc = start;
94  // warn that a change has been made at this place
95  bool change_made = true;
96  // to make sure we don't go off edge...
97  ValueLoc * heap_end = (&(_heap[0])) + _heap.size();
98 
99  // now work our way up the heap
100  while(change_made) {
101  ValueLoc * here = &(_heap[loc]);
102  change_made = false;
103 
104  // if we were pointing to start, then we must re-initialise things
105  if (here->minloc == start) {
106  here->minloc = here; change_made = true;
107  }
108 
109  // now compare current location to children (at 2*loc+1, 2*loc+2)
110  ValueLoc * child = &(_heap[2*loc+1]);
111  if (child < heap_end && child->minloc->value < here->minloc->value ) {
112  here->minloc = child->minloc;
113  change_made = true;}
114  child++;
115  if (child < heap_end && child->minloc->value < here->minloc->value ) {
116  here->minloc = child->minloc;
117  change_made = true;}
118 
119  // then we move up (loc ->(loc-1)/2) if there's anywhere to go
120  if (loc == 0) {break;}
121  loc = (loc-1)/2;
122  }
123 
124 }
125 
126 FASTJET_END_NAMESPACE
127