Use Flatbush RTree for Spatial Joins

[jagill]

There are memory problems with spatial joins, potentially due to the
large number of objects created by the JTS RTree. A Flatbush is an
RTree with a low memory footprint, minimal heap allocation, and fast
lookup. As the name would suggest, it is very flat -- the tree itself
is just a DoubleBigArray, with a couple other small fields. Lookup is
also pretty fast: the main difficulty is some index ninjitsu to handle
parent/child relationships.

This PR is branched off of #13235, which contains formatting and
other changes that are good independent of this one.

== RELEASE NOTES ==
General
--------
* Use Hilbert Packed RTree (Flatbush) instead of the JTS STR RTree for spatial joins.  The Flatbush is a flattened tree that uses vastly fewer heap allocations, and has faster build and query times.

[jagill]

cc @viczhang861 @mbasmanova

[jagill]

cc @jwass

[mbasmanova]

@jagill James, thanks for replacing JTS R-Tree with the one that doesn't create so many objects. Where can I read about the Flatbush you are using? Is https://github.com/mourner/flatbush/ the right place?

[mbasmanova]

Ignore order of rows in spatial join tests

Looks good, but I would split it into two commits: (1) reformat the code; (2) switch to assertOperatorEqualsIgnoreOrder

[mbasmanova]

Add HasExtent interface

Looks good.

[mbasmanova]

@OverRide

[mbasmanova]

nit: I prefer getEstimatedSize or getEstimatedSizeInBytes because often we can't know exactly how much memory an object is using and return an estimate or approximation.

[mbasmanova]

Flushing comments as I ran out of time.

Implement Flatbush Hilbert RTree

Commit message is indented (most likely unintentionally) - would you fix that?

[mbasmanova]

FLATBUSH_INSTANCE_SIZE -> INSTANCE_SIZE by convention

[mbasmanova]

• Flatbush.calculateLevelIndices -> calculateLevelIndices
• calculateLevelIndices -> calculateLevelOffsets to avoid using multiple terms for the same thing

[mbasmanova]

nit: this.buildTree() -> buildTree()

[mbasmanova]

• levelIndices -> offsets
• new ArrayList() -> new ArrayList<>()

        List<Integer> offsets = new ArrayList<>();

[mbasmanova]

Integer.valueOf(0) -> 0

offsets.add(0);

[mbasmanova]

this method is not used

[mbasmanova]

this.items -> items

[mbasmanova]

this variable is not used until next loop; move it there

[mbasmanova]

childLevel -> level

[jagill]

@mbasmanova For further reading on the flatbush, the mourner/flatbush repo is a good start; it also links to the wikipedia page which describes the general class/algorith: https://en.wikipedia.org/wiki/Hilbert_R-tree#Packed_Hilbert_R-trees

Note in this case we only sort the leaf notes, and rely on the hierarchy of the Hilbert curve to keep all parent nodes roughly in the right place.

[mbasmanova]

@jagill Finished reading Implement Flatbush Hilbert RTree commit. Overall looks good. Here are some comments.

[mbasmanova]

