将三维numpy数组传递给C

Question

我正在为我的Python程序写一个C扩展，主要是为了提高速度，但在传递一个三维的numpy数组时遇到了一些奇怪的问题。传递二维数组时一切正常，但我觉得在处理指针时，三维数组的部分我搞错了。不过更奇怪的是，如果我直接传递一个三维数组，程序就会崩溃，出现总线错误。而如果我在Python中先创建一个二维数组，然后再用三维数组覆盖它，这就能正常工作。如果变量一开始是一个空数组，然后变成三维数组，程序就会崩溃，出现段错误。这到底是怎么回事呢？

另外，有人能帮我让三维数组正常工作吗？还是我干脆放弃，传递一个二维数组，然后自己在里面调整一下？

这是我的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。这在处理一维和二维数组时没问题，但三维数组就不行了。所以我把Python中l3的定义改成了浮点数，现在一切都正常了（非常感谢Bi Rico）。

但是现在，A行和B行又出现了奇怪的行为！如果两行都注释掉，程序能正常工作。如果B行在但A行被注释掉，程序也能正常运行；如果两行都取消注释，程序也能正常运行。但是如果A行在而B行被注释掉，我又遇到了那个可怕的总线错误。我真的希望以后能避免这些错误，有人知道为什么Python变量的声明会有这样的影响吗？

编辑2：虽然这些错误看起来很疯狂，但它们都是因为我传入的三维numpy数组。如果我只传入一维或二维数组，程序就会按预期运行，其他Python变量的操作也不会有影响。这让我觉得问题可能出在Python的引用计数上。在C代码中，三维数组的引用计数减少得比应该的要多，当那个函数返回时，Python尝试清理对象，并试图删除一个空指针。这只是我的猜测，我尝试过Py_INCREF();所有我能想到的东西，但都没有效果。我想我还是会使用二维数组，然后在C中调整它。

Answer 1

根据这个链接的内容：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数组。

Answer 2

与其把数据转换成C语言风格的数组，我通常直接用 PyArray_GETPTR 来访问numpy数组里的元素（详细信息可以查看 这个链接）。

举个例子，如果你想访问一个三维的numpy数组中的某个元素，类型是double，可以这样写：double elem=*((double *)PyArray_GETPTR3(list3_obj,i,j,k))。

在你的应用中，你可以用 PyArray_NDIM 来检测每个数组的维度，然后用合适的 PyArray_GETPTR 版本来访问元素。

Answer 3

我之前在评论中提到过这个，但我希望详细说明一下，让它更清楚。

当你在C语言中使用numpy数组时，明确数组的类型是很重要的。具体来说，你似乎把指针声明为 double ***list3，但在你的Python代码中创建 l3 的方式，会得到一个类型为 npy_intp 的数组（我想是这样）。你可以通过在创建数组时明确指定类型来解决这个问题。

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语言的速度提升。