briarpig: code: bheap

bheap

I wrote about binary min heap priority queues earlier this month. Material below is a verbatim copy of what appeared in the log. The only reason I put a copy here is to clarify license and copyright.

purpose

The purpose of cy_event and cy_heap is to support a simple kind of priority queue for events delivered under a virtual clock.

12apr09 ¶ tastier tomatoes

event simulation ¶

Below I present some code I plan to use for event simulation. I need something like this at work, and it seems trivial to me, but last week my manager pulled me aside and asked "have you written an event simulator before?" and then looked at me with incredulity, no matter what I said. I understand someone at work wants to download an open source event simulator and go to town. How about this one? All the code below using a cy_ prefix has a BSD license. Normally I'd never open source anything with any relation to a work task. But in this case, it was made clear to me no work task would arise related to event simulation, because it sounds like a waste of my time. For some reason, my time is hoarded like something too valuable to waste. (Because it "sounds" complex, folks assume something complex is involved. But that's not so. Tell me if you think this code is hard.)

I wrote all this today (Sunday afternoon) instead of writing more fiction. It was also kinda fun, so I don't miss writing the fiction, which can wait.

However, before the code, first a little introduction. I wrote a discrete event simulator in college because a professor asked for one. It wasn't a class assignment. I was one of the few students she could ask for bespoke code and expect to get something useful. This was a long time ago, but still during a period in which my memory was excellent. (I still recall most things from that period just as clearly—I think.)

To model virtual time is surprisingly easy. All you need to do is queue up events to be delivered, then deliver a next event whose time of delivery is earliest. When this event is delivered, you set your virtual clock to use the time of this event. All I can say about this is: it was obvious to me. (Lots of things are.)

The code below is designed around the idea just keeping a priority queue in the form of a binary heap is sufficient. I poked around in Wikipedia for references and add these links to a comment, so folks reading the code later can go find a clue:

/* SEE ALSO: * http://en.wikipedia.org/wiki/Heapsort * http://en.wikipedia.org/wiki/Binary_heap * http://en.wikipedia.org/wiki/Heap_(data_structure) * http://en.wikipedia.org/wiki/Heapsort * http://en.wikipedia.org/wiki/Priority_queue * http://en.wikipedia.org/wiki/Discrete_event_simulation */

All this code is in C, or something very near C, because some clients will be anti C++. If I write it in C++, there's too big a chance it will be spat upon. I don't enjoy writing C as much as C++, but I no longer particularly notice any more when C is my only choice, because I transliterate C++ to C as I write the code, hand-rolling vtables as needed. I understand this kind of C code disturbs a few folks. (Too bad.)

Here's a standard macro to inform the compiler a branch is not expected to be taken:

#define cy_unlikely(x) __builtin_expect((x)!=0,0)

And here's one of my zillions of hand-rolled pseudo random number generators. (This one comes from the hash demo.)

static inline unsigned long cy_u64hash32(unsigned long long n) { /* note: n is 64-bit ONLY so n * 48271 cannot */ /* overflow 32 bits; note: largest prime in */ /* 32-bits: 4294967291 == 0xfffffffb. */ return (unsigned long) ((n * 48271) % 4294967291ULL); };

Later stages of simulation will use floating point guassian distributions based on the Box Muller transformation shown in the rand demo, but that's not used in code below. However, it should be obvious that pseudo random time delays will involve a gaussian means and standard deviations. (I guess I first learned about use of gaussian distributions in discrete event simulation when I wrote the code in school many years ago.)

Folks make it hard to use standard iovecs, so:

typedef struct cy_iovec_ { /* clone of iovec */ void* iov_base; size_t iov_len; } cy_iovec_t;

Below is an event to put in priority queues. It's small because it assumes we need nothing more than virtual time of delivery and an object. The "object" part is code (function pointer) and state (opaque data context). These are values meant to be copied in and out.

typedef void (*cy_event_cb)(void* cx, uint64_t when); typedef struct cy_event_ { /* function vtable */ uint64_t ev_when; /* virtual microseconds */ cy_event_cb ev_cb; /* event callback */ void* ev_cx; /* event context */ } cy_event; static inline void /* convenience inline method */ cy_event_cb_do(const cy_event* e) { /* run callback */ if (e->ev_cb) { (*e->ev_cb)(e->ev_cx, e->ev_when); } }