if (child % degree == 0) { might be clearer and make the comment unnecessary

[mbasmanova]

What are your thoughts re:

if (child % degree == 0) { might be clearer and make the comment unnecessary

?

[mbasmanova]

if (numChildren % degree > 0) { might be clearer

[mbasmanova]

What are your thoughts re:

if (readOffset % degree > 0) { might be clearer and make the comment unnecessary

?

[mbasmanova]

capacity -> numChildren to match buildTree?

[mbasmanova]

Let's not add this dependency. This interface is trivial. Let's just make our own.

    public interface ItemVisitor<T>
    {
        void visitItem(T item);
    }

[mbasmanova]

This method is used only in TestFlatbush. I'd move it there and use results::add:

    private static List<Rectangle> findIntersections(Flatbush<Rectangle> flatbush, Rectangle rectangle)
    {
        List<Rectangle> results = new ArrayList<>();
        flatbush.findIntersections(rectangle, results::add);
        return results;
    }

[mbasmanova]

• Don't abbreviate.
• I think a better structure would to be have helper method take two lists of rectangles and have the tests compose these lists as they see fit.

    @Test
    public void testRandomRectanglesSmall()
    {
        Random random = new Random(42);
        assertRectangles(makeRectangles(random, 1_000), makeRectangles(random, 100));
    }

    @Test
    public void testRandomRectanglesMedium()
    {
        Random random = new Random(123);
        assertRectangles(makeRectangles(random, 10_000), makeRectangles(random, 1_000));
    }

    @Test
    public void testRandomRectanglesLarge()
    {
        Random random = new Random(321);
        assertRectangles(makeRectangles(random, 50_000), makeRectangles(random, 5_000));
    }

    private static void assertRectangles(List<Rectangle> buildRectanges, List<Rectangle> probeRectanges)
    {
        Flatbush<Rectangle> rtree = new Flatbush<>(buildRectanges.toArray(new Rectangle[] {}));
        for (Rectangle probe : probeRectanges) {
            List<Rectangle> expected = buildRectanges.stream()
                    .filter(rectangle -> rectangle.intersects(probe))
                    .collect(toImmutableList());
            assertEqualsNoOrder(findIntersections(rtree, probe), expected);
        }
    }

    private static void assertEqualsNoOrder(List<Rectangle> actual, List<Rectangle> expected)
    {
        List<Rectangle> actualSorted = actual.stream().sorted(RECTANGLE_COMPARATOR).collect(toImmutableList());
        List<Rectangle> expectedSorted = expected.stream().sorted(RECTANGLE_COMPARATOR).collect(toImmutableList());

        assertEquals(actualSorted, expectedSorted);
    }

[mbasmanova]

move this to the previous line:

Flatbush<OGCGeometryWrapper> rtree = new Flatbush<>(new OGCGeometryWrapper[] {pointX, pointY, pointZ, pointW});

[mbasmanova]

use assertEqualsNoOrder helper function

        assertEqualsNoOrder(findIntersections(rtree, octagonA.getExtent()), ImmutableList.of(pointX, pointY), Comparator.naturalOrder());
        assertEqualsNoOrder(findIntersections(rtree, octagonB.getExtent()), ImmutableList.of(pointY, pointZ), Comparator.naturalOrder());

    private static <T> void assertEqualsNoOrder(List<T> actual, List<T> expected, Comparator<T> comparator)
    {
        List<T> actualSorted = actual.stream().sorted(comparator).collect(toImmutableList());
        List<T> expectedSorted = expected.stream().sorted(comparator).collect(toImmutableList());

        assertEquals(actualSorted, expectedSorted);
    }

[mbasmanova]

• Don't abbreviate
• rng -> random, minVal -> min, maxVal -> max
• makeRandomRectangle -> makeRectangle; make it private
• makeRandomRectangles -> makeRectanges

public static List<Rectangle> makeRectangles(Random random, int count)
    {
        return IntStream.range(0, count)
                .mapToObj(i -> makeRectangle(random))
                .collect(toImmutableList());
    }

    private static double randomDouble(Random random, double min, double max)
    {
        return min + random.nextDouble() * (max - min);
    }

    /*
     * Make a random rectangle "near" origin of size "about" size.
     */
    private static Rectangle makeRectangle(Random random)
    {
        double minX = randomDouble(random, -100, 100);
        double minY = randomDouble(random, -100, 100);
        double sizeX = randomDouble(random, 0.0, 10);
        double sizeY = randomDouble(random, 0.0, 10);
        return new Rectangle(minX, minY, minX + sizeX, minY + sizeY);
    }

[mbasmanova]

This is a very generic name. Perhaps, RectangleTestUtils?

[mbasmanova]

Use Flatbush instead of JTS STR RTree for spatial joins

This is awesome!!!

I'd change commit title to Use Flatbush R-Tree in spatial join

[mbasmanova]

static import + inline: rtree.findIntersections(getExtent(probeGeometry), item -> {

[jagill]

We can't inline, because it's used in the closure (for testReferencePoint).

[mbasmanova]

could you move this class in a separate commit?

[mbasmanova]

Use Flatbush instead of JTS STR RTree for spatial joins

Very nice. Thank you!

[mbasmanova]

do not abbreviate: queryRect -> query or rectangle

[mbasmanova]

Thanks for adding benchmarks. Some of them show very high error bounds. Would you adjust the settings on the benchmarks to bring errors down?

293.261 ± 41.055
444.502 ± 86.900

[mbasmanova]

move main method to the end of the class; ditto other benchmarks

[mbasmanova]

don't abbreviate: numBuildRects -> numRectangles and buildRects -> rectangles

[mbasmanova]

don't abbreviate

[mbasmanova]

don't abbreviate throughout this commit

[mbasmanova]

CC: @pettyjamesm

[jagill]

cc @mourner: This PR ports your Flatbush to Java for use with Presto. Just thought you should know, and happy to hear any feedback you have!

[jagill]

cc @rawrunprotected We used your log(n) 2D hilbert curve algorithm here. Just wanted to let you know, and to see if you had any feedback! We had to port to Java signed ints (wrapped in a long), because Java doesn't have unsigned ints :(.

[rawrunprotected]

cc @rawrunprotected We used your log(n) 2D hilbert curve algorithm here. Just wanted to let you know, and to see if you had any feedback! We had to port to Java signed ints (wrapped in a long), because Java doesn't have unsigned ints :(.

Thanks for letting me know, awesome to see my stuff being used :)

I haven't touched Java in a decade, but I think the code should work fine with just plain ints and without explicit masking to 16bit intermediates. Bitwise operations don't care about signedness, and the code already uses logical right shifts (>>>), so unsigned vs signed doesn't really matter at all. The only caveat is that you'll want to zero-extend the result at the very end to long in order to allow using it as a non-negative offset. That can be done with a quick pair of shifts. See attached, but I only did some very superficial testing ;)

hilbert.txt

[mourner]

@jagill this is awesome, thanks a lot for letting me know! Can you share the benchmark results here to get an idea for the improvement?

@rawrunprotected thanks so much for chiming in here and for the amazing original implementation! In the JS version (which this PR ports), I use it like this, although it could probably be simplified further because I limit the input to 16-bit uints (plenty for spatial indexing needs) so that the resulting hilbert values can be stored in a uint32 array.

[jagill]

@rawrunprotected Thanks! That little change sped up build times about 10%. Prior to using (the first version of) your code, the building bottleneck was the Hilbert index calculation, using the bit twiddling massively sped up the sorting.

@mourner I've included some benchmarks in the description of the commit that adds them: c1d6b4d . In general, both querying and building are about 4x faster when we get to 1M geometries. Critically, it has a better scaling factor than the JTS STR RTree.

[mbasmanova]

Add HasExtent interface

Looks good.

[mbasmanova]

Implement Flatbush Hilbert RTree

Looks good.

[mbasmanova]

Use Flatbush instead of JTS STR RTree for spatial joins

Looks good.

[mbasmanova]

Ignore order of rows in spatial join tests

Thanks.

[mbasmanova]

Use Flatbush RTree for Spatial Joins

Looks good.

[mbasmanova]

@jagill Looks good to me. This is a very nice and desirable improvement for the spatial joins. Thanks for researching and implementing such an elegant fix.

[dr-jts]

@jagill is it inefficient to compute the hilbert index dynamically like this? Won't there be many computations of the same index?

[jagill]

The index is calculated on average O(ln N) times per entry, as you said. I originally used davidmoten/hilbert-curve, and (via profiling) index calculation did indeed dominate construction CPU. However, this new calculation is very efficient (all bit twiddling), and so now it's only a few % of construction CPU.

[dr-jts]

Ok, makes sense, and certainly simpler than maintaining a pre-computed value.

[dr-jts]

Does this logic compute "query contains node", or "query intersects node"?