Coordinate compression

The GeoNames database has 45,227 cities with ten thousand people or more. Each one's latitude and longitude are recorded to a thousandth of a degree, and if you set up a table “how many people live at this point” with a cell for every possible pair of coordinates, it would have 180,001 × 360,001 cells, sixty-five billion, almost all empty. There are only 45 thousand cities. Below are all of them. Select a rectangle, and the scene counts how many cities and people are in it. Then press “Compressed”.

The map is redrawn. Each city moves from its real latitude to that latitude's rank among all distinct latitudes, and from its longitude to the longitude's rank. The cities have 32,407 distinct latitudes and 38,493 longitudes, and each gets a strip of the same width. So strips swell where there are many cities anywhere in the world: the longitudes from Gibraltar to the Volga, Hindustan, the east of Asia, the latitudes from the tropics to the fifty-fifth parallel. The Pacific, where for thousands of kilometres there is not a single longitude with a city, shrinks into narrow strips at the left and right edges of the map, the north beyond the Arctic Circle and the south beyond the fortieth parallel into thin strips too, Antarctica into a line along the bottom edge. Yet the answer about a rectangle does not change: in the rectangle from the Atlantic to the Dnieper both maps have 10,061 cities and 470,893,364 people. The numbers are different, the order is the same.

The idea

Many algorithms care not about the values themselves but only about which of them is larger. Then each value can be replaced by its place in the sorted list of distinct values: the smallest becomes zero, the next one, and so on. The ranks run consecutively from 0 to m − 1, where m is the number of distinct values, and ordinary arrays, Fenwick trees and segment trees can be set up over them, however huge the original numbers.

This is coordinate compression. The word “coordinates” is conventional: anything that can be compared is compressed, dates, sums in tenge, map node IDs, strings. In statistics the same thing is called a rank transform, in SQL the function DENSE_RANK, in pandas factorize.

The recipe

Three steps. Collect all values that will occur, both in the data and in the queries, if the queries are known in advance. Sort and drop duplicates. The rank of a value is its index in the resulting list; turning a rank back into a value is just indexing the list.

A rank can be found in two ways. A dictionary “value → rank” answers in O(1) on average, but knows only the values present at compression time. Binary search in the sorted list, bisect_left in Python, sort.SearchInts in Go, works in O(log m) and suits any number: it answers how many distinct values are smaller than a given one. For a value from the list that is its rank, and for a foreign one, the rank of the place where it would go. The last line of the listing checks that both ways agree on all 45,227 latitudes.

Along the way the listing reports that Almaty at latitude 43.252 gets rank 25,008 of 32,407, and that 33,996 cities of the world, 75.2%, lie south of it. A rank in a sorted list is the same as a place in a ranking, and compression answers questions about shares and places at once.

In numpy one function does all three steps: np.unique(a, return_inverse=True) returns both the list of distinct values and the rank of each element. In pandas, pd.factorize(a, sort=True); without sort=True ranks are handed out in order of first appearance, an order unrelated to the values, and ranks can no longer be compared with each other.

Why the answers do not change

Let x₀ < x₁ < … < x_m−1 be the distinct values, and id(v) the rank of the value v, that is, the i with x_i = v. For values not in the list, define two functions: lb(a), the number of distinct values smaller than a, and ub(b), the number of distinct values no greater than b. These are bisect_left and bisect_right.

The main consequence follows from the first part: any algorithm that looks at values only through the comparisons “less” and “equal” makes exactly the same steps on ranks and gives the same answer. Sorting, searching, counting inversions, the longest increasing subsequence, a Fenwick tree over values are all like that. The second part translates a query with bounds not present in the data into the language of ranks: the rectangle “latitude from 42.252 to 44.252” becomes “latitude ranks from lb(42,252) to ub(44,252) − 1”. The scene shows this translation under the map.

What the theorem does not promise is visible from its wording: it says nothing about differences and sums of values. Ranks 7 and 8 may stand for the numbers 10 and 11, or for 10 and 10⁹. An algorithm that needs distances, lengths, areas, “is it within a kilometre”, will go wrong on ranks. We will come back to this in the section on the trap.

Two dimensions: one compression is not enough

Back to rectangles. If both axes are compressed, the table “how many people per cell” shrinks from 6.5·10¹⁰ cells to 32,407 × 38,493 ≈ 1.25·10⁹. Fifty times smaller, but at eight bytes per cell that is still ten gigabytes, almost all zeros. In two dimensions compression removes empty rows and columns, but not empty cells.

