Binary heap

A box holds a hundred and fifty discs. Each flies in a straight line and bounces off the walls and its neighbours, and the program needs to know one thing: what happens first. Moving all the discs a little at a time and checking whether they have run into each other is slow and inaccurate: with a large step the discs pass through one another. So the program computes in advance the moment of impact for neighbouring discs, and for each disc the moment it reaches a wall, and writes all this into a calendar. Then, over and over, it takes the earliest entry from the calendar, and it needs nothing else from the calendar. The calendar is a binary heap, and its top is shown next to the box.

At t = 10 the scene stops: by then the program has taken 34,103 entries from the heap, and 22,279 of them, almost two thirds, turned out stale. The forecast “these two discs will collide at 4.317” stops being true as soon as one of them hits a third disc earlier. The program did not hunt down and delete such entries; it threw each one away when it surfaced at the top by itself. There were 11,824 real events: 3061 collisions, 937 bounces off the walls and 7826 crossings of grid lines (why there is a grid, the section with the code explains). The heap grew to 1600 entries, and all its operations cost 679 thousand comparisons. Had the calendar been an ordinary list searched for the nearest entry by a scan, it would have needed 46.7 million comparisons, almost 69 times more.

Watch the histogram too. At the start all discs have the same speed, 0.5, and there is a single bar. By t ≈ 2 the impacts mix the speeds, and the histogram settles on the curve James Clerk Maxwell derived in 1860, or rather on its version for a plane. There are no random numbers in the program after the start: the distribution comes out of the mechanics of elastic impacts alone.

The idea

A binary heap is a complete binary tree in which each node is no greater than its children. Hence the smallest element is always at the root. The tree lies in an array without a single pointer: the children of cell i are cells 2i + 1 and 2i + 2, the parent is cell ⌊(i − 1)/2⌋. Adding an element or extracting the smallest means repairing the order along one path from the root to a leaf, and such a path is never longer than log₂ n.

The condition “a parent is no greater than its children” is called the heap property. It is surprisingly weak. Siblings are not ordered among themselves, neighbouring branches have no say over each other, and finding an arbitrary value in a heap is as hard as in a sack. But the vertical order, keys never decrease from the root to any leaf, gives what a calendar needs: the earliest entry is on top. A sorted array and a search tree keep the full order and pay for it: an insertion into the middle of an array shifts the tail, a search tree keeps pointers and balances itself with rotations. A heap does not need the full order and does not pay for it.

Completeness of the tree is the second condition: levels are filled from the top down, and the last one from left to right without gaps. So a tree of n nodes is always as low as possible, has no way to degenerate into a stick, and the cell numbers define the links by themselves. If keys have to be extracted from largest to smallest, the comparison is reversed. In Python the heapq module long supported only a heap with the minimum on top, and the largest was extracted by pushing keys with a minus sign; the functions heappush_max, heappop_max and their relatives appeared in version 3.14. A gentler introduction to the heap, through a hospital emergency room, is in the course chapter “Who is next”; here: the structure, the proofs, and what is needed in practice.

Sifting

All heap operations come down to two movements. Up: an element that turned out smaller than its parent is swapped with it until it meets a parent no greater than itself or becomes the root. Down: an element that turned out greater than one of its children is swapped with the smaller child, so that the new parent is no greater than the other child too, until the order is restored or the tree runs out.

• Insert: put at the end of the array (the tree stays complete) and sift up.
• Extract the minimum: take the root, put the last element of the array in its place and sift it down.
• Change a key: write the new one and sift: up if the key decreased, down if it grew.
• Delete from the middle: put the last element into the freed place and sift it up or down, whichever way it does not fit.

The last two operations require knowing which cell the element is in. The heap does not know this: elements move with every swap. So next to the array a second, inverse one is kept, “where it lies”, and updated with every swap. Below are a server's timers: a name and in how many seconds to fire. This is the same heap as in the program timers in the section with the code.

Try cancelling a timer from the middle: the last one in the array takes its place, and whether it goes up or down depends on which branch it was in before. Switch the heap to 4-ary: there will be fewer levels, but on each level four children have to be compared. Then untick “sift after a change” and make some bottom timer the earliest of all. The edge to it turns red, and “Fire” will release the timers out of order: a heap does not check its property, it only maintains it.

Why it is correct, and why a logarithm

There is a third operation that runs faster than one might expect: building a heap from a ready array. Adding the elements one at a time gives n times O(log n). Robert Floyd proposed in 1964 to do it the other way round: sift down every node that has children, starting from the last one and ending with the root. When a node's turn comes, both its subtrees are already heaps, which is the initial condition of the claim above.