The binary min heap used as a priority queue below is slightly more complex because it has "virtual" methods for memory allocation so the event array can be dynamically resized if you provide function pointers to alloc and free memory. (However, note I haven't tested the reserve() method even though I wrote it carefully. I only used statically allocated memory in testing.)

typedef struct cy_bheap_ cy_bheap; /* forward for methods */ typedef void* (*cy_bh_malloc_fn)(cy_bheap* h, size_t size); typedef void (*cy_bh_free_fn)(cy_bheap* h, void *ptr);

Those are malloc() and free() virtual methods.

/* index zero is never used in this binary heap formulation: */ #define cy_bheap_ROOT 1 /* root member always at index ONE */

I use the binary heap definition ignoring the zero array element, so a left child and right child are defined as follows:

/* arg i MUST NOT BE ZERO in left() and right(): */ static inline unsigned /* index of left child of i */ cy_bheap_left(unsigned i) { return 2 * i; } static inline unsigned /* index of right child of i */ cy_bheap_right(unsigned i) { return (2 * i)+1; }

At the root (least) member is always at index 1:

/* index zero is never used in this binary heap formulation: */ #define cy_bheap_ROOT 1 /* root member always at index ONE */

Here's a few more useful inlines related to root and children. (Sometimes I don't bother using these inlines because they obfuscate code. But I like having them around to define rules clearly.)

static inline int /* left child of j, or -1 if it's missing */ cy_bheap_lch(const cy_bheap* h, unsigned j) { unsigned l = 2 * j; if (l <= h->bh_size) return (int) l; /* left is in bounds with a context */ else return -1; } static inline int /* right child of j, or -1 if it's missing */ cy_bheap_rch(const cy_bheap* h, unsigned j) { unsigned r = (2 * j)+1; if (r <= h->bh_size) return (int) r; /* right is in bounds with a context */ else return -1; } static inline unsigned /* index of parent for child at i */ cy_bheap_parent(unsigned i) { return i / 2; } extern cy_event cy_event_oob; /* out-of-bounds default event */ static inline cy_event /* least event now in heap */ cy_bheap_root(const cy_bheap* h, unsigned i) { if (h->bh_events && h->bh_size) { /* not empty? */ return h->bh_events[cy_bheap_ROOT]; } return cy_event_oob; /* default for out-of-bounds */ }

The binary heap object looks like this:

#define cy_bheap_TAG 0x42684550 /*'BhEP'*/ struct cy_bheap_ { /* binary 'heap' (for priority queue)*/ /* bh_now tracks the usecs in each least event removed: */ uint64_t bh_now; /* virtual clock for this heap*/ cy_bh_malloc_fn bh_malloc; /* malloc() */ cy_bh_free_fn bh_free; /* free() */ void* bh_cx; /* if malloc() and free() need */ uint32_t bh_tag; /* must equal TAG */ cy_event* bh_events; /* array of events */ /* size is the index of the last USED array slot: */ uint32_t bh_size; /* count of USED events */ /* REQUIRED: capacity >= size+1 (zero index is unused) */ uint32_t bh_capacity; /* physical array length */ };

Whose constructor and destructor are:

/*you can either alloc all events up front, or dynamic alloc: */ static inline void cy_bheap_ctor(cy_bheap* h, cy_event* v, unsigned capacity) { h->bh_now = 0; /* beginning of time: zero microseconds */ h->bh_malloc = 0; /* replace if you need auto-grow */ h->bh_free = 0; /* replace if you need auto-grow */ h->bh_cx = 0; /* replace if methods above need a context */ h->bh_tag = cy_bheap_TAG; /* magic signature */ h->bh_events = v; /* vector of starting empty space*/ h->bh_size = 0; /* no members in heap at first */ h->bh_capacity = (v)? capacity : 0; /*starting events (if any)*/ } static inline void cy_bheap_dtor(cy_bheap* h) { void* events = h->bh_events; cy_bh_free_fn fn = h->bh_free; /* nil if not needed */ if (cy_bheap_TAG == h->bh_tag) { /* heap seems valid? */ h->bh_tag = 0xdeadbeef; if (fn) { /* we have a method to deallocate events? */ h->bh_events = 0; /* never do this again */ h->bh_free = 0; (*fn)(h, events); /* free array space if desired */ } } else { /* remove this block if you want no err checking */ cy_printf("bh_tag=0x%lx not TAG=0x%lx", (long) h->bh_tag, (long) cy_bheap_TAG); assert(cy_bheap_TAG == h->bh_tag); /* die */ } }

bheap api ¶

The primary api looks like this:

