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