数组下标错误：无效类型‘float[int]’及传递变量到scipy.weave.inline

Question

我最近在玩Scipy的内联工具（通过weave），觉得挺有意思，但遇到了一些麻烦。我的C语言有点生疏，感觉可能漏掉了什么简单的东西。

下面这个函数是用来处理一个三维的float32类型的numpy数组。我正在使用一大堆网格化的气象数据，但这个函数应该适用于任何三维数组。它的作用是对每个点j,k（也就是说，如果i是时间轴，j和k是经纬度，那么我就是在对每个网格点的时间进行平均）在i轴上计算算术平均值。

我希望我的代码能够做到这一点，并且避免numpy中的NaN值（我相信isnan()在内联C/C++中是可以用的...？）。但是，不管代码是否能正常工作，我在编译时总是遇到一些错误，比如：

tools.py: In function ‘PyObject* compiled_func(PyObject*, PyObject*)’:
tools.py:93:45: error: invalid types ‘float[int]’ for array subscript
tools.py:95:51: error: invalid types ‘float[int]’ for array subscript
tools.py: In function ‘PyObject* compiled_func(PyObject*, PyObject*)’:
tools.py:93:45: error: invalid types ‘float[int]’ for array subscript
tools.py:95:51: error: invalid types ‘float[int]’ for array subscript

我觉得我在声明和初始化方面做得还不错，所以可能是某些东西没有以我想的方式传递给weave？如果有人能帮我解决这个问题，我会非常感激。以下是这个函数的代码：

from scipy.weave import inline

def foo(x):
    xi = np.shape(x)[0]
    xj = np.shape(x)[1]
    xk = np.shape(x)[2]
    code = """
           #line 87 "tools.py"
           int n;
           float out[xj][xk];
           for (int k = 0; k < xk; k++) {
               for (int j = 0; j < xj; j++) {
                   n = 0;
                   for (int i = 0; i < xi; i++) {
                       if (!isnan(x[i][j][k])) {
                           n += 1;
                           out[j][k] += x[i][j][k];
                       }
                   }
                   out[j][k] = out[j][k]/n;
               }
           }
           return_val = out;
           """
    awesomeness = inline(code, ['x', 'xi', 'xj', 'xk'], compiler = 'gcc')
    return(awesomeness)

Answer 1

C语言不支持像out[xj][xk]这样的动态数组大小。你要么得把数组的大小写死，也就是固定大小，要么就得用malloc这样的函数，或者使用Weave支持的其他方法来动态分配数据。

Answer 2

你可以先在Python中创建一个叫做out的数组，然后把它传给C++。在C++中，你可以通过Nx[0]、Nx[1]和Nx[2]来获取x的形状。你还可以使用为数组定义的宏来访问它的元素。例如，X3(k,j,i)在C++中和Python里的x[k,j,i]是一样的，而OUT2(j,i)在C++中和out[j,i]在Python里也是一样的。你可以查看自动生成的C++代码，了解可以用来操作数组的变量和宏。要获取C++代码的文件夹：

from scipy import weave
print weave.catalog.default_dir()

我的编译器不支持isnan()这个函数，所以我用tmp==tmp来检查。

# -*- coding: utf-8 -*-
import scipy.weave as weave
import numpy as np

def foo(x):
    out = np.zeros(x.shape[1:])
    code = """
    int i,j,k,n;
    for(i=0;i<Nx[2];i++)
    {
        for(j=0;j<Nx[1];j++)
        {
            n = 0;
            for(k=0;k<Nx[0];k++)
            {
                double tmp = X3(k,j,i);
                if(tmp == tmp) // if isnan() is not available
                {
                    OUT2(j,i) += tmp;
                    n++;
                }
            }
            OUT2(j,i) /= n;
        }
    }
    """
    weave.inline(code, ["x","out"], headers=["<math.h>"], compiler="gcc")
    return out

np.random.seed(0)    
x = np.random.rand(3,4,5)
x[0,0,0] = np.nan
mx = np.ma.array(x, mask=np.isnan(x))
avg1 = foo(x)
avg2 = np.ma.average(mx, axis=0)
print np.all(avg1 == avg2)

你还可以使用blitz转换器来在C++中访问数组。想了解更多细节，可以去谷歌搜索：weave.converters.blitz