#ifdef cy_bheap_PRINTING void cy_bheap_print(const cy_bheap* h, cy_sink* s); #endif /*cy_bheap_PRINTING*/ void cy_bheap_test(uint32_t seed); /* needs printing code */ /* reserve() tries to make sure capacity is at least min: */ int /* zero: success; nonzero: errno (probably no memory) */ cy_bheap_reserve(cy_bheap* h, unsigned minCapacity); /* insert() appends event, then bubbles upward */ int /* zero: success; nonzero: errno (probably no memory) */ cy_bheap_insert(cy_bheap* h, cy_event e); /* add copy of e */ /* append() is more efficient than N adds for large N */ int /* zero: success; nonzero: errno (probably no memory) */ cy_bheap_append(cy_bheap* h, const cy_event* v, unsigned len); /* side effect of remove(): bh_now = max(bh_now, ev_when) */ cy_event /* the root of the tree, now removed from the heap*/ cy_bheap_remove(cy_bheap* h); /* event with least ev_when*/ /* delay means: insert(h, cy_event(h->bh_now+usecs, cb, cx)) */ int /* queue future event at now + usecs (delayed by usecs) */ cy_bheap_delay(cy_bheap* h, unsigned usecs, cy_event_cb cb, void* cx) { cy_event e; e.ev_when = h->bh_now + usecs; e.ev_cb = cb; e.ev_cx = cx; return cy_bheap_insert(h, e); }

bheap algorithm ¶

I added a helpful note on algorithms:

/* BASIC ALGORITHMS: remove and insert, as follows: * REMOVE: * cy_event outputEv = bh_events[cy_bheap_ROOT]; * unsigned i = bh_size--; * bh_events[cy_bheap_ROOT] = bh_events[i]; // can't be too low * cy_bheap_siftdown(this, 1); // mv down while any child is less * INSERT: * cy_bheap_reserve(bh_size+1) // ensure room exists * unsigned i = ++bh_size; // new last event in array * bh_events[i] = inputEv; // can't be too hi starting at bottom * cy_bheap_siftup(this, i); // mv child up while parent is more */

bheap reserve ¶

Here's how you reserve capacity. Note this hasn't been tested or even run. If it looks debugged, that just reflects how I think carefully. Yes, I know this sort of thing is hard for some folks. But it's old hat to me.

/* reserve() tries to make sure capacity is at least min: */ int /* zero: success; nonzero: errno (probably no memory) */ cy_bheap_reserve(cy_bheap* h, unsigned minCapacity) { assert(cy_bheap_TAG == h->bh_tag); if (minCapacity > 1024*1024) /* over a million?*/ minCapacity = 1024*1024; /* that's enough */ if (h->bh_capacity < minCapacity) { /* need more? */ cy_bh_malloc_fn fn = h->bh_malloc; if (fn) { /* we can alloc memory? */ cy_event* old = h->bh_events; /* old content */ unsigned have = h->bh_size + 1; /* unused zero */ unsigned want = minCapacity + 1024; /* more */ unsigned cmin = 4 * (h->bh_capacity/3); /* + 1/3 */ if (want < cmin) /* less than one third increase? */ want = cmin; /* geometric growth when bigger */ unsigned newvol = sizeof(cy_event) * want; /* bytes */ cy_event* e = (cy_event*) (*fn)(h, newvol); /* alloc */ if (e) { /* good allocation? */ if (old) { /* anything to copy or free? */ uint8_t* p = (uint8_t*) e; /* byte pointer */ unsigned oldvol = have * sizeof(cy_event); memcpy(p, old, oldvol); /* copy old events */ unsigned extra = newvol - oldvol; /* bytes */ memset(p + oldvol, 0, extra); /* zero new */ assert(0 != h->bh_free); /* required */ (*h->bh_free)(h, old); /* reclaim */ } else { /* all zero with nothing to copy */ memset(e, 0, newvol); } h->bh_events = e; h->bh_capacity = want; return 0; /* bigger capacity is good to go */ } } return ENOMEM; /* cannot alloc or failed alloc */ } return 0; /* capacity already big enough */ }

bheap insert ¶

Below a single new event is inserted.

