Did you use the rectangle-contains-point benchmark? I see that as a 'worst case' as by far most GEOS operations would be more heavy.

Yes, the contains operation is indeed one of the relatively fast operations.

I think that the two forloop approach should be done completely in C

OK, I tried something similar to your suggestion.

The current function (with expanded macro's) for contains (YY_b) looks like:

static void YY_b_func(char **args, npy_intp *dimensions,
                      npy_intp *steps, void *data)
{
    FuncGEOS_YY_b *func = (FuncGEOS_YY_b *)data;
    void *context_handle = geos_context[0];
    GEOSGeometry *in1, *in2;

    // BINARY_LOOP
    char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];
    npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];
    npy_intp n = dimensions[0];
    npy_intp i;
    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1) {
        //INPUT_YY;
        if (!get_geom(*(GeometryObject **)ip1, &in1)) { return; }
        if (!get_geom(*(GeometryObject **)ip2, &in2)) { return; }

        npy_bool ret = func(context_handle, in1, in2);

        //OUTPUT_b;
        if (! ((ret == 0) || (ret == 1))) {
            RAISE_ILLEGAL_GEOS;
            return;
        }
        *(npy_bool *)op1 = ret;
    }
}

If I change this to

static void YY_b_func(char **args, npy_intp *dimensions,
                      npy_intp *steps, void *data)
{
    FuncGEOS_YY_b *func = (FuncGEOS_YY_b *)data;
    void *context_handle = geos_context[0];
    GEOSGeometry *in1, *in2;

    // BINARY_LOOP
    char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];
    npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];
    npy_intp n = dimensions[0];
    npy_intp i;

    GEOSGeometry **geoms1 = malloc(sizeof(void *) * dimensions[0]);
    GEOSGeometry **geoms2 = malloc(sizeof(void *) * dimensions[0]);

    for(i = 0; i < n; i++, ip1 += is1, ip2 += is2) {
        if (!get_geom(*(GeometryObject **)ip1, &in1)) { return; }
        if (!get_geom(*(GeometryObject **)ip2, &in2)) { return; }
        geoms1[i] = in1;
        geoms2[i] = in2;
    }

    NPY_BEGIN_THREADS_DEF;
    NPY_BEGIN_THREADS;
    for(i = 0; i < n; i++, op1 += os1) {
        //INPUT_YY;
        in1 = geoms1[i];
        in2 = geoms2[i];

        npy_bool ret = func(context_handle, in1, in2);

        //OUTPUT_b;
        if (! ((ret == 0) || (ret == 1))) {
            RAISE_ILLEGAL_GEOS;
            return;
        }
        *(npy_bool *)op1 = ret;
    }
    NPY_END_THREADS;
    if (geoms1 != NULL) { free(geoms1); }
    if (geoms2 != NULL) { free(geoms2); }
}

where in a first loop over the inputs I create GEOSGeometry ** arrays, and then in the second loop over the output I do the actual GEOS operation in a part that releases the GIL.
This still seems to work correctly (I only suppose that the RAISE_ILLEGAL_GEOS is not safe to do inside the nogil part, that needs some further handling). And it gives (for the polygon-contains-points example) a 2x à 3x speed-up when done in parallel.

Does the above code looks more or less OK? (I would need to convert it again into macro's then)

Release the GIL where possible #7

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions