1 /* Heaps and priority queues.
2 * See TH Cormen, CE Leiserson, and RL Rivest, _Introduction to Algorithms_, MIT Press, 1999.
3 *
4 * Contents:
5 * 1. The <ESL_HEAP> object: creation, access.
6 * 2. Rest of the API: inserting, extracting values.
7 * 3. Debugging, development.
8 * 4. Internal functions.
9 * 5. Unit tests.
10 * 6. Test driver.
11 */
12 #include "esl_config.h"
13
14 #include <stdlib.h>
15 #include <stdio.h>
16
17 #include "easel.h"
18 #include "esl_heap.h"
19
20 static int heap_grow(ESL_HEAP *hp);
21 static void iheapify (ESL_HEAP *hp, int idx);
22
23 /*****************************************************************
24 * 1. The <ESL_HEAP> object.
25 *****************************************************************/
26
27 /* Function: esl_heap_ICreate()
28 * Synopsis: Create a heap for storing integers.
29 *
30 * Purpose: Create a heap for storing integers. <maxormin> is
31 * <eslHEAP_MIN> or <eslHEAP_MAX>; it states whether
32 * minimum or maximum values are sorted to the top of the
33 * heap.
34 *
35 * Args: maxormin : <eslHEAP_MIN | eslHEAP_MAX>
36 *
37 * Returns: a pointer to the new heap.
38 *
39 * Throws: <NULL> on allocation failure.
40 */
41 ESL_HEAP *
esl_heap_ICreate(int maxormin)42 esl_heap_ICreate(int maxormin)
43 {
44 ESL_HEAP *hp = NULL;
45 int status;
46
47 ESL_DASSERT1(( maxormin == eslHEAP_MIN || maxormin == eslHEAP_MAX));
48
49 ESL_ALLOC(hp, sizeof(ESL_HEAP));
50 hp->idata = NULL;
51 hp->n = 0;
52 hp->maxormin = maxormin;
53
54 ESL_ALLOC(hp->idata, sizeof(int) * eslHEAP_INITALLOC);
55 hp->nalloc = eslHEAP_INITALLOC;
56
57 return hp;
58
59 ERROR:
60 esl_heap_Destroy(hp);
61 return NULL;
62 }
63
64 /* Function: esl_heap_GetCount()
65 * Synopsis: Returns the number of items in the heap.
66 */
67 int
esl_heap_GetCount(ESL_HEAP * hp)68 esl_heap_GetCount(ESL_HEAP *hp)
69 {
70 return hp->n;
71 }
72
73 /* Function: esl_heap_IGetTopVal()
74 * Synopsis: Peeks at and returns the best (topmost) value in the heap.
75 *
76 * Purpose: Peek at the best (topmost) value in heap <hp> and
77 * return it. The heap is unaffected. If the heap is
78 * empty, return 0.
79 */
80 int
esl_heap_IGetTopVal(ESL_HEAP * hp)81 esl_heap_IGetTopVal(ESL_HEAP *hp)
82 {
83 return (hp->n ? hp->idata[0] : 0);
84 }
85
86
87 /* Function: esl_heap_Reuse()
88 * Synopsis: Reuse a heap.
89 *
90 * Purpose: As an alternative to destroy'ing an old heap and
91 * create'ing a new one, empty this heap and reinitialize
92 * it, as if it is a freshly created heap of the same
93 * data type and same <maxormin>.
94 *
95 * Returns: <eslOK>
96 */
97 int
esl_heap_Reuse(ESL_HEAP * hp)98 esl_heap_Reuse(ESL_HEAP *hp)
99 {
100 hp->n = 0;
101 return eslOK;
102 }
103
104
105 /* Function: esl_heap_Destroy()
106 * Synopsis: Free a heap.
107 *
108 * Purpose: Destroys heap <hp>, of any data type.
109 *
110 * Returns: (void)
111 *
112 * Throws: (no abnormal error conditions)
113 */
114 void
esl_heap_Destroy(ESL_HEAP * hp)115 esl_heap_Destroy(ESL_HEAP *hp)
116 {
117 if (hp)
118 {
119 if (hp->idata) free(hp->idata);
120 free (hp);
121 }
122 }
123
124 /*****************************************************************
125 * 2. Rest of the API: inserting, extracting values
126 *****************************************************************/
127
128 /* Function: esl_heap_IInsert()
129 * Synopsis: Insert a value into the heap.
130 *
131 * Purpose: Insert value <val> into heap <hp>.
132 *
133 * Returns: <eslOK> on success.
134 *
135 * Throws: <eslEMEM> on allocation failure.
136 */
137 int
esl_heap_IInsert(ESL_HEAP * hp,int val)138 esl_heap_IInsert(ESL_HEAP *hp, int val)
139 {
140 int idx, parentidx;
141 int status;
142
143 if (hp->n == hp->nalloc && (status = heap_grow(hp)) != eslOK) return status;
144
145 hp->n++;
146 idx = hp->n - 1;
147 while (idx > 0 && (hp->maxormin == eslHEAP_MIN ? hp->idata[ESL_HEAP_PARENT(idx)] > val : hp->idata[ESL_HEAP_PARENT(idx)] < val))
148 {
149 parentidx = ESL_HEAP_PARENT(idx);
150 hp->idata[idx] = hp->idata[parentidx];
151 idx = parentidx;
152 }
153 hp->idata[idx] = val;
154 return eslOK;
155 }
156
157
158 /* Function: esl_heap_IExtractTop()
159 * Synopsis: Extract the top value from the heap.
160 *
161 * Purpose: Extract the best (topmost) value from heap <hp>.
162 * Delete it from the heap. Return it in <*opt_val>.
163 *
164 * To simply delete the topmost value (without retrieving
165 * its value), pass <NULL> for <opt_val>.
166
167 * If the heap is empty, return <eslEOD>, and
168 * <*opt_val> is 0.
169 *
170 * Returns: <eslOK> on success, and <*opt_val> is the extracted
171 * topmost value.
172 *
173 * <eslEOD> if the heap is empty; <*opt_val> is 0.
174 *
175 * Throws: (no abnormal error conditions)
176 */
177 int
esl_heap_IExtractTop(ESL_HEAP * hp,int * opt_val)178 esl_heap_IExtractTop(ESL_HEAP *hp, int *opt_val)
179 {
180 int bestval;
181
182 if (hp->n == 0) { *opt_val = 0; return eslEOD; }
183
184 bestval = hp->idata[0];
185
186 hp->idata[0] = hp->idata[hp->n-1];
187 hp->n--;
188 iheapify(hp, 0);
189
190 if (opt_val) *opt_val = bestval;
191 return eslOK;
192 }
193
194
195 /*****************************************************************
196 * 3. Debugging, development.
197 *****************************************************************/
198
199 int
esl_heap_Validate(ESL_HEAP * hp,char * errbuf)200 esl_heap_Validate(ESL_HEAP *hp, char *errbuf)
201 {
202 int idx, lidx, ridx;
203
204 for (idx = 0; idx < hp->n; idx++)
205 {
206 lidx = ESL_HEAP_LEFT(idx);
207 ridx = lidx+1;
208 if (lidx < hp->n && ( hp->maxormin == eslHEAP_MIN ? hp->idata[lidx] < hp->idata[idx] : hp->idata[lidx] > hp->idata[idx] ))
209 ESL_FAIL(eslFAIL, errbuf, "at %d (value %d): left child %d (value %d) is better", idx, hp->idata[idx], lidx, hp->idata[lidx]);
210 if (ridx < hp->n && ( hp->maxormin == eslHEAP_MIN ? hp->idata[ridx] < hp->idata[idx] : hp->idata[ridx] > hp->idata[idx] ))
211 ESL_FAIL(eslFAIL, errbuf, "at %d (value %d): right child %d (value %d) is better", idx, hp->idata[idx], ridx, hp->idata[ridx]);
212 }
213 return eslOK;
214 }
215
216 /*****************************************************************
217 * 4. Internal functions
218 *****************************************************************/
219
220 static int
heap_grow(ESL_HEAP * hp)221 heap_grow(ESL_HEAP *hp)
222 {
223 int status;
224
225 if (hp->idata) {
226 ESL_REALLOC(hp->idata, sizeof(int) * (hp->nalloc*2));
227 hp->nalloc += hp->nalloc;
228 }
229 return eslOK;
230
231 ERROR:
232 return eslEMEM;
233 }
234
235
236 static void
iheapify(ESL_HEAP * hp,int idx)237 iheapify(ESL_HEAP *hp, int idx)
238 {
239 int bestidx = idx;
240 int leftidx, rightidx;
241
242 while (1) /* while loop avoids recursive heapify call */
243 {
244 leftidx = ESL_HEAP_LEFT(idx);
245 rightidx = leftidx+1;
246 if (leftidx < hp->n && (hp->maxormin == eslHEAP_MIN ? hp->idata[leftidx] < hp->idata[idx] : hp->idata[leftidx] > hp->idata[idx]) ) bestidx = leftidx;
247 if (rightidx < hp->n && (hp->maxormin == eslHEAP_MIN ? hp->idata[rightidx] < hp->idata[bestidx] : hp->idata[rightidx] > hp->idata[bestidx]) ) bestidx = rightidx;
248 if (bestidx == idx) break; /* nothing needed to be changed: either because <idx> satisfies heap property, or because it has no children */
249 ESL_SWAP(hp->idata[idx], hp->idata[bestidx], int);
250 idx = bestidx;
251 }
252 }
253
254
255 /*****************************************************************
256 * 5. Unit tests
257 *****************************************************************/
258 #ifdef eslHEAP_TESTDRIVE
259 #include "esl_random.h"
260
261 static void
utest_sorting(ESL_RANDOMNESS * rng)262 utest_sorting(ESL_RANDOMNESS *rng)
263 {
264 char *msg = "utest_sorting():: unit test failure";
265 ESL_HEAP *hp = NULL;
266 int *val = NULL;
267 int nv = 1 + esl_rnd_Roll(rng, 10000);
268 char errbuf[eslERRBUFSIZE];
269 int i,n2,x;
270
271 /* Create an array of numbers 1..nv in randomized order. */
272 if ( (val = malloc(sizeof(int) * nv)) == NULL) esl_fatal("utest_sorting():: allocation failed");
273 for (i = 0; i < nv; i++) val[i] = i+1;
274 for (n2 = nv; n2 > 1; n2--)
275 { /* a compact Fisher-Yates shuffle. Can't put the Roll() into the ESL_SWAP(), because it's a macro: avoid double evaluation */
276 i = esl_rnd_Roll(rng, n2);
277 ESL_SWAP( val[i], val[n2-1], int);
278 }
279
280 /* Add those numbers (in their randomized order) to a min heap */
281 if ( (hp = esl_heap_ICreate(eslHEAP_MIN)) == NULL) esl_fatal(msg);
282 for (i = 0; i < nv; i++)
283 if (esl_heap_IInsert(hp, val[i]) != eslOK) esl_fatal(msg);
284 if (esl_heap_Validate(hp, errbuf) != eslOK) esl_fatal("utest: heap validation fails: %s", errbuf);
285
286 /* Now if we pull numbers off the heap, they'll come off in sorted order, 1..nv */
287 for (i = 1; i <= nv; i++)
288 {
289 if (esl_heap_IExtractTop(hp, &x) != eslOK) esl_fatal(msg);
290 if (x != i) esl_fatal(msg);
291 if (hp->n != nv-i) esl_fatal(msg);
292 }
293
294 esl_heap_Destroy(hp);
295 free(val);
296 }
297
298 #endif /*eslHEAP_TESTDRIVE*/
299
300 /*****************************************************************
301 * 6. Test driver
302 *****************************************************************/
303 #ifdef eslHEAP_TESTDRIVE
304
305 #include "easel.h"
306 #include "esl_getopts.h"
307
308 static ESL_OPTIONS options[] = {
309 /* name type default env range togs reqs incomp help docgrp */
310 {"-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show help and usage", 0},
311 {"-s", eslARG_INT, "0", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0},
312 { 0,0,0,0,0,0,0,0,0,0},
313 };
314 static char usage[] = "[-options]";
315 static char banner[] = "test driver for ESL_HEAP: heaps and priority queues";
316
317 int
main(int argc,char ** argv)318 main(int argc, char **argv)
319 {
320 ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 0, argc, argv, banner, usage);
321 ESL_RANDOMNESS *rng = esl_randomness_Create(esl_opt_GetInteger(go, "-s"));
322
323 fprintf(stderr, "## %s\n", argv[0]);
324 fprintf(stderr, "# rng seed = %" PRIu32 "\n", esl_randomness_GetSeed(rng));
325
326 utest_sorting(rng);
327
328 fprintf(stderr, "# status = ok\n");
329
330 esl_getopts_Destroy(go);
331 esl_randomness_Destroy(rng);
332 return 0;
333 }
334
335 #endif /*eslHEAP_TESTDRIVE*/
336