10

何か新しいことを学ぶために、現在 C で numpy.mean() 関数を再実装しようとしています。3D 配列を取り、軸 0 に沿った要素の平均を含む 2D 配列を返す必要があります。すべての値の平均ですが、新しい配列を Python に返す方法がよくわかりません。

これまでの私のコード:

#include <Python.h>
#include <numpy/arrayobject.h>

// Actual magic here:
static PyObject*
myexts_std(PyObject *self, PyObject *args)
{
    PyArrayObject *input=NULL;
    int i, j, k, x, y, z, dims[2];
    double out = 0.0; 

    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &input))
        return NULL;

    x = input->dimensions[0];
    y = input->dimensions[1];
    z = input->dimensions[2];

    for(k=0;k<z;k++){
        for(j=0;j<y;j++){
            for(i=0;i < x; i++){
                out += *(double*)(input->data + i*input->strides[0] 
+j*input->strides[1] + k*input->strides[2]);
            }
        }
    }
    out /= x*y*z;
    return Py_BuildValue("f", out);
}

// Methods table - this defines the interface to python by mapping names to
// c-functions    
static PyMethodDef myextsMethods[] = {
    {"std", myexts_std, METH_VARARGS,
        "Calculate the standard deviation pixelwise."},
    {NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC initmyexts(void)
{
    (void) Py_InitModule("myexts", myextsMethods);
    import_array();
}

私がこれまでに理解していること (間違っている場合は訂正してください) はPyArrayObject、出力となる新しい を作成する必要があるということです (おそらくPyArray_FromDims? を使用)。次に、この配列のメモリへのアドレスの配列が必要で、データを入力します。これについてどうすればいいですか？

編集：

ポインターについてさらに読んだ後 (ここ: http://pw1.netcom.com/~tjensen/ptr/pointers.htm )、私が目指していたものを達成しました。ここで別の疑問が生じます: numpy.mean() の元の実装はどこにありますか? Pythonの操作が私のバージョンよりもはるかに高速であることを確認したいと思います。醜いループを回避すると思います。

これが私の解決策です：

static PyObject*
myexts_std(PyObject *self, PyObject *args)
{
    PyArrayObject *input=NULL, *output=NULL; // will be pointer to actual numpy array ?
    int i, j, k, x, y, z, dims[2]; // array dimensions ?
    double *out = NULL;
    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &input))
        return NULL;

    x = input->dimensions[0];
    y = dims[0] = input->dimensions[1];
    z = dims[1] = input->dimensions[2];
    output = PyArray_FromDims(2, dims, PyArray_DOUBLE);    
    for(k=0;k<z;k++){
        for(j=0;j<y;j++){
            out = output->data + j*output->strides[0] + k*output->strides[1];
            *out = 0;
            for(i=0;i < x; i++){
                *out += *(double*)(input->data + i*input->strides[0] +j*input->strides[1] + k*input->strides[2]);
            }
            *out /= x;
        }
    }
    return PyArray_Return(output);
}
4

1 に答える 1

3

The Numpy API has a function PyArray_Mean that accomplishes what you're trying to do without the "ugly looping" ;).

static PyObject *func1(PyObject *self, PyObject *args) {
    PyArrayObject *X, *meanX;
    int axis;

    PyArg_ParseTuple(args, "O!i", &PyArray_Type, &X, &axis);
    meanX = (PyArrayObject *) PyArray_Mean(X, axis, NPY_DOUBLE, NULL);

    return PyArray_Return(meanX);
}
于 2013-09-16T16:57:00.810 に答える