So one axis is compressed, and the other is swept. Sort the cities from west to east and add their populations one at a time into a Fenwick tree whose cells are latitude ranks. When the tree holds all cities west of some longitude, the sum over latitude ranks from lo to hi − 1 is the population of the strip “west of this longitude, between these latitudes”. A rectangle is the difference of two such strips: up to its eastern edge and up to its western edge. Each rectangle sets two moments of the sweep, and if these moments are sorted together with the cities, one pass is enough for all rectangles at once. The tree is 32,407 cells of eight bytes, a quarter of a megabyte.

rectangle                            latitude            longitude  towns      people
a degree around Almaty       42.252 ..  44.252     75.912 ..  77.912     10     2217955
the Nile delta               29.800 ..  31.600     29.800 ..  32.400    137    36672070
the island of Java           -8.800 ..  -5.800    105.100 .. 114.600    243    51941160
the whole world             -90.000 ..  90.000   -180.000 .. 180.000  45227  4210365101

A two-by-two-degree square around Almaty holds 10 cities and 2.2 million people, a rectangle around the island of Java 243 cities and almost 52 million. All 45,227 cities, according to GeoNames, are home to 4.2 billion people, about half of humanity. The numbers are checked by brute force: at the end the listing compares each answer with a direct count.

A sweep with deferred questions is a special case of a general technique: if all queries are known in advance, they can be asked not in the order they came but in the order convenient for answering. The sweep itself is a sweep line; other techniques of the same kind are reading queries from the end in the article on disjoint set union and Mo's algorithm in square root decomposition.

Code

Both listings above are written the way people write in real life. In Python compression is sorted(set(values)) and a dictionary from enumerate; in Go, slices.Sort, then slices.Compact, which throws out neighbouring duplicates, and a map. Turning query bounds into ranks is bisect_left and bisect_right in Python, sort.SearchInts(ys, a) and sort.SearchInts(ys, b + 1) in Go: for integers “the first one greater than b” is the same as “the first one no less than b + 1”.

Two small things people get wrong most often. First: everything that will be compared must be compressed, both the data and the query bounds, if an array is set up over them. Second: a dictionary does not know the rank of a new value that was absent at compression time, and binary search gives the rank of a neighbour. If something is later written at that rank, the write lands in someone else's cell.

What it costs

Compression is a sort, O(n log n), and a pass to drop duplicates, O(n). The rank of each of the n values is another O(n) with a dictionary or O(n log m) with binary search. Memory is the list of distinct values and an array of ranks, O(n). It pays off where a structure over values is needed next: a Fenwick tree over ranks weighs O(m), not O(span of the numbers), and every operation on it costs O(log m).

But it does not pay off at once. One rectangle is faster to count by brute force: 45,227 checks against sorting all the cities. The measurement below is q random rectangles on the same cities, two ways. The axes are logarithmic.

Brute force grows in exact proportion to q. The sweep at first stays flat: almost all its time is sorting and the pass over the cities, and rectangles add almost nothing. In Python the curves cross at about 35 rectangles, and in Go at about 250: brute force in Go spends less than a nanosecond per city and is harder to beat. Beyond that the gap grows with q: at a thousand rectangles the sweep in Python is 23 times faster than brute force, at four thousand in Go 13 times. This is the general rule of data preparation: it costs about one or two passes over the data and pays off when the prepared data are walked many times.

The trap: vanishing gaps

A classic competition problem. Posters are pasted one after another along a fence, each new one covering the old ones. The poster ends are numbers up to a billion. How many posters are visible, at least a little? Painting the fence metre by metre is impossible: it is too long. Compressing the poster ends and painting an array of ranks is possible. But look at what comes of it.

Three posters from the listing: [1, 10], and on top of it [1, 4] and [6, 10]. All three are visible on the fence: the fifth metre of the first one is left showing. Compress the ends: the values 1, 4, 6, 10 get ranks 0, 1, 2, 3, the second poster occupies cells 0–1, the third 2–3, and the first is covered entirely. The fifth metre between ranks 1 and 2 simply does not exist: compression preserved the order of the ends, but not what lay between them. This is exactly what the theorem did not promise.

The cure is to compress not points but elementary segments. A poster [l, r] occupies the metres from l to r inclusive, that is, the half-open interval [l, r + 1). Cut the fence at all points l and r + 1. Each piece between neighbouring cuts is either entirely under a poster or entirely outside it: inside a piece no poster begins or ends. So a piece can be treated as one cell, and painting the cells reproduces the painting of the fence exactly. In the example the cuts are 1, 5, 6, 11, there are three cells, and the middle one, [5, 6), is that fifth metre.

