A SpatialPrefixTree based on the Hilbert Curve and variable grid sizes [LUCENE-4922]

[asfimport]

My wish-list for an ideal SpatialPrefixTree has these properties:

• Hilbert Curve ordering
• Variable grid size per level (ex: 256 at the top, 64 at the bottom, 16 for all in-between)
• Compact binary encoding (so-called "Morton number")
• Works for geodetic (i.e. lat & lon) and non-geodetic

Some bonus wishes for use in geospatial:

• Use an equal-area projection such that each cell has an equal area to all others at the same level.
• When advancing a grid level, if a cell's width is less than half its height. then divide it as 4 vertically stacked instead of 2 by 2. The point is to avoid super-skinny cells which occurs towards the poles and degrades performance.

All of this requires some basic performance benchmarks to measure the effects of these characteristics.

---

Migrated from LUCENE-4922 by David Smiley (@dsmiley), 2 votes, updated Aug 22 2014
Attachments: HilbertConverter.zip, LUCENE-4922.patch

[asfimport]

Hatim Hakeel (migrated from JIRA)

Hi David,

I believe this would require to implement a class that extends the SpatialPrefixTree class. It should support all porperties mentioned. We get it to perform Hilbert space filling curve based spatial data structure traversal in making range queries etc.
Apart from performance benchmarks to measure effects of the new feature, are the bonus wishes for use in geospatial to be considered as optimizations given availability of time?

Thanks,
Hatim.

[asfimport]

David Smiley (@dsmiley) (migrated from JIRA)

Thanks for your interest Hatim!

Your assessment of just needing to extend SpatialPrefixTree is mostly correct. It also involves implementing Cell (formerly known as Node) which is tightly related. There are lots of potential optimizations here. As part of this task, modifying SpatialPrefixTree & Cell itself is absolutely okay if we find better (more optimized) ways of doing things, though we need to retain backwards compatibility with the two existing SpatialPrefixTree implementations – geohash & quad. If retaining that backwards compatibility turns out to be too constraining on making this SPT be even better, then we could make SPT an interface. But the general API of SPT & Cell should be the same as it is now.

What will help in this task is to know that there are extensive tests of the PrefixTree based spatial search strategies, such that if supplying a new SPT causes a failure, then you know you've got a bug ;-). But there are few unit level tests of the SPT itself, and so knowing precisely what is broken could be hard.

You are correct that the bonus wishes are purely optional optimizations, time permitting. Very optional; in general I want to see that without these optimizations that the SPT is as good as it can reasonably be. It may turn out that these geo specific optimizations have extra runtime costs (such as calculating the sin() a lot or complexity that make it not worthwhile; I don't know.

[asfimport]

Hatim Hakeel (migrated from JIRA)

Hi David,

Thanks for the feedback. I have drafted a proposal based on your comments and stuff I figured out by digging into the background.

https://google-melange.appspot.com/gsoc/proposal/review/google/gsoc2013/hatimhakeel/1

Hope it satisfies the requirements of the project.

Thanks,
Hatim.

[asfimport]

David Smiley (@dsmiley) (migrated from JIRA)

Thanks for the proposal Hatim! That's a good proposal. I enjoyed reading some of the research papers you referenced. I'm excited to think that at the end of the summer, there's going to be an awesome optimized PrefixTree; and I'm a tad jealous that a lucky student is going to build it instead of me :-)

p.s. I independently developed my own outline of tasks to do, which I'll share with the student that is accepted (which is likely you but theoretically could be another student if there is interest from others)

[asfimport]

Hatim Hakeel (migrated from JIRA)

Hi David,

Thanks for the compliments :). I'm eagerly looking forward to the announcement of accepted proposals with fingers crossed.

[asfimport]

John Berryman (migrated from JIRA)

In case you find it useful, this python function converts from an x,y coordinate (vals[0],vals[1]) on a grid from mins and maxs values to a point on the Hilbert curve in base4