void /* swap child at j w/ parent that's bigger, then recurse */ cy_bheap_siftup(cy_bheap* h, unsigned j) { cy_event* v = h->bh_events; uint64_t when = v[j].ev_when; /* key to move upward */ while (j > cy_bheap_ROOT) { /* j has a parent? */ unsigned p = j / 2; /* parent */ if (v[p].ev_when > when) { /* bigger parent?*/ cy_event t = v[j]; v[j] = v[p]; v[p] = t; /* swap */ j = p; /* recurse on parent */ } else { /* done when parent is not bigger */ break; } } } /* insert() appends event, then bubbles upward */ int /* zero: success; nonzero: errno (probably no memory) */ cy_bheap_insert(cy_bheap* h, cy_event e) { /* add copy of e */ unsigned i = h->bh_size + 1; /* where to insert at first */ assert(cy_bheap_TAG == h->bh_tag); if (cy_unlikely(i >= h->bh_capacity)) { /* need slots? */ int err = cy_bheap_reserve(h, i + 32); if (err) /* ENOMEM? */ return err; /* cannot insert */ } assert(i < h->bh_capacity); /* must have slot to use */ h->bh_size = i; /* alloc first unused slot */ h->bh_events[i] = e; cy_bheap_siftup(h, i); /* move upward while parent is less */ return 0; }

bheap remove ¶

Removing the least event (the root) works like this:

void /* swap event at j with smaller child, then recurse */ cy_bheap_siftdown(cy_bheap* h, unsigned j) { cy_event* v = h->bh_events; uint64_t when = v[j].ev_when; /* key to move downward */ unsigned last = h->bh_size; /* last used event */ unsigned parent = last / 2; /* parent of last event */ while (j <= parent) { /* might have children? */ unsigned l = 2 * j; /* left child */ unsigned r = l + 1; /* right child */ unsigned c = l; /* assume left is least child */ if (r <= last && v[r].ev_when < v[c].ev_when) { c = r; /* use r as index of least child */ } if (when > v[c].ev_when) { /* smaller child?*/ cy_event t = v[j]; v[j] = v[c]; v[c] = t; /* swap */ j = c; /* recurse on least child */ } else { /* done when neither child is smaller */ break; } } } /* side effect of remove(): bh_now = max(bh_now, ev_when) */ cy_event /* the root of the tree, now removed from the heap*/ cy_bheap_remove(cy_bheap* h) { /* event with least ev_when*/ cy_event outEv; cy_event* v = h->bh_events; unsigned i = h->bh_size; /* used heap array elems */ if (v && i) { /* have anything to remove? */ assert(cy_bheap_TAG == h->bh_tag); assert(h->bh_capacity > i); /* must exceed size */ outEv = v[cy_bheap_ROOT]; /* least event in heap */ if (h->bh_now < outEv.ev_when) { h->bh_now = outEv.ev_when; } v[cy_bheap_ROOT] = v[i]; /* last event becomes root */ h->bh_size = i - 1; /* heap is one smaller */ cy_bheap_siftdown(h, cy_bheap_ROOT); } else { memset(&outEv, 0, sizeof(cy_event)); /* nothing */ } return outEv; }

bheap append ¶

You append an array of unordered events like this:

/* append() is more efficient than N adds for large N */ int /* zero: success; nonzero: errno (probably no memory) */ cy_bheap_append(cy_bheap* h, const cy_event* vec, unsigned len) { assert(cy_bheap_TAG == h->bh_tag); unsigned sz = h->bh_size + len; /* resulting size */ if (cy_unlikely(sz >= h->bh_capacity)) { /* need slots? */ int err = cy_bheap_reserve(h, sz + 32); if (err) /* ENOMEM? */ return err; /* cannot insert */ } assert(sz < h->bh_capacity); /* must have slots to use */ cy_event* v = h->bh_events; /* event array (maybe new) */ unsigned addvol = sizeof(cy_event) * len; /* bytes to add */ unsigned prefix = h->bh_size + 1; /* include unused zero */ memcpy(v + prefix, vec, addvol); /* append events */ h->bh_size = sz; /* count all new events as added */ unsigned parent = sz / 2; /* parent of last event */ /* note: all events at indexes > sz / 2 have no children */ for (unsigned i = parent; i >= cy_bheap_ROOT; --i) { cy_bheap_siftdown(h, i); /* down while child is more */ } return 0; }

bheap testing ¶

Now, you might find part of my test below a little disturbing. I ported my C++ out stream class from the C++ (in the out demo) to C, so I could test the binary heap using a C based stream api.

I only did this so folks who dislike C++ can print the binary heaps using C only. And perhaps I might find that C stream useful again. However, the C code for the stream is quite long—it took me an hour or two write. I'll show that code last, after test code and resulting output.

