(Numpy C API) Iterating over a single array: NpyIter vs for loop (with PyArray_DATA)

Question

I am writing some C extension code for a python module. The function I want to write is (in python)

output = 1./(1. + input)

where input is a numpy array of any shape.

Originally I was using NpyIter_MultiNew:

static PyObject *
helper_calc1(PyObject *self, PyObject *args){

    PyObject * input;
    PyObject * output = NULL;

    if (!PyArg_ParseTuple(args, "O", &input)){
        return NULL;
    }

    // -- input -----------------------------------------------
    PyArrayObject * in_arr;

    in_arr = (PyArrayObject *) PyArray_FROM_OTF(input, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (in_arr == NULL){
        goto fail;
    }

    // -- set up iterator -------------------------------------
    PyArrayObject * op[2];
    npy_uint32 op_flags[2];
    npy_uint32 flags;

    op[0] = in_arr;
    op_flags[0] = NPY_ITER_READONLY;

    op[1] = NULL;
    op_flags[1] = NPY_ITER_WRITEONLY | NPY_ITER_ALLOCATE;

    flags = NPY_ITER_EXTERNAL_LOOP | NPY_ITER_BUFFERED | NPY_ITER_GROWINNER;

    NpyIter * iter = NpyIter_MultiNew(2, op,
                            flags, NPY_KEEPORDER, NPY_NO_CASTING,
                            op_flags, NULL);

    if (iter == NULL){
        goto fail;
    };

    NpyIter_IterNextFunc * iternext = NpyIter_GetIterNext(iter, NULL);
    if (iternext == NULL){
        NpyIter_Deallocate(iter);
        goto fail;
    };

    // -- iterate ---------------------------------------------
    npy_intp count;
    char ** dataptr = NpyIter_GetDataPtrArray(iter);
    npy_intp * strideptr = NpyIter_GetInnerStrideArray(iter);
    npy_intp * innersizeptr = NpyIter_GetInnerLoopSizePtr(iter);

    do {
        count = *innersizeptr;

        while (count--){

            *(double *) dataptr[1] = 1. / (1. + *(double *)dataptr[0]);

            dataptr[0] += strideptr[0];
            dataptr[1] += strideptr[1];
        }

    } while (iternext(iter));

    output = NpyIter_GetOperandArray(iter)[1];

    if (NpyIter_Deallocate(iter) != NPY_SUCCEED){
        goto fail;
    }

    Py_DECREF(in_arr);

    return output;

    fail:
        Py_XDECREF(in_arr);
        Py_XDECREF(output);
        return NULL;
}

However, since this is just a single array (i.e. I don't need to be concerned about broadcasting multiple arrays), Is there any reason I can't iterate over the data myself using, PyArray_DATA, a for loop and the array size?

static PyObject *
helper_calc2(PyObject *self, PyObject *args){

    PyObject * input;
    PyObject * output = NULL;

    if (!PyArg_ParseTuple(args, "O", & in)){
        return NULL;
    }

    // -- input -----------------------------------------------
    PyArrayObject * in_arr;

    in_arr = (PyArrayObject *) PyArray_FROM_OTF(input, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
    if (in_arr == NULL){
        Py_XDECREF(in_arr);
        return NULL;
    }

    int ndim = PyArray_NDIM(in_arr);
    npy_intp * shape = PyArray_DIMS(in_arr);
    int size = (int) PyArray_SIZE(in_arr);

    double * in_data = (double *) PyArray_DATA(in_arr);

    output =  PyArray_SimpleNew(ndim, shape, NPY_DOUBLE);
    double * out_data = (double *) PyArray_DATA((PyArrayObject *) output);

    for (int i = 0; i < size; i++){
        out_data[i] = 1. / (1. + in_data[i]);
    }

    Py_DECREF(in_arr);
    return output;

fail:
    Py_XDECREF(in_arr);
    Py_XDECREF(output);
    return NULL;
}

This second version runs more quickly and the code is shorter.

Are there any dangers I need to look out for when using,PyArray_DATA with a for loop instead of NpyIter_MultiNew?

From the PyArray_DATA documentation:

If you have not guaranteed a contiguous and/or aligned array then be sure you understand how to access the data in the array to avoid memory and/or alignment problems.

But I believe that this is taken care of by PyArray_FROM_OTF with the NPY_ARRAY_IN_ARRAY flag.

score 0 · Answer 1 · answered Feb 12 '21 at 00:05

You should be ok in this case. Like you mentioned NPY_ARRAY_IN_ARRAY with PyArray_FROM_OTF takes care of that issue and since you are casting the data pointer to the right type, you should be fine. However, in general, as you seem to be aware of, you can't use this approach if you are directly accepting a NumPy array from Python code.

Happy C extension writing.

(Numpy C API) Iterating over a single array: NpyIter vs for loop (with PyArray_DATA)

1 Answers1