def hilbert_point(vals,depth,mins=[],maxs=[]):
    """
    cell_names

    3 ---- 2
          |    
    0 ---- 1


    second generation

    33    30 -- 23 -- 22
    |     |           |
    32 -- 31    20 -- 21
                |
    01 -- 02    13 -- 12
    |     |           |
    00    03 -- 10 -- 11


    """
    if(mins):
        # then convert vals to some place in 1-by-1 grid at the origin
        vals[0] = (vals[0]-mins[0])/(maxs[0]-mins[0])
        vals[1] = (vals[1]-mins[1])/(maxs[1]-mins[1])
    
    if(depth == 0):
        return ''

    if( vals[0]<=0.5 and vals[1]<=0.5 ):
        cell_name = '0'
        x = vals[1]*2
        y = vals[0]*2
    elif( vals[0]>=0.5 and vals[1]<=0.5 ):
        cell_name = '1'
        x = (vals[0] - 0.5)*2
        y = vals[1]*2
    elif( vals[0]>=0.5 and vals[1]>=0.5 ):
        cell_name = '2'
        x = (vals[0] - 0.5)*2
        y = (vals[1] - 0.5)*2
    elif( vals[0]<=0.5 and vals[1]>=0.5 ):
        cell_name = '3'
        x = 1 - (vals[1] - 0.5)*2
        y = 1 - vals[0]*2

    return cell_name + hilbert_point([x,y],depth-1)

Try hilbert_point([0.125435,0.854320], 8, [0,0], [1,1]). It should be 33133010.

[asfimport]

David Smiley (@dsmiley) (migrated from JIRA)

Cool John. I suggest that any future code-samples be attached to the issue.

BTW another thing about the implementation of this grid I want to explore is the efficacy of working with integers and fixed precision (configured), versus floating point math.

[asfimport]

John Berryman (migrated from JIRA)

Hmmm... integer representation huh. Well here's a thought then:

As a first got at this idea, let's define something like a geohash where x are interleaved, but here's how we do it. At the top level, number squares from 0 to 3.

   0 --- 1
   |     |
   2 --- 3

At the next level, number things similarly,

    00 -- 01 -- 10 -- 11
    |     |     |     |
    02 -- 03 -- 12 -- 13
    |     |     |     |
    20 -- 21 -- 30 -- 31
    |     |     |     |
    22 -- 23 -- 32 -- 33

Even though this looks like the hilbert thing I did above, notice that this is actually the Z-ordering which is a little easier to compute.

In this case, the first two bits encodes which of the four big boxes the point is in, the next two bits encodes which of the four sub boxes the point is in, etc. So for example [0.375, 0.625] would be encoded to a depth of 2 by "03" which can be stored in half a byte.

Got it? So... now since we have the original point encoded in z-ordering, we can create a new hilbert_point algorithm that takes a byte array representing the z-ordering encoding of a point rather than a 2-vector of doubles. And the code looks much the same except that instead of the "val[0]*2" etc. we're actually just iterating through the byte array 2 bits at a time which is effectively the same as multiplying by 2.

This would make for some exquisitely indecipherable byte-munging code. But would it ultimately be more efficient? It largely depends upon how complex the Z-ordering encoding is. What do you think?

[asfimport]

John Berryman (migrated from JIRA)

This is a Java example of converting from a set of x,y values (within specified bounds) to a HilbertOrdered value.

The strategy is to

1. Convert the x,y value so doubles between 0 and 1.
2. Convert these doubles to large integers with max value 256^numBytes (user specifies numBytes). Note: there's probably a way to do these two last steps simultaneously and regain some precision. Note: numBytes must be <= 7 for now.
3. Interleave the bits so that you get a z-order value (note, I did this the "obvious way". In the code I've pointed to a website with a much more efficient method - so called morten numbers.
4. Convert the z-ordered value to a hilbert-ordered value.

How to do this last step really deserves a big whiteboard session - it's an inherently visual discussion. However, as a clue to what's happening in the code:

• There are only 4 shapes that compose the Hilbert curve. In the code I've called them "D,U,n, and C" because of the way they look. These are the states of a state machine.
• I convert from z to hilbert 2 bits at a time.
• On the first iteration I assume that I'm in the "D" state. In this simplistic case, I convert from 2 z-ordered bits to 2 hilbert-ordered bits based upon a lookup table that goes with the D state. I replace the z-ordered bits with the hilbert-ordered bits.
• I then check for which state I should go to next based upon a different lookup table that goes with the D state. It directs me to another state.
• I then get the next 2 bits from the byte array and repeat this method until I'm out of bits.

I've spot checked the input/output and it looks good (you'll see where I've done this in the code). No tests!

Also. This method could be slower than expected because I'm doing all operations on 2 bits at a time. As is, the method in the python code might even be faster because (correct me if I'm wrong) a double multiply can take place in one clock cycle.

