GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
minmaxheap.h
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * MODULE: iostream
4  *
5 
6  * COPYRIGHT (C) 2007 Laura Toma
7  *
8  *
9 
10  * Iostream is a library that implements streams, external memory
11  * sorting on streams, and an external memory priority queue on
12  * streams. These are the fundamental components used in external
13  * memory algorithms.
14 
15  * Credits: The library was developed by Laura Toma. The kernel of
16  * class STREAM is based on the similar class existent in the GPL TPIE
17  * project developed at Duke University. The sorting and priority
18  * queue have been developed by Laura Toma based on communications
19  * with Rajiv Wickremesinghe. The library was developed as part of
20  * porting Terraflow to GRASS in 2001. PEARL upgrades in 2003 by
21  * Rajiv Wickremesinghe as part of the Terracost project.
22 
23  *
24  * This program is free software; you can redistribute it and/or modify
25  * it under the terms of the GNU General Public License as published by
26  * the Free Software Foundation; either version 2 of the License, or
27  * (at your option) any later version.
28  *
29 
30  * This program is distributed in the hope that it will be useful,
31  * but WITHOUT ANY WARRANTY; without even the implied warranty of
32  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
33  * General Public License for more details. *
34  * **************************************************************************/
35 
36 
37 #ifndef _MINMAXHEAP_H
38 #define _MINMAXHEAP_H
39 
40 #include <stdio.h>
41 #include <assert.h>
42 #include <stdlib.h>
43 #include <math.h>
44 
45 #ifdef log2
46 #undef log2
47 #endif
48 
49 #include <sstream>
50 
51 #include "mm_utils.h"
52 #include "ami_config.h" //for SAVE_MEMORY flag
53 /* this flag is set if we are stingy on memory; in that case 'reset'
54  deletes the pq memory and the subsequent insert reallocates it; if
55  the operation following a reset is not an insert or an operation
56  which does not touch the array A, behaviour is unpredictable (core
57  dump probably) */
58 
59 
60 
61 
62 
63 /*****************************************************************
64  *****************************************************************
65  *****************************************************************
66 
67 Priority queue templated on a single type (assumed to be a class with
68 getPriority() and getValue() implemented);
69 
70 Supported operations: min, extract_min, insert, max, extract_max in
71 O(lg n)
72 
73 *****************************************************************
74 *****************************************************************
75 *****************************************************************/
76 
77 #undef XXX
78 #define XXX if(0)
79 
80 
81 #define MY_LOG_DEBUG_ID(x) //inhibit debug printing
82 //#define MY_LOG_DEBUG_ID(x) LOG_DEBUG_ID(x)
83 
84 typedef unsigned int HeapIndex;
85 
86 
87 template <class T>
89 protected:
91  HeapIndex lastindex; // last used position (0 unused) (?)
92  T *A;
93 
94 protected:
95  /* couple of memory mgt functions to keep things consistent */
96  static T *allocateHeap(HeapIndex n);
97  static void freeHeap(T *);
98 
99 public:
100  BasicMinMaxHeap(HeapIndex size) : maxsize(size) {
101  char str[100];
102  sprintf(str, "BasicMinMaxHeap: allocate %ld\n",
103  (long)((size+1)*sizeof(T)));
104  // MEMORY_LOG(str);
105 
106  lastindex = 0;
107  MY_LOG_DEBUG_ID("minmaxheap: allocation");
108  A = allocateHeap(maxsize);
109  };
110 
111  virtual ~BasicMinMaxHeap(void) {
112  MY_LOG_DEBUG_ID("minmaxheap: deallocation");
113  freeHeap(A);
114  };
115 
116  bool empty(void) const { return size() == 0; };
117  HeapIndex size() const {
118  assert(A || !lastindex);
119  return lastindex;
120  };
121 
122  T get(HeapIndex i) const { assert(i <= size()); return A[i]; }
123 
124  //build a heap from an array of elements;
125  //if size > maxsize, insert first maxsize elements from array;
126  //return nb of elements that did not fit;
127 
128  void insert(const T& elt);
129 
130  bool min(T& elt) const ;
131  bool extract_min(T& elt);
132  bool max(T& elt) const;
133  bool extract_max(T& elt);
134  //extract all elts with min key, add them and return their sum
135  bool extract_all_min(T& elt);
136 
137  void reset();
138  void clear(); /* mark all the data as deleted, but don't do free */
139 
140  void destructiveVerify();
141 
142  void verify();
143 
144  void print() const;
145  void print_range() const;
146  friend ostream& operator<<(ostream& s, const BasicMinMaxHeap<T> &pq) {
147  HeapIndex i;
148  s << "[";
149  for(i = 1; i <= pq.size(); i++) {
150  s << " " << pq.get(i);
151  }
152  s << "]";
153  return s;
154  }
155 
156 protected:
157  virtual void grow()=0;
158 
159 private:
160  long log2(long n) const;
161  int isOnMaxLevel(HeapIndex i) const { return (log2(i) % 2); };
162  int isOnMinLevel(HeapIndex i) const { return !isOnMaxLevel(i); };
163 
164  HeapIndex leftChild(HeapIndex i) const { return 2*i; };
165  HeapIndex rightChild(HeapIndex i) const { return 2*i + 1; };
166  int hasRightChild(HeapIndex i) const { return (rightChild(i) <= size()); };
167  int hasRightChild(HeapIndex i, HeapIndex *c) const { return ((*c=rightChild(i)) <= size()); };
168  HeapIndex parent(HeapIndex i) const { return (i/2); };
169  HeapIndex grandparent(HeapIndex i) const { return (i/4); };
170  int hasChildren(HeapIndex i) const { return (2*i) <= size(); }; // 1 or more
171  void swap(HeapIndex a, HeapIndex b);
172 
173  T leftChildValue(HeapIndex i) const;
174  T rightChildValue(HeapIndex i) const;
175  HeapIndex smallestChild(HeapIndex i) const;
176  HeapIndex smallestChildGrandchild(HeapIndex i) const;
177  HeapIndex largestChild(HeapIndex i) const;
178  HeapIndex largestChildGrandchild(HeapIndex i) const;
179  int isGrandchildOf(HeapIndex i, HeapIndex m) const;
180 
181  void trickleDownMin(HeapIndex i);
182  void trickleDownMax(HeapIndex i);
183  void trickleDown(HeapIndex i);
184 
185  void bubbleUp(HeapIndex i);
186  void bubbleUpMin(HeapIndex i);
187  void bubbleUpMax(HeapIndex i);
188 };
189 
190 
191 // index 0 is invalid
192 // index <= size
193 
194 // ----------------------------------------------------------------------
195 
196 template <class T>
197 long BasicMinMaxHeap<T>::log2(long n) const {
198  long i=-1;
199  // let log2(0)==-1
200  while(n) {
201  n = n >> 1;
202  i++;
203  }
204  return i;
205 }
206 
207 
208 // ----------------------------------------------------------------------
209 
210 template <class T>
212  T tmp;
213  tmp = A[a];
214  A[a] = A[b];
215  A[b] = tmp;
216 }
217 
218 
219 // ----------------------------------------------------------------------
220 
221 // child must exist
222 template <class T>
224  HeapIndex p = leftChild(i);
225  assert(p <= size());
226  return A[p];
227 }
228 
229 // ----------------------------------------------------------------------
230 
231 // child must exist
232 template <class T>
234  HeapIndex p = rightChild(i);
235  assert(p <= size());
236  return A[p];
237 }
238 
239 
240 // ----------------------------------------------------------------------
241 
242 // returns index of the smallest of children of node
243 // it is an error to call this function if node has no children
244 template <class T>
246  assert(hasChildren(i));
247  if(hasRightChild(i) && (leftChildValue(i) > rightChildValue(i))) {
248  return rightChild(i);
249  } else {
250  return leftChild(i);
251  }
252 }
253 
254 // ----------------------------------------------------------------------
255 
256 template <class T>
258  assert(hasChildren(i));
259  if(hasRightChild(i) && (leftChildValue(i) < rightChildValue(i))) {
260  return rightChild(i);
261  } else {
262  return leftChild(i);
263  }
264 }
265 
266 // ----------------------------------------------------------------------
267 
268 // error to call on node without children
269 template <class T>
271  HeapIndex p,q;
272  HeapIndex minpos = 0;
273 
274  assert(hasChildren(i));
275 
276  p = leftChild(i);
277  if(hasChildren(p)) {
278  q = smallestChild(p);
279  if(A[p] > A[q]) p = q;
280  }
281  // p is smallest of left child, its grandchildren
282  minpos = p;
283 
284  if(hasRightChild(i,&p)) {
285  //p = rightChild(i);
286  if(hasChildren(p)) {
287  q = smallestChild(p);
288  if(A[p] > A[q]) p = q;
289  }
290  // p is smallest of right child, its grandchildren
291  if(A[p] < A[minpos]) minpos = p;
292  }
293  return minpos;
294 }
295 
296 // ----------------------------------------------------------------------
297 
298 template <class T>
300  HeapIndex p,q;
301  HeapIndex maxpos = 0;
302 
303  assert(hasChildren(i));
304 
305  p = leftChild(i);
306  if(hasChildren(p)) {
307  q = largestChild(p);
308  if(A[p] < A[q]) p = q;
309  }
310  // p is smallest of left child, its grandchildren
311  maxpos = p;
312 
313  if(hasRightChild(i,&p)) {
314  //p = rightChild(i);
315  if(hasChildren(p)) {
316  q = largestChild(p);
317  if(A[p] < A[q]) p = q;
318  }
319  // p is smallest of right child, its grandchildren
320  if(A[p] > A[maxpos]) maxpos = p;
321  }
322  return maxpos;
323 }
324 
325 // ----------------------------------------------------------------------
326 
327 // this is pretty loose - only to differentiate between child and grandchild
328 template <class T>
330  return (m >= i*4);
331 }
332 
333 // ----------------------------------------------------------------------
334 
335 template <class T>
337  HeapIndex m;
338  bool done = false;
339 
340  while (!done) {
341 
342  if (!hasChildren(i)) {
343  done = true;
344  return;
345  }
346  m = smallestChildGrandchild(i);
347  if(isGrandchildOf(i, m)) {
348  if(A[m] < A[i]) {
349  swap(i, m);
350  if(A[m] > A[parent(m)]) {
351  swap(m, parent(m));
352  }
353  //trickleDownMin(m);
354  i = m;
355  } else {
356  done = true;
357  }
358  } else {
359  if(A[m] < A[i]) {
360  swap(i, m);
361  }
362  done = true;
363  }
364  }//while
365 }
366 
367 // ----------------------------------------------------------------------
368 
369 // unverified
370 template <class T>
372  HeapIndex m;
373  bool done = false;
374 
375  while (!done) {
376  if(!hasChildren(i)) {
377  done = true;
378  return;
379  }
380 
381  m = largestChildGrandchild(i);
382  if(isGrandchildOf(i, m)) {
383  if(A[m] > A[i]) {
384  swap(i, m);
385  if(A[m] < A[parent(m)]) {
386  swap(m, parent(m));
387  }
388  //trickleDownMax(m);
389  i = m;
390  } else {
391  done = true;
392  }
393  } else {
394  if(A[m] > A[i]) {
395  swap(i, m);
396  }
397  done = true;
398  }
399  } //while
400 }
401 
402 
403 // ----------------------------------------------------------------------
404 
405 
406 template <class T>
408  if(isOnMinLevel(i)) {
409  trickleDownMin(i);
410  } else {
411  trickleDownMax(i);
412  }
413 }
414 
415 // ----------------------------------------------------------------------
416 template <class T>
418  HeapIndex m;
419  m = parent(i);
420 
421  if(isOnMinLevel(i)) {
422  if (m && (A[i] > A[m])) {
423  swap(i, m);
424  bubbleUpMax(m);
425  } else {
426  bubbleUpMin(i);
427  }
428  } else {
429  if (m && (A[i] < A[m])) {
430  swap(i, m);
431  bubbleUpMin(m);
432  } else {
433  bubbleUpMax(i);
434  }
435  }
436 }
437 
438 
439 // ----------------------------------------------------------------------
440 template <class T>
442  HeapIndex m;
443  m = grandparent(i);
444 
445  while (m && (A[i] < A[m])) {
446  swap(i,m);
447  //bubbleUpMin(m);
448  i = m;
449  m = grandparent(i);
450 
451  }
452 }
453 
454 
455 
456 // ----------------------------------------------------------------------
457 template <class T>
459  HeapIndex m;
460  m = grandparent(i);
461 
462  while(m && (A[i] > A[m])) {
463  swap(i,m);
464  //bubbleUpMax(m);
465  i=m;
466  m = grandparent(i);
467  }
468 }
469 
470 
471 #if(0)
472 // ----------------------------------------------------------------------
473 template <class T>
474 void BasicMinMaxHeap<T>::print_rajiv() const {
475  HeapIndex i;
476  ostrstream *ostr = new ostrstream();
477 
478  *ostr << "[1]";
479  for(i=1; i<=size(); i++) {
480  *ostr << " " << A[i];
481  if(ostr->pcount() > 70) {
482  cout << ostr->str() << endl;
483  delete ostr;
484  ostr = new ostrstream();
485  *ostr << "[" << i << "]";
486  }
487  }
488  cout << ostr->str() << endl;
489 }
490 #endif
491 
492 
493 
494 // ----------------------------------------------------------------------
495 template <class T>
497  cout << "[";
498  for (unsigned int i=1; i<=size(); i++) {
499  cout << A[i].getPriority() <<",";
500  }
501  cout << "]" << endl;
502 }
503 
504 // ----------------------------------------------------------------------
505 template <class T>
507  cout << "[";
508  T a, b;
509  min(a);
510  max(b);
511  if (size()) {
512  cout << a.getPriority() << ".."
513  << b.getPriority();
514  }
515  cout << " (" << size() << ")]";
516 }
517 
518 
519 // ----------------------------------------------------------------------
520 template <class T>
521 void BasicMinMaxHeap<T>::insert(const T& elt) {
522 #ifdef SAVE_MEMORY
523  if (!A) {
524  MY_LOG_DEBUG_ID("minmaxheap: re-allocation");
526  }
527 #endif
528 
529  if(lastindex == maxsize) grow();
530 
531  XXX cerr << "insert: " << elt << endl;
532 
533  lastindex++;
534  A[lastindex] = elt;
535  bubbleUp(lastindex);
536 }
537 
538 // ----------------------------------------------------------------------
539 template <class T>
541 
542  assert(A);
543 
544  if(lastindex == 0) return false;
545 
546  elt = A[1];
547  A[1] = A[lastindex];
548  lastindex--;
549  trickleDown(1);
550 
551  return true;
552 }
553 
554 // ----------------------------------------------------------------------
555 //extract all elts with min key, add them and return their sum
556 template <class T>
558  T next_elt;
559  bool done = false;
560 
561  //extract first elt
562  if (!extract_min(elt)) {
563  return false;
564  } else {
565  while (!done) {
566  //peek at the next min elt to see if matches
567  if ((!min(next_elt)) ||
568  !(next_elt.getPriority() == elt.getPriority())) {
569  done = true;
570  } else {
571  extract_min(next_elt);
572  elt = elt + next_elt;
573  }
574  }
575  }
576  return true;
577 }
578 
579 // ----------------------------------------------------------------------
580 template <class T>
582 
583  assert(A);
584 
585  HeapIndex p; // max
586  if(lastindex == 0) return false;
587 
588  if(hasChildren(1)) {
589  p = largestChild(1);
590  } else {
591  p = 1;
592  }
593  elt = A[p];
594  A[p] = A[lastindex];
595  lastindex--;
596  trickleDown(p);
597 
598  return true;
599 }
600 
601 // ----------------------------------------------------------------------
602 template <class T>
603 bool BasicMinMaxHeap<T>::min(T& elt) const {
604 
605  assert(A);
606 
607  if(lastindex == 0) return false;
608 
609  elt = A[1];
610  return true;
611 }
612 
613 // ----------------------------------------------------------------------
614 template <class T>
615 bool BasicMinMaxHeap<T>::max(T& elt) const {
616 
617  assert(A);
618 
619  HeapIndex p; // max
620  if(lastindex == 0) return false;
621 
622  if(hasChildren(1)) {
623  p = largestChild(1);
624  } else {
625  p = 1;
626  }
627  elt = A[p];
628  return true;
629 }
630 
631 
632 
633 // ----------------------------------------------------------------------
634 //free memory if SAVE_MEMORY is set
635 template <class T>
637 #ifdef SAVE_MEMORY
638  assert(empty());
639  MY_LOG_DEBUG_ID("minmaxheap: deallocation");
640  freeHeap(A);
641  A = NULL;
642 #endif
643 }
644 
645 // ----------------------------------------------------------------------
646 
647 template <class T>
648 void
650  lastindex = 0;
651 }
652 
653 // ----------------------------------------------------------------------
654 template <class T>
655 T *
657  T *p;
658 #ifdef USE_LARGEMEM
659  p = (T*)LargeMemory::alloc(sizeof(T) * (n+1));
660 #else
661  p = new T[n+1];
662 #endif
663  return p;
664 }
665 
666 // ----------------------------------------------------------------------
667 template <class T>
668 void
670  if (p) {
671 #ifdef USE_LARGEMEM
673 #else
674  delete [] p;
675 #endif
676  }
677 }
678 
679 
680 // ----------------------------------------------------------------------
681 
682 template <class T>
683 void
685  HeapIndex n = size();
686  T val, prev;
687  bool ok;
688 
689  if(!n) return;
690 
691  XXX print();
692 
693  /* make sure that min works */
694  extract_min(prev);
695  for(HeapIndex i=1; i<n; i++) {
696  ok = min(val);
697  assert(ok);
698  XXX cerr << i << ": " << val << endl;
699  if(val.getPriority() < prev.getPriority()) { // oops!
700  print();
701  cerr << "n=" << n << endl;
702  cerr << "val=" << val << endl;
703  cerr << "prev=" << prev << endl;
704  cerr << "looks like minmaxheap.min is broken!!" << endl;
705  assert(0);
706  return;
707  }
708  prev = val;
709  ok = extract_min(val);
710  assert(ok);
711  assert(prev == val);
712  }
713 }
714 
715 
716 // ----------------------------------------------------------------------
717 
718 template <class T>
719 void
721  long n = size();
722  T *dup;
723 
724  if(!n) return;
725 
726  dup = allocateHeap(maxsize);
727  for(HeapIndex i=0; i<n+1; i++) {
728  dup[i] = A[i];
729  }
731  freeHeap(A);
732  /* restore the heap */
733  A = dup;
734  lastindex = n;
735 }
736 
737 
738 // ----------------------------------------------------------------------
739 // ----------------------------------------------------------------------
740 
741 template <class T>
742 class MinMaxHeap : public BasicMinMaxHeap<T> {
743 public:
745  virtual ~MinMaxHeap() {};
746  bool full(void) const { return this->size() >= this->maxsize; };
747  HeapIndex get_maxsize() const { return this->maxsize; };
748  HeapIndex fill(T* arr, HeapIndex n);
749 
750 protected:
751  virtual void grow() { fprintf(stderr, "MinMaxHeap::grow: not implemented\n"); assert(0); exit(1); };
752 };
753 
754 // ----------------------------------------------------------------------
755 //build a heap from an array of elements;
756 //if size > maxsize, insert first maxsize elements from array;
757 //return nb of elements that did not fit;
758 template <class T>
760  HeapIndex i;
761  //heap must be empty
762  assert(this->size()==0);
763  for (i = 0; !full() && i<n; i++) {
764  this->insert(arr[i]);
765  }
766  if (i < n) {
767  assert(i == this->maxsize);
768  return n - i;
769  } else {
770  return 0;
771  }
772 }
773 
774 
775 
776 #define MMHEAP_INITIAL_SIZE 1024
777 
778 template <class T>
780 public:
783  virtual ~UnboundedMinMaxHeap() {};
784 protected:
785  virtual void grow();
786 };
787 
788 template <class T>
790  T *old = this->A;
791  this->maxsize *= 2;
792 
793  assert(this->maxsize > 0);
794 
795  if(old) {
796  HeapIndex n = this->size();
797  this->A = this->allocateHeap(this->maxsize); /* allocate a new array */
798  /* copy over the old values */
799  assert(this->maxsize > n);
800  for(HeapIndex i=0; i<=n; i++) { /* why extra value? -RW */
801  this->A[i] = old[i];
802  }
803  this->freeHeap(old); /* free up old storage */
804  }
805 
806 }
807 
808 
809 #endif
MinMaxHeap(HeapIndex size)
Definition: minmaxheap.h:744
#define MMHEAP_INITIAL_SIZE
Definition: minmaxheap.h:776
virtual void grow()=0
Definition: lz4.c:487
bool full(void) const
Definition: minmaxheap.h:746
HeapIndex maxsize
Definition: minmaxheap.h:90
static T * allocateHeap(HeapIndex n)
Definition: minmaxheap.h:656
BasicMinMaxHeap(HeapIndex size)
Definition: minmaxheap.h:100
void insert(const T &elt)
Definition: minmaxheap.h:521
UnboundedMinMaxHeap(HeapIndex size)
Definition: minmaxheap.h:782
#define XXX
Definition: minmaxheap.h:78
bool max(T &elt) const
Definition: minmaxheap.h:615
void free(void *)
HeapIndex lastindex
Definition: minmaxheap.h:91
#define NULL
Definition: ccmath.h:32
bool extract_min(T &elt)
Definition: minmaxheap.h:540
static void freeHeap(T *)
Definition: minmaxheap.h:669
virtual void grow()
Definition: minmaxheap.h:789
#define assert(condition)
Definition: lz4.c:324
double b
Definition: r_raster.c:39
#define MY_LOG_DEBUG_ID(x)
Definition: minmaxheap.h:81
void print_range() const
Definition: minmaxheap.h:506
virtual ~BasicMinMaxHeap(void)
Definition: minmaxheap.h:111
bool extract_all_min(T &elt)
Definition: minmaxheap.h:557
bool extract_max(T &elt)
Definition: minmaxheap.h:581
bool empty(void) const
Definition: minmaxheap.h:116
virtual ~UnboundedMinMaxHeap()
Definition: minmaxheap.h:783
HeapIndex fill(T *arr, HeapIndex n)
Definition: minmaxheap.h:759
HeapIndex size() const
Definition: minmaxheap.h:117
void destructiveVerify()
Definition: minmaxheap.h:684
void print() const
Definition: minmaxheap.h:496
virtual ~MinMaxHeap()
Definition: minmaxheap.h:745
HeapIndex get_maxsize() const
Definition: minmaxheap.h:747
bool min(T &elt) const
Definition: minmaxheap.h:603
unsigned int HeapIndex
Definition: minmaxheap.h:84
virtual void grow()
Definition: minmaxheap.h:751