## 2010-06-10

### You are Here

I am a terrible person: I read code, I think to myself "I can do better" and I set out thus, out of pride. I did just that with an OCaml implementation of the geohash geographical encoding. Even if the code I read was perfectly workable, I found it unidiomatic and somewhat bloated. I wrote two versions, the first one compact but impenetrable, the second one from first principles: what would be an executable description of the algorithm? To wit:

```let decode =
(narrow (-90., 90.) *** narrow (-180., 180.)) % swap % split
% List.concat % List.map (to_bits 5 % decode_base32) % explode
```

Because of (or despite!) the point-free, combinator-heavy style I intended this to be read as a pipeline or as a recipe, from right to left (it is an applicative-order pipeline), much as the algorithm description does:

1. Explode the code string into a list of characters
2. For each character, decode it as a base-32 digit decomposing it into a list of 5 bits
3. Flatten the list of lists of bits into a big list of bits
4. Split the list of bits into a pair of lists for the even bits and the odd bits
5. Swap both components of the pair, since the coding uses the even bits for the longitude
6. Narrow each list independently into a coordinate starting from the corresponding interval

The result is a (latitude, longitude) pair. Supposing that every step is invertible, the algorithm for encoding a coordinate into a geohash is:

```let encode nchars =
let bits = 5 * nchars in
let lo = bits / 2 in
let hi = bits - lo in
implode % List.map (encode_base32 % of_bits) % group 5 % join
% swap % (expand lo (-90., 90.) *** expand hi (-180., 180.))
```

Again, this is easy to follow as a recipe or as a succint description:

1. Given the desired hash length, compute the number of bits for the latitude and for the longitude
2. Expand each coordinate independently into a list of bits describing its position in the corresponding interval
3. Swap both components of the pair
4. Join both lists by alternating even and odd elements
5. Group the bits in the list into sublists by taking 5 bits at a time
6. For each group, convert the list of 5 bits into an integer and encode it as a base-32 digit
7. Implode the list of characters into a code string

Now I have a number of inverse or near-inverse pairs of functions to write, namely:

• `explode` and `implode`
• `decode_base32` and `encode_base32`
• `to_bits` and `of_bits`
• `List.concat` and `group`
• `split` and `join`
• `narrow` and `expand`

Let's start with the (well-known to Haskellers) combinators (though they pronounce `%` as `.`):

```let ( \$ ) f x = f x
let ( % ) f g x = f (g x)
```

The first is the application operator and the second is the composition operator. Next are some operators on pairs borrowed from `Control.Arrow`:

```let first f (x, y) = (f x, y)
let second f (x, y) = (x, f y)
let ( *** ) f g (x, y) = (f x, g y)
let swap (x, y) = (y, x)
```

They allow applying a function to one or both members of a pair, and to manipulate its components. Note that `swap``swap``id`. This completes the functional scaffolding. The first and more imperative pair of functions is `explode` and `implode`:

```let explode s =
let res = ref [] in
for i = String.length s - 1 downto 0 do
res := String.unsafe_get s i :: !res
done;
!res

let implode l =
let b = Buffer.create 10 in
Buffer.contents b
```