That said, the methods here can be extended to operate at the processor word size.

What's more, I'm working with lookup tables. I suspect that you could use magic numbers and bitmasks and all that jazz and make something that took up very little space or time.

Also, I couldn't find them now, but there exist known efficient algorithms for doing this conversion. (I guess this weekend I just felt like making my own.:-/) I've run across algorithms that are even for higher dimension Hilbert Curves

[asfimport]

John Berryman (migrated from JIRA)

I just attached code called "HilbertConverter". This is a Java example of converting from a set of x,y values (within specified bounds) to a HilbertOrdered value.

The strategy is to

1. Convert the x,y value so doubles between 0 and 1.
2. Convert these doubles to large integers with max value 256^numBytes (user specifies numBytes). Note: there's probably a way to do these two last steps simultaneously and regain some precision. Note: numBytes must be <= 7 for now.
3. Interleave the bits so that you get a z-order value (note, I did this the "obvious way". In the code I've pointed to a website with a much more efficient method - so called morten numbers.
4. Convert the z-ordered value to a hilbert-ordered value.

How to do this last step really deserves a big whiteboard session - it's an inherently visual discussion. However, as a clue to what's happening in the code:

• There are only 4 shapes that compose the Hilbert curve. In the code I've called them "D,U,n, and C" because of the way they look. These are the states of a state machine.
• I convert from z to hilbert 2 bits at a time.
• On the first iteration I assume that I'm in the "D" state. In this simplistic case, I convert from 2 z-ordered bits to 2 hilbert-ordered bits based upon a lookup table that goes with the D state. I replace the z-ordered bits with the hilbert-ordered bits.
• I then check for which state I should go to next based upon a different lookup table that goes with the D state. It directs me to another state.
• I then get the next 2 bits from the byte array and repeat this method until I'm out of bits.

I've spot checked the input/output and it looks good (you'll see where I've done this in the code). No tests!

Also. This method could be slower than expected because I'm doing all operations on 2 bits at a time. As is, the method in the python code might even be faster because (correct me if I'm wrong) a double multiply can take place in one clock cycle.

That said, the methods here can be extended to operate at the processor word size.

What's more, I'm working with lookup tables. I suspect that you could use magic numbers and bitmasks and all that jazz and make something that took up very little space or time.

Also, I couldn't find them now, but there exist known efficient algorithms for doing this conversion. (I guess this weekend I just felt like making my own.:-/) I've run across algorithms that are even for higher dimension Hilbert Curves

[asfimport]

David Smiley (@dsmiley) (migrated from JIRA)

Thanks for your interest in this John. You've gotten the ball rolling. I confess I don't do any byte or bit manipulation work so it'll take some time to fully digest what's going on in this code which uses a lot of it. More comments would have helped. It might be interesting to look at how Lucene Trie floating point fields work. Thanks for the reference to efficient ways to interleave bits; I'll re-post here for convenient reference: http://www-graphics.stanford.edu/\~seander/bithacks.html#InterleaveBMN
(there are 2 approaches there). I'm not sure when I'll have significant time to contribute more directly other than add to the conversation.

[asfimport]

John Berryman (migrated from JIRA)

That's fine. Pull me into the conversations you have with Hatim and I can probably help things move along faster.

[asfimport]

Varun V Shenoy (migrated from JIRA)

I am Varun and I will be really interested to pursue this as my GSOC2014 project. Any documentation, references, papers etc on this topic will be extremely helpful.

[asfimport]

Hatim Hakeel (migrated from JIRA)

Hi,
I have provided a link to the proposal I prepared last year for the benefit of anyone looking for related background material like documentation, papers etc. Hope it would be helpful.

http://hatimhakeel.blogspot.com/2014/03/this-is-my-gsoc-2013-project-proposal.html

[asfimport]

Varun V Shenoy (migrated from JIRA)

Thanks Hatim, I have been going through your proposal, and it helped me a lot. I am getting really excited.