> On Nov 28, 2020, at 1:13 PM, Paul Ramsey <pram...@cleverelephant.ca> wrote:
> 
> 
> 
>> On Nov 28, 2020, at 1:11 PM, Paul Ramsey <pram...@cleverelephant.ca> wrote:
>> 
>> 
>> 
>>> On Nov 28, 2020, at 1:05 PM, Joris Van den Bossche 
>>> <jorisvandenboss...@gmail.com> wrote:
>>> 
>>> 
>>> On Sat, 28 Nov 2020 at 21:53, Paul Ramsey <pram...@cleverelephant.ca> wrote:
>>> 
>>> 
>>>> On Nov 28, 2020, at 12:44 PM, Joris Van den Bossche 
>>>> <jorisvandenboss...@gmail.com> wrote:
>>>> 
>>>> Thanks for trying to reproduce it in C/C++. One obvious difference that I 
>>>> can spot is that we use an integer for the "item" that gets inserted, and 
>>>> not the geometry itself, but I wouldn't expect that to influence the 
>>>> result. 
>>>> Although, trying to update your test case to do that, the test fails. But 
>>>> that might also be an issue on my side due to my limited C++ experience 
>>>> (it already fails on the "geoms.size()" check):
>>> 
>>> Nope still not seeing it... one small mistake in your query
>>> 
>>> Also not if you try my original code *with* the mistake? 
>> 
>> Well, with the mistake the code is saying "here's a null pointer" (since 
>> that's what casting 0 to a void* will get you) index this. I'll see! It 
>> shouldn't really matter from an index point-of-view, it should still index 
>> it and return a null pointer out the back end.
> 
> And indeed it does not work. So, something to examine... why can we not index 
> with a null item?

Because I was using a non-null item pointer as a cue that the node was a leaf 
node, so it wasn't picking up null items as leaves. 
So, I've pushed that fix, does that change your test suite all that much?

P


> 
> P
> 
>> 
>>> Because your version indeed passes for me as well, but I *think* the 
>>> version I wrote resembles more closely the PyGEOS code (so it might be an 
>>> issue in our C code on how we use the tree).
>>> 
>>> , trying to cast the int to a void, instead of passing in a the pointer to 
>>> the address, here's one that works. Unfortunately that leaves us no closer 
>>> to knowing why the SimpleSTRtree is unhappy in the python context. I fear I 
>>> may just have to revert the CAPI to the old tree.
>>> 
>>> 
>>> // querying tree with box
>>> template<>
>>> template<>
>>> void object::test<9>
>>> ()
>>> {
>>>   GEOSSTRtree* tree = GEOSSTRtree_create(10);
>>> 
>>>   GEOSGeometry* g = GEOSGeomFromWKT("POINT (2 3)");
>>>   int payload = 876;
>>>   GEOSSTRtree_insert(tree, g, &payload);
>>> 
>>>   GEOSGeometry* q = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 
>>> 0))");
>>> 
>>>   typedef std::vector<int*> IList;
>>>   IList items;
>>>   ensure_equals(items.size(), 0);
>>>   GEOSSTRtree_query(
>>>       tree,
>>>       q,
>>>       [](void* item, void* userdata) {
>>>           IList* items = (IList*)userdata;
>>>           items->push_back((int*)item);
>>>       },
>>>       &items);
>>> 
>>>   ensure_equals(items.size(), 1);
>>> 
>>>   ensure_equals(*(items[0]), payload);
>>> 
>>>   GEOSGeom_destroy(q);
>>>   GEOSGeom_destroy(g);
>>>   GEOSSTRtree_destroy(tree);
>>> }
>>> 
>>> 
>>> 
>>> 
>>>> 
>>>> --- a/tests/unit/capi/GEOSSTRtreeTest.cpp
>>>> +++ b/tests/unit/capi/GEOSSTRtreeTest.cpp
>>>> @@ -268,10 +268,11 @@ void object::test<8>
>>>> {
>>>>    GEOSSTRtree* tree = GEOSSTRtree_create(10);
>>>>    GEOSGeometry* g = GEOSGeomFromWKT("POINT (2 3)");
>>>> -    GEOSSTRtree_insert(tree, g, g);
>>>> +    int idx = 0;
>>>> +    GEOSSTRtree_insert(tree, g, (void*)idx);
>>>>    GEOSGeometry* q = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 
>>>> 0))");
>>>> 
>>>> -    typedef std::vector<GEOSGeometry*> GList;
>>>> +    typedef std::vector<int> GList;
>>>>    GList geoms;
>>>>    ensure_equals(geoms.size(), 0);
>>>>    GEOSSTRtree_query(
>>>> @@ -279,23 +280,16 @@ void object::test<8>
>>>>        q,
>>>>        [](void* item, void* userdata) {
>>>>            GList* geoms = (GList*)userdata;
>>>> -            geoms->push_back((GEOSGeometry*)item);
>>>> +            geoms->push_back(*((int *)item));
>>>>        },
>>>>        &geoms);
>>>> 
>>>>    ensure_equals(geoms.size(), 1);
>>>> -    const GEOSCoordSequence* seq = GEOSGeom_getCoordSeq(geoms[0]);
>>>> -
>>>> -    double x = -1;
>>>> -    double y = -1;
>>>> -    GEOSCoordSeq_getXY(seq,  0, &x, &y);
>>>> -    ensure_equals(x, 2.0);
>>>> -    ensure_equals(y, 3.0);
>>>> +    ensure_equals(geoms.at(0), 0);
>>>> 
>>>> 
>>>> On Sat, 28 Nov 2020 at 20:55, Paul Ramsey <pram...@cleverelephant.ca> 
>>>> wrote:
>>>> 
>>>> 
>>>>> On Nov 28, 2020, at 8:11 AM, Joris Van den Bossche 
>>>>> <jorisvandenboss...@gmail.com> wrote:
>>>>> 
>>>>> Hi all,
>>>>> 
>>>>> On the CI of PyGEOS we have a build testing against GEOS master, and 
>>>>> somewhere in the last 4 days, a lot of the STRtree tests started failing 
>>>>> (see eg https://github.com/pygeos/pygeos/runs/1465460418#step:9:86). 
>>>>> Looking at the commits of the last days, this might be related to the 
>>>>> SimpleSTRtree work?
>>>>> 
>>>>> A small (python) example of a tree consisting of a single point, which no 
>>>>> longer is returned when querying the tree with a big polygon that 
>>>>> certainly contains the point:
>>>>> 
>>>>> Using released version of GEOS:
>>>>> 
>>>>>>>> import pygeos
>>>>>>>> pygeos.geos_version
>>>>> (3, 8, 1)
>>>>>>>> point = pygeos.Geometry("POINT (2 3)") 
>>>>>>>> tree = pygeos.STRtree([point]) 
>>>>>>>> tree.query(pygeos.box(0, 0, 10, 10)) 
>>>>> array([0])
>>>>> 
>>>>> This is correctly returning the index of the single point. But when 
>>>>> running with the latest GEOS master, the query doesn't find any point of 
>>>>> the tree:
>>>>> 
>>>>>>>> import pygeos
>>>>>>>> pygeos.geos_version
>>>>> (3, 9, 0)
>>>>>>>> point = pygeos.Geometry("POINT (2 3)")
>>>>>>>> tree = pygeos.STRtree([point])
>>>>>>>> tree.query(pygeos.box(0, 0, 10, 10))
>>>>> array([], dtype=int64)
>>>>> 
>>>>> Are there changes expected in how the STRtree C API functions or required 
>>>>> changes in user code? Or maybe we are using it in some 
>>>>> incorrect/unexpected way? (code is at 
>>>>> https://github.com/pygeos/pygeos/blob/master/src/strtree.c)
>>>> 
>>>> There are changes, I don't think you're mis-using anything. I swapped the 
>>>> CAPI to use the SimpleSTRtree, figuring it would be good to share the 
>>>> performance win with downstream. However, I can swap it back to the 
>>>> original STRtree if this remains a problem.
>>>> 
>>>> One thing I noticed when trying to construct GEOS envelopes directly was 
>>>> that annoyingly they were xmin xmax, ymin ymax, but I doubt that would be 
>>>> a problem in your pre-existing working test. 
>>>> 
>>>> I just reconstructed your test in the GEOS CAPI suite, and it works as one 
>>>> would expect. (Namely, it finds the one point.) So I'm not sure why your 
>>>> test is getting different results.
>>>> 
>>>> 
>>>> // querying tree with box
>>>> template<>
>>>> template<>
>>>> void object::test<8>
>>>> ()
>>>> {
>>>>   GEOSSTRtree* tree = GEOSSTRtree_create(10);
>>>>   GEOSGeometry* g = GEOSGeomFromWKT("POINT (2 3)");
>>>>   GEOSSTRtree_insert(tree, g, g);
>>>>   GEOSGeometry* q = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 
>>>> 0))");
>>>> 
>>>>   typedef std::vector<GEOSGeometry*> GList;
>>>>   GList geoms;
>>>>   ensure_equals(geoms.size(), 0);
>>>>   GEOSSTRtree_query(
>>>>       tree,
>>>>       q,
>>>>       [](void* item, void* userdata) {
>>>>           GList* geoms = (GList*)userdata;
>>>>           geoms->push_back((GEOSGeometry*)item);
>>>>       },
>>>>       &geoms);
>>>> 
>>>>   ensure_equals(geoms.size(), 1);
>>>>   const GEOSCoordSequence* seq = GEOSGeom_getCoordSeq(geoms[0]);
>>>> 
>>>>   double x = -1;
>>>>   double y = -1;
>>>>   GEOSCoordSeq_getXY(seq,  0, &x, &y);
>>>>   ensure_equals(x, 2.0);
>>>>   ensure_equals(y, 3.0);
>>>> 
>>>>   GEOSGeom_destroy(q);
>>>>   GEOSGeom_destroy(g);
>>>>   GEOSSTRtree_destroy(tree);
>>>> }
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>>>>> 
>>>>> Best,
>>>>> Joris
>>>>> 
>>>>> 
>>>>> On Wed, 25 Nov 2020 at 00:44, Paul Ramsey <pram...@cleverelephant.ca> 
>>>>> wrote:
>>>>> Hey all, just truing up what's underway and nearly there...
>>>>> 
>>>>> - Am I right that Z coordinates are nearly done? What's the status there?
>>>>> 
>>>>> I've been trying to address some performance issues, with some success 
>>>>> and some ... other things. 
>>>>> 
>>>>> The success is the SimpleSTRtree, which is just the current STRtree but 
>>>>> without the inheritance structure and with the nodes stored all next to 
>>>>> each other in contiguous memory for more locality. For at least one use 
>>>>> case I've seen 20% speed-ups on overlays, using the SimpleSTRtree in 
>>>>> place of the STRtree inside the MCIndexNoder. I have not seen any 
>>>>> slow-downs. I have pushed the SimpleSTRtree into master.
>>>>> 
>>>>> While I have implemented the nearestNeighbor() functionality on the 
>>>>> SimpleSTRtree, I haven't hooked it up to anything yet. It could go into 
>>>>> the IndexedFacetDistance, if anyone is super enthusiastic about it. From 
>>>>> there it would affect searching in PreparedGeometry of various sorts.
>>>>> 
>>>>> I also tried using a similar trick with the MonotoneChainBuilder that 
>>>>> sits inside the MCIndexNoder, replacing individual heap allocations with 
>>>>> slabs by putting objects onto a std::deque, and incidentally stripping 
>>>>> out some book-keeping. While that seems to pick up about 3-5% speedwise, 
>>>>> unfortunately something about my implementation is incorrect (and in a 
>>>>> wonderfully subtle way) as it fails testing on some platforms (not mine). 
>>>>> https://github.com/pramsey/geos/tree/monotone-chain-builder
>>>>> 
>>>>> I've put that work to the side for now.
>>>>> 
>>>>> All the performance talk is mostly because JTS still runs a lot faster 
>>>>> than GEOS for some bulk processing. My current test is a big union of 
>>>>> watershed boundaries, about 6MB of data, which takes about 20s under GEOS 
>>>>> and about 25% of that under JTS. It's a big gap, and in theory the two 
>>>>> code bases are pretty aligned right now. Same overlayNG engine, etc. So I 
>>>>> figure there has to be a big implementation ball of performance hiding 
>>>>> under the covers somewhere. No luck thus far.
>>>>> 
>>>>> I think we're close, looking forward to release :)
>>>>> 
>>>>> P
>>>>> _______________________________________________
>>>>> geos-devel mailing list
>>>>> geos-devel@lists.osgeo.org
>>>>> https://lists.osgeo.org/mailman/listinfo/geos-devel
>>>>> _______________________________________________
>>>>> geos-devel mailing list
>>>>> geos-devel@lists.osgeo.org
>>>>> https://lists.osgeo.org/mailman/listinfo/geos-devel
>>>> 
>>>> _______________________________________________
>>>> geos-devel mailing list
>>>> geos-devel@lists.osgeo.org
>>>> https://lists.osgeo.org/mailman/listinfo/geos-devel
>>>> _______________________________________________
>>>> geos-devel mailing list
>>>> geos-devel@lists.osgeo.org
>>>> https://lists.osgeo.org/mailman/listinfo/geos-devel
>>> 
>>> _______________________________________________
>>> geos-devel mailing list
>>> geos-devel@lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/geos-devel
>>> _______________________________________________
>>> geos-devel mailing list
>>> geos-devel@lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/geos-devel

_______________________________________________
geos-devel mailing list
geos-devel@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/geos-devel

Reply via email to