Release the GIL for geometry-creating operations

[caspervdw]

This addresses #155 , now only for the geometry-to-geometry operations (like envelope and convex_hull)

The same pattern could be repeated for all other geometry-creating operations.

[jorisvandenbossche]

Cool. I assume the python wrappers still need to be decorated with multithreading_enabled ?

Now, I am trying this out, and I don't see any speedup. I actually don't see usage of multiple threads on my laptop (I was testing with make_valid, both with dask as with stdlib threading). Which seems to point that there is still somewhere the GIL being reacquired in the loop, but don't directly see where that would be.
(but would need to do a more elaborate profile to be sure where it is spending time)

[jorisvandenbossche]

.. ahum, most stupid error ever ;) I forgot to restart my python kernel after building this branch

Now I see a good speedup. For make_valid (an expensive operation), I see around 3x speedup by running in parallel on my laptop

Also no noticeable slowdown in the single threaded setting compared to master.

[caspervdw]

.. ahum, most stupid error ever ;) I forgot to restart my python kernel after building this branch

Now I see a good speedup. For make_valid (an expensive operation), I see around 3x speedup by running in parallel on my laptop

Also no noticeable slowdown in the single threaded setting compared to master.

Good! Thanks for checking :) I'll proceed with the other geometry-returning functions. And indeed, I should add the decorators as well.

[caspervdw]

@jorisvandenbossche I am running into issues with the set operations. If you do this two-step procedure (create GEOSGeometry without GIL, create PyObject with GIL), the reduce functions do not work, as they expect the full PyObject to be created in the ufunc loop.

[brendan-ward]

@caspervdw what do you think about moving this up below line 345, with a return immediately after? Something like:

GEOSGeometry **geom_arr = malloc(sizeof(void *) * dimensions[0]);
if (geom_arr == NULL) {
    errstate = PGERR_NO_MALLOC;
    return
}

Then the "happy path" becomes simpler:

GEOS_INIT_THREADS;
BINARY_LOOP { ... }
GEOS_FINISH_THREADS;

if (errstate == PGERR_SUCCESS) {
    geom_arr_to_npy(geom_arr, args[2], steps[2], dimensions[0]);
}
free(geom_arr);

[caspervdw]

Thanks for the suggestion! The errstate variable doesn't exist outside of the GEOS_INIT and GEOS_FINISH, but I'll just set the python errorstate directly and return instead.

[jorisvandenbossche]

I am running into issues with the set operations. If you do this two-step procedure (create GEOSGeometry without GIL, create PyObject with GIL), the reduce functions do not work, as they expect the full PyObject to be created in the ufunc loop.

Hmm. So you mean that they expect that at the end of each loop iteration, there is a PyObject that can then be used as starting point for the next iteration?

Basically, a reduce operation looks like (so where the result gets reused as input)

r = identity
for i in range(len(A)):
  r = func(r, A[i])
return r

and currently, we have the OUTPUT_Y inside the loop (which converts the GEOSGeometry pointer to a PyObject and assigns this to op1).
And I assume that numpy then sets the input pointers and strides in such a way that op1 gets used as input again? (just guessing here, to try to understand what is exactly going on).

So then I was wondering: can't we solve it with adding yet another loop? First a loop unpacking the input objects into an array pointers. Then the loop with the actual operation only dealing with GEOSGeometry pointers, and then afterwards the loop wrapping in python objects.
But I suppose this actually doesn't help at all since when doing a reduce, you don't have two input arrays that can be unwrapped (do you know what ip1 and ip2 point to when it is called from reduce ?)

[jorisvandenbossche]

And I assume that numpy then sets the input pointers and strides in such a way that op1 gets used as input again? (just guessing here, to try to understand what is exactly going on).

That indeed seems to be the case:

https://github.com/numpy/numpy/blob/a5d021a1b6f439a19812926bc4d796ef5f346c44/numpy/core/src/umath/ufunc_object.c#L3500-L3508

            /* Turn the two items into three for the inner loop */
            dataptrs_copy[0] = dataptrs[0];
            dataptrs_copy[1] = dataptrs[1];
            dataptrs_copy[2] = dataptrs[0];
            strides_copy[0] = strides[0];
            strides_copy[1] = strides[1];
            strides_copy[2] = strides[0];
            innerloop(dataptrs_copy, &count,
                        strides_copy, innerloopdata);

