Передача 3-мерного массива numpy на C

Question

Я пишу расширение C для моей программы Python для скорости и запускаю какое-то очень странное поведение, пытающееся перейти в 3-мерный массив numpy. Он работает с 2-мерным массивом, но я уверен, что я что-то прикрутил указателями, пытаясь заставить его работать с 3-м измерением. Но вот странная часть. Если я просто перехожу в трехмерный массив, он сбой с Ошибка шины. Если (в Python) я сначала создаю свою переменную как 2D-массив, а затем перезаписываю ее с помощью 3D-массива, работает отлично. Если переменная является пустым массивом сначала, а затем 3D-массив, он сбой с Seg Fault. Как это может случиться?

Кроме того, может ли кто-нибудь помочь мне получить 3D-массив? Или я должен просто сдаться и пройти в 2D-массив и изменить его сам?

Вот мой C код:

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]); 
} 

А вот мой Python код, который вызывает функции выше:

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) 

Так что, если я закомментируйте обе линии А и В, он падает с Ошибка шины. Если строка A есть, но строка B закомментирована, она работает правильно, без ошибок. Если строка B есть, но строка A закомментирована, она печатает правильные номера, но затем Seg faults. Наконец, если обе строки присутствуют, они также печатают правильные цифры, а затем Seg faults. Что, черт возьми, здесь происходит?

EDIT: Ok. Вау. Поэтому я использовал int в Python, но набрал их double в C. И это отлично работало с 1D и 2D массивами. Но не 3D. Поэтому я изменил определение Python l3 на float, и теперь все работает фантастически (Большое спасибо Bi Rico).

Но теперь более странное поведение с линиями A & B! Теперь, если обе строки закомментированы, программа работает. Если присутствует строка B, но A закомментирован, она работает, а если обе раскоментированы. Но если строка A присутствует и B закомментирован, я снова получаю эту фантастическую ошибку шины. Я бы очень хотел избежать этого в будущем, так и есть ли у кого-нибудь понять, почему объявление переменной Python может иметь такое влияние?

EDIT 2:. Ну, как безумные, как эти ошибки, они все из-за 3-мерный массив Numpy я прохожу в случае, если я прохожу только в 1- или 2-D массивов, он ведет себя как ожидается, и манипуляции с другими переменными Python ничего не делают. Это заставляет меня думать, что проблема лежит где-то в подсчете ссылок Python. В C-коде счетчик ссылок уменьшается больше, чем нужно для трехмерных массивов, и когда эта функция возвращает Python пытается очистить объекты и пытается удалить указатель NULL. Это только моя догадка, и я пытался сделать Py_INCREF(); все, что я мог придумать безрезультатно. Я думаю, я просто использовать 2D-массив и видоизменение его в С.

Answer 1

Я уже упоминал об этом в комментарии, но, надеюсь, его немного промыть поможет сделать его более понятным.

Когда вы работаете с массивами numpy в C, хорошо говорить о типизации ваших массивов. В частности, похоже, что вы указали свои указатели как double ***list3, но они, как вы создаете l3 в своем коде на python, вы получите массив с dtype npy_intp (я думаю). Вы можете исправить это, явно используя dtype при создании своих массивов.

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) 

Еще одно замечание, потому как питон работает почти невозможно для «линии А» и «В линии», чтобы иметь никакого влияния на код C, что так всегда. Я знаю, что это, похоже, противоречит вашему эмпирическому опыту, но я уверен в этом.

Я немного менее уверен в этом, но, основываясь на моем опыте с C, ошибки шины и segfaults не детерминированы. Они зависят от распределения памяти, выравнивания и адресов. В какой-то ситуации код, кажется, работает нормально 10 раз и не работает на 11-м запуске, хотя ничего не изменилось.

Считаете ли вы использование cython? Я знаю, что это не вариант для всех, но если это вариант, вы можете получить почти ускорение уровня C, используя typed memoryviews.

Answer 2

Согласно http://docs.scipy.org/doc/numpy/reference/c-api.array.html?highlight=pyarray_ascarray#PyArray_AsCArray:

Примечание Моделирование массива C-стиля не является полным для 2-й и 3- d массивов. Например, моделируемые массивы указателей не могут быть переданы подпрограммам, ожидающим конкретные, статически заданные массивы с 2-го и 3-мерным массивами. Чтобы перейти к функциям, требующим такого рода входов, вы должны статически определить требуемый массив и скопировать данные.

Я думаю, что это означает, что PyArray_AsCArray возвращает блок памяти с данными в нем в порядке C. Однако для доступа к этим данным требуется дополнительная информация (см. http://www.phy225.dept.shef.ac.uk/mediawiki/index.php/Arrays,_dynamic_array_allocation). Это можно достичь, зная размеры заблаговременно, объявив массив, а затем скопировав данные в нужном порядке. Однако я подозреваю, что более общий случай более полезен: вы не знаете размеры до тех пор, пока они не будут возвращены. Я думаю, что следующий код создаст необходимую C-указательную структуру C, чтобы разрешить данные.

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); 
} 

Я не мог найти определение для npy_intp, выше предполагает, что это то же самое, как int. Если это не так, вам нужно будет преобразовать dim2 и dim3 в int массивы перед выполнением кода.

Answer 3

Вместо преобразования в массив c-style я обычно обращаюсь к элементам массива numpy напрямую, используя PyArray_GETPTR (см. http://docs.scipy.org/doc/numpy/reference/c-api.array.html#data-access).

Например, для доступа к элементу трехмерного массива numpy типа double use double elem=*((double *)PyArray_GETPTR3(list3_obj,i,j,k)).

Для вашего приложения вы можете определить правильное количество измерений для каждого массива с помощью PyArray_NDIM, затем обращаться к элементам с использованием соответствующей версии PyArray_GETPTR.

Answer 4

Существовал ошибка в Numpy C-API, который должен быть исправлена:

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