Passing 3-dimensional numpy array to C

Question

I'm writing a C extension to my Python program for speed purposes, and running into some very strange behaviour trying to pass in a 3-dimensional numpy array. It works with a 2-dimensional array, but I'm sure I'm screwing something up with the pointers trying to get it to work with the 3rd dimension. But here's the weird part. If I just pass in a 3-D array, it crashes with a Bus Error. If (in Python) I create my variable as a 2D array first, and then overwrite it with a 3D array, it works perfectly. If the variable is an empty array first and then a 3D array, it crashes with a Seg Fault. How can that possibly happen?

Also, can anyone help me get a 3D array working? Or should I just give up and pass in a 2D array and reshape it myself?

Here's my C code:

static PyObject* func(PyObject* self, PyObject* args) {
  PyObject *list2_obj;
  PyObject *list3_obj;
  if (!PyArg_ParseTuple(args, "OO", &list2_obj, &list3_obj))
    return NULL;

  double **list2;
  double ***list3;

  //Create C arrays from numpy objects:
  int typenum = NPY_DOUBLE;
  PyArray_Descr *descr;
  descr = PyArray_DescrFromType(typenum);
  npy_intp dims[3];
  if (PyArray_AsCArray(&list2_obj, (void **)&list2, dims, 2, descr) < 0 || PyArray_AsCArray(&list3_obj, (void ***)&list3, dims, 3, descr) < 0) {
    PyErr_SetString(PyExc_TypeError, "error converting to c array");
    return NULL;
  }
  printf("2D: %f, 3D: %f.\n", list2[3][1], list3[1][0][2]);
}

And here's my Python code that calls the above function:

import cmod, numpy
l2 = numpy.array([[1.0,2.0,3.0], [4.0,5.0,6.0], [7.0,8.0,9.0], [3.0, 5.0, 0.0]])

l3 = numpy.array([[2,7, 1], [6, 3, 9], [1, 10, 13], [4, 2, 6]])  # Line A
l3 = numpy.array([])                                             # Line B

l3 = numpy.array([[[2,7, 1, 11], [6, 3, 9, 12]],
                 [[1, 10, 13, 15], [4, 2, 6, 2]]])

cmod.func(l2, l3)

So, if I comment out both Line A and B, it crashes with a Bus error. If Line A is there, but Line B is commented out, it runs correctly with no errors. If Line B is there but Line A is commented out, it prints the correct numbers but then Seg faults. Finally if both lines are present it also prints the correct numbers and then Seg faults. What in the hell is going on here?

EDIT: Ok. Wow. So I was using int in Python but calling them double in C. And that was working fine with 1D and 2D arrays. But not 3D. So I changed the Python definition of l3 to be floats, and now it all works fantastically (Thank you very much Bi Rico).

But now, more strange behaviour with Lines A & B! Now if both Lines are commented out, the program works. If Line B is present but A is commented out, it works, and ditto if both are uncommented. But if Line A is present and B is commented out, I get that fantastic Bus error again. I'd really like to avoid these in the future, so does anyone have any clue why the declaration of a Python variable can have this kind of impact?

EDIT 2: Well, as insane as these errors are, they're all due to the 3-dimensional numpy array I pass in. If I only pass in 1- or 2-D arrays, it behaves as expected, and manipulation of the other Python variables does nothing. This leads me to believe that the problem lies somewhere in Python's reference counting. In the C-code the reference count is decreased more than it should for the 3-D arrays, and when that function returns Python tries to clean up objects, and attempts to delete a NULL pointer. This is just my guess, and I've tried to Py_INCREF(); everything I could think of to no avail. I guess I'll just be using a 2D array and reshaping it in C.

Circe

Are you sure that (void **) is correct, shouldn't you just pass in a (void*)? – Gabriel

My C sucks but... Isn't your expression in the if short-circuited if the first call to PyArray_AsCArray suceeds? It may very well be that the second call, i.e. the one for list3, is never made. – Lyndy

@Gabriel I'm not positive that (void **) is correct, but (void*) causes a Bus error. @Lyndy No, that function returns negative values only if it fails, most likely if the malloc it calls fails. – Circe

@Gabriel Ok... now that I took Bi Rico's advice and tried python floats, both a single and double (or triple) star seem to work fine. Any idea which is best/right? – Circe

Getting the following error: "error: ‘NPY_DOUBLE’ undeclared (first use in this function)". – Penicillium

Solution to the error above: Include "numpy/arrayobject.h". Ref: mail.python.org/pipermail/scipy-user/2011-October/030801.html – Penicillium

Answer 1

I already mentioned this in a comment, but I hope flushing it out a little helps make it more clear.

When you're working with numpy arrays in C it's good to be explicit about the typing of your arrays. Specifically it looks like you're declaring your pointers as double ***list3, but they way you're creating l3 in your python code you'll get an array with dtype npy_intp (I think). You can fix this by explicitly using the dtype when creating your arrays.

import cmod, numpy
l2 = numpy.array([[1.0,2.0,3.0],
                  [4.0,5.0,6.0],
                  [7.0,8.0,9.0],
                  [3.0, 5.0, 0.0]], dtype="double")

l3 = numpy.array([[[2,7, 1, 11], [6, 3, 9, 12]],
                  [[1, 10, 13, 15], [4, 2, 6, 2]]], dtype="double")

cmod.func(l2, l3)

Another note, because of the way python works it's nearly impossible for "line A" and "line B" to have any effect on the C code what so ever. I know that this seems to conflict with your empirical experience, but I'm pretty sure on this point.

