## Celebrities die 2.7182 at a time

The claim that celebrities die in threes is usually dismissed as the result of the human propensity to see patterns where there are none. But celebrities don’t die at regularly spaced intervals either. It would be very weird if a celebrity predictably died on the 14th of every month. And once you deviate from a regularly spaced pattern, some amount of clustering is inevitable. Can we make this more precise?

Rather than trying to define exactly what constitutes a celebrity, I’ll simply assume that they die at a fixed rate and that they do so independently of each other (The Day the Music Died notwithstanding). It follows that celebrity deaths is a Poisson process with intensity $\lambda$ where $\lambda$ is the number of deaths that occur in some fixed time period.

As an example, suppose we define celebrityhood in such a way that twelve celebrities die each year on average. Then $\lambda = 12/\text{year}$, and because the time between events in a Poisson process is exponentially distributed with parameter $\lambda$, the average time between two deaths is $1/\lambda$ = 1/12th year, or one month.

What does it mean for celebrities to die $n$ at a time? We will simply say that two celebrities die together if the period between their deaths is shorter than expected. If the celebrity death rate is 12/year, then two celebrities died together if their deaths were less than one month apart. Similarly, three celebrities died together if the period between death 1 and death 2 and the period between death 2 and death 3 were both shorter than a month. In general, $k$ celebrities died together if the $k - 1$ periods between their deaths were all shorter than expected.

Here is a diagram of 10 years worth of randomly generated deaths with 12 deaths per year and clusters as defined above highlighted:

Average cluster size
Suppose a celebrity has just died after a longer than average wait. This death will start a new cluster, and we want to figure out what the size of it is.

In a Poisson process the waiting time between two events is exponentially distributed with parameter $\lambda$, so it can be modelled with a stochastic variable $W \sim Exp(\lambda)$. The cluster size itself is modelled with another stochastic variable, $C$, whose distribution is derived as follows.

The cluster size will be 1 when the waiting time for the next death is larger than or equal to the average (which is $1/\lambda$ for the exponential distribution):

$\text{P}(C = 1) = \text{P}(W > 1/\lambda)$

The probability that the cluster will have size 2 is the same as the probability that the next waiting time is shorter than average and the next one after that is longer:

$\text{P}(C = 2) = \text{P}(W \le 1/\lambda)\cdot \text{P}(W > 1/\lambda)$

For size three, it’s the probability that the next two waiting times are shorter and the third one longer:

$\text{P}(C = 3) = \text{P}(W \le 1/\lambda)^2\cdot \text{P}(W > 1/\lambda)$

In general, the probability that the next cluster will be size $k$ is:

$\text{P}(C = k) = \text{P}(W \le 1/\lambda)^{k - 1}\cdot \text{P}(W > 1/\lambda)$

So what’s the average size of a Celebrity Death Cluster? The expected value of $C$ is given by:

$\displaystyle \text{E}[C] = \sum_{k=1}^\infty k \cdot \text{P}(C = k) = \sum_{k=1}^\infty k\cdot \text{P}(W \le 1/\lambda)^{k - 1}\cdot \text{P}(W > 1/\lambda)$

Plugging in the distribution function for the exponential distribution, we get:

\begin{align*} \text{E}[C] &= \sum_{k=1}^\infty k \cdot (1 - e^{- \lambda \cdot (1/\lambda) })^{k - 1} \cdot (1 - (1 - e^{- \lambda \cdot (1 / \lambda)}))\\ &= \sum_{k=1}^\infty k \cdot (1 - e^{- 1})^{k - 1} \cdot e^{-1} \end{align*}

It’s not hard to show that this infinite series has sum $e$ (Hint: Use the fact that $k x^{k - 1}$ is the derivative of $x^k$), so on average, celebrities die 2.7182 at a time.

## The Radix Heap

The Radix Heap is a priority queue that has better caching behavior than the well-known binary heap, but also two restrictions: (a) that all the keys in the heap are integers and (b) that you can never insert a new item that is smaller than all the other items currently in the heap.

These restrictions are not that severe. The Radix Heap still works in many algorithms that use heaps as a subroutine: Dijkstra’s shortest-path algorithm, Prim’s minimum spanning tree algorithm, various sweepline algorithms in computational geometry.

Here is how it works. If we assume that the keys are 32 bit integers, the radix heap will have 33 buckets, each one containing a list of items. We also maintain one global value last_deleted, which is initially MIN_INT and otherwise contains the last value extracted from the queue.

The invariant is this:

The items in bucket $k$ differ from last_deleted in bit $k - 1$, but not in bit $k$ or higher. The items in bucket 0 are equal to last_deleted.

For example, if we compare an item from bucket 10 to last_deleted, we will find that bits 31–10 are equal, bit 9 is different, and bits 8–0 may or may not be different.

Here is an example of a radix heap where the last extracted value was 7:

