Segfault when registering gufunc with signature "()->(4)" in C (with numpy 1.15)

[jorisvandenbossche]

In PyGEOS, we are writing (g)ufuncs in C, and are running into a problem with a generalized ufunc on older numpy versions: pygeos/pygeos#69

We are defining a generalized ufunc with a ()->(4) signature. With numpy 1.17 this works fine, but with older numpy (tested 1.15 locally, and on travis also 1.13) it is failing.
It doesn't fail when running, but already segfaults when registering the ufunc with that signature (to be more precise: it segfaults on importing the python package that on import loads the extension module that registers the ufuncs).

We are doing something like (the full diff can be seen in the issue linked above):

ufunc = PyUFunc_FromFuncAndDataAndSignature(
    bounds_funcs, null_data, bounds_dtypes, 1, 1, 1, PyUFunc_None,
    "bounds", "", 0, "()->(4)");
PyDict_SetItemString(d, "bounds", ufunc)

and there doesn't seem to be anything wrong with the inner functions and dtypes (bounds_funcs and bounds_dtypes), because if I change the signature from "()->(4)" to "()->()", it actually "works" (it no longer segfaults on registration, and it runs, only giving a wrong (partly) result).

I know this is a rather vague issue description, and in addition only happens for older numpy (but still, as a dependent library, we would like to support a range of numpy versions). But:

• Was there a known problems with signatures like "()->(4)"? Or was this not a supported signature?
• Would there be a workaround to still do something like this with older numpy versions like 1.15?

[WarrenWeckesser]

The ability to specify gufunc signatures with fixed-size dimensions was added in NumPy 1.16: https://numpy.org/devdocs/release/1.16.0-notes.html#generalized-ufunc-signatures-now-allow-fixed-size-dimensions

I don't think there is a simple way to work around this for older versions of NumPy, but maybe @mattip, @mvhk, or another developer involved with the implementation of NEP 20 will have a recommendation.

[jorisvandenbossche]

Thanks for the answer! That at least explains our segfault ;)

And from reading the description of the NEP, I might have spotted a workaround (emphasize mine):

Fixed-size dimensions. ... E.g., the signature of a function that converts a polar angle to a two-dimensional cartesian unit vector would currently have to be ()->(n), with there being no way to indicate that n has to equal 2. Indeed, this signature is particularly annoying since without putting in an output argument, the current gufunc wrapper code fails because it cannot determine n.

So we can actually use ()->(n) and then provide an empty array as output argument with the proper shape (original shape + (4,)).
Since we still wrap the "raw" ufuncs in python functions, this should be a sufficient workaround for our use case (so for end users we can hide the need to provide this out inside that function).

From a quick test this seems to work (now compiled with ()->(n) and using numpy 1.15):

In [1]: import pygeos 

# create array with 2 points
In [6]: a = pygeos.points([1, 2], [3, 4])   

In [7]: a      
Out[7]: 
array([<pygeos.Geometry POINT (1 3)>, <pygeos.Geometry POINT (2 4)>],
      dtype=object)

In [8]: pygeos.bounds(a) 
...
ValueError: bounds: Output operand 0 has core dimension 0 unspecified, with gufunc signature ()->(n)

In [10]: out = np.empty(a.shape + (4,), dtype='float64')  

In [11]: pygeos.bounds(a, out=out)  
Out[11]: 
array([[1., 3., 1., 3.],
       [2., 4., 2., 4.]])

[jorisvandenbossche]

OK, closing this, as it seems solved for us. Thanks again for the pointer!

[WarrenWeckesser]

@jorisvandenbossche, thanks for the update. Glad to hear you worked out a solution.