fastjet 2.4.5
|
00001 //STARTHEADER 00002 // $Id: SearchTree.hh 3157 2013-07-11 16:02:51Z soyez $ 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 __FASTJET_SEARCHTREE_HH__ 00033 #define __FASTJET_SEARCHTREE_HH__ 00034 00035 #include<vector> 00036 #include<cassert> 00037 #include<cstddef> 00038 #include "fastjet/internal/base.hh" 00039 00040 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00041 00042 00043 //====================================================================== 00048 template<class T> class SearchTree { 00049 public: 00050 00051 class Node; 00052 class circulator; 00053 class const_circulator; 00054 00056 SearchTree(const std::vector<T> & init); 00057 00060 SearchTree(const std::vector<T> & init, unsigned int max_size); 00061 00063 void remove(unsigned node_index); 00064 void remove(typename SearchTree::Node * node); 00065 void remove(typename SearchTree::circulator & circ); 00066 00069 //Node * insert(const T & value); 00070 circulator insert(const T & value); 00071 00072 const Node & operator[](int i) const {return _nodes[i];}; 00073 00075 unsigned int size() const {return _nodes.size() - _available_nodes.size();} 00076 00078 void verify_structure(); 00079 void verify_structure_linear() const; 00080 void verify_structure_recursive(const Node * , const Node * , const Node * ) const; 00081 00083 void print_elements(); 00084 00085 // tracking the depth may have some speed overhead -- so leave it 00086 // out for the time being... 00087 #ifdef TRACK_DEPTH 00088 00089 inline unsigned int max_depth() const {return _max_depth;}; 00090 #else 00091 inline unsigned int max_depth() const {return 0;}; 00092 #endif 00093 00094 int loc(const Node * node) const ; 00095 00097 Node * _find_predecessor(const Node *); 00099 Node * _find_successor(const Node *); 00100 00101 const Node & operator[](unsigned int i) const {return _nodes[i];}; 00102 00105 const_circulator somewhere() const; 00106 circulator somewhere(); 00107 00108 private: 00109 00110 void _initialize(const std::vector<T> & init); 00111 00112 std::vector<Node> _nodes; 00113 std::vector<Node *> _available_nodes; 00114 Node * _top_node; 00115 unsigned int _n_removes; 00116 00117 00122 void _do_initial_connections(unsigned int this_one, unsigned int scale, 00123 unsigned int left_edge, unsigned int right_edge, 00124 unsigned int depth); 00125 00126 00127 #ifdef TRACK_DEPTH 00128 unsigned int _max_depth; 00129 #endif 00130 00131 }; 00132 00133 00134 //====================================================================== 00135 template<class T> class SearchTree<T>::Node{ 00136 public: 00137 Node() {}; 00138 00139 00141 bool treelinks_null() const { 00142 return ((parent==0) && (left==0) && (right==0));}; 00143 00145 inline void nullify_treelinks() { 00146 parent = NULL; 00147 left = NULL; 00148 right = NULL; 00149 }; 00150 00153 void reset_parents_link_to_me(Node * XX); 00154 00155 T value; 00156 Node * left; 00157 Node * right; 00158 Node * parent; 00159 Node * successor; 00160 Node * predecessor; 00161 }; 00162 00163 //---------------------------------------------------------------------- 00164 template<class T> void SearchTree<T>::Node::reset_parents_link_to_me(typename SearchTree<T>::Node * XX) { 00165 if (parent == NULL) {return;} 00166 if (parent->right == this) {parent->right = XX;} 00167 else {parent->left = XX;} 00168 } 00169 00170 00171 00172 //====================================================================== 00173 template<class T> class SearchTree<T>::circulator{ 00174 public: 00175 00176 // so that it can access out _node object; 00177 // note: "class U" needed for clang (v1.1 branches/release_27) compilation 00178 template<class U> friend class SearchTree<U>::const_circulator; 00179 friend class SearchTree<T>; 00180 00181 circulator() : _node(NULL) {} 00182 00183 circulator(Node * node) : _node(node) {} 00184 00185 const T * operator->() const {return &(_node->value);} 00186 T * operator->() {return &(_node->value);} 00187 const T & operator*() const {return _node->value;} 00188 T & operator*() {return _node->value;} 00189 00191 circulator & operator++() { 00192 _node = _node->successor; 00193 return *this;} 00194 00197 circulator operator++(int) { 00198 circulator tmp = *this; 00199 _node = _node->successor; 00200 return tmp;} 00201 00203 circulator & operator--() { 00204 _node = _node->predecessor; 00205 return *this;} 00206 00209 circulator operator--(int) { 00210 circulator tmp = *this; 00211 _node = _node->predecessor; 00212 return tmp;} 00213 00215 circulator next() const { 00216 return circulator(_node->successor);} 00217 00219 circulator previous() const { 00220 return circulator(_node->predecessor);} 00221 00222 bool operator!=(const circulator & other) const {return other._node != _node;} 00223 bool operator==(const circulator & other) const {return other._node == _node;} 00224 00225 private: 00226 Node * _node; 00227 }; 00228 00229 00230 //====================================================================== 00231 template<class T> class SearchTree<T>::const_circulator{ 00232 public: 00233 00234 const_circulator() : _node(NULL) {} 00235 00236 const_circulator(const Node * node) : _node(node) {} 00237 const_circulator(const circulator & circ) :_node(circ._node) {} 00238 00239 const T * operator->() {return &(_node->value);} 00240 const T & operator*() const {return _node->value;} 00241 00243 const_circulator & operator++() { 00244 _node = _node->successor; 00245 return *this;} 00246 00249 const_circulator operator++(int) { 00250 const_circulator tmp = *this; 00251 _node = _node->successor; 00252 return tmp;} 00253 00254 00256 const_circulator & operator--() { 00257 _node = _node->predecessor; 00258 return *this;} 00259 00262 const_circulator operator--(int) { 00263 const_circulator tmp = *this; 00264 _node = _node->predecessor; 00265 return tmp;} 00266 00268 const_circulator next() const { 00269 return const_circulator(_node->successor);} 00270 00272 const_circulator previous() const { 00273 return const_circulator(_node->predecessor);} 00274 00275 00276 00277 bool operator!=(const const_circulator & other) const {return other._node != _node;} 00278 bool operator==(const const_circulator & other) const {return other._node == _node;} 00279 00280 private: 00281 const Node * _node; 00282 }; 00283 00284 00285 00286 00287 //---------------------------------------------------------------------- 00290 template<class T> SearchTree<T>::SearchTree(const std::vector<T> & init, 00291 unsigned int max_size) : 00292 _nodes(max_size) { 00293 00294 _available_nodes.reserve(max_size); 00295 _available_nodes.resize(max_size - init.size()); 00296 for (unsigned int i = init.size(); i < max_size; i++) { 00297 _available_nodes[i-init.size()] = &(_nodes[i]); 00298 } 00299 00300 _initialize(init); 00301 } 00302 00303 //---------------------------------------------------------------------- 00305 template<class T> SearchTree<T>::SearchTree(const std::vector<T> & init) : 00306 _nodes(init.size()), _available_nodes(0) { 00307 00308 // reserve space for the list of available nodes 00309 _available_nodes.reserve(init.size()); 00310 _initialize(init); 00311 } 00312 00313 //---------------------------------------------------------------------- 00315 template<class T> void SearchTree<T>::_initialize(const std::vector<T> & init) { 00316 00317 _n_removes = 0; 00318 unsigned n = init.size(); 00319 assert(n>=1); 00320 00321 // reserve space for the list of available nodes 00322 //_available_nodes.reserve(); 00323 00324 #ifdef TRACK_DEPTH 00325 _max_depth = 0; 00326 #endif 00327 00328 00329 // validate the input 00330 for (unsigned int i = 1; i<n; i++) { 00331 assert(!(init[i] < init[i-1])); 00332 } 00333 00334 // now initialise the vector; link neighbours in the sequence 00335 for(unsigned int i = 0; i < n; i++) { 00336 _nodes[i].value = init[i]; 00337 _nodes[i].predecessor = (& (_nodes[i])) - 1; 00338 _nodes[i].successor = (& (_nodes[i])) + 1; 00339 _nodes[i].nullify_treelinks(); 00340 } 00341 // make a loop structure so that we can circulate... 00342 _nodes[0].predecessor = (& (_nodes[n-1])); 00343 _nodes[n-1].successor = (& (_nodes[0])); 00344 00345 // now label the rest of the nodes 00346 unsigned int scale = (n+1)/2; 00347 unsigned int top = std::min(n-1,scale); 00348 _nodes[top].parent = NULL; 00349 _top_node = &(_nodes[top]); 00350 _do_initial_connections(top, scale, 0, n, 0); 00351 00352 // make sure things are sensible... 00353 //verify_structure(); 00354 } 00355 00356 00357 00358 //---------------------------------------------------------------------- 00359 template<class T> inline int SearchTree<T>::loc(const Node * node) const {return node == NULL? 00360 -999 : node - &(_nodes[0]);} 00361 00362 00363 //---------------------------------------------------------------------- 00366 template<class T> void SearchTree<T>::_do_initial_connections( 00367 unsigned int this_one, 00368 unsigned int scale, 00369 unsigned int left_edge, 00370 unsigned int right_edge, 00371 unsigned int depth 00372 ) { 00373 00374 #ifdef TRACK_DEPTH 00375 // keep track of tree depth for checking things stay reasonable... 00376 _max_depth = max(depth, _max_depth); 00377 #endif 00378 00379 //std::cout << this_one << " "<< scale<< std::endl; 00380 unsigned int ref_new_scale = (scale+1)/2; 00381 00382 // work through children to our left 00383 unsigned new_scale = ref_new_scale; 00384 bool did_child = false; 00385 while(true) { 00386 int left = this_one - new_scale; // be careful here to use signed int... 00387 // if there is something unitialised to our left, link to it 00388 if (left >= static_cast<int>(left_edge) 00389 && _nodes[left].treelinks_null() ) { 00390 _nodes[left].parent = &(_nodes[this_one]); 00391 _nodes[this_one].left = &(_nodes[left]); 00392 // create connections between left_edge and this_one 00393 _do_initial_connections(left, new_scale, left_edge, this_one, depth+1); 00394 did_child = true; 00395 break; 00396 } 00397 // reduce the scale so as to try again 00398 unsigned int old_new_scale = new_scale; 00399 new_scale = (old_new_scale + 1)/2; 00400 // unless we've reached end of tree 00401 if (new_scale == old_new_scale) break; 00402 } 00403 if (!did_child) {_nodes[this_one].left = NULL;} 00404 00405 00406 // work through children to our right 00407 new_scale = ref_new_scale; 00408 did_child = false; 00409 while(true) { 00410 unsigned int right = this_one + new_scale; 00411 if (right < right_edge && _nodes[right].treelinks_null()) { 00412 _nodes[right].parent = &(_nodes[this_one]); 00413 _nodes[this_one].right = &(_nodes[right]); 00414 // create connections between this_one+1 and right_edge 00415 _do_initial_connections(right, new_scale, this_one+1,right_edge,depth+1); 00416 did_child = true; 00417 break; 00418 } 00419 // reduce the scale so as to try again 00420 unsigned int old_new_scale = new_scale; 00421 new_scale = (old_new_scale + 1)/2; 00422 // unless we've reached end of tree 00423 if (new_scale == old_new_scale) break; 00424 } 00425 if (!did_child) {_nodes[this_one].right = NULL;} 00426 00427 } 00428 00429 00430 00431 //---------------------------------------------------------------------- 00432 template<class T> void SearchTree<T>::remove(unsigned int node_index) { 00433 remove(&(_nodes[node_index])); 00434 } 00435 00436 //---------------------------------------------------------------------- 00437 template<class T> void SearchTree<T>::remove(circulator & circ) { 00438 remove(circ._node); 00439 } 00440 00441 //---------------------------------------------------------------------- 00442 // Useful reference for this: 00443 // http://en.wikipedia.org/wiki/Binary_search_tree#Deletion 00444 template<class T> void SearchTree<T>::remove(typename SearchTree<T>::Node * node) { 00445 00446 // we don't remove things from the tree if we've reached the last 00447 // elements... (is this wise?) 00448 assert(size() > 1); // switch this to throw...? 00449 assert(!node->treelinks_null()); 00450 00451 // deal with relinking predecessor and successor 00452 node->predecessor->successor = node->successor; 00453 node->successor->predecessor = node->predecessor; 00454 00455 if (node->left == NULL && node->right == NULL) { 00456 // node has no children, so remove it by nullifying the pointer 00457 // from the parent 00458 node->reset_parents_link_to_me(NULL); 00459 00460 } else if (node->left != NULL && node->right == NULL){ 00461 // make parent point to my child 00462 node->reset_parents_link_to_me(node->left); 00463 // and child to parent 00464 node->left->parent = node->parent; 00465 // sort out the top node... 00466 if (_top_node == node) {_top_node = node->left;} 00467 00468 } else if (node->left == NULL && node->right != NULL){ 00469 // make parent point to my child 00470 node->reset_parents_link_to_me(node->right); 00471 // and child to parent 00472 node->right->parent = node->parent; 00473 // sort out the top node... 00474 if (_top_node == node) {_top_node = node->right;} 00475 00476 } else { 00477 // we have two children; we will put a replacement in our place 00478 Node * replacement; 00479 //SearchTree<T>::Node * replacements_child; 00480 // chose predecessor or successor (one, then other, then first, etc...) 00481 bool use_predecessor = (_n_removes % 2 == 1); 00482 if (use_predecessor) { 00483 // Option 1: put predecessor in our place, and have its parent 00484 // point to its left child (as a predecessor it has no right child) 00485 replacement = node->predecessor; 00486 assert(replacement->right == NULL); // guaranteed if it's our predecessor 00487 // we have to be careful of replacing certain links when the 00488 // replacement is this node's child 00489 if (replacement != node->left) { 00490 if (replacement->left != NULL) { 00491 replacement->left->parent = replacement->parent;} 00492 replacement->reset_parents_link_to_me(replacement->left); 00493 replacement->left = node->left; 00494 } 00495 replacement->parent = node->parent; 00496 replacement->right = node->right; 00497 } else { 00498 // Option 2: put successor in our place, and have its parent 00499 // point to its right child (as a successor it has no left child) 00500 replacement = node->successor; 00501 assert(replacement->left == NULL); // guaranteed if it's our successor 00502 if (replacement != node->right) { 00503 if (replacement->right != NULL) { 00504 replacement->right->parent = replacement->parent;} 00505 replacement->reset_parents_link_to_me(replacement->right); 00506 replacement->right = node->right; 00507 } 00508 replacement->parent = node->parent; 00509 replacement->left = node->left; 00510 } 00511 node->reset_parents_link_to_me(replacement); 00512 00513 // make sure node's original children now point to the replacement 00514 if (node->left != replacement) {node->left->parent = replacement;} 00515 if (node->right != replacement) {node->right->parent = replacement;} 00516 00517 // sort out the top node... 00518 if (_top_node == node) {_top_node = replacement;} 00519 } 00520 00521 // make sure we leave something nice and clean... 00522 node->nullify_treelinks(); 00523 node->predecessor = NULL; 00524 node->successor = NULL; 00525 00526 // for bookkeeping (and choosing whether to use pred. or succ.) 00527 _n_removes++; 00528 // for when we next need access to a free node... 00529 _available_nodes.push_back(node); 00530 } 00531 00532 00533 //---------------------------------------------------------------------- 00534 //template<class T> typename SearchTree<T>::Node * SearchTree<T>::insert(const T & value) { 00535 00536 //---------------------------------------------------------------------- 00537 template<class T> typename SearchTree<T>::circulator SearchTree<T>::insert(const T & value) { 00538 // make sure we don't exceed allowed number of nodes... 00539 assert(_available_nodes.size() > 0); 00540 00541 Node * node = _available_nodes.back(); 00542 _available_nodes.pop_back(); 00543 node->value = value; 00544 00545 Node * location = _top_node; 00546 Node * old_location = NULL; 00547 bool on_left = true; // (init not needed -- but soothes g++4) 00548 // work through tree until we reach its end 00549 #ifdef TRACK_DEPTH 00550 unsigned int depth = 0; 00551 #endif 00552 while(location != NULL) { 00553 #ifdef TRACK_DEPTH 00554 depth++; 00555 #endif 00556 old_location = location; 00557 on_left = value < location->value; 00558 if (on_left) {location = location->left;} 00559 else {location = location->right;} 00560 } 00561 #ifdef TRACK_DEPTH 00562 _max_depth = max(depth, _max_depth); 00563 #endif 00564 // now create tree links 00565 node->parent = old_location; 00566 if (on_left) {node->parent->left = node;} 00567 else {node->parent->right = node;} 00568 node->left = NULL; 00569 node->right = NULL; 00570 // and create predecessor / successor links 00571 node->predecessor = _find_predecessor(node); 00572 if (node->predecessor != NULL) { 00573 // it exists, so make use of its info (will include a cyclic case, 00574 // when successor is round the bend) 00575 node->successor = node->predecessor->successor; 00576 node->predecessor->successor = node; 00577 node->successor->predecessor = node; 00578 } else { 00579 // deal with case when we are left-most edge of tree (then successor 00580 // will exist...) 00581 node->successor = _find_successor(node); 00582 assert(node->successor != NULL); // can only happen if we're sole element 00583 // (but not allowed, since tree size>=1) 00584 node->predecessor = node->successor->predecessor; 00585 node->successor->predecessor = node; 00586 node->predecessor->successor = node; 00587 } 00588 00589 return circulator(node); 00590 } 00591 00592 00593 //---------------------------------------------------------------------- 00594 template<class T> void SearchTree<T>::verify_structure() { 00595 00596 // do a check running through all elements 00597 verify_structure_linear(); 00598 00599 // do a recursive check down tree from top 00600 00601 // first establish the extremities 00602 const Node * left_limit = _top_node; 00603 while (left_limit->left != NULL) {left_limit = left_limit->left;} 00604 const Node * right_limit = _top_node; 00605 while (right_limit->right != NULL) {right_limit = right_limit->right;} 00606 00607 // then actually do recursion 00608 verify_structure_recursive(_top_node, left_limit, right_limit); 00609 } 00610 00611 00612 //---------------------------------------------------------------------- 00613 template<class T> void SearchTree<T>::verify_structure_recursive( 00614 const typename SearchTree<T>::Node * element, 00615 const typename SearchTree<T>::Node * left_limit, 00616 const typename SearchTree<T>::Node * right_limit) const { 00617 00618 assert(!(element->value < left_limit->value)); 00619 assert(!(right_limit->value < element->value)); 00620 00621 const Node * left = element->left; 00622 if (left != NULL) { 00623 assert(!(element->value < left->value)); 00624 if (left != left_limit) { 00625 // recurse down the tree with this element as the right-hand limit 00626 verify_structure_recursive(left, left_limit, element);} 00627 } 00628 00629 const Node * right = element->right; 00630 if (right != NULL) { 00631 assert(!(right->value < element->value)); 00632 if (right != right_limit) { 00633 // recurse down the tree with this element as the left-hand limit 00634 verify_structure_recursive(right, element, right_limit);} 00635 } 00636 } 00637 00638 //---------------------------------------------------------------------- 00639 template<class T> void SearchTree<T>::verify_structure_linear() const { 00640 00641 //print_elements(); 00642 00643 unsigned n_top = 0; 00644 unsigned n_null = 0; 00645 for(unsigned i = 0; i < _nodes.size(); i++) { 00646 const typename SearchTree<T>::Node * node = &(_nodes[i]); 00647 // make sure node is defined 00648 if (node->treelinks_null()) {n_null++; continue;} 00649 00650 // make sure of the number of "top" nodes 00651 if (node->parent == NULL) { 00652 n_top++; 00653 //assert(node->left != NULL); 00654 //assert(node->right != NULL); 00655 } else { 00656 // make sure that I am a child of my parent... 00657 //assert((node->parent->left == node) || (node->parent->right == node)); 00658 assert((node->parent->left == node) ^ (node->parent->right == node)); 00659 } 00660 00661 // when there is a left child make sure it's value is ordered 00662 // (note use of !(b<a), to allow for a<=b while using just the < 00663 // operator) 00664 if (node->left != NULL) { 00665 assert(!(node->value < node->left->value ));} 00666 00667 // when there is a right child make sure it's value is ordered 00668 if (node->right != NULL) { 00669 assert(!(node->right->value < node->value ));} 00670 00671 } 00672 assert(n_top == 1 || (n_top == 0 && size() <= 1) ); 00673 assert(n_null == _available_nodes.size() || 00674 (n_null == _available_nodes.size() + 1 && size() == 1)); 00675 } 00676 00677 00678 //---------------------------------------------------------------------- 00679 template<class T> typename SearchTree<T>::Node * SearchTree<T>::_find_predecessor(const typename SearchTree<T>::Node * node) { 00680 00681 typename SearchTree<T>::Node * newnode; 00682 if (node->left != NULL) { 00683 // go down left, and then down right as far as possible. 00684 newnode = node->left; 00685 while(newnode->right != NULL) {newnode = newnode->right;} 00686 return newnode; 00687 } else { 00688 const typename SearchTree<T>::Node * lastnode = node; 00689 newnode = node->parent; 00690 // go up the tree as long as we're going right (when we go left then 00691 // we've found something smaller, so stop) 00692 while(newnode != NULL) { 00693 if (newnode->right == lastnode) {return newnode;} 00694 lastnode = newnode; 00695 newnode = newnode->parent; 00696 } 00697 return newnode; 00698 } 00699 } 00700 00701 00702 //---------------------------------------------------------------------- 00703 template<class T> typename SearchTree<T>::Node * SearchTree<T>::_find_successor(const typename SearchTree<T>::Node * node) { 00704 00705 typename SearchTree<T>::Node * newnode; 00706 if (node->right != NULL) { 00707 // go down right, and then down left as far as possible. 00708 newnode = node->right; 00709 while(newnode->left != NULL) {newnode = newnode->left;} 00710 return newnode; 00711 } else { 00712 const typename SearchTree<T>::Node * lastnode = node; 00713 newnode = node->parent; 00714 // go up the tree as long as we're going left (when we go right then 00715 // we've found something larger, so stop) 00716 while(newnode != NULL) { 00717 if (newnode->left == lastnode) {return newnode;} 00718 lastnode = newnode; 00719 newnode = newnode->parent; 00720 } 00721 return newnode; 00722 } 00723 } 00724 00725 00726 //---------------------------------------------------------------------- 00727 // print out all the elements for visual checking... 00728 template<class T> void SearchTree<T>::print_elements() { 00729 typename SearchTree<T>::Node * base_node = &(_nodes[0]); 00730 typename SearchTree<T>::Node * node = base_node; 00731 00732 int n = _nodes.size(); 00733 for(; node - base_node < n ; node++) { 00734 printf("%4d parent:%4d left:%4d right:%4d pred:%4d succ:%4d value:%10.6f\n",loc(node), loc(node->parent), loc(node->left), loc(node->right), loc(node->predecessor),loc(node->successor),node->value); 00735 } 00736 } 00737 00738 //---------------------------------------------------------------------- 00739 template<class T> typename SearchTree<T>::circulator SearchTree<T>::somewhere() { 00740 return circulator(_top_node); 00741 } 00742 00743 00744 //---------------------------------------------------------------------- 00745 template<class T> typename SearchTree<T>::const_circulator SearchTree<T>::somewhere() const { 00746 return const_circulator(_top_node); 00747 } 00748 00749 00750 FASTJET_END_NAMESPACE 00751 00752 #endif // __FASTJET_SEARCHTREE_HH__