[caspervdw]

And I assume that numpy then sets the input pointers and strides in such a way that op1 gets used as input again? (just guessing here, to try to understand what is exactly going on).

That indeed seems to be the case:

That is indeed the case. I added

printf("Args: %d, %d, %d Steps: %d, %d, %d\n", args[0], args[1], args[2], steps[0], steps[1], steps[2]);

into the ufunc and I got

>>> pygeos.lib.intersection.reduce(arr)
Args: 30428160, 26932968, 30428160 Steps: 0, 8, 0

I presume that the output array is filled with the first geometry in the array, else this would never work ;)

So then I was wondering: can't we solve it with adding yet another loop? First a loop unpacking the input objects into an array pointers. Then the loop with the actual operation only dealing with GEOSGeometry pointers, and then afterwards the loop wrapping in python objects.

I think this is our safest option. This might be tricky to get right though.

[caspervdw]

@jorisvandenbossche I tried 2 approaches, one with the temporary array, and one with some custom logic inside the inner loop. The second one turned out to be much cleaner than the first.

[brendan-ward]

This is shaping up nicely - thanks for your work on this @caspervdw

How does one test the use of multiple threads? I don't see that our test cases here have anything that covers that.

[caspervdw]

Thanks for your review @brendan-ward , I gave especially the reduce function some additional attention.

[jorisvandenbossche]

I was initially going to say: we can also simply "punt" on having a no-gil reduce version if it turns out to be too difficult / complex (it's not the most important thing to have).
But it seems you already found a relatively clean solution in the meantime! ;) Taking a look now

[jorisvandenbossche]

Super! It adds quite some complexity, but being able to release the GIL in those cases will also be really worth it IMO

Apart from the minor inline comments, I am wondering if some repeated snippets could be factored out in functions / templates, but that's certainly something that could also be done in a follow-up if we want this.

[jorisvandenbossche]

Should we call destroy_geom_arr when the operation was not succesful?
For example an error could have happened somewhere in the middle of the array, and then the first set of geometries still have to be destroyed before raising an error?

[caspervdw]

That is done inside the loop, whenever an error is set (lines 283 and 296)

[jorisvandenbossche]

Suggested change

	* second loop: out[0] = func(in1[0], in2[2]) [ = func(out[0], in2[1]) ]
	* second loop: out[0] = func(in1[0], in2[2]) [ = func(out[0], in2[2]) ]

[jorisvandenbossche]

I don't fully understand this comment. In what kind of cases is in1 owned by python? If it is the original input geometry?

[jorisvandenbossche]

What is the meaning of this dimension actually? Or more specifically, why is it larger than in case of reduce ?

[caspervdw]

Thanks @jorisvandenbossche, I tried to clarify the docstrings a bit.

Also, this PR is now tested for memory leaks using valgrind (see #159). It passes the tests.

[caspervdw]

@jorisvandenbossche and @brendan-ward , I think I addressed the issues. OK to merge?

[brendan-ward]

Thanks for your work on this @caspervdw - a couple minor comments.

[brendan-ward]

Should this be updated to use GEOS_INIT_THREADS / GEOS_FINISH_THREADS? If not, it might be good to add a comment indicating why it is not compatible with multiple threads.

[caspervdw]

Thanks for spotting that omission, I added it. Also, the creation functions like points need attention still, but I'd like to leave that for later.

Care to take a look at the new snap_func?

[brendan-ward]

This looks good - thanks @caspervdw !

I added a few follow on issues:

• TEST: Add multithreading tests for geometry-creating operations #164 for adding tests
• ENH: Consolidate repeated snippets used in ufuncs.c #165 for consolidating repeated snippets (I think they can be done after this PR)
• ENH: add documentation / examples for using multiple threads with geometry-creating operations #166 for documentation / example

[jorisvandenbossche]

Thanks a lot @caspervdw !

[jorisvandenbossche]

(BTW, I recently experimented with those multithreading capabilities again a bit in dask-geopandas -> https://github.com/jsignell/dask-geopandas/)

[caspervdw]

Thanks! @jorisvandenbossche maybe we should have a hangouts session sometime about this dask integration. I have this work project at github.com/nens/dask-geomodeling that also parallelizes geopandas.

[jorisvandenbossche]

@caspervdw for some reason I only notice your comment now, but: certainly interested in a hangout chat about dask integration!