Distances on Directed Graphs in R
mpadge.tngl.sh
1/* File fheap.c - Fibonacci Heap
2 * ----------------------------------------------------------------------------
3 * Mark Padgham, adapted from code by Shane Saunders
4 */
5#include <cmath>
6#if defined(FHEAP_DUMP) && FHEAP_DUMP > 0
7#include <cstdio>
8#endif
9#include "fheap.h"
10
11/*--- FHeap (public methods) ------------------------------------------------*/
12
13/* --- Constructor ---
14 * Creates an FHeap object capable of holding up to $n$ items.
15 */
16FHeap::FHeap(size_t n)
17{
18#if FHEAP_DUMP
19 Rcpp::Rcout << "new, ";
20#endif
21 maxTrees = 1 + static_cast <size_t>(1.44 *
22 log(static_cast <double> (n)) / log (2.0));
23 maxNodes = n;
24
25 trees = new FHeapNode *[maxTrees];
26 for(size_t i = 0; i < maxTrees; i++) trees[i] = std::nullptr_t ();
27
28 nodes = new FHeapNode *[n];
29 for(size_t i = 0; i < n; i++) nodes[i] = std::nullptr_t ();
30
31 itemCount = 0;
32
33 /* The $treeSum$ of the heap helps to keep track of the maximum rank while
34 * nodes are inserted or deleted.
35 */
36 treeSum = 0;
37
38 /* For experimental purposes, we keep a count of the number of key
39 * comparisons.
40 */
41 compCount = 0;
42
43#if FHEAP_DUMP
44 Rcpp::Rcout << "new-exited, ";
45#endif
46}
47
48/* --- Destructor ---
49*/
50FHeap::~FHeap()
51{
52#if FHEAP_DUMP
53 Rcpp::Rcout << "delete, ";
54#endif
55
56 for(size_t i = 0; i < maxNodes; i++) delete nodes[i];
57 delete [] nodes;
58 delete [] trees;
59
60#if FHEAP_DUMP
61 Rcpp::Rcout << "delete-exited, ";
62#endif
63}
64
65/* --- insert() ---
66 * Inserts an item $item$ with associated key $k$ into the heap.
67 */
68void FHeap::insert(size_t item, double k)
69{
70 FHeapNode *newNode;
71
72#if FHEAP_DUMP
73 Rcpp::Rcout << "insert, ";
74#endif
75
76 /* create an initialise the new node */
77 newNode = new FHeapNode;
78 newNode->child = std::nullptr_t ();
79 newNode->left = newNode->right = newNode;
80 newNode->rank = 0;
81 newNode->item = item;
82 newNode->key = k;
83
84 /* maintain a pointer to $item$'s new node in the heap */
85 nodes[item] = newNode;
86
87 /* meld the new node into the heap */
88 meld(newNode);
89
90 /* update the heaps node count */
91 itemCount++;
92
93#if FHEAP_DUMP
94 Rcpp::Rcout << "insert-exited, ";
95#endif
96}
97
98/* --- deleteMin() ---
99 * Deletes and returns the minimum item from the heap.
100 */
101size_t FHeap::deleteMin()
102{
103 FHeapNode *minNode, *child, *next;
104 double k, k2;
105 size_t r, v, item;
106
107#if FHEAP_DUMP
108 Rcpp::Rcout << "deleteMin, ";
109#endif
110
111 /* First we determine the maximum rank in the heap. */
112 v = treeSum;
113 // MP Note: r was an int formerly intialised at -1, but it's used as a
114 // direct array index below, so really should be unsigned! The r-- at end
115 // fixes that.
116 r = 0;
117 while(v) {
118 v = v >> 1;
119 r++;
120 };
121 r--;
122
123 /* Now determine which root node is the minimum. */
124 minNode = trees[r];
125 k = minNode->key;
126 while(r > 0) {
127 r--;
128 next = trees[r];
129 if(next) {
130 if((k2 = next->key) < k) {
131 k = k2;
132 minNode = next;
133 }
134 compCount++;
135 }
136 }
137
138 /* We remove the minimum node from the heap but keep a pointer to it. */
139 r = minNode->rank;
140 trees[r] = std::nullptr_t ();
141 treeSum -= (1 << r);
142
143 child = minNode->child;
144 if(child) meld(child);
145
146 /* Record the vertex no of the old minimum node before deleting it. */
147 item = minNode->item;
148 nodes[item] = std::nullptr_t ();
149 delete minNode;
150 itemCount--;
151
152#if FHEAP_DUMP
153 Rcpp::Rcout << "deleteMin-exited, ";
154#endif
155
156 return item;
157}
158
159/* --- decreaseKey() ---
160 * Decreases the key used for item $item$ to the value newValue. It is left
161 * for the user to ensure that newValue is in-fact less than the current value
162 */
163void FHeap::decreaseKey(size_t item, double newValue)
164{
165 FHeapNode *cutNode, *parent, *newRoots, *r, *l;
166 size_t prevRank;
167
168#if FHEAP_DUMP
169 Rcpp::Rcout << "decreaseKey on vn = " << item << ", ";
170#endif
171
172 /* Obtain a pointer to the decreased node and its parent then decrease the
173 * nodes key.
174 */
175 cutNode = nodes[item];
176 parent = cutNode->parent;
177 cutNode->key = newValue;
178
179 /* No reinsertion occurs if the node changed was a root. */
180 if(!parent) {
181#if FHEAP_DUMP
182 Rcpp::Rcout << "decreaseKey-exited, ";
183#endif
184 return;
185 }
186
187 /* Update the left and right pointers of cutNode and its two neighbouring
188 * nodes.
189 */
190 l = cutNode->left;
191 r = cutNode->right;
192 l->right = r;
193 r->left = l;
194 cutNode->left = cutNode->right = cutNode;
195
196 /* Initially the list of new roots contains only one node. */
197 newRoots = cutNode;
198
199 /* While there is a parent node that is marked a cascading cut occurs. */
200 while(parent && parent->marked) {
201
202 /* Decrease the rank of cutNode's parent and update its child pointer.
203 */
204 parent->rank--;
205 if(parent->rank) {
206 if(parent->child == cutNode) parent->child = r;
207 }
208 else {
209 parent->child = std::nullptr_t ();
210 }
211
212 /* Update the cutNode and parent pointers to the parent. */
213 cutNode = parent;
214 parent = cutNode->parent;
215
216 /* Update the left and right pointers of cutNodes two neighbouring
217 * nodes.
218 */
219 l = cutNode->left;
220 r = cutNode->right;
221 l->right = r;
222 r->left = l;
223
224 /* Add cutNode to the list of nodes to be reinserted as new roots. */
225 l = newRoots->left;
226 newRoots->left = l->right = cutNode;
227 cutNode->left = l;
228 cutNode->right = newRoots;
229 newRoots = cutNode;
230 }
231
232 /* If the root node is being relocated then update the trees[] array.
233 * Otherwise mark the parent of the last node cut.
234 */
235 if(!parent) {
236 prevRank = cutNode->rank + 1;
237 trees[prevRank] = std::nullptr_t ();
238 treeSum -= (1 << prevRank);
239 }
240 else {
241 /* Decrease the rank of cutNode's parent an update its child pointer.
242 */
243 parent->rank--;
244 if(parent->rank) {
245 if(parent->child == cutNode) parent->child = r;
246 }
247 else {
248 parent->child = std::nullptr_t ();
249 }
250
251 parent->marked = 1;
252 }
253
254 /* Meld the new roots into the heap. */
255 meld(newRoots);
256
257#if FHEAP_DUMP
258 Rcpp::Rcout << "decreaseKey-exited, ";
259#endif
260}
261
262/*--- FHeap (private methods) -----------------------------------------------*/
263
264/* --- meld() ---
265 * melds the linked list of trees pointed to by $treeList$ into the heap.
266 */
267void FHeap::meld(FHeapNode *treeList)
268{
269 FHeapNode *first, *next, *nodePtr, *newRoot, *temp, *temp2, *lc, *rc;
270 size_t r;
271
272#if FHEAP_DUMP
273 Rcpp::Rcout << "meld: ";
274#endif
275
276 /* We meld each tree in the circularly linked list back into the root level
277 * of the heap. Each node in the linked list is the root node of a tree.
278 * The circularly linked list uses the sibling pointers of nodes. This
279 * makes melding of the child nodes from a deleteMin operation simple.
280 */
281 nodePtr = first = treeList;
282
283 do {
284
285#if FHEAP_DUMP
286 Rcpp::Rcout << nodePtr->item << ", ";
287#endif
288
289 /* Keep a pointer to the next node and remove sibling and parent links
290 * from the current node. nodePtr points to the current node.
291 */
292 next = nodePtr->right;
293 nodePtr->right = nodePtr->left = nodePtr;
294 nodePtr->parent = std::nullptr_t ();
295
296 /* We merge the current node, nodePtr, by inserting it into the
297 * root level of the heap.
298 */
299 newRoot = nodePtr;
300 r = nodePtr->rank;
301
302 /* This loop inserts the new root into the heap, possibly restructuring
303 * the heap to ensure that only one tree for each degree exists.
304 */
305 do {
306
307 /* Check if there is already a tree of degree r in the heap.
308 * If there is then we need to link it with newRoot so it will be
309 * reinserted into a new place in the heap.
310 */
311 if((temp = trees[r])) {
312
313 /* temp will be linked to newRoot and relocated so we no
314 * longer will have a tree of degree r.
315 */
316 trees[r] = std::nullptr_t ();
317 treeSum -= (1 << r);
318
319 /* Swap temp and newRoot if necessary so that newRoot always
320 * points to the root node which has the smaller key of the
321 * two.
322 */
323 if(temp->key < newRoot->key) {
324 temp2 = newRoot;
325 newRoot = temp;
326 temp = temp2;
327 }
328 compCount++;
329
330 /* Link temp with newRoot, making sure that sibling pointers
331 * get updated if rank is greater than 0. Also, increase r for
332 * the next pass through the loop since the rank of new has
333 * increased.
334 */
335 if(r++ > 0) {
336 rc = newRoot->child;
337 lc = rc->left;
338 temp->left = lc;
339 temp->right = rc;
340 lc->right = rc->left = temp;
341 }
342 newRoot->child = temp;
343 newRoot->rank = r;
344 temp->parent = newRoot;
345 temp->marked = 0;
346 }
347 /* Otherwise if there is not a tree of degree r in the heap we
348 * allow newRoot, which possibly carries moved trees in the heap,
349 * to be a tree of degree r in the heap.
350 */
351 else {
352
353 trees[r] = newRoot;
354 treeSum += (1 << r);;
355
356 /* NOTE: Because newRoot is now a root we ensure it is
357 * marked.
358 */
359 newRoot->marked = 1;
360 }
361
362 /* Note that temp will be NULL if and only if there was not a tree
363 * of degree r.
364 */
365 } while(temp);
366
367 nodePtr = next;
368
369 } while(nodePtr != first);
370
371#if FHEAP_DUMP
372 Rcpp::Rcout << "meld-exited, ";
373#endif
374}
375
376/*--- FHeap (debugging) -----------------------------------------------------*/
377
378/* --- dumpNodes() ---
379 * Recursively print a text representation of the node subtree starting at the
380 * node poined to by $node$, using an indentationlevel of $level$.
381 */
382void FHeap::dumpNodes(FHeapNode *node, size_t level)
383{
384#if FHEAP_DUMP
385 FHeapNode *childNode, *partner;
386 size_t childCount;
387
388 /* Print leading whitespace for this level. */
389 for(size_t i = 0; i < level; i++)
390 Rcpp::Rcout << " ";
391
392 Rcpp::Rcout << node->item << "(" << node->key << ")[" << node->rank <<
393 "]" << std::endl;
394
395 if((childNode = node->child)) {
396 childNode = node->child->right;
397
398 childCount = 0;
399
400 do {
401 dumpNodes(childNode, level+1);
402 if(childNode->dim > node->dim) {
403 for(size_t i = 0; i < level+1; i++) Rcpp::Rcout << " ";
404 throw std::runtime_error ("error(dim)");
405 }
406 if(childNode->parent != node) {
407 for(size_t i = 0; i < level+1; i++) Rcpp::Rcout << " ";
408 throw std::runtime_error ("error(parent)");
409 }
410 childNode = childNode->right;
411 childCount++;
412 } while(childNode != node->child->right);
413
414 if(childCount != node->dim) {
415 for(size_t i = 0; i < level; i++) Rcpp::Rcout << " ";
416 throw std::runtime_error ("error(childCount)");
417 }
418 }
419 else {
420 if(node->dim != 0) {
421 for(size_t i = 0; i < level; i++) Rcpp::Rcout << " ";
422 throw std::runtime_error ("error(dim)");
423 }
424 }
425#endif
426}
427
428/* --- dump() ---
429 * Print a text representation of the heap to the standard output.
430 */
431void FHeap::dump() const
432{
433#if FHEAP_DUMP
434 FHeapNode *node;
435
436 Rcpp::Rcout << std::endl;
437 Rcpp::Rcout << "treeSum = " << treeSum << std::endl;
438 Rcpp::Rcout << "array entries 0..maxTrees =";
439 for(size_t i = 0; i < maxTrees; i++) {
440 bool tempval = trees [i] ? 1 : 0;
441 Rcpp::Rcout << " " << tempval;
442 }
443 Rcpp::Rcout << std::endl << std::endl;
444 for(size_t i = 0; i < maxTrees; i++) {
445 if((node = trees[i])) {
446 Rcpp::Rcout << "tree " << i;
447 dumpNodes(node, 0);
448 Rcpp::Rcout << std::endl;
449 }
450 }
451 Rcpp::Rcout.flush ();
452#endif
453}
454
455
456/*---------------------------------------------------------------------------*/