Half the nodes are leaves with nowhere to sink; a quarter go down at most one level; and only the root may get the whole path. The sum ⌊n/2⌋ + ⌊n/4⌋ + … can, by the way, be computed exactly: it equals n minus the number of ones in the binary notation of n, which is Legendre's formula for the power of two in n!. For ten timers (1010 in binary) the sum of heights is 8, the comparisons are at most 16, and the program below makes 12.

How to cancel a timer

A textbook heap can insert and extract. A real calendar needs more: a request timeout is usually cancelled because the answer came earlier; a connection with traffic on it postpones its idle timer every time; in the gas every impact invalidates the forecasts of two discs. An element has to be removed from the middle of the heap, and the heap does not know where it is. There are two ways to cope, and both are used in standard libraries.

Mark and forget. The entry is left alone but marked as invalid; when it surfaces at the top, it is silently thrown away. This is lazy deletion. In the gas no mark is even needed: each entry remembers how many times its discs had changed course at the time of the forecast, and if the counter has grown since, the forecast is stale. The calendar in Dijkstra's algorithm is built the same way: a crossroads to which a shorter path has been found is simply pushed into the heap again. The price is memory and extra extractions: in the gas two thirds of the heap's work goes on garbage. So lazy heaps are cleaned from time to time:

• in asyncio a cancelled timer is only marked; the event loop throws such timers off the top of the heap, and when the heap holds more than a hundred timers and more than half of them are cancelled, it builds a new heap from the live ones, in O(n), with heapify;
• in the Go runtime a stopped timer becomes a “zombie” and stays in its processor's heap; zombies are cleared along the way, and if there are more than a quarter of them, the whole heap is gone through;
• the heapq documentation recommends the same recipe: keep a dictionary “task → entry”, on cancellation replace the task in the entry with a REMOVED marker, and on a priority change push a new entry.

Know where it lies. Each element remembers the number of its cell, and every swap updates it. Then cancelling and moving cost O(log n): the last element into the freed place, and sift. In Go the example from the container/heap documentation is built this way: an element has an index field, and the functions heap.Fix and heap.Remove take a cell number. In Java ScheduledThreadPoolExecutor keeps delayed tasks this way: each task records its index in the heap, and, as a comment in the source says, this cuts removal from O(n) to O(log n). By default, though, a cancelled task still lies in the queue until its time comes; immediate removal is a separate setting, available since Java 7.

Code

First the timer heap from the scene. In Go it is built on container/heap: the type implements the five methods of the interface, and Swap also updates the index field. In Python heapq has no way to learn where an element has moved, so the class is written from scratch, and it repeats step by step what container/heap does, down to the order of comparisons. Hence the same output and the same comparison counts.

Moving timer H from 55 seconds to 2 cost three comparisons, one per level as H rose to the root. Cancelling E cost one: the last timer, B31, took E's place, had nowhere to go down, and did not need to go up. The line slots is that inverse array, “where it lies”.

Now the gas from the first scene. The calendar is heapq in Python and container/heap in Go, and an entry is a tuple (time, number, kind of event, discs, their move counters). The entry number is needed so that with equal times the comparison does not reach the other fields. A disc's position is stored as of its last event: between events a disc flies in a straight line, and there is no point recomputing all discs at every step.

That leaves the question of whom to predict collisions with. If with everyone, every impact means going through all the discs; that is how it was computed in 1959. The program divides the box into 10 × 10 cells, no smaller than a disc's diameter: only discs in neighbouring cells can collide. A disc crossing a grid line is also a calendar event: in the new cell the disc has new neighbours, and forecasts have to be made with them. Hence 7826 crossings out of 11,824 events: cheap events that save going through the whole box.

Note the function mul in the Go version. The language specification allows the compiler to fuse a multiplication and an addition into one FMA instruction, which rounds once instead of twice; on ARM processors the compiler does so. The explicit conversion float64(a * b) forbids it. Without it the same program on a Mac with an ARM processor takes 34,081 entries from the heap by t = 10 instead of 34,103 and counts 3037 collisions instead of 3061. The difference starts in the last bit of one number, but the gas is chaotic: every impact amplifies the discrepancy, and after a hundred impacts the history is different. The energy is conserved to the twelfth digit in both versions; physics does not care which of the two histories took place.

What it costs

Memory is the array itself, O(n), with no pointers and no bookkeeping fields; if moving and cancelling are needed, the inverse array “where it lies” is added. The last line of the table is the binary heap's weak spot: two heaps cannot be merged faster than by building a new one. If merging is needed often, there are pointer-based heaps, leftist, binomial, Fibonacci, which merge in O(log n) and faster, but in ordinary work they are slower than an array.

