Cython overhead of numpy arrays

Question

I wrote a simple interpolation function in cython, to be called (a lot) from other cython code. One of the parameters is a numpy array:

@cython.boundscheck(False) 
@cython.cdivision(True)
@cython.wraparound(False) 
cdef double interpU_cython(double kX,double kY,int iX,int iY, int iTheta,int nbX,int nbY,np.ndarray[double, ndim=3] u, double outVal):
    cdef double uPt, u0, u1
    if (iX >= 0 and iY >= 0 and iX < nbX-1 and iY < nbY-1):
        u0 = u[iX,iY,iTheta] + (u[iX+1,iY,iTheta]-u[iX,iY,iTheta]) * kX
        u1 = u[iX,iY+1,iTheta] + (u[iX+1,iY+1,iTheta]-u[iX,iY+1,iTheta]) * kX
        uPt = u0 + (u1-u0) * kY
    else:
        uPt = outVal
    return uPt

I checked the python calls with cython -a, and it looks like the function call relies on several python calls:

+01: cimport cython
 02: cimport numpy as np
+03: import numpy as np
 04: 
 05: @cython.boundscheck(False)
 06: @cython.cdivision(True)
 07: @cython.wraparound(False)
+08: cdef double interpU_cython(double kX,double kY,int iX,int iY, int iTheta,int nbX,int nbY,np.ndarray[double, ndim=3] u, double outVal):
static double __pyx_f_10FSM_cython_interpU_cython(double __pyx_v_kX, double __pyx_v_kY, int __pyx_v_iX, int __pyx_v_iY, int __pyx_v_iTheta, int __pyx_v_nbX, int __pyx_v_nbY, PyArrayObject *__pyx_v_u, double __pyx_v_outVal) {
  double __pyx_v_uPt;
  double __pyx_v_u0;
  double __pyx_v_u1;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_u;
  __Pyx_Buffer __pyx_pybuffer_u;
  double __pyx_r;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("interpU_cython", 0);
  __pyx_pybuffer_u.pybuffer.buf = NULL;
  __pyx_pybuffer_u.refcount = 0;
  __pyx_pybuffernd_u.data = NULL;
  __pyx_pybuffernd_u.rcbuffer = &__pyx_pybuffer_u;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_u.rcbuffer->pybuffer, (PyObject*)__pyx_v_u, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 8, __pyx_L1_error)
  }
  __pyx_pybuffernd_u.diminfo[0].strides = __pyx_pybuffernd_u.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_u.diminfo[0].shape = __pyx_pybuffernd_u.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_u.diminfo[1].strides = __pyx_pybuffernd_u.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_u.diminfo[1].shape = __pyx_pybuffernd_u.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_u.diminfo[2].strides = __pyx_pybuffernd_u.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_u.diminfo[2].shape = __pyx_pybuffernd_u.rcbuffer->pybuffer.shape[2];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_u.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_WriteUnraisable("FSM_cython.interpU_cython", __pyx_clineno, __pyx_lineno, __pyx_filename, 0, 0);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_u.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
 09:     cdef double uPt, u0, u1
+10:     if (iX >= 0 and iY >= 0 and iX < nbX-1 and iY < nbY-1):
+11:         u0 = u[iX,iY,iTheta] + (u[iX+1,iY,iTheta]-u[iX,iY,iTheta]) * kX
+12:         u1 = u[iX,iY+1,iTheta] + (u[iX+1,iY+1,iTheta]-u[iX,iY+1,iTheta]) * kX
+13:         uPt = u0 + (u1-u0) * kY
 14:     else:
+15:         uPt = outVal
+16:     return uPt

Is there an efficient way to pass and use numpy arrays with no significant overhead, or should I just use c arrays for everything in the c compiled parts of the code?

Is it considered safe to use the pointer to the numpy array first element and add array size in the parameters to use it as a one-dimensional array?

Thanks

Answer 1

Have a look at https://jakevdp.github.io/blog/2012/08/08/memoryview-benchmarks/

This blog post compares several possibilities to use NumPy arrays in Cython.

The short answer is that you should use typed memoryviews that are declared as double[:,:,:] u instead of np.ndarray[double, ndim=3] u. Docs: http://docs.cython.org/en/latest/src/userguide/memoryviews.html

Edit: you can query the shape of a memoryview