As an example, consider the item 13 in bucket 4. The bit pattern of 7 is 0111 and the bit pattern of 13 is 1101, so the highest bit that is different is bit number 3. Therefore the item 13 belongs in bucket $3 + 1 = 4$. Buckets 1, 2, and 3 are empty, but that’s because a number that differs from 7 in bits 0, 1, or 2 would be smaller than 7 and so isn’t allowed in the heap according to restriction (b).

Operations
When a new item is inserted, it has to be added to the correct bucket. How can we compute the bucket number? We have to find the highest bit where the new item differs from last_deleted. This is easily done by XORing them together and then finding the highest bit in the result. Adding one then gives the bucket number:

bucket_no = highest_bit (new_element XOR last_deleted) + 1


where highest_bit(x) is a function that returns the highest set bit of x, or $-1$ if x is 0.

Inserting the item clearly preserves the invariant because the new item will be in the correct bucket, and last_deleted didn’t change, so all the existing items are still in the right place.

Extracting the minimum involves first finding the minimal item by walking the lowest-numbered non-empty bucket and finding the minimal item in that bucket. Then that item is deleted and last_deleted is updated. Then the bucket is walked again and all the items are redistributed into new buckets according to the new last_deleted item.

The extracted item will be the minimal one in the data structure because we picked the minimal item in the redistributed bucket, and all the buckets with lower numbers are empty. And if there were a smaller item in one of the buckets with higher numbers, it would be differing from last_deleted in one of the more significant bits, say bit $k$. But since the items in the redistributed bucket are equal to last_deleted in bit $k$, the hypothetical smaller item would then have to also be smaller than last_deleted, which it can’t be because of restriction (b) mentioned in the introduction. Note that this argument also works for two-complement signed integers.

We have to be sure this doesn’t violate the invariant. First note that all the items that are being redistributed will satisfy the invariant because they are simply being inserted. The items in a bucket with a higher number $k$ were all different from the old last_deleted in the $(k-1)$th bit. This bit must then necessarily also be different from the $(k-1)$th bit in the new last_deleted, because if it weren’t, the new last_deleted would itself have belonged in bucket $k$. And finally, since the bucket being redistributed is the lowest-numbered non-empty one, there can’t be any items in a bucket with a lower number. So the invariant still holds.

In the example above, if we extract the two ‘7’s from bucket 0 and the ‘8’ from bucket 4, the new heap will look like this:

Notice that bucket 4, where the ‘8’ came from, is now empty.

Performance
Inserting into the radix heap takes constant time because all we have to do is add the new item to a list. Determining the highest set bit can be done in constant time with an instruction such as bsr.

The performance of extraction is dominated by the redistribution of items. When a bucket is redistributed, it ends up being empty. To see why, remember that all the items are different from last_deleted in the $(k - 1)$th bit. Because the new last_deleted comes from bucket $k$, the items are now all equal to last_deleted in the $(k - 1)th$ bit. Hence they will all be redistributed to a lower-numbered bucket.

Now consider the life-cycle of a single element. In the worst case it starts out being added to bucket 31 and every time it is redistributed, it moves to a lower-numbered bucket. When it reaches bucket 0, it will be next in line for extraction. It follows that the maximum number of redistributions that an element can experience is 31.

Since a redistribution takes constant time per element distributed, and since an element will only be redistributed $d$ times, where $d$ is the number of bits in the element, it follows that the amortized time complexity of extraction is $O(d)$. In practice we will often do better though, because most items will not move through all the buckets.

Caching performance
Some descriptions of the radix heap recommend implementing the buckets as doubly linked lists, but that would be a mistake because linked lists have terrible cache locality. It is better to implement them as dynamically growing arrays. If you do that, the top of the buckets will tend to be hot which means the per-item number of cache misses during redistribution of a bucket will tend to be $O(1/B)$, where $B$ is the number of integers in a cache line. This means the amortized cache-miss complexity of extraction will be closer to $O(d/B)$ than to $O(d)$.

In a regular binary heap, both insertion and extraction require $\Theta(\log n)$ swaps in the worst case, and each swap (except for those very close to the top of the heap) will cause a cache miss.

In other words, if $d = \Theta(\log n)$, extraction from a radix heap will tend to generate $\Theta(\log n / B)$ cache misses, where a binary heap will require $\Theta(\log n)$.

## Porter/Duff Compositing and Blend Modes

In the Porter/Duff compositing algebra, images are equipped with an alpha channel that determines on a per-pixel basis whether the image is there or not. When the alpha channel is 1, the image is fully there, when it is 0, the image isn’t there at all, and when it is in between, the image is partially there. In other words, the alpha channel describes the shape of the image, it does not describe opacity. The way to think of images with an alpha channel is as irregularly shaped pieces of cardboard, not as colored glass. Consider these two images:

When we combine them, each pixel of the result can be divided into four regions:

One region where only the source is present, one where only the destination is present, one where both are present, and one where neither is present.

By deciding on what happens in each of the four regions, various effects can be generated. For example, if the destination-only region is treated as blank, the source-only region is filled with the source color, and the ‘both’ region is filled with the destination color like this:

