1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2003-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard street, Cambridge, MA 02139 USA
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 #include "igraph_types.h"
25 #include "igraph_types_internal.h"
26 #include "igraph_memory.h"
27 #include "igraph_error.h"
28 #include "igraph_math.h"
29 #include "config.h"
30
31 #include <assert.h>
32 #include <string.h> /* memcpy & co. */
33 #include <stdlib.h>
34
35 #define PARENT(x) (((x)+1)/2-1)
36 #define LEFTCHILD(x) (((x)+1)*2-1)
37 #define RIGHTCHILD(x) (((x)+1)*2)
38
39 /**
40 * \ingroup indheap
41 * \brief Initializes an indexed heap (constructor).
42 *
43 * @return Error code:
44 * - <b>IGRAPH_ENOMEM</b>: out of memory
45 */
46
igraph_indheap_init(igraph_indheap_t * h,long int alloc_size)47 int igraph_indheap_init (igraph_indheap_t* h, long int alloc_size) {
48 if (alloc_size <= 0 ) {
49 alloc_size = 1;
50 }
51 h->stor_begin = igraph_Calloc(alloc_size, igraph_real_t);
52 if (h->stor_begin == 0) {
53 h->index_begin = 0;
54 IGRAPH_ERROR("indheap init failed", IGRAPH_ENOMEM);
55 }
56 h->index_begin = igraph_Calloc(alloc_size, long int);
57 if (h->index_begin == 0) {
58 igraph_Free(h->stor_begin);
59 h->stor_begin = 0;
60 IGRAPH_ERROR("indheap init failed", IGRAPH_ENOMEM);
61 }
62
63 h->stor_end = h->stor_begin + alloc_size;
64 h->end = h->stor_begin;
65 h->destroy = 1;
66
67 return 0;
68 }
69
igraph_indheap_clear(igraph_indheap_t * h)70 int igraph_indheap_clear(igraph_indheap_t *h) {
71 h->end = h->stor_begin;
72 return 0;
73 }
74
75 /**
76 * \ingroup indheap
77 * \brief Initializes and build an indexed heap from a C array (constructor).
78 *
79 * @return Error code:
80 * - <b>IGRAPH_ENOMEM</b>: out of memory
81 */
82
igraph_indheap_init_array(igraph_indheap_t * h,igraph_real_t * data,long int len)83 int igraph_indheap_init_array (igraph_indheap_t *h, igraph_real_t* data, long int len) {
84 long int i;
85
86 h->stor_begin = igraph_Calloc(len, igraph_real_t);
87 if (h->stor_begin == 0) {
88 h->index_begin = 0;
89 IGRAPH_ERROR("indheap init from array failed", IGRAPH_ENOMEM);
90 }
91 h->index_begin = igraph_Calloc(len, long int);
92 if (h->index_begin == 0) {
93 igraph_Free(h->stor_begin);
94 h->stor_begin = 0;
95 IGRAPH_ERROR("indheap init from array failed", IGRAPH_ENOMEM);
96 }
97 h->stor_end = h->stor_begin + len;
98 h->end = h->stor_end;
99 h->destroy = 1;
100
101 memcpy(h->stor_begin, data, (size_t) len * sizeof(igraph_real_t));
102 for (i = 0; i < len; i++) {
103 h->index_begin[i] = i + 1;
104 }
105
106 igraph_indheap_i_build (h, 0);
107
108 return 0;
109 }
110
111 /**
112 * \ingroup indheap
113 * \brief Destroys an initialized indexed heap.
114 */
115
igraph_indheap_destroy(igraph_indheap_t * h)116 void igraph_indheap_destroy (igraph_indheap_t* h) {
117 assert(h != 0);
118 if (h->destroy) {
119 if (h->stor_begin != 0) {
120 igraph_Free(h->stor_begin);
121 h->stor_begin = 0;
122 }
123 if (h->index_begin != 0) {
124 igraph_Free(h->index_begin);
125 h->index_begin = 0;
126 }
127 }
128 }
129
130 /**
131 * \ingroup indheap
132 * \brief Checks whether a heap is empty.
133 */
134
igraph_indheap_empty(igraph_indheap_t * h)135 igraph_bool_t igraph_indheap_empty (igraph_indheap_t* h) {
136 assert(h != 0);
137 assert(h->stor_begin != 0);
138 return h->stor_begin == h->end;
139 }
140
141 /**
142 * \ingroup indheap
143 * \brief Adds an element to an indexed heap.
144 */
145
igraph_indheap_push(igraph_indheap_t * h,igraph_real_t elem)146 int igraph_indheap_push (igraph_indheap_t* h, igraph_real_t elem) {
147 assert(h != 0);
148 assert(h->stor_begin != 0);
149
150 /* full, allocate more storage */
151 if (h->stor_end == h->end) {
152 long int new_size = igraph_indheap_size(h) * 2;
153 if (new_size == 0) {
154 new_size = 1;
155 }
156 IGRAPH_CHECK(igraph_indheap_reserve(h, new_size));
157 }
158
159 *(h->end) = elem;
160 h->end += 1;
161 *(h->index_begin + igraph_indheap_size(h) - 1) = igraph_indheap_size(h) - 1;
162
163 /* maintain indheap */
164 igraph_indheap_i_shift_up(h, igraph_indheap_size(h) - 1);
165
166 return 0;
167 }
168
169 /**
170 * \ingroup indheap
171 * \brief Adds an element to an indexed heap with a given index.
172 */
173
igraph_indheap_push_with_index(igraph_indheap_t * h,long int idx,igraph_real_t elem)174 int igraph_indheap_push_with_index(igraph_indheap_t* h, long int idx, igraph_real_t elem) {
175 assert(h != 0);
176 assert(h->stor_begin != 0);
177
178 /* full, allocate more storage */
179 if (h->stor_end == h->end) {
180 long int new_size = igraph_indheap_size(h) * 2;
181 if (new_size == 0) {
182 new_size = 1;
183 }
184 IGRAPH_CHECK(igraph_indheap_reserve(h, new_size));
185 }
186
187 *(h->end) = elem;
188 h->end += 1;
189 *(h->index_begin + igraph_indheap_size(h) - 1) = idx;
190
191 /* maintain indheap */
192 igraph_indheap_i_shift_up(h, igraph_indheap_size(h) - 1);
193
194 return 0;
195 }
196
197 /**
198 * \ingroup indheap
199 * \brief Modifies an element in an indexed heap.
200 */
201
igraph_indheap_modify(igraph_indheap_t * h,long int idx,igraph_real_t elem)202 int igraph_indheap_modify(igraph_indheap_t* h, long int idx, igraph_real_t elem) {
203 long int i, n;
204
205 assert(h != 0);
206 assert(h->stor_begin != 0);
207
208 n = igraph_indheap_size(h);
209 for (i = 0; i < n; i++)
210 if (h->index_begin[i] == idx) {
211 h->stor_begin[i] = elem;
212 break;
213 }
214
215 if (i == n) {
216 return 0;
217 }
218
219 /* maintain indheap */
220 igraph_indheap_i_build(h, 0);
221
222 return 0;
223 }
224
225 /**
226 * \ingroup indheap
227 * \brief Returns the largest element in an indexed heap.
228 */
229
igraph_indheap_max(igraph_indheap_t * h)230 igraph_real_t igraph_indheap_max (igraph_indheap_t* h) {
231 assert(h != NULL);
232 assert(h->stor_begin != NULL);
233 assert(h->stor_begin != h->end);
234
235 return h->stor_begin[0];
236 }
237
238 /**
239 * \ingroup indheap
240 * \brief Removes the largest element from an indexed heap.
241 */
242
igraph_indheap_delete_max(igraph_indheap_t * h)243 igraph_real_t igraph_indheap_delete_max(igraph_indheap_t* h) {
244 igraph_real_t tmp;
245
246 assert(h != NULL);
247 assert(h->stor_begin != NULL);
248
249 tmp = h->stor_begin[0];
250 igraph_indheap_i_switch(h, 0, igraph_indheap_size(h) - 1);
251 h->end -= 1;
252 igraph_indheap_i_sink(h, 0);
253
254 return tmp;
255 }
256
257 /**
258 * \ingroup indheap
259 * \brief Gives the number of elements in an indexed heap.
260 */
261
igraph_indheap_size(igraph_indheap_t * h)262 long int igraph_indheap_size (igraph_indheap_t* h) {
263 assert(h != 0);
264 assert(h->stor_begin != 0);
265 return h->end - h->stor_begin;
266 }
267
268 /**
269 * \ingroup indheap
270 * \brief Reserves more memory for an indexed heap.
271 *
272 * @return Error code:
273 * - <b>IGRAPH_ENOMEM</b>: out of memory
274 */
275
igraph_indheap_reserve(igraph_indheap_t * h,long int size)276 int igraph_indheap_reserve (igraph_indheap_t* h, long int size) {
277 long int actual_size = igraph_indheap_size(h);
278 igraph_real_t *tmp1;
279 long int *tmp2;
280 assert(h != 0);
281 assert(h->stor_begin != 0);
282
283 if (size <= actual_size) {
284 return 0;
285 }
286
287 tmp1 = igraph_Calloc(size, igraph_real_t);
288 if (tmp1 == 0) {
289 IGRAPH_ERROR("indheap reserve failed", IGRAPH_ENOMEM);
290 }
291 IGRAPH_FINALLY(igraph_free, tmp1);
292 tmp2 = igraph_Calloc(size, long int);
293 if (tmp2 == 0) {
294 IGRAPH_ERROR("indheap reserve failed", IGRAPH_ENOMEM);
295 }
296 IGRAPH_FINALLY(igraph_free, tmp2);
297 memcpy(tmp1, h->stor_begin, (size_t) actual_size * sizeof(igraph_real_t));
298 memcpy(tmp2, h->index_begin, (size_t) actual_size * sizeof(long int));
299 igraph_Free(h->stor_begin);
300 igraph_Free(h->index_begin);
301
302 h->stor_begin = tmp1;
303 h->index_begin = tmp2;
304 h->stor_end = h->stor_begin + size;
305 h->end = h->stor_begin + actual_size;
306
307 IGRAPH_FINALLY_CLEAN(2);
308 return 0;
309 }
310
311 /**
312 * \ingroup indheap
313 * \brief Returns the index of the largest element in an indexed heap.
314 */
315
igraph_indheap_max_index(igraph_indheap_t * h)316 long int igraph_indheap_max_index(igraph_indheap_t *h) {
317 assert(h != 0);
318 assert(h->stor_begin != 0);
319 return h->index_begin[0];
320 }
321
322 /**
323 * \ingroup indheap
324 * \brief Builds an indexed heap, this function should not be called
325 * directly.
326 */
327
igraph_indheap_i_build(igraph_indheap_t * h,long int head)328 void igraph_indheap_i_build(igraph_indheap_t* h, long int head) {
329
330 long int size = igraph_indheap_size(h);
331 if (RIGHTCHILD(head) < size) {
332 /* both subtrees */
333 igraph_indheap_i_build(h, LEFTCHILD(head) );
334 igraph_indheap_i_build(h, RIGHTCHILD(head));
335 igraph_indheap_i_sink(h, head);
336 } else if (LEFTCHILD(head) < size) {
337 /* only left */
338 igraph_indheap_i_build(h, LEFTCHILD(head));
339 igraph_indheap_i_sink(h, head);
340 } else {
341 /* none */
342 }
343 }
344
345 /**
346 * \ingroup indheap
347 * \brief Moves an element up in the heap, don't call this function
348 * directly.
349 */
350
igraph_indheap_i_shift_up(igraph_indheap_t * h,long int elem)351 void igraph_indheap_i_shift_up(igraph_indheap_t *h, long int elem) {
352
353 if (elem == 0 || h->stor_begin[elem] < h->stor_begin[PARENT(elem)]) {
354 /* at the top */
355 } else {
356 igraph_indheap_i_switch(h, elem, PARENT(elem));
357 igraph_indheap_i_shift_up(h, PARENT(elem));
358 }
359 }
360
361 /**
362 * \ingroup indheap
363 * \brief Moves an element down in the heap, don't call this function
364 * directly.
365 */
366
igraph_indheap_i_sink(igraph_indheap_t * h,long int head)367 void igraph_indheap_i_sink(igraph_indheap_t* h, long int head) {
368
369 long int size = igraph_indheap_size(h);
370 if (LEFTCHILD(head) >= size) {
371 /* no subtrees */
372 } else if (RIGHTCHILD(head) == size ||
373 h->stor_begin[LEFTCHILD(head)] >= h->stor_begin[RIGHTCHILD(head)]) {
374 /* sink to the left if needed */
375 if (h->stor_begin[head] < h->stor_begin[LEFTCHILD(head)]) {
376 igraph_indheap_i_switch(h, head, LEFTCHILD(head));
377 igraph_indheap_i_sink(h, LEFTCHILD(head));
378 }
379 } else {
380 /* sink to the right */
381 if (h->stor_begin[head] < h->stor_begin[RIGHTCHILD(head)]) {
382 igraph_indheap_i_switch(h, head, RIGHTCHILD(head));
383 igraph_indheap_i_sink(h, RIGHTCHILD(head));
384 }
385 }
386 }
387
388 /**
389 * \ingroup indheap
390 * \brief Switches two elements in a heap, don't call this function
391 * directly.
392 */
393
igraph_indheap_i_switch(igraph_indheap_t * h,long int e1,long int e2)394 void igraph_indheap_i_switch(igraph_indheap_t* h, long int e1, long int e2) {
395 if (e1 != e2) {
396 igraph_real_t tmp = h->stor_begin[e1];
397 h->stor_begin[e1] = h->stor_begin[e2];
398 h->stor_begin[e2] = tmp;
399
400 tmp = h->index_begin[e1];
401 h->index_begin[e1] = h->index_begin[e2];
402 h->index_begin[e2] = (long int) tmp;
403 }
404 }
405
406
407 /**
408 * \ingroup doubleindheap
409 * \brief Initializes an empty doubly indexed heap object (constructor).
410 *
411 * @return Error code:
412 * - <b>IGRAPH_ENOMEM</b>: out of memory
413 */
414
igraph_d_indheap_init(igraph_d_indheap_t * h,long int alloc_size)415 int igraph_d_indheap_init (igraph_d_indheap_t* h, long int alloc_size) {
416 if (alloc_size <= 0 ) {
417 alloc_size = 1;
418 }
419 h->stor_begin = igraph_Calloc(alloc_size, igraph_real_t);
420 if (h->stor_begin == 0) {
421 h->index_begin = 0;
422 h->index2_begin = 0;
423 IGRAPH_ERROR("d_indheap init failed", IGRAPH_ENOMEM);
424 }
425 h->stor_end = h->stor_begin + alloc_size;
426 h->end = h->stor_begin;
427 h->destroy = 1;
428 h->index_begin = igraph_Calloc(alloc_size, long int);
429 if (h->index_begin == 0) {
430 igraph_Free(h->stor_begin);
431 h->stor_begin = 0;
432 h->index2_begin = 0;
433 IGRAPH_ERROR("d_indheap init failed", IGRAPH_ENOMEM);
434 }
435 h->index2_begin = igraph_Calloc(alloc_size, long int);
436 if (h->index2_begin == 0) {
437 igraph_Free(h->stor_begin);
438 igraph_Free(h->index_begin);
439 h->stor_begin = 0;
440 h->index_begin = 0;
441 IGRAPH_ERROR("d_indheap init failed", IGRAPH_ENOMEM);
442 }
443
444 return 0;
445 }
446
447 /**
448 * \ingroup doubleindheap
449 * \brief Destroys an initialized doubly indexed heap object.
450 */
451
igraph_d_indheap_destroy(igraph_d_indheap_t * h)452 void igraph_d_indheap_destroy (igraph_d_indheap_t* h) {
453 assert(h != 0);
454 if (h->destroy) {
455 if (h->stor_begin != 0) {
456 igraph_Free(h->stor_begin);
457 h->stor_begin = 0;
458 }
459 if (h->index_begin != 0) {
460 igraph_Free(h->index_begin);
461 h->index_begin = 0;
462 }
463 if (h->index2_begin != 0) {
464 igraph_Free(h->index2_begin);
465 h->index2_begin = 0;
466 }
467 }
468 }
469
470 /**
471 * \ingroup doubleindheap
472 * \brief Decides whether a heap is empty.
473 */
474
igraph_d_indheap_empty(igraph_d_indheap_t * h)475 igraph_bool_t igraph_d_indheap_empty (igraph_d_indheap_t* h) {
476 assert(h != 0);
477 assert(h->stor_begin != 0);
478 return h->stor_begin == h->end;
479 }
480
481 /**
482 * \ingroup doubleindheap
483 * \brief Adds an element to the heap.
484 */
485
igraph_d_indheap_push(igraph_d_indheap_t * h,igraph_real_t elem,long int idx,long int idx2)486 int igraph_d_indheap_push (igraph_d_indheap_t* h, igraph_real_t elem,
487 long int idx, long int idx2) {
488 assert(h != 0);
489 assert(h->stor_begin != 0);
490
491 /* full, allocate more storage */
492 if (h->stor_end == h->end) {
493 long int new_size = igraph_d_indheap_size(h) * 2;
494 if (new_size == 0) {
495 new_size = 1;
496 }
497 IGRAPH_CHECK(igraph_d_indheap_reserve(h, new_size));
498 }
499
500 *(h->end) = elem;
501 h->end += 1;
502 *(h->index_begin + igraph_d_indheap_size(h) - 1) = idx ;
503 *(h->index2_begin + igraph_d_indheap_size(h) - 1) = idx2 ;
504
505 /* maintain d_indheap */
506 igraph_d_indheap_i_shift_up(h, igraph_d_indheap_size(h) - 1);
507
508 return 0;
509 }
510
511 /**
512 * \ingroup doubleindheap
513 * \brief Returns the largest element in the heap.
514 */
515
igraph_d_indheap_max(igraph_d_indheap_t * h)516 igraph_real_t igraph_d_indheap_max (igraph_d_indheap_t* h) {
517 assert(h != NULL);
518 assert(h->stor_begin != NULL);
519 assert(h->stor_begin != h->end);
520
521 return h->stor_begin[0];
522 }
523
524 /**
525 * \ingroup doubleindheap
526 * \brief Removes the largest element from the heap.
527 */
528
igraph_d_indheap_delete_max(igraph_d_indheap_t * h)529 igraph_real_t igraph_d_indheap_delete_max(igraph_d_indheap_t* h) {
530 igraph_real_t tmp;
531
532 assert(h != NULL);
533 assert(h->stor_begin != NULL);
534
535 tmp = h->stor_begin[0];
536 igraph_d_indheap_i_switch(h, 0, igraph_d_indheap_size(h) - 1);
537 h->end -= 1;
538 igraph_d_indheap_i_sink(h, 0);
539
540 return tmp;
541 }
542
543 /**
544 * \ingroup doubleindheap
545 * \brief Gives the number of elements in the heap.
546 */
547
igraph_d_indheap_size(igraph_d_indheap_t * h)548 long int igraph_d_indheap_size (igraph_d_indheap_t* h) {
549 assert(h != 0);
550 assert(h->stor_begin != 0);
551 return h->end - h->stor_begin;
552 }
553
554 /**
555 * \ingroup doubleindheap
556 * \brief Allocates memory for a heap.
557 *
558 * @return Error code:
559 * - <b>IGRAPH_ENOMEM</b>: out of memory
560 */
561
igraph_d_indheap_reserve(igraph_d_indheap_t * h,long int size)562 int igraph_d_indheap_reserve (igraph_d_indheap_t* h, long int size) {
563 long int actual_size = igraph_d_indheap_size(h);
564 igraph_real_t *tmp1;
565 long int *tmp2, *tmp3;
566 assert(h != 0);
567 assert(h->stor_begin != 0);
568
569 if (size <= actual_size) {
570 return 0;
571 }
572
573 tmp1 = igraph_Calloc(size, igraph_real_t);
574 if (tmp1 == 0) {
575 IGRAPH_ERROR("d_indheap reserve failed", IGRAPH_ENOMEM);
576 }
577 IGRAPH_FINALLY(igraph_free, tmp1);
578 tmp2 = igraph_Calloc(size, long int);
579 if (tmp2 == 0) {
580 IGRAPH_ERROR("d_indheap reserve failed", IGRAPH_ENOMEM);
581 }
582 IGRAPH_FINALLY(igraph_free, tmp2);
583 tmp3 = igraph_Calloc(size, long int);
584 if (tmp3 == 0) {
585 IGRAPH_ERROR("d_indheap reserve failed", IGRAPH_ENOMEM);
586 }
587 IGRAPH_FINALLY(igraph_free, tmp3);
588
589 memcpy(tmp1, h->stor_begin, (size_t) actual_size * sizeof(igraph_real_t));
590 memcpy(tmp2, h->index_begin, (size_t) actual_size * sizeof(long int));
591 memcpy(tmp3, h->index2_begin, (size_t) actual_size * sizeof(long int));
592 igraph_Free(h->stor_begin);
593 igraph_Free(h->index_begin);
594 igraph_Free(h->index2_begin);
595
596 h->stor_begin = tmp1;
597 h->stor_end = h->stor_begin + size;
598 h->end = h->stor_begin + actual_size;
599 h->index_begin = tmp2;
600 h->index2_begin = tmp3;
601
602 IGRAPH_FINALLY_CLEAN(3);
603 return 0;
604 }
605
606 /**
607 * \ingroup doubleindheap
608 * \brief Gives the indices of the maximal element in the heap.
609 */
610
igraph_d_indheap_max_index(igraph_d_indheap_t * h,long int * idx,long int * idx2)611 void igraph_d_indheap_max_index(igraph_d_indheap_t *h, long int *idx, long int *idx2) {
612 assert(h != 0);
613 assert(h->stor_begin != 0);
614 (*idx) = h->index_begin[0];
615 (*idx2) = h->index2_begin[0];
616 }
617
618 /**
619 * \ingroup doubleindheap
620 * \brief Builds the heap, don't call it directly.
621 */
622
igraph_d_indheap_i_build(igraph_d_indheap_t * h,long int head)623 void igraph_d_indheap_i_build(igraph_d_indheap_t* h, long int head) {
624
625 long int size = igraph_d_indheap_size(h);
626 if (RIGHTCHILD(head) < size) {
627 /* both subtrees */
628 igraph_d_indheap_i_build(h, LEFTCHILD(head) );
629 igraph_d_indheap_i_build(h, RIGHTCHILD(head));
630 igraph_d_indheap_i_sink(h, head);
631 } else if (LEFTCHILD(head) < size) {
632 /* only left */
633 igraph_d_indheap_i_build(h, LEFTCHILD(head));
634 igraph_d_indheap_i_sink(h, head);
635 } else {
636 /* none */
637 }
638 }
639
640 /**
641 * \ingroup doubleindheap
642 * \brief Moves an element up in the heap, don't call it directly.
643 */
644
igraph_d_indheap_i_shift_up(igraph_d_indheap_t * h,long int elem)645 void igraph_d_indheap_i_shift_up(igraph_d_indheap_t *h, long int elem) {
646
647 if (elem == 0 || h->stor_begin[elem] < h->stor_begin[PARENT(elem)]) {
648 /* at the top */
649 } else {
650 igraph_d_indheap_i_switch(h, elem, PARENT(elem));
651 igraph_d_indheap_i_shift_up(h, PARENT(elem));
652 }
653 }
654
655 /**
656 * \ingroup doubleindheap
657 * \brief Moves an element down in the heap, don't call it directly.
658 */
659
igraph_d_indheap_i_sink(igraph_d_indheap_t * h,long int head)660 void igraph_d_indheap_i_sink(igraph_d_indheap_t* h, long int head) {
661
662 long int size = igraph_d_indheap_size(h);
663 if (LEFTCHILD(head) >= size) {
664 /* no subtrees */
665 } else if (RIGHTCHILD(head) == size ||
666 h->stor_begin[LEFTCHILD(head)] >= h->stor_begin[RIGHTCHILD(head)]) {
667 /* sink to the left if needed */
668 if (h->stor_begin[head] < h->stor_begin[LEFTCHILD(head)]) {
669 igraph_d_indheap_i_switch(h, head, LEFTCHILD(head));
670 igraph_d_indheap_i_sink(h, LEFTCHILD(head));
671 }
672 } else {
673 /* sink to the right */
674 if (h->stor_begin[head] < h->stor_begin[RIGHTCHILD(head)]) {
675 igraph_d_indheap_i_switch(h, head, RIGHTCHILD(head));
676 igraph_d_indheap_i_sink(h, RIGHTCHILD(head));
677 }
678 }
679 }
680
681 /**
682 * \ingroup doubleindheap
683 * \brief Switches two elements in the heap, don't call it directly.
684 */
685
igraph_d_indheap_i_switch(igraph_d_indheap_t * h,long int e1,long int e2)686 void igraph_d_indheap_i_switch(igraph_d_indheap_t* h, long int e1, long int e2) {
687 if (e1 != e2) {
688 long int tmpi;
689 igraph_real_t tmp = h->stor_begin[e1];
690 h->stor_begin[e1] = h->stor_begin[e2];
691 h->stor_begin[e2] = tmp;
692
693 tmpi = h->index_begin[e1];
694 h->index_begin[e1] = h->index_begin[e2];
695 h->index_begin[e2] = tmpi;
696
697 tmpi = h->index2_begin[e1];
698 h->index2_begin[e1] = h->index2_begin[e2];
699 h->index2_begin[e2] = tmpi;
700 }
701 }
702
703 /*************************************************/
704
705 #undef PARENT
706 #undef LEFTCHILD
707 #undef RIGHTCHILD
708 #define PARENT(x) ((x)/2)
709 #define LEFTCHILD(x) ((x)*2+1)
710 #define RIGHTCHILD(x) ((x)*2)
711 #define INACTIVE IGRAPH_INFINITY
712 #define UNDEFINED 0.0
713 #define INDEXINC 1
714
igraph_i_cutheap_switch(igraph_i_cutheap_t * ch,long int hidx1,long int hidx2)715 void igraph_i_cutheap_switch(igraph_i_cutheap_t *ch,
716 long int hidx1, long int hidx2) {
717 if (hidx1 != hidx2) {
718 long int idx1 = (long int) VECTOR(ch->index)[hidx1];
719 long int idx2 = (long int) VECTOR(ch->index)[hidx2];
720
721 igraph_real_t tmp = VECTOR(ch->heap)[hidx1];
722 VECTOR(ch->heap)[hidx1] = VECTOR(ch->heap)[hidx2];
723 VECTOR(ch->heap)[hidx2] = tmp;
724
725 VECTOR(ch->index)[hidx1] = idx2;
726 VECTOR(ch->index)[hidx2] = idx1;
727
728 VECTOR(ch->hptr)[idx1] = hidx2 + INDEXINC;
729 VECTOR(ch->hptr)[idx2] = hidx1 + INDEXINC;
730 }
731 }
732
igraph_i_cutheap_sink(igraph_i_cutheap_t * ch,long int hidx)733 void igraph_i_cutheap_sink(igraph_i_cutheap_t *ch, long int hidx) {
734 long int size = igraph_vector_size(&ch->heap);
735 if (LEFTCHILD(hidx) >= size) {
736 /* leaf node */
737 } else if (RIGHTCHILD(hidx) == size ||
738 VECTOR(ch->heap)[LEFTCHILD(hidx)] >=
739 VECTOR(ch->heap)[RIGHTCHILD(hidx)]) {
740 /* sink to the left if needed */
741 if (VECTOR(ch->heap)[hidx] < VECTOR(ch->heap)[LEFTCHILD(hidx)]) {
742 igraph_i_cutheap_switch(ch, hidx, LEFTCHILD(hidx));
743 igraph_i_cutheap_sink(ch, LEFTCHILD(hidx));
744 }
745 } else {
746 /* sink to the right */
747 if (VECTOR(ch->heap)[hidx] < VECTOR(ch->heap)[RIGHTCHILD(hidx)]) {
748 igraph_i_cutheap_switch(ch, hidx, RIGHTCHILD(hidx));
749 igraph_i_cutheap_sink(ch, RIGHTCHILD(hidx));
750 }
751 }
752 }
753
igraph_i_cutheap_shift_up(igraph_i_cutheap_t * ch,long int hidx)754 void igraph_i_cutheap_shift_up(igraph_i_cutheap_t *ch, long int hidx) {
755 if (hidx == 0 || VECTOR(ch->heap)[hidx] < VECTOR(ch->heap)[PARENT(hidx)]) {
756 /* at the top */
757 } else {
758 igraph_i_cutheap_switch(ch, hidx, PARENT(hidx));
759 igraph_i_cutheap_shift_up(ch, PARENT(hidx));
760 }
761 }
762
igraph_i_cutheap_init(igraph_i_cutheap_t * ch,igraph_integer_t nodes)763 int igraph_i_cutheap_init(igraph_i_cutheap_t *ch, igraph_integer_t nodes) {
764 ch->dnodes = nodes;
765 IGRAPH_VECTOR_INIT_FINALLY(&ch->heap, nodes); /* all zero */
766 IGRAPH_CHECK(igraph_vector_init_seq(&ch->index, 0, nodes - 1));
767 IGRAPH_FINALLY(igraph_vector_destroy, &ch->index);
768 IGRAPH_CHECK(igraph_vector_init_seq(&ch->hptr, INDEXINC, nodes + INDEXINC - 1));
769 IGRAPH_FINALLY_CLEAN(2);
770 return 0;
771 }
772
igraph_i_cutheap_destroy(igraph_i_cutheap_t * ch)773 void igraph_i_cutheap_destroy(igraph_i_cutheap_t *ch) {
774 igraph_vector_destroy(&ch->hptr);
775 igraph_vector_destroy(&ch->index);
776 igraph_vector_destroy(&ch->heap);
777 }
778
igraph_i_cutheap_empty(igraph_i_cutheap_t * ch)779 igraph_bool_t igraph_i_cutheap_empty(igraph_i_cutheap_t *ch) {
780 return igraph_vector_empty(&ch->heap);
781 }
782
783 /* Number of active vertices */
784
igraph_i_cutheap_active_size(igraph_i_cutheap_t * ch)785 igraph_integer_t igraph_i_cutheap_active_size(igraph_i_cutheap_t *ch) {
786 return (igraph_integer_t) igraph_vector_size(&ch->heap);
787 }
788
789 /* Number of all (defined) vertices */
790
igraph_i_cutheap_size(igraph_i_cutheap_t * ch)791 igraph_integer_t igraph_i_cutheap_size(igraph_i_cutheap_t *ch) {
792 return (igraph_integer_t) (ch->dnodes);
793 }
794
igraph_i_cutheap_maxvalue(igraph_i_cutheap_t * ch)795 igraph_real_t igraph_i_cutheap_maxvalue(igraph_i_cutheap_t *ch) {
796 return VECTOR(ch->heap)[0];
797 }
798
igraph_i_cutheap_popmax(igraph_i_cutheap_t * ch)799 igraph_integer_t igraph_i_cutheap_popmax(igraph_i_cutheap_t *ch) {
800 long int size = igraph_vector_size(&ch->heap);
801 igraph_integer_t maxindex = (igraph_integer_t) VECTOR(ch->index)[0];
802 /* put the last element to the top */
803 igraph_i_cutheap_switch(ch, 0, size - 1);
804 /* remove the last element */
805 VECTOR(ch->hptr)[(long int) igraph_vector_tail(&ch->index)] = INACTIVE;
806 igraph_vector_pop_back(&ch->heap);
807 igraph_vector_pop_back(&ch->index);
808 igraph_i_cutheap_sink(ch, 0);
809
810 return maxindex;
811 }
812
813 /* Update the value of an active vertex, if not active it will be ignored */
814
igraph_i_cutheap_update(igraph_i_cutheap_t * ch,igraph_integer_t index,igraph_real_t add)815 int igraph_i_cutheap_update(igraph_i_cutheap_t *ch, igraph_integer_t index,
816 igraph_real_t add) {
817 igraph_real_t hidx = VECTOR(ch->hptr)[(long int)index];
818 if (hidx != INACTIVE && hidx != UNDEFINED) {
819 long int hidx2 = (long int) (hidx - INDEXINC);
820 /* printf("updating vertex %li, heap index %li\n", (long int) index, hidx2); */
821 VECTOR(ch->heap)[hidx2] += add;
822 igraph_i_cutheap_sink(ch, hidx2);
823 igraph_i_cutheap_shift_up(ch, hidx2);
824 }
825 return 0;
826 }
827
828 /* Reset the value of all vertices to zero and make them active */
829
igraph_i_cutheap_reset_undefine(igraph_i_cutheap_t * ch,long int vertex)830 int igraph_i_cutheap_reset_undefine(igraph_i_cutheap_t *ch, long int vertex) {
831 long int i, j, n = igraph_vector_size(&ch->hptr);
832 /* undefine */
833 VECTOR(ch->hptr)[vertex] = UNDEFINED;
834 ch->dnodes -= 1;
835
836 IGRAPH_CHECK(igraph_vector_resize(&ch->heap, ch->dnodes));
837 igraph_vector_null(&ch->heap);
838
839 IGRAPH_CHECK(igraph_vector_resize(&ch->index, ch->dnodes));
840
841 j = 0;
842 for (i = 0; i < n; i++) {
843 if (VECTOR(ch->hptr)[i] != UNDEFINED) {
844 VECTOR(ch->index)[j] = i;
845 VECTOR(ch->hptr)[i] = j + INDEXINC;
846 j++;
847 }
848 }
849
850 return 0;
851 }
852
853 /* -------------------------------------------------- */
854 /* Two-way indexed heap */
855 /* -------------------------------------------------- */
856
857 #undef PARENT
858 #undef LEFTCHILD
859 #undef RIGHTCHILD
860 #define PARENT(x) (((x)+1)/2-1)
861 #define LEFTCHILD(x) (((x)+1)*2-1)
862 #define RIGHTCHILD(x) (((x)+1)*2)
863
864 /* This is a smart indexed heap. In addition to the "normal" indexed heap
865 it allows to access every element through its index in O(1) time.
866 In other words, for this heap the indexing operation is O(1), the
867 normal heap does this in O(n) time.... */
868
igraph_i_2wheap_switch(igraph_2wheap_t * h,long int e1,long int e2)869 void igraph_i_2wheap_switch(igraph_2wheap_t *h,
870 long int e1, long int e2) {
871 if (e1 != e2) {
872 long int tmp1, tmp2;
873 igraph_real_t tmp3 = VECTOR(h->data)[e1];
874 VECTOR(h->data)[e1] = VECTOR(h->data)[e2];
875 VECTOR(h->data)[e2] = tmp3;
876
877 tmp1 = VECTOR(h->index)[e1];
878 tmp2 = VECTOR(h->index)[e2];
879
880 VECTOR(h->index2)[tmp1] = e2 + 2;
881 VECTOR(h->index2)[tmp2] = e1 + 2;
882
883 VECTOR(h->index)[e1] = tmp2;
884 VECTOR(h->index)[e2] = tmp1;
885 }
886 }
887
igraph_i_2wheap_shift_up(igraph_2wheap_t * h,long int elem)888 void igraph_i_2wheap_shift_up(igraph_2wheap_t *h,
889 long int elem) {
890 if (elem == 0 || VECTOR(h->data)[elem] < VECTOR(h->data)[PARENT(elem)]) {
891 /* at the top */
892 } else {
893 igraph_i_2wheap_switch(h, elem, PARENT(elem));
894 igraph_i_2wheap_shift_up(h, PARENT(elem));
895 }
896 }
897
igraph_i_2wheap_sink(igraph_2wheap_t * h,long int head)898 void igraph_i_2wheap_sink(igraph_2wheap_t *h,
899 long int head) {
900 long int size = igraph_2wheap_size(h);
901 if (LEFTCHILD(head) >= size) {
902 /* no subtrees */
903 } else if (RIGHTCHILD(head) == size ||
904 VECTOR(h->data)[LEFTCHILD(head)] >= VECTOR(h->data)[RIGHTCHILD(head)]) {
905 /* sink to the left if needed */
906 if (VECTOR(h->data)[head] < VECTOR(h->data)[LEFTCHILD(head)]) {
907 igraph_i_2wheap_switch(h, head, LEFTCHILD(head));
908 igraph_i_2wheap_sink(h, LEFTCHILD(head));
909 }
910 } else {
911 /* sink to the right */
912 if (VECTOR(h->data)[head] < VECTOR(h->data)[RIGHTCHILD(head)]) {
913 igraph_i_2wheap_switch(h, head, RIGHTCHILD(head));
914 igraph_i_2wheap_sink(h, RIGHTCHILD(head));
915 }
916 }
917 }
918
919 /* ------------------ */
920 /* These are public */
921 /* ------------------ */
922
igraph_2wheap_init(igraph_2wheap_t * h,long int size)923 int igraph_2wheap_init(igraph_2wheap_t *h, long int size) {
924 h->size = size;
925 /* We start with the biggest */
926 IGRAPH_CHECK(igraph_vector_long_init(&h->index2, size));
927 IGRAPH_FINALLY(igraph_vector_long_destroy, &h->index2);
928 IGRAPH_VECTOR_INIT_FINALLY(&h->data, 0);
929 IGRAPH_CHECK(igraph_vector_long_init(&h->index, 0));
930 /* IGRAPH_FINALLY(igraph_vector_long_destroy, &h->index); */
931
932 IGRAPH_FINALLY_CLEAN(2);
933 return 0;
934 }
935
igraph_2wheap_destroy(igraph_2wheap_t * h)936 void igraph_2wheap_destroy(igraph_2wheap_t *h) {
937 igraph_vector_destroy(&h->data);
938 igraph_vector_long_destroy(&h->index);
939 igraph_vector_long_destroy(&h->index2);
940 }
941
igraph_2wheap_clear(igraph_2wheap_t * h)942 int igraph_2wheap_clear(igraph_2wheap_t *h) {
943 igraph_vector_clear(&h->data);
944 igraph_vector_long_clear(&h->index);
945 igraph_vector_long_null(&h->index2);
946 return 0;
947 }
948
igraph_2wheap_empty(const igraph_2wheap_t * h)949 igraph_bool_t igraph_2wheap_empty(const igraph_2wheap_t *h) {
950 return igraph_vector_empty(&h->data);
951 }
952
igraph_2wheap_push_with_index(igraph_2wheap_t * h,long int idx,igraph_real_t elem)953 int igraph_2wheap_push_with_index(igraph_2wheap_t *h,
954 long int idx, igraph_real_t elem) {
955
956 /* printf("-> %.2g [%li]\n", elem, idx); */
957
958 long int size = igraph_vector_size(&h->data);
959 IGRAPH_CHECK(igraph_vector_push_back(&h->data, elem));
960 IGRAPH_CHECK(igraph_vector_long_push_back(&h->index, idx));
961 VECTOR(h->index2)[idx] = size + 2;
962
963 /* maintain heap */
964 igraph_i_2wheap_shift_up(h, size);
965 return 0;
966 }
967
igraph_2wheap_size(const igraph_2wheap_t * h)968 long int igraph_2wheap_size(const igraph_2wheap_t *h) {
969 return igraph_vector_size(&h->data);
970 }
971
igraph_2wheap_max_size(const igraph_2wheap_t * h)972 long int igraph_2wheap_max_size(const igraph_2wheap_t *h) {
973 return h->size;
974 }
975
igraph_2wheap_max(const igraph_2wheap_t * h)976 igraph_real_t igraph_2wheap_max(const igraph_2wheap_t *h) {
977 return VECTOR(h->data)[0];
978 }
979
igraph_2wheap_max_index(const igraph_2wheap_t * h)980 long int igraph_2wheap_max_index(const igraph_2wheap_t *h) {
981 return VECTOR(h->index)[0];
982 }
983
igraph_2wheap_has_elem(const igraph_2wheap_t * h,long int idx)984 igraph_bool_t igraph_2wheap_has_elem(const igraph_2wheap_t *h, long int idx) {
985 return VECTOR(h->index2)[idx] != 0;
986 }
987
igraph_2wheap_has_active(const igraph_2wheap_t * h,long int idx)988 igraph_bool_t igraph_2wheap_has_active(const igraph_2wheap_t *h, long int idx) {
989 return VECTOR(h->index2)[idx] > 1;
990 }
991
igraph_2wheap_get(const igraph_2wheap_t * h,long int idx)992 igraph_real_t igraph_2wheap_get(const igraph_2wheap_t *h, long int idx) {
993 long int i = VECTOR(h->index2)[idx] - 2;
994 return VECTOR(h->data)[i];
995 }
996
igraph_2wheap_delete_max(igraph_2wheap_t * h)997 igraph_real_t igraph_2wheap_delete_max(igraph_2wheap_t *h) {
998
999 igraph_real_t tmp = VECTOR(h->data)[0];
1000 long int tmpidx = VECTOR(h->index)[0];
1001 igraph_i_2wheap_switch(h, 0, igraph_2wheap_size(h) - 1);
1002 igraph_vector_pop_back(&h->data);
1003 igraph_vector_long_pop_back(&h->index);
1004 VECTOR(h->index2)[tmpidx] = 0;
1005 igraph_i_2wheap_sink(h, 0);
1006
1007 /* printf("<-max %.2g\n", tmp); */
1008
1009 return tmp;
1010 }
1011
igraph_2wheap_deactivate_max(igraph_2wheap_t * h)1012 igraph_real_t igraph_2wheap_deactivate_max(igraph_2wheap_t *h) {
1013
1014 igraph_real_t tmp = VECTOR(h->data)[0];
1015 long int tmpidx = VECTOR(h->index)[0];
1016 igraph_i_2wheap_switch(h, 0, igraph_2wheap_size(h) - 1);
1017 igraph_vector_pop_back(&h->data);
1018 igraph_vector_long_pop_back(&h->index);
1019 VECTOR(h->index2)[tmpidx] = 1;
1020 igraph_i_2wheap_sink(h, 0);
1021
1022 return tmp;
1023 }
1024
igraph_2wheap_delete_max_index(igraph_2wheap_t * h,long int * idx)1025 igraph_real_t igraph_2wheap_delete_max_index(igraph_2wheap_t *h, long int *idx) {
1026
1027 igraph_real_t tmp = VECTOR(h->data)[0];
1028 long int tmpidx = VECTOR(h->index)[0];
1029 igraph_i_2wheap_switch(h, 0, igraph_2wheap_size(h) - 1);
1030 igraph_vector_pop_back(&h->data);
1031 igraph_vector_long_pop_back(&h->index);
1032 VECTOR(h->index2)[tmpidx] = 0;
1033 igraph_i_2wheap_sink(h, 0);
1034
1035 if (idx) {
1036 *idx = tmpidx;
1037 }
1038 return tmp;
1039 }
1040
igraph_2wheap_modify(igraph_2wheap_t * h,long int idx,igraph_real_t elem)1041 int igraph_2wheap_modify(igraph_2wheap_t *h, long int idx, igraph_real_t elem) {
1042
1043 long int pos = VECTOR(h->index2)[idx] - 2;
1044
1045 /* printf("-- %.2g -> %.2g\n", VECTOR(h->data)[pos], elem); */
1046
1047 VECTOR(h->data)[pos] = elem;
1048 igraph_i_2wheap_sink(h, pos);
1049 igraph_i_2wheap_shift_up(h, pos);
1050
1051 return 0;
1052 }
1053
1054 /* Check that the heap is in a consistent state */
1055
igraph_2wheap_check(igraph_2wheap_t * h)1056 int igraph_2wheap_check(igraph_2wheap_t *h) {
1057 long int size = igraph_2wheap_size(h);
1058 long int i;
1059 igraph_bool_t error = 0;
1060
1061 /* Check the heap property */
1062 for (i = 0; i < size; i++) {
1063 if (LEFTCHILD(i) >= size) {
1064 break;
1065 }
1066 if (VECTOR(h->data)[LEFTCHILD(i)] > VECTOR(h->data)[i]) {
1067 error = 1; break;
1068 }
1069 if (RIGHTCHILD(i) >= size) {
1070 break;
1071 }
1072 if (VECTOR(h->data)[RIGHTCHILD(i)] > VECTOR(h->data)[i]) {
1073 error = 1; break;
1074 }
1075 }
1076
1077 if (error) {
1078 IGRAPH_ERROR("Inconsistent heap", IGRAPH_EINTERNAL);
1079 }
1080
1081 return 0;
1082 }
1083