I'm a little less sure about this, but based on my experience with C, bus-errors and segfaults are not deterministic. They depend on memory allocation, alignment, and addresses. In some situation code seems to run fine 10 times, and fails on the 11th run even though nothing has changed.

Have you considered using cython? I know it's not an option for everyone, but if it is an option you could get nearly C level speedups using typed memoryviews.

Trawick

The next time I have a need to write a C extension I'm fairly sure that I'll spend the time to learn cython. And yes, everything I know about Python and C says that there should be no way that "Line A and B" could possibly impact the C program, as each time L2 is declared it gets a new memory address. But they absolutely are for me, and that's a major reason I started this question. I could paste the whole files if someone else wants to try on their system, as I'd love to get to the bottom of this. – Circe

Answer 2

Rather than converting to a c-style array, I usually access numpy array elements directly using PyArray_GETPTR (see https://numpy.org/doc/stable/reference/c-api/array.html#data-access).

For instance, to access an element of a 3-dimensional numpy array of type double use double elem=*((double *)PyArray_GETPTR3(list3_obj,i,j,k)).

For your application, you could detect the correct number of dimensions for each array using PyArray_NDIM, then access elements using the appropriate version of PyArray_GETPTR.

Crummy

I wanted to convert to a regular C array because I assumed that it would be faster. I also assumed it would be simpler, but that was clearly wrong... – Circe

Any idea if this is slower or faster? – Marisamariscal

This is a little slower but not much. The GETPTR is a simple macro, you can see how it's defined at numpy's github under numpy/core/numpy/include/numpy/ndarrayobject.h. Say we have PyArrayObject *myArray and nb = number of bytes per element. This returns a pointer to address myArray->data + (i * myArray->strides[0] / (bytes per element)) + (j * myArray->strides[1] / (bytes per element)) + (k * myArray->strides[2] / (bytes per element)). There is a small performance penalty to that as opposed to incrementing a pointer that is fairly negligible but could add up if done inside an inner loop. – Hereunto

Answer 3

According to http://docs.scipy.org/doc/numpy/reference/c-api.array.html?highlight=pyarray_ascarray#PyArray_AsCArray:

Note The simulation of a C-style array is not complete for 2-d and 3-d arrays. For example, the simulated arrays of pointers cannot be passed to subroutines expecting specific, statically-defined 2-d and 3-d arrays. To pass to functions requiring those kind of inputs, you must statically define the required array and copy data.

I think that this means that PyArray_AsCArray returns a block of memory with the data in it in C order. However, to access that data, more information is needed (see http://www.phy225.dept.shef.ac.uk/mediawiki/index.php/Arrays,_dynamic_array_allocation). This can either be achieved by knowing the dimensions ahead of time, declaring an array, and then copying the data in in the right order. However, I suspect that more general case is more useful: you don't know the dimensions until they are returned. I think that the following code will create the necessary C pointer framework to allow the data to be addressed.

static PyObject* func(PyObject* self, PyObject* args) {
    PyObject *list2_obj;
    PyObject *list3_obj;
    if (!PyArg_ParseTuple(args, "OO", &list2_obj, &list3_obj)) return NULL;

    double **list2;
    double ***list3;

    // For the final version
    double **final_array2;
    double **final_array2;

    // For loops
    int i,j;

    //Create C arrays from numpy objects:
    int typenum = NPY_DOUBLE;
    PyArray_Descr *descr;
    descr = PyArray_DescrFromType(typenum);

    // One per array coming back ...
    npy_intp dims2[2];
    npy_intp dims3[3];

    if (PyArray_AsCArray(&list2_obj, (void **)&list2, dims2, 2, descr) < 0 || PyArray_AsCArray(&list3_obj, (void ***)&list3, dims3, 3, descr) < 0) {
        PyErr_SetString(PyExc_TypeError, "error converting to c array");
        return NULL;
    }

    // Create the pointer arrays needed to access the data

    // 2D array
    final_array2 = calloc(dim2[0], sizeof(double *));
    for (i=0; i<dim[0]; i++) final_array2[i] = list2 + dim2[1]*sizeof(double);

    // 2D array
    final_array3    = calloc(dim3[0], sizeof(double **));
    final_array3[0] = calloc(dim3[0]*dim3[1], sizeof(double *));
    for (i=0; i<dim[0]; i++) {
         final_array3[i] = list2 + dim3[1]*sizeof(double *);
         for (j=0; j<dim[1]; j++) {
             final_array[i][j] = final_array[i] + dim3[2]*sizeof(double);
         }
    }

    printf("2D: %f, 3D: %f.\n", final_array2[3][1], final_array3[1][0][2]);
    // Do stuff with the arrays

    // When ready to complete, free the array access stuff
    free(final_array2);

    free(final_array3[0]);
    free(final_array3);

    // I would guess you also need to free the stuff allocated by PyArray_AsCArray, if so:
    free(list2);
    free(list3);
}

I couldn't find a definition for npy_intp, the above assumes it is the same as int. If it isn't you will need to convert dim2 and dim3 into int arrays before doing the code.

Retsina

Not sure about the downvoter. You're right about just creating the pointer, but the calls to PyArray_AsCArray() do the malloc for me. I'm not great in C, so I don't really know why I need to (void **)&list2, but the program crashes with a Bus error if I don't. – Circe

-1: Your answer is incorrect, because the OP does not need to allocate memory for the arrays. read the function definition: docs.scipy.org/doc/numpy-1.3.x/reference/… – Bolger

@Bolger Thanks, I've completely rewritten the answer to cope with this scenario, hopefully now correct. – Retsina

Answer 4

There was a bug in the numpy C-API, that should be fixed now:

https://github.com/numpy/numpy/pull/5314

Doornail