将3维numpy数组传递给C

Question

我正在为我的Python程序编写C扩展以提高速度，并且在尝试传递一个三维numpy数组时遇到了一些非常奇怪的行为。它可以使用二维数组，但我确信在尝试使用第三个维度时出了问题。但是这里有一个奇怪的部分，如果我只传入一个三维数组，它会崩溃并显示总线错误。 如果（在Python中）我首先将变量创建为2D数组，然后用3D数组覆盖它，它完美地工作。如果变量首先为空数组，然后是3D数组，则会崩溃并显示段错误。这怎么可能发生？

此外，有人能帮助我使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两行注释掉，程序会出现总线错误。如果只有A行而B行被注释掉，它可以正常运行没有错误。如果只有B行而A行被注释掉，它会输出正确的数字但随后出现段错误。最后，如果两行都存在，它也会输出正确的数字并然后发生段错误。这到底是怎么回事？

编辑：好的，哇，所以我在Python中使用了int，但在C中调用它们时使用了double。对于1D和2D数组，这很好用，但不适用于3D。所以我改变了Python的l3定义为浮点数，现在一切都工作得很好（非常感谢Bi Rico）。

但现在，使用A和B两行出现了更奇怪的行为！如果两行都被注释掉，程序就可以正常工作。如果B行存在而A行被注释掉，它也可以正常工作，如果两行都未被注释掉，则同样如此。但是如果A行存在但B行被注释掉，我又得到了那个奇妙的总线错误。我真的很想避免这种问题，那么有人知道为什么Python变量的声明会产生这种影响吗？

编辑2：嗯，尽管这些错误很疯狂，但它们都是由我传递的3维numpy数组引起的。如果只传递1或2维数组，则行为符合预期，并且对其他Python变量的操作没有任何影响。这使我认为问题出在Python的引用计数上。在C代码中，3D数组的引用计数减少得比应该多，当函数返回时Python试图清除对象，并尝试删除一个空指针。这只是我的猜测，我已经尝试过尽可能多地Py_INCREF();，但无济于事。我想我只能使用2D数组并在C中重新调整形状。

- DaveTheScientist

---

你确定 (void **) 是正确的吗？难道不应该只传入 (void*) 吗？ - seberg

我的 C 程序很烂，但是……如果第一次调用 PyArray_AsCArray 成功，你在 if 中的表达式不是会短路吗？很可能第二次调用，也就是 list3 的那个，根本没有执行。 - Jaime

@seberg 我不能确定(void **)是否正确，但使用(void*)会导致总线错误。 @Jaime 不，该函数仅在失败时返回负值，最有可能是它调用的malloc失败了。 - DaveTheScientist

@seberg 好的...现在我采纳了 Bi Rico 的建议并尝试了 Python 浮点数，无论是单星号还是双星号（或三星号）都能正常工作。你有什么想法哪个更好/正确？ - DaveTheScientist

出现以下错误："error: ‘NPY_DOUBLE’ undeclared (first use in this function)"。 - vineeshvs

以上错误的解决方案：包含“numpy/arrayobject.h”。参考：https://mail.python.org/pipermail/scipy-user/2011-October/030801.html - vineeshvs

Answer 1

与转换为c风格数组不同，我通常直接使用PyArray_GETPTR访问numpy数组元素（请参阅https://numpy.org/doc/stable/reference/c-api/array.html#data-access）。

例如，要访问类型为double的三维numpy数组的元素，请使用double elem=*((double *)PyArray_GETPTR3(list3_obj,i,j,k))。

对于您的应用程序，您可以使用PyArray_NDIM检测每个数组的正确维数，然后使用适当版本的PyArray_GETPTR访问元素。

- user1149913

---

我想转换成常规的C数组，因为我认为它会更快。我还以为它会更简单，但显然是错误的... - DaveTheScientist

这个是慢还是快，你有什么想法吗？ - Mike Vella

这会慢一点，但不是很明显。GETPTR是一个简单的宏，你可以在numpy的github下的numpy/core/numpy/include/numpy/ndarrayobject.h中看到它的定义。假设我们有PyArrayObject *myArray和nb =每个元素的字节数。这将返回一个指针，指向地址myArray->data + (i * myArray->strides[0] / (bytes per element)) + (j * myArray->strides[1] / (bytes per element)) + (k * myArray->strides[2] / (bytes per element))。与在内部循环中执行指针递增相比，这存在一定的性能损失，但非常小，不过如果在内部循环中执行多次，则可能会累加。 - HappyDog

Answer 2

我之前在评论中提到过这一点，但希望进一步阐明。

当您在C语言中使用numpy数组时，最好明确您的数组的类型。具体来说，看起来您正在将指针声明为double ***list3，但是通过在Python代码中创建l3，您将获得一个带有dtypenpy_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)

另外需要注意的是，由于Python的工作方式，"A行"和"B行"几乎不可能对C代码产生任何影响。我知道这似乎与你的经验相矛盾，但我非常确定这一点。

根据我的C语言经验，总线错误和段错误并不是确定性的。它们取决于内存分配、对齐和地址。在某些情况下，即使没有任何变化，代码似乎跑了10次都很好，但第11次会出现错误。

您是否考虑过使用Cython？我知道它并非适用于每个人，但如果可以使用，您可以使用类型化的内存视图来获得接近C级别的加速。

- Bi Rico

---

下次我需要编写C扩展时，我相当确定我会花时间学习Cython。是的，根据我对Python和C的了解，"Line A and B"不应该对C程序产生任何影响，因为每次声明L2时都会获得一个新的内存地址。但是对我来说它们确实有影响，这是我提出这个问题的一个重要原因。如果其他人愿意在他们的系统上尝试，我可以粘贴整个文件，因为我很想弄清楚其中的原因。 - DaveTheScientist

Answer 3

现在numpy C-API中的一个错误已经修复：

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

- Hensing

Answer 4

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

注意：对于二维和三维数组，C风格数组的模拟不完整。例如，指针的模拟数组不能传递给期望特定的静态定义的二维和三维数组的子例程。要传递给需要这些类型输入的函数，必须静态定义所需的数组并复制数据。

我认为这意味着PyArray_AsCArray返回一个按C顺序排列的数据块内存。但是，要访问该数据，需要更多信息（请参见http://www.phy225.dept.shef.ac.uk/mediawiki/index.php/Arrays,_dynamic_array_allocation）。可以通过提前知道维度，声明数组，然后按正确顺序将数据复制到其中来实现此目的。然而，我怀疑更一般的情况更有用：您不知道维度直到它们被返回。我认为以下代码将创建必要的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 数组后再执行代码。

- Neil Townsend

---

不确定downvoter怎么想的。你关于只创建指针的做法是对的，但是对PyArray_AsCArray()的调用可以为我分配内存。我在C语言方面并不擅长，所以我不知道为什么需要(void **)&list2，但是如果我不这样做，程序就会因为总线错误而崩溃。 - DaveTheScientist

-1：您的答案不正确，因为 OP 不需要为数组分配内存。请阅读函数定义：http://docs.scipy.org/doc/numpy-1.3.x/reference/c-api.array.html#PyArray_AsCArray - meyumer

@meyumer 谢谢，我已经完全重写了答案以应对这种情况，希望现在是正确的。 - Neil Townsend