The effect is as if the destination image is trimmed to match the source image, and then held up in front of it:

The Porter/Duff operator that does this is called “Dest Atop”.

There are twelve of these operators, each one characterized by its behavior in the three regions: source, destination and both. The ‘neither’ region is always blank. The source and destination regions can either be blank or filled with the source or destination colors respectively.

The formula for the operators is a linear combination of the contents of the four regions, where the weights are the areas of each region:

$A_\text{src} \cdot [s] + A_\text{dest} \cdot [d] + A_\text{both} \cdot [b]$

Where $[s]$ is either 0 or the color of the source pixel, $[d]$ either 0 or the color of the destination pixel, and $[b]$ is either 0, the color of the source pixel, or the color of the destination pixel. With the alpha channel being interpreted as coverage, the areas are given by these formulas:

$A_\text{src} = \alpha_\text{s} \cdot (1 - \alpha_\text{d})$
$A_\text{dst} = \alpha_\text{d} \cdot (1 - \alpha_\text{s})$
$A_\text{both} = \alpha_\text{s} \cdot \alpha_\text{d}$

The alpha channel of the result is computed in a similar way:

$A_\text{src} \cdot [\text{as}] + A_\text{dest} \cdot [\text{ad}] + A_\text{both} \cdot [\text{ab}]$

where $[\text{as}]$ and $[\text{ad}]$ are either 0 or 1 depending on whether the source and destination regions are present, and where $[\text{ab}]$ is 0 when the ‘both’ region is blank, and 1 otherwise.

Here is a table of all the Porter/Duff operators:

 $[\text{s}]$ $[\text{d}]$ $[\text{b}]$ Src $s$ $0$ s Atop $0$ $d$ s Over $s$ $d$ s In $0$ $0$ s Out $s$ $0$ $0$ Dest $0$ $d$ d DestAtop $s$ $0$ d DestOver $s$ $d$ d DestIn $0$ $0$ d DestOut $0$ $d$ $0$ Clear $0$ $0$ $0$ Xor $s$ $d$ $0$

And here is how they look:

Despite being referred to as alpha blending and despite alpha often being used to model opacity, in concept Porter/Duff is not a way to blend the source and destination shapes. It is way to overlay, combine and trim them as if they were pieces of cardboard. The only places where source and destination pixels are actually blended is where the antialiased edges meet.

Blending
Photoshop and the Gimp have a concept of layers which are images stacked on top of each other. In Porter/Duff, stacking images on top of each other is done with the “Over” operator, which is also what Photoshop/Gimp use by default to composite layers:

Conceptually, two pieces of cardboard are held up with one in front of the other. Neither shape is trimmed, and in places where both are present, only the top layer is visible.

A layer in these programs also has an associated Blend Mode which can be used to modify what happens in places where both are visible. For example, the ‘Color Dodge’ blend mode computes a mix of source and destination according to this formula:

$\begin{equation*} B(s,d)= \begin{cases} 0 & \text{if $$d=0$$,} \\ 1 & \text{if $$d \ge (1 - s)$$,} \\ d / (1 - s) & \text{otherwise} \end{cases} \end{equation*}$

The result is this:

Unlike with the regular Over operator, in this case there is a substantial chunk of the output where the result is actually a mix of the source and destination.

Layers in Photoshop and Gimp are not tailored to each other (except for layer masks, which we will ignore here), so the compositing of the layer stack is done with the source-only and destination-only region set to source and destination respectively. However, there is nothing in principle stopping us from setting the source-only and destination-only regions to blank, but keeping the blend mode in the ‘both’ region, so that tailoring could be supported alongside blending. For example, we could set the ‘source’ region to blank, the ‘destination’ region to the destination color, and the ‘both’ region to ColorDodge:

Here are the four combinations that involve a ColorDodge blend mode:

In this model the original twelve Porter/Duff operators can be viewed as the results of three simple blend modes:

 Source: $B(s, d) = s$ Dest: $B(s, d) = d$ Zero: $B(s, d) = 0$

In this generalization of Porter/Duff the blend mode is chosen from a large set of formulas, and each formula gives rise to four new compositing operators characterized by whether the source and destination are blank or contain the corresponding pixel color.

Here is a table of the operators that are generated by various blend modes:

The general formula is still an area weighted average:

$A_\text{src} \cdot [s] + A_\text{dest} \cdot [d] + A_\text{both}\cdot B(s, d)$

where [s] and [d] are the source and destination colors respectively or 0, but where $B(s, d)$ is no longer restricted to one of $0$, $s$, and $d$, but can instead be chosen from a large set of formulas.

The output of the alpha channel is the same as before:

$A_\text{src} \cdot [\text{as}] + A_\text{dest} \cdot [\text{ad}] + A_\text{both} \cdot [\text{ab}]$

except that [ab] is now determined by the blend mode. For the Zero blend mode there is no coverage in the both region, so [ab] is 0; for most others, there is full coverage, so [ab] is 1.