C Function in Python: Return an Array and Variables

Question

I just got transited into Python from Matlab and I know that calling C function is different than Matlab mex. With lacks of proper documentation, I have scourge the web for days now and couldn't find the solution to my simple problem.

Basically I want to call a C-Function in Python, input 2 integers and an 2D array, do some calculations, and return a 2D Array. I would also like to output some other variables too (this might require the use of structure). I know this are very basic stuffs but if someone can help me, I will greatly appreciated.

So what I am hoping for is the equivalent of this in matlab! Thx!!!

[Nxy,outArray] = Function(Nx,Ny,inArray)

Code for setup.py

from distutils.core import setup, Extension
import numpy.distutils.misc_util

setup(
    ext_modules=[Extension("myfunc", ["myfunc.c"])],
    include_dirs=numpy.distutils.misc_util.get_numpy_include_dirs(),
)

Code for myfunc.c

static char module_docstring[] =
    "This function does some calculations...";
static char Run_docstring[] =
    "Run what ever algorithm there is!";

static PyObject *Run(PyObject *self, PyObject *args)
{
    int i, j, Nx, Ny;
    PyObject *Data;

    /* Parse the input tuple */
    if (!PyArg_ParseTuple(args, "iiO", &Nx, &Ny, &Data))   // Data is a 2D array  
        return NULL;

    PyObject *array = PyArray_FROM_OTF(Data, NPY_DOUBLE, NPY_IN_ARRAY); // Interpret as numpy array
    double *newData = (double*)PyArray_DATA(array); // Pointers to the data as C-types

    double outData[Ny][Nx]; // Creating output 2D Array
    int outCount;

    // Calculations
    outCount = Nx*Ny;

    for (i=0; i<Nx; i++){
        for (j=0; i<Ny; j++){
           outData[j][i] = sqrt(Data[j][i]) + sqrt(outCount);
        }
    }

    // Free memory used in PyObject
    Py_DECREF(array);

    // Return output Data
    PyObject *ret = Py_BuildValue("i", outCount);
    return ret, PyArray_Return(outData);   

}

static PyMethodDef module_methods[] = {
    {"Run", Run, METH_VARARGS, Run_docstring},
    {NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC initmyfunc(void)
{
    PyObject *m = Py_InitModule3("myfunc", module_methods, module_docstring);
    if (m == NULL)
        return;
    import_array();
}

Answer 1

It is possible to call optimised C functions from Python using Cython.

In this particular case, we can for instance create a myfunc.pyx file,

import numpy as np
cimport numpy as np
from libc.math cimport sqrt

cpdef tuple myfunc(int Nx, int Ny, double[:,::1] inArray):

        cdef double [:,::1] outData = np.zeros((Nx, Ny))
        cdef int i,j, res

        with nogil:
            for i in range(Nx):
                for j in range(Ny):
                    outData[i, j] = sqrt(inArray[i,j]) + sqrt(<double> Nx*Ny)
        res = 0 # not sure how res is computed

        return res, outData.base

that can be compiled with the following setup.py,

from distutils.core import setup, Extension
import numpy as np
from Cython.Distutils import build_ext

setup(
    ext_modules=[Extension("myfunc", ["myfunc.pyx"])],
    cmdclass = {'build_ext': build_ext},
    include_dirs=[np.get_include()])

using

$ python setup.py build_ext --inplace

This generates and compiles the myfunc.c. The resulting Python module can be then used as follows,

from myfunc import myfunc
import numpy as np
Nx, Ny = 2, 2
inArray =  np.ones((Nx,Ny))
res, outArray = myfunc(Ny,Ny, inArray)
print(outArray)
# which would return
[[ 3.  3.]
[ 3.  3.]]

Note, that it in this case it is not necessary to pass the array dimensions Nx, Ny to the function as they can be accessed through inArray.shape in Cython.

Please refer to the Cython documentation regarding Numpy for further optimisation details.