将 3 维 numpy 数组传递给 C

Question

为了提高速度，我正在为我的 Python 程序编写一个 C 扩展，并且在尝试传入 3 维 numpy 数组时遇到了一些非常奇怪的行为。它适用于二维数组，但我确信我正在用指针试图让它与第 3 维一起工作。但这是奇怪的部分。如果我只是传入一个 3-D 数组，它会因 Bus Error 而崩溃。如果（在 Python 中）我首先将变量创建为 2D 数组，然后用 3D 数组覆盖它，它可以完美运行。如果变量先是空数组，然后是 3D 数组，则会因 Seg Fault 而崩溃。这怎么可能发生？

另外，谁能帮我让 3D 阵列正常工作？还是我应该放弃并传入一个二维数组并自己重塑它？

这是我的 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 行被注释掉，它会打印正确的数字，但随后会出现 Seg 错误。最后，如果两行都存在，它还会打印正确的数字，然后是 Seg 错误。这到底是怎么回事？

编辑：好的。哇。所以我在 Python 中使用int，但在 C 中调用它们为double。这适用于一维和二维数组。但不是 3D。所以我将 l3 的 Python 定义更改为浮点数，现在一切正常（非常感谢 Bi Rico）。

但是现在，A 行和 B 行出现了更奇怪的行为！现在，如果两行都被注释掉，程序就可以工作了。如果 B 行存在，但 A 被注释掉，它可以工作，如果两者都未注释，则同上。但是，如果 A 行存在并且 B 被注释掉，我会再次收到那个奇妙的总线错误。我真的很想在以后避免这些，所以有人知道为什么 Python 变量的声明会产生这种影响吗？

编辑 2： 好吧，尽管这些错误很疯狂，但它们都是由于我传入的 3 维 numpy 数组。如果我只传入 1 维或 2 维数组，它的行为与预期一样，并且对其他 Python 变量的操作什么也不做。这让我相信问题出在 Python 的引用计数中。在 C 代码中，引用计数比 3-D 数组减少的要多，当该函数返回时，Python 会尝试清理对象，并尝试删除 NULL 指针。这只是我的猜测，我已经尝试Py_INCREF(); 我能想到的一切都无济于事。我想我只会使用二维数组并用 C 对其进行整形。

Answer 1

我通常使用PyArray_GETPTR（参见http://docs.scipy.org/doc/numpy/reference/c-api.array.html#data-access）直接访问numpy 数组元素，而不是转换为c 样式的数组。

例如，要访问 double 类型的 3 维 numpy 数组的元素 double elem=*((double *)PyArray_GETPTR3(list3_obj,i,j,k)).

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

Answer 2

我已经在评论中提到了这一点，但我希望将其刷新一下有助于使其更清楚。

当您在 C 中使用 numpy 数组时，最好明确说明数组的类型。具体来说，您似乎将指针声明为double ***list3，但是您在python代码中创建l3的方式您将获得一个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)

另一个注意事项，由于 python 的工作方式，“line A”和“line B”几乎不可能对 C 代码产生任何影响。我知道这似乎与您的经验经验相冲突，但我很确定这一点。

对此我不太确定，但根据我在 C 语言方面的经验，总线错误和段错误不是确定性的。它们取决于内存分配、对齐和地址。在某些情况下，代码似乎可以正常运行 10 次，但在第 11 次运行时失败，即使没有任何变化。

您是否考虑过使用cython？我知道这不是每个人的选择，但如果它是一个选项，您可以使用 typed memoryviews 获得接近 C 级的加速。

Answer 3

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

注意对于 2-d 和 3-d 数组，C 样式数组的模拟是不完整的。例如，模拟的指针数组不能传递给期望特定的、静态定义的 2-d 和 3-d 数组的子程序。要传递给需要这些输入的函数，您必须静态定义所需的数组并复制数据。

我认为这意味着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 数组。

Answer 4

numpy C-API 中有一个错误，现在应该修复：

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