void cy_bheap_test(uint32_t seed) { /* needs printing code */ unsigned i = 0; cy_event finiteVec[ 2048 ]; cy_event inboundVec[ 1024 ]; /* for vector append */ cy_event e; cy_bheap minheap; char buf[ 2048 + 4 ]; cy_fdsink sout; cy_iovec_t bufiov; bufiov.iov_base = buf; bufiov.iov_len = 2048; cy_fdsink_ctor(&sout, bufiov, STDOUT_FILENO); cy_bheap_ctor(&minheap, finiteVec, 2048); for ( ; i < 16; ++i) { seed = cy_u64hash32(seed); e.ev_when = seed; cy_bheap_insert(&minheap, e); } cy_bheap_print(&minheap, &sout.f_sink); cy_sink_n(&sout.f_sink); for (i = 0; i < 8; ++i) { seed = cy_u64hash32(seed); inboundVec[i].ev_when = seed; cy_sink_f(&sout.f_sink, " %08lx", seed); } cy_sink_n(&sout.f_sink); cy_bheap_append(&minheap, inboundVec, 8); cy_bheap_print(&minheap, &sout.f_sink); while (minheap.bh_size) { (void) cy_bheap_remove(&minheap); cy_bheap_print(&minheap, &sout.f_sink); } }

The print method and the fdsink stream api appear in the next section. Here's the output on stdout:

bheap print ¶

This printing code is recursive. (Duh.) If you don't like the stream api, too bad. Note code detecting inverted time order between parents and children. It happened and I needed to debug it, so I added this code.

#ifdef cy_bheap_PRINTING void cy_bheap_show(const cy_bheap* h, cy_sink* s, int j) { int inverted = (j < 0); /* negative index? */ if (inverted) /* note inverted order with parent */ j = -j; /* only use positive indexes */ uint64_t here = h->bh_events[j].ev_when; int r = cy_bheap_rch(h, j); /* right: -1 if missing */ if (r > 0) { /* has a right child too? */ if (here > h->bh_events[r].ev_when) { /* inverted?? */ r = -r; /* show inversion */ assert(here <= h->bh_events[r].ev_when); /* die */ } cy_sink_t(s); /* tab then newline indent */ cy_bheap_show(h, s, r); cy_sink_u(s); /* untab */ } if (inverted) /* show inversion of order with parent: */ cy_sink_nf(s, "%u @@ %08llx ", j, here); else cy_sink_nf(s, "%u: %08llx ", j, here); int l = cy_bheap_lch(h, j); /* left: -1 if missing */ if (l > 0) { /* has a left child? */ if (here > h->bh_events[l].ev_when) { /* inverted?? */ l = -l; /* show inversion */ assert(here <= h->bh_events[l].ev_when); /* die */ } cy_sink_t(s); /* tab then newline indent */ cy_bheap_show(h, s, l); cy_sink_u(s); /* untab */ } } void cy_bheap_print(const cy_bheap* h, cy_sink* s) { cy_sink_f(s, "<bheap now=%08llx size=%d cap=%d>", (long long) h->bh_now, (int) h->bh_size, (int) h->bh_capacity); if (h->bh_size) { cy_bheap_show(h, s, cy_bheap_ROOT); } cy_sink_nf(s, "</bheap>"); cy_sink_n(s); sink_flush_do(s); } #endif /*cy_bheap_PRINTING*/

menu

Here's a menu of pages on cy code.

vector - std::vector clone
bheap - binary min heap
label - ascync descriptors
misc - bunch of basic utils
pool - C-based vats and pools
deck - scatter/gather suite
sink - C-based out stream api
row - a new deck rewrite

license

See license and copyright for code here. For more context, see the cy page.

comments

Compared to a thorn demo, I explain cy code less: I care little whether folks use or grasp cy source. But since I aim to get ideas across, I reveal a point to code constructs so you see intentions.

If you ask: What was this for? That's the only question I address: why a thing was done. If you what to know how code works or what loose ends remain, figure it out.

color coding

Library source code appears appears in amber (orange/brown):

amber is_source(code* c);

Source .cpp code appears in red:

void cy_logf(int, const char* f, ...) { char temp[ 2048 + 4 ]; va_list args; va_start(args,f); vsnprintf(temp, 2048, f, args); va_end(args); temp[2048] = 0; printf("%s\n", temp); }

Sample test code is purple:

o << "# purple=test green=stdout" << cy_newl;

Printed output on stdout is green:

# purple=test green=stdout

I know these aren't the best color cues. (Amber and green might appear the same hue to color blind folks. I have excellent color discrimination myself.)