Asymptotics do not tell everything. Below is a measurement of the model used to test event calendars: take the nearest event and schedule a new one a little later, many times in a row, with a constant number n of events in the calendar. Three calendars: a heap; a sorted array, where the nearest lies at the end and a new one is inserted in its place; and an unordered array, where the nearest is found by a scan.

The heap barely notices the growth: from a hundred to a million events the time in Go doubles, in Python also doubles, although there are three times as many levels. The scan grows linearly and on large calendars loses to the heap thousands of times over. The sorted array is the most interesting. Its insertion is linear too, but it is a single memory shift, which the processor does very fast. At a hundred events it beats the heap twofold in Go and slightly in Python, at a thousand it runs almost level, and only beyond that does linear growth take over: five times slower than the heap at ten thousand, hundreds of times at a million. A calendar of a dozen timers does not need a heap; a calendar for a gas of thousands of discs cannot do without one.

Where it is used

The heap is the standard priority queue in almost every language: heapq in Python, container/heap in Go, PriorityQueue in Java (its source says outright: a balanced binary heap, the children of cell n in cells 2n + 1 and 2(n + 1)), std::priority_queue in C++. Dijkstra, Prim and A* rest on it, as do Huffman coding, picking the k best from a stream, merging sorted files in an external sort (the problem on regional catalogues takes it apart), and heapsort. And all event calendars: the timers of asyncio and the Go runtime, Java's delayed tasks, event-driven models from a queue at a bank to molecular dynamics.

The binary heap is a special case of the d-ary heap, where each node has d children: the parent of cell i is ⌊(i − 1)/d⌋, the children are di + 1, …, di + d. There are log_d n levels: sifting up gets cheaper, and sifting down costs d comparisons per level. The Go runtime keeps timers in a 4-ary heap; in the source this is the constant timerHeapN = 4. A large comparison of priority queues carried out in 2014 by Daniel Larkin, Siddhartha Sen and Robert Tarjan showed that the running time of heaps is most closely tied to first-level cache misses: in sorting the fastest was a plain 4-ary heap in an array, in Dijkstra's algorithm on a nationwide road network the 4-ary heap was also among the best, and the Fibonacci heap about three times slower. The “4-ary” switch in the second scene shows where the gain comes from: fewer levels, and neighbouring children lie next to each other in memory.

Alder and Wainwright's list

In 1957 Berni Alder and Thomas Wainwright published the computation with which molecular dynamics begins: a machine traced the motion of a few hundred hard spheres, and the result indicated that even spheres with no attraction to each other pass from a liquid into a crystal. In a 1959 paper they described the method. For each pair of spheres the moment of collision is computed, all spheres are moved to the earliest one, and the two that collided exchange velocities. After that there is no need to recompute everything: it is enough to check the pairs involving one of the two spheres, 2N − 3 of them. Forecasts involving those two spheres are struck out, since, as the paper explains, they were made on the mistaken assumption that the spheres would not collide. This is the same lazy staleness as in our first scene, only the calendar was a list.

The UNIVAC managed about a hundred collisions an hour for a hundred spheres; the IBM 704 about two thousand an hour for a hundred spheres and about five hundred for five hundred. As a guard against machine faults the authors checked the total energy, which must stay constant; the same line is in our program's output. There too they proposed, for the future, dividing the box into cells and checking only neighbours. The JavaScript version of our program, running in the first scene, computes about two hundred thousand collisions a second on an ordinary laptop.

The heap itself appeared in 1964: in June Communications of the ACM published J. W. J. Williams's note “Algorithm 232: Heapsort”, where the heap was described as a data structure in its own right, and in December Robert Floyd printed “Treesort 3” with bottom-up building in linear time. Williams's story is in the course chapter. Ten years later the heap began to be improved again, this time for graphs: in the algorithms of Dijkstra and Prim the main operation is decreasing a key, and Donald Johnson proposed d-ary heaps in 1975, in which this operation is cheaper; the title of the paper says it plainly: “Priority queues with update and finding minimum spanning trees”. In 1984 Michael Fredman and Robert Tarjan invented Fibonacci heaps, where decreasing a key costs O(1) amortized over a sequence of operations, and Dijkstra got the bound O(E + V log V). On paper this is better, but in practice, as the measurement thirty years later showed, the array was faster again.

Problems

In the first problem timers are extended and cancelled, and one of the tests is the real USGS earthquake catalogue. In the second there are two heaps, and both have to be deleted from, lazily; the data are 85 years of Almaty weather.