The first traverses the string from right to left to accumulate each character on the head of the resulting list; the second uses a string buffer to build the resulting string a character at a time. I use `unsafe_get` because, well, these are as low-level functions as they come, so squeezing a bit extra from them doesn't seem inappropriate. Each is the other's inverse, or in more abstract terms the pair `(explode, implode)` witnesses the monoid isomorphism between OCaml strings and lists of characters (that's why Haskell identifies both datatypes).

Aside: I won't attempt any proofs to avoid making a long post even longer. Some are easy, most are tedious, and the preceding is probably difficult. I've verified a few of the inverse pairs but not all of them, so I don't consider this code verified by any stretch of the word. I am, however, confident that the code works.

Now `decode_base32`and `encode_base32` are also fairly imperative, but just because I chose to use an array as the map between characters and their base-32 encodings:

```let codetable = [|
'0'; '1'; '2'; '3'; '4'; '5'; '6'; '7';
'8'; '9'; 'b'; 'c'; 'd'; 'e'; 'f'; 'g';
'h'; 'j'; 'k'; 'm'; 'n'; 'p'; 'q'; 'r';
's'; 't'; 'u'; 'v'; 'w'; 'x'; 'y'; 'z'
|]
```

This means that for encoding a 5-bit number indexing into the array is enough, but decoding a base-32 character requires a lookup. Since the array is sorted by character, a binary search is a good choice:

```let binsearch ?(compare=Pervasives.compare) v e =
let i = ref 0
and j = ref (Array.length v) in
if !j = 0 then raise Not_found else
while !i <> !j - 1 do
let m = (!i + !j) / 2 in
if compare v.(m) e <= 0
then i := m
else j := m
done;
if v.(!i) = e then !i else raise Not_found

let encode_base32 i = codetable.(i land 31)
and decode_base32 c = binsearch codetable (Char.lowercase c)
```

(This binary search is completely general). Now `decode_base32``encode_base32``id` (for the restricted domain of 5-bit integers) if and only if `binsearch v v.(i)` = i for i in range of v. Similarly, `encode_base32``decode_base32``id` (again, for the restricted domain of base-32 characters) if and only if `v.(binsearch v x)` = x. The conditions on `binsearch` amount to it being correct (and would make for excelent unit or random testing).

The next inverse pair is `to_bits` and `of_bits`:

```let iota =
let rec go l i =
if i <= 0 then l else let i = i - 1 in go (i :: l) i
in go []

let to_bits k n = List.rev_map (fun i -> (n land (1 lsl i)) != 0) \$ iota k

let of_bits = List.fold_left (fun n b -> 2 * n + if b then 1 else 0) 0
```

The first maps a list of bit positions (given by `iota`, the initial natural interval) into the corresponding bits by direct testing, most-significant (highest) bit first; the second computes the binary value of the list of bits by a Horner iteration.

The following functions are very general functions on lists. Any function f : α `list`α `list list` such that `List.concat`f`id` is called a partition. One such function is `group`:

```let group k l =
let rec go gs gss i = function
| []           -> List.rev (if gs = [] then gss else List.rev gs :: gss)
| l when i = 0 -> go [] (List.rev gs :: gss) k l
| x :: xs      -> go (x :: gs) gss (i - 1) xs
in go [] [] k l
```

It collects up to k elements in the stack gs, pushing complete groups into the stack gss. The last group might be incomplete and it is added to the top of the stack before returning it. A function that I will need later is the the converse of `List.fold_right`, appropriately called `unfold`:

```let rec unfold f e n =
if n <= 0 then [] else
let (x, e') = f e in
x :: unfold f e' (n - 1)
```

`unfold` generates a list of the given length n from a seed e and a function f that computes the head and the next seed. Now `split` and `join` are another very functional couple:

```let cons x xs = x :: xs

let split l =
let rec even = function
| [] -> [], []
| x :: xs -> first (cons x) \$ odd xs
and odd = function
| [] -> [], []
| x :: xs -> second (cons x) \$ even xs
in even l

let rec join (xs, ys) = match xs with
| [] -> ys
| x :: xs -> x :: join (ys, xs)
```

In general, `split``join` is not the identity on pairs of lists (think of `([], l)` for nonempty l), but indeed `join``split``id`, and a proof could proceed by induction. Now all these functions are completely general (except perhaps for the base-32 encoding and decoding) and have little if anything to do with geocaching. The last pair, `narrow` and `expand`, are the meat of the algorithm:

```let avg (x, y) = (x +. y) /. 2.

let bisect (lo, hi as interval) b =
let md = avg interval in
if b then (md, hi) else (lo, md)

let narrow interval = avg % List.fold_left bisect interval
```

`narrow` repeatedly `bisect`s the interval, keeping the upper or lower half depending on the next bit in the input list, finally returning the midpoint of the final interval. Its inverse is `expand`:

```let divide x (lo, hi as interval) =
let md = avg interval in
if x > md
then true , (md, hi)
else false, (lo, md)

let expand bits interval x = unfold (divide x) interval bits
```

It `unfold`s each bit in turn, by `divide`ing the interval and finding which half contains x. These functions are inverses in the sense that, for any interval i = (min, max) and list l, `expand (List.length l) i (narrow i l)` = l and conversely, for any n ≥ 0 and minxmax:

```|`narrow i (expand n i x)` - x| ≤ (max - min)2-n-1
```

These are the crucial properties that make `encode` the true inverse of `decode`. The astute reader might have noticed the resemblance of geohashing with arithmetic coding: in a sense, a geohash is the base-32 encoding of the coordinates lossly compressed with an uniform probability model. One way to improve the precision of the encoding would be by allocating a predefined weight to each hemispherical quadrant, by taking into account that more landmass is present in the north-eastern quadrant of the globe. This would require splitting the intervals not along the middle but, say, the first third, thus representing more precise locations with the same number of bits in continental Europe and the Far East, then in North America, then in Africa and Oceania, and lastly in South America.

Jonathan said...

The usual (subrange of) character-as-int indexed lookup table is probably faster, and definitely easier to prove correct than binary search.

I don't mind that only in ASCII/Unicode/etc. is the array actually sorted, though.

Matías Giovannini said...

@Jonathan, indeed, such an inverse mapping can be built as:

let invtable =
let t = Array.make 256 (-1) in
(Array.iteri (fun i c -> t.(Char.code c) <- i) codetable; t)

which directly satisfies the conditions on encode_int32 and decode_int32.