On two thousand random posters on a fence a billion millimetres long, naive compression counted 642 visible posters, and correct compression 708. Sixty-six posters are visible only in gaps between other posters' ends, and naive compression lost them. Such an answer does not crash and looks plausible, which is why this mistake lives long.

The same reasoning gives a rule for problems where lengths matter: keep the values themselves next to the ranks. Cell i is the segment from x_i to x_i+1, and its length x_i+1 − x_i has not gone anywhere, it just lies in another array. This is how the length of a union of segments and the area of a union of rectangles are computed: go through the cells and add up the lengths of the covered ones.

Three smaller traps.

• Equal values. Compression glues equal values into one rank. If an algorithm needs strict inequalities, “strictly increasing subsequence”, “strictly less”, that is what is needed. If equal elements must be told apart, pairs (value, index) are compressed.
• Fractional numbers. They can be compressed too, but equality of fractions is treacherous: 0.1 + 0.2 is not equal to 0.3, and two “identical” values computed differently get different ranks. Before compression fractions are rounded to the needed precision or turned into integers, as the coordinates of cities are turned into thousandths of a degree here.
• Values that arrive later. Ranks are handed out once for the whole list. A new value between x₅ and x₆ has nowhere to go: everything larger would have to be renumbered.

When values arrive one at a time

The last trap is the technique's main limitation: compression requires knowing all values in advance. In competition problems that is almost always so: the input is read whole, and so are all queries. In a live system where events come as a stream, people manage differently. You can keep a structure not over ranks but over the values themselves, but create segment tree nodes only where something lies: such an “implicit” tree on the range [0, 10⁹] spends about 30 nodes per insertion. You can take a balanced search tree that keeps the values themselves and can count how many of them are smaller than a given one. Both ways are slower than an array, but they do not require knowing the future.

Columnar databases know the same difficulty. They keep a column's repeated values, city names for example, as codes from a dictionary, and if the dictionary preserves order, sorting and searching can be done directly on the codes. Carsten Binnig, Stefan Hildenbrand and Franz Färber, in a 2009 paper on such dictionaries, note the same trouble: fixed-length order-preserving codes are hard to extend with new values.

Where it is used

• Maps and routes. OpenStreetMap node IDs have long exceeded 2³²; a comment in the source of the OSRM router says so, and an OSM node ID takes 64 bits there. Internally OSRM numbers nodes consecutively with 32-bit numbers, and the graphs for path finding are kept in arrays over these numbers. An array with a cell for every original ID, and there are more than four billion of them, would take tens of gigabytes.
• Data and databases. np.unique(…, return_inverse=True), pd.factorize, pandas categorical columns, the window function DENSE_RANK in SQL, dictionary encoding of columns in columnar stores: all of these compress values into ranks.
• Statistics. Spearman's rank correlation is an ordinary correlation computed on the ranks of values instead of the values themselves. That is why outliers do not throw it off: a billionaire in the sample gets a rank one higher than the neighbour's, not a thousand times the income. Histogram equalisation of a photograph works similarly: a pixel's brightness is replaced by the share of pixels no brighter than it.
• Competitions. A limit of “coordinates up to 10⁹, points up to 2·10⁵” is almost always a hint: a Fenwick or segment tree is built over ranks. This is how inversions are counted (Fenwick tree), the longest increasing subsequence in O(n log n), the area of a union of rectangles (sweep line).

History

The technique has no author: it has been used ever since arrays were indexed by something other than a loop counter. But its best-known early record is tied to Klee's problem. In 1977 Victor Klee asked in a journal note whether the length of a union of n segments can be found faster than in O(n log n). In the same year Jon Bentley, in unpublished notes “Solutions to Klee's rectangle problems”, solved the two-dimensional version, the area of a union of rectangles, with a sweep line. For it Bentley devised the segment tree, and the segment tree in its classic form is built not over coordinates but over elementary segments between sorted ends, those very pieces that saved the posters in the section on the trap. Michael Fredman and Bruce Weide showed in 1978 that for segments nothing faster than O(n log n) is possible in the natural model of computation.

Problems

The first is rectangles on a map, as in the listing, but with tens of thousands of queries, so brute force will not fit in the time. The second is posters on a fence; the most natural solution loses the gaps, and the tests will notice.