To learn something new, I'm currently trying to redefine the numpy.mean () function in C. It needs to take a 3D array and return a 2D array with the average value of the elements along the 0 axis. I control to calculate the average of all the values, but I don’t know how I will return the new array in Python. From the sources that I read, I think that this is due to serious juggling with pointers, etc., with which I am not very familiar (but ready to become one).
My code is:
#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(); }
What I understand so far (and please correct me if I'm wrong), I need to create a new PyArrayObject that will be my output (maybe with PyArray_FromDims?). Then I need an array of addresses in the memory of this array and fill it with data. How can i do this?
EDIT:
After reading pointers in more detail (here: http://pw1.netcom.com/~tjensen/ptr/pointers.htm ), I achieved what I was aiming for. Another question arises: where to find the original implementation of numpy.mean ()? I would like to see how this happens, that the python operation is much faster than my version. I guess he escapes the ugly cycle.
Here is my solution:
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); }