GRASS Programmer's Manual
6.4.2(2012)
|
00001 /**************************************************************************** 00002 * 00003 * MODULE: iostream 00004 * 00005 * COPYRIGHT (C) 2007 Laura Toma 00006 * 00007 * This program is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU General Public License as published by 00009 * the Free Software Foundation; either version 2 of the License, or 00010 * (at your option) any later version. 00011 * 00012 * This program is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU General Public License for more details. 00016 * 00017 *****************************************************************************/ 00018 00019 #ifndef _MINMAXHEAP_H 00020 #define _MINMAXHEAP_H 00021 00022 #include <stdio.h> 00023 #include <assert.h> 00024 #include <stdlib.h> 00025 #include <math.h> 00026 00027 #ifdef log2 00028 #undef log2 00029 #endif 00030 00031 #include <sstream> 00032 using namespace std; 00033 00034 #include "mm_utils.h" 00035 #include "ami_config.h" //for SAVE_MEMORY flag 00036 /* this flag is set if we are stingy on memory; in that case 'reset' 00037 deletes the pq memory and the subsequent insert reallocates it; if 00038 the operation following a reset is not an insert or an operation 00039 which does not touch the array A, behaviour is unpredictable (core 00040 dump probably) */ 00041 00042 00043 00044 00045 00046 /***************************************************************** 00047 ***************************************************************** 00048 ***************************************************************** 00049 00050 Priority queue templated on a single type (assumed to be a class with 00051 getPriority() and getValue() implemented); 00052 00053 Supported operations: min, extract_min, insert, max, extract_max in 00054 O(lg n) 00055 00056 ***************************************************************** 00057 ***************************************************************** 00058 *****************************************************************/ 00059 00060 #undef XXX 00061 #define XXX if(0) 00062 00063 00064 #define MY_LOG_DEBUG_ID(x) //inhibit debug printing 00065 //#define MY_LOG_DEBUG_ID(x) LOG_DEBUG_ID(x) 00066 00067 typedef unsigned int HeapIndex; 00068 00069 00070 template <class T> 00071 class BasicMinMaxHeap { 00072 protected: 00073 HeapIndex maxsize; 00074 HeapIndex lastindex; // last used position (0 unused) (?) 00075 T *A; 00076 00077 protected: 00078 /* couple of memory mgt functions to keep things consistent */ 00079 static T *allocateHeap(HeapIndex n); 00080 static void freeHeap(T *); 00081 00082 public: 00083 BasicMinMaxHeap(HeapIndex size) : maxsize(size) { 00084 char str[100]; 00085 sprintf(str, "BasicMinMaxHeap: allocate %ld\n", 00086 (long)((size+1)*sizeof(T))); 00087 // MEMORY_LOG(str); 00088 00089 lastindex = 0; 00090 MY_LOG_DEBUG_ID("minmaxheap: allocation"); 00091 A = allocateHeap(maxsize); 00092 }; 00093 00094 virtual ~BasicMinMaxHeap(void) { 00095 MY_LOG_DEBUG_ID("minmaxheap: deallocation"); 00096 freeHeap(A); 00097 }; 00098 00099 bool empty(void) const { return size() == 0; }; 00100 HeapIndex size() const { 00101 assert(A || !lastindex); 00102 return lastindex; 00103 }; 00104 00105 T get(HeapIndex i) const { assert(i <= size()); return A[i]; } 00106 00107 //build a heap from an array of elements; 00108 //if size > maxsize, insert first maxsize elements from array; 00109 //return nb of elements that did not fit; 00110 00111 void insert(const T& elt); 00112 00113 bool min(T& elt) const ; 00114 bool extract_min(T& elt); 00115 bool max(T& elt) const; 00116 bool extract_max(T& elt); 00117 //extract all elts with min key, add them and return their sum 00118 bool extract_all_min(T& elt); 00119 00120 void reset(); 00121 void clear(); /* mark all the data as deleted, but don't do free */ 00122 00123 void destructiveVerify(); 00124 00125 void verify(); 00126 00127 void print() const; 00128 void print_range() const; 00129 friend ostream& operator<<(ostream& s, const BasicMinMaxHeap<T> &pq) { 00130 HeapIndex i; 00131 s << "["; 00132 for(i = 1; i <= pq.size(); i++) { 00133 s << " " << pq.get(i); 00134 } 00135 s << "]"; 00136 return s; 00137 } 00138 00139 protected: 00140 virtual void grow()=0; 00141 00142 private: 00143 long log2(long n) const; 00144 int isOnMaxLevel(HeapIndex i) const { return (log2(i) % 2); }; 00145 int isOnMinLevel(HeapIndex i) const { return !isOnMaxLevel(i); }; 00146 00147 HeapIndex leftChild(HeapIndex i) const { return 2*i; }; 00148 HeapIndex rightChild(HeapIndex i) const { return 2*i + 1; }; 00149 int hasRightChild(HeapIndex i) const { return (rightChild(i) <= size()); }; 00150 int hasRightChild(HeapIndex i, HeapIndex *c) const { return ((*c=rightChild(i)) <= size()); }; 00151 HeapIndex parent(HeapIndex i) const { return (i/2); }; 00152 HeapIndex grandparent(HeapIndex i) const { return (i/4); }; 00153 int hasChildren(HeapIndex i) const { return (2*i) <= size(); }; // 1 or more 00154 void swap(HeapIndex a, HeapIndex b); 00155 00156 T leftChildValue(HeapIndex i) const; 00157 T rightChildValue(HeapIndex i) const; 00158 HeapIndex smallestChild(HeapIndex i) const; 00159 HeapIndex smallestChildGrandchild(HeapIndex i) const; 00160 HeapIndex largestChild(HeapIndex i) const; 00161 HeapIndex largestChildGrandchild(HeapIndex i) const; 00162 int isGrandchildOf(HeapIndex i, HeapIndex m) const; 00163 00164 void trickleDownMin(HeapIndex i); 00165 void trickleDownMax(HeapIndex i); 00166 void trickleDown(HeapIndex i); 00167 00168 void bubbleUp(HeapIndex i); 00169 void bubbleUpMin(HeapIndex i); 00170 void bubbleUpMax(HeapIndex i); 00171 }; 00172 00173 00174 // index 0 is invalid 00175 // index <= size 00176 00177 // ---------------------------------------------------------------------- 00178 00179 template <class T> 00180 long BasicMinMaxHeap<T>::log2(long n) const { 00181 long i=-1; 00182 // let log2(0)==-1 00183 while(n) { 00184 n = n >> 1; 00185 i++; 00186 } 00187 return i; 00188 } 00189 00190 00191 // ---------------------------------------------------------------------- 00192 00193 template <class T> 00194 void BasicMinMaxHeap<T>::swap(HeapIndex a, HeapIndex b) { 00195 T tmp; 00196 tmp = A[a]; 00197 A[a] = A[b]; 00198 A[b] = tmp; 00199 } 00200 00201 00202 // ---------------------------------------------------------------------- 00203 00204 // child must exist 00205 template <class T> 00206 T BasicMinMaxHeap<T>::leftChildValue(HeapIndex i) const { 00207 HeapIndex p = leftChild(i); 00208 assert(p <= size()); 00209 return A[p]; 00210 } 00211 00212 // ---------------------------------------------------------------------- 00213 00214 // child must exist 00215 template <class T> 00216 T BasicMinMaxHeap<T>::rightChildValue(HeapIndex i) const { 00217 HeapIndex p = rightChild(i); 00218 assert(p <= size()); 00219 return A[p]; 00220 } 00221 00222 00223 // ---------------------------------------------------------------------- 00224 00225 // returns index of the smallest of children of node 00226 // it is an error to call this function if node has no children 00227 template <class T> 00228 HeapIndex BasicMinMaxHeap<T>::smallestChild(HeapIndex i) const { 00229 assert(hasChildren(i)); 00230 if(hasRightChild(i) && (leftChildValue(i) > rightChildValue(i))) { 00231 return rightChild(i); 00232 } else { 00233 return leftChild(i); 00234 } 00235 } 00236 00237 // ---------------------------------------------------------------------- 00238 00239 template <class T> 00240 HeapIndex BasicMinMaxHeap<T>::largestChild(HeapIndex i) const { 00241 assert(hasChildren(i)); 00242 if(hasRightChild(i) && (leftChildValue(i) < rightChildValue(i))) { 00243 return rightChild(i); 00244 } else { 00245 return leftChild(i); 00246 } 00247 } 00248 00249 // ---------------------------------------------------------------------- 00250 00251 // error to call on node without children 00252 template <class T> 00253 HeapIndex BasicMinMaxHeap<T>::smallestChildGrandchild(HeapIndex i) const { 00254 HeapIndex p,q; 00255 HeapIndex minpos = 0; 00256 00257 assert(hasChildren(i)); 00258 00259 p = leftChild(i); 00260 if(hasChildren(p)) { 00261 q = smallestChild(p); 00262 if(A[p] > A[q]) p = q; 00263 } 00264 // p is smallest of left child, its grandchildren 00265 minpos = p; 00266 00267 if(hasRightChild(i,&p)) { 00268 //p = rightChild(i); 00269 if(hasChildren(p)) { 00270 q = smallestChild(p); 00271 if(A[p] > A[q]) p = q; 00272 } 00273 // p is smallest of right child, its grandchildren 00274 if(A[p] < A[minpos]) minpos = p; 00275 } 00276 return minpos; 00277 } 00278 00279 // ---------------------------------------------------------------------- 00280 00281 template <class T> 00282 HeapIndex BasicMinMaxHeap<T>::largestChildGrandchild(HeapIndex i) const { 00283 HeapIndex p,q; 00284 HeapIndex maxpos = 0; 00285 00286 assert(hasChildren(i)); 00287 00288 p = leftChild(i); 00289 if(hasChildren(p)) { 00290 q = largestChild(p); 00291 if(A[p] < A[q]) p = q; 00292 } 00293 // p is smallest of left child, its grandchildren 00294 maxpos = p; 00295 00296 if(hasRightChild(i,&p)) { 00297 //p = rightChild(i); 00298 if(hasChildren(p)) { 00299 q = largestChild(p); 00300 if(A[p] < A[q]) p = q; 00301 } 00302 // p is smallest of right child, its grandchildren 00303 if(A[p] > A[maxpos]) maxpos = p; 00304 } 00305 return maxpos; 00306 } 00307 00308 // ---------------------------------------------------------------------- 00309 00310 // this is pretty loose - only to differentiate between child and grandchild 00311 template <class T> 00312 int BasicMinMaxHeap<T>::isGrandchildOf(HeapIndex i, HeapIndex m) const { 00313 return (m >= i*4); 00314 } 00315 00316 // ---------------------------------------------------------------------- 00317 00318 template <class T> 00319 void BasicMinMaxHeap<T>::trickleDownMin(HeapIndex i) { 00320 HeapIndex m; 00321 bool done = false; 00322 00323 while (!done) { 00324 00325 if (!hasChildren(i)) { 00326 done = true; 00327 return; 00328 } 00329 m = smallestChildGrandchild(i); 00330 if(isGrandchildOf(i, m)) { 00331 if(A[m] < A[i]) { 00332 swap(i, m); 00333 if(A[m] > A[parent(m)]) { 00334 swap(m, parent(m)); 00335 } 00336 //trickleDownMin(m); 00337 i = m; 00338 } else { 00339 done = true; 00340 } 00341 } else { 00342 if(A[m] < A[i]) { 00343 swap(i, m); 00344 } 00345 done = true; 00346 } 00347 }//while 00348 } 00349 00350 // ---------------------------------------------------------------------- 00351 00352 // unverified 00353 template <class T> 00354 void BasicMinMaxHeap<T>::trickleDownMax(HeapIndex i) { 00355 HeapIndex m; 00356 bool done = false; 00357 00358 while (!done) { 00359 if(!hasChildren(i)) { 00360 done = true; 00361 return; 00362 } 00363 00364 m = largestChildGrandchild(i); 00365 if(isGrandchildOf(i, m)) { 00366 if(A[m] > A[i]) { 00367 swap(i, m); 00368 if(A[m] < A[parent(m)]) { 00369 swap(m, parent(m)); 00370 } 00371 //trickleDownMax(m); 00372 i = m; 00373 } else { 00374 done = true; 00375 } 00376 } else { 00377 if(A[m] > A[i]) { 00378 swap(i, m); 00379 } 00380 done = true; 00381 } 00382 } //while 00383 } 00384 00385 00386 // ---------------------------------------------------------------------- 00387 00388 00389 template <class T> 00390 void BasicMinMaxHeap<T>::trickleDown(HeapIndex i) { 00391 if(isOnMinLevel(i)) { 00392 trickleDownMin(i); 00393 } else { 00394 trickleDownMax(i); 00395 } 00396 } 00397 00398 // ---------------------------------------------------------------------- 00399 template <class T> 00400 void BasicMinMaxHeap<T>::bubbleUp(HeapIndex i) { 00401 HeapIndex m; 00402 m = parent(i); 00403 00404 if(isOnMinLevel(i)) { 00405 if (m && (A[i] > A[m])) { 00406 swap(i, m); 00407 bubbleUpMax(m); 00408 } else { 00409 bubbleUpMin(i); 00410 } 00411 } else { 00412 if (m && (A[i] < A[m])) { 00413 swap(i, m); 00414 bubbleUpMin(m); 00415 } else { 00416 bubbleUpMax(i); 00417 } 00418 } 00419 } 00420 00421 00422 // ---------------------------------------------------------------------- 00423 template <class T> 00424 void BasicMinMaxHeap<T>::bubbleUpMin(HeapIndex i) { 00425 HeapIndex m; 00426 m = grandparent(i); 00427 00428 while (m && (A[i] < A[m])) { 00429 swap(i,m); 00430 //bubbleUpMin(m); 00431 i = m; 00432 m = grandparent(i); 00433 00434 } 00435 } 00436 00437 00438 00439 // ---------------------------------------------------------------------- 00440 template <class T> 00441 void BasicMinMaxHeap<T>::bubbleUpMax(HeapIndex i) { 00442 HeapIndex m; 00443 m = grandparent(i); 00444 00445 while(m && (A[i] > A[m])) { 00446 swap(i,m); 00447 //bubbleUpMax(m); 00448 i=m; 00449 m = grandparent(i); 00450 } 00451 } 00452 00453 00454 #if(0) 00455 // ---------------------------------------------------------------------- 00456 template <class T> 00457 void BasicMinMaxHeap<T>::print_rajiv() const { 00458 HeapIndex i; 00459 ostrstream *ostr = new ostrstream(); 00460 00461 *ostr << "[1]"; 00462 for(i=1; i<=size(); i++) { 00463 *ostr << " " << A[i]; 00464 if(ostr->pcount() > 70) { 00465 cout << ostr->str() << endl; 00466 delete ostr; 00467 ostr = new ostrstream(); 00468 *ostr << "[" << i << "]"; 00469 } 00470 } 00471 cout << ostr->str() << endl; 00472 } 00473 #endif 00474 00475 00476 00477 // ---------------------------------------------------------------------- 00478 template <class T> 00479 void BasicMinMaxHeap<T>::print() const { 00480 cout << "["; 00481 for (unsigned int i=1; i<=size(); i++) { 00482 cout << A[i].getPriority() <<","; 00483 } 00484 cout << "]" << endl; 00485 } 00486 00487 // ---------------------------------------------------------------------- 00488 template <class T> 00489 void BasicMinMaxHeap<T>::print_range() const { 00490 cout << "["; 00491 T a, b; 00492 min(a); 00493 max(b); 00494 if (size) { 00495 cout << a.getPriority() << ".." 00496 << b.getPriority(); 00497 } 00498 cout << " (" << size() << ")]"; 00499 } 00500 00501 00502 // ---------------------------------------------------------------------- 00503 template <class T> 00504 void BasicMinMaxHeap<T>::insert(const T& elt) { 00505 #ifdef SAVE_MEMORY 00506 if (!A) { 00507 MY_LOG_DEBUG_ID("minmaxheap: re-allocation"); 00508 A = allocateHeap(maxsize); 00509 } 00510 #endif 00511 00512 if(lastindex == maxsize) grow(); 00513 00514 XXX cerr << "insert: " << elt << endl; 00515 00516 lastindex++; 00517 A[lastindex] = elt; 00518 bubbleUp(lastindex); 00519 } 00520 00521 // ---------------------------------------------------------------------- 00522 template <class T> 00523 bool BasicMinMaxHeap<T>::extract_min(T& elt) { 00524 00525 assert(A); 00526 00527 if(lastindex == 0) return false; 00528 00529 elt = A[1]; 00530 A[1] = A[lastindex]; 00531 lastindex--; 00532 trickleDown(1); 00533 00534 return true; 00535 } 00536 00537 // ---------------------------------------------------------------------- 00538 //extract all elts with min key, add them and return their sum 00539 template <class T> 00540 bool BasicMinMaxHeap<T>::extract_all_min(T& elt) { 00541 T next_elt; 00542 bool done = false; 00543 00544 //extract first elt 00545 if (!extract_min(elt)) { 00546 return false; 00547 } else { 00548 while (!done) { 00549 //peek at the next min elt to see if matches 00550 if ((!min(next_elt)) || 00551 !(next_elt.getPriority() == elt.getPriority())) { 00552 done = true; 00553 } else { 00554 extract_min(next_elt); 00555 elt = elt + next_elt; 00556 } 00557 } 00558 } 00559 return true; 00560 } 00561 00562 // ---------------------------------------------------------------------- 00563 template <class T> 00564 bool BasicMinMaxHeap<T>::extract_max(T& elt) { 00565 00566 assert(A); 00567 00568 HeapIndex p; // max 00569 if(lastindex == 0) return false; 00570 00571 if(hasChildren(1)) { 00572 p = largestChild(1); 00573 } else { 00574 p = 1; 00575 } 00576 elt = A[p]; 00577 A[p] = A[lastindex]; 00578 lastindex--; 00579 trickleDown(p); 00580 00581 return true; 00582 } 00583 00584 // ---------------------------------------------------------------------- 00585 template <class T> 00586 bool BasicMinMaxHeap<T>::min(T& elt) const { 00587 00588 assert(A); 00589 00590 if(lastindex == 0) return false; 00591 00592 elt = A[1]; 00593 return true; 00594 } 00595 00596 // ---------------------------------------------------------------------- 00597 template <class T> 00598 bool BasicMinMaxHeap<T>::max(T& elt) const { 00599 00600 assert(A); 00601 00602 HeapIndex p; // max 00603 if(lastindex == 0) return false; 00604 00605 if(hasChildren(1)) { 00606 p = largestChild(1); 00607 } else { 00608 p = 1; 00609 } 00610 elt = A[p]; 00611 return true; 00612 } 00613 00614 00615 00616 // ---------------------------------------------------------------------- 00617 //free memory if SAVE_MEMORY is set 00618 template <class T> 00619 void BasicMinMaxHeap<T>::reset() { 00620 #ifdef SAVE_MEMORY 00621 assert(empty()); 00622 MY_LOG_DEBUG_ID("minmaxheap: deallocation"); 00623 freeHeap(A); 00624 A = NULL; 00625 #endif 00626 } 00627 00628 // ---------------------------------------------------------------------- 00629 00630 template <class T> 00631 void 00632 BasicMinMaxHeap<T>::clear() { 00633 lastindex = 0; 00634 } 00635 00636 // ---------------------------------------------------------------------- 00637 template <class T> 00638 T * 00639 BasicMinMaxHeap<T>::allocateHeap(HeapIndex n) { 00640 T *p; 00641 #ifdef USE_LARGEMEM 00642 p = (T*)LargeMemory::alloc(sizeof(T) * (n+1)); 00643 #else 00644 p = new T[n+1]; 00645 #endif 00646 return p; 00647 } 00648 00649 // ---------------------------------------------------------------------- 00650 template <class T> 00651 void 00652 BasicMinMaxHeap<T>::freeHeap(T *p) { 00653 if (p) { 00654 #ifdef USE_LARGEMEM 00655 LargeMemory::free(p); 00656 #else 00657 delete [] p; 00658 #endif 00659 } 00660 } 00661 00662 00663 // ---------------------------------------------------------------------- 00664 00665 template <class T> 00666 void 00667 BasicMinMaxHeap<T>::destructiveVerify() { 00668 HeapIndex n = size(); 00669 T val, prev; 00670 bool ok; 00671 00672 if(!n) return; 00673 00674 XXX print(); 00675 00676 /* make sure that min works */ 00677 extract_min(prev); 00678 for(HeapIndex i=1; i<n; i++) { 00679 ok = min(val); 00680 assert(ok); 00681 XXX cerr << i << ": " << val << endl; 00682 if(val.getPriority() < prev.getPriority()) { // oops! 00683 print(); 00684 cerr << "n=" << n << endl; 00685 cerr << "val=" << val << endl; 00686 cerr << "prev=" << prev << endl; 00687 cerr << "looks like minmaxheap.min is broken!!" << endl; 00688 assert(0); 00689 return; 00690 } 00691 prev = val; 00692 ok = extract_min(val); 00693 assert(ok); 00694 assert(prev == val); 00695 } 00696 } 00697 00698 00699 // ---------------------------------------------------------------------- 00700 00701 template <class T> 00702 void 00703 BasicMinMaxHeap<T>::verify() { 00704 long n = size(); 00705 T *dup; 00706 00707 if(!n) return; 00708 00709 dup = allocateHeap(maxsize); 00710 for(HeapIndex i=0; i<n+1; i++) { 00711 dup[i] = A[i]; 00712 } 00713 destructiveVerify(); 00714 freeHeap(A); 00715 /* restore the heap */ 00716 A = dup; 00717 lastindex = n; 00718 } 00719 00720 00721 // ---------------------------------------------------------------------- 00722 // ---------------------------------------------------------------------- 00723 00724 template <class T> 00725 class MinMaxHeap : public BasicMinMaxHeap<T> { 00726 public: 00727 MinMaxHeap(HeapIndex size) : BasicMinMaxHeap<T>(size) {}; 00728 virtual ~MinMaxHeap() {}; 00729 bool full(void) const { return this->size() >= this->maxsize; }; 00730 HeapIndex get_maxsize() const { return this->maxsize; }; 00731 HeapIndex fill(T* arr, HeapIndex n); 00732 00733 protected: 00734 virtual void grow() { fprintf(stderr, "MinMaxHeap::grow: not implemented\n"); assert(0); exit(1); }; 00735 }; 00736 00737 // ---------------------------------------------------------------------- 00738 //build a heap from an array of elements; 00739 //if size > maxsize, insert first maxsize elements from array; 00740 //return nb of elements that did not fit; 00741 template <class T> 00742 HeapIndex MinMaxHeap<T>::fill(T* arr, HeapIndex n) { 00743 HeapIndex i; 00744 //heap must be empty 00745 assert(this->size()==0); 00746 for (i = 0; !full() && i<n; i++) { 00747 insert(arr[i]); 00748 } 00749 if (i < n) { 00750 assert(i == this->maxsize); 00751 return n - i; 00752 } else { 00753 return 0; 00754 } 00755 } 00756 00757 00758 00759 #define MMHEAP_INITIAL_SIZE 1024 00760 00761 template <class T> 00762 class UnboundedMinMaxHeap : public BasicMinMaxHeap<T> { 00763 public: 00764 UnboundedMinMaxHeap() : BasicMinMaxHeap<T>(MMHEAP_INITIAL_SIZE) {}; 00765 UnboundedMinMaxHeap(HeapIndex size) : BasicMinMaxHeap<T>(size) {}; 00766 virtual ~UnboundedMinMaxHeap() {}; 00767 protected: 00768 virtual void grow(); 00769 }; 00770 00771 template <class T> 00772 void UnboundedMinMaxHeap<T>::grow() { 00773 T *old = this->A; 00774 this->maxsize *= 2; 00775 00776 assert(this->maxsize > 0); 00777 00778 if(old) { 00779 HeapIndex n = this->size(); 00780 this->A = allocateHeap(this->maxsize); /* allocate a new array */ 00781 /* copy over the old values */ 00782 assert(this->maxsize > n); 00783 for(HeapIndex i=0; i<=n; i++) { /* why extra value? -RW */ 00784 this->A[i] = old[i]; 00785 } 00786 freeHeap(old); /* free up old storage */ 00787 } 00788 00789 } 00790 00791 00792 #endif