0

I have a simple numerical function for cython. Following the numpy tutorial, I got

import numpy as np
from cython.view cimport array as cvarray


cpdef double[:, ::1] invert_xyz(double[:, ::1] xyz, double[:] center):
    """
    Inversion operation on `xyz` with `center` as inversion center.

    :param xyz:    Nx3 coordinate array
    """
    cdef Py_ssize_t i, n
    n = xyz.shape[0]
    # see https://stackoverflow.com/questions/18462785/what-is-the-recommended-way-of-allocating-memory-for-a-typed-memory-view
    got = cvarray(shape=(n, 3), itemsize=sizeof(double), format="d")
    cdef double[:, ::1] mv = got

    for i in range(n):
        mv[i, 0] = 2 * center[0] - xyz[i, 0]
        mv[i, 1] = 2 * center[1] - xyz[i, 1]
        mv[i, 2] = 2 * center[2] - xyz[i, 2]
    return mv

But the lines inside the for loop still contain python interactions (they are yellow when compiled with annotation). For example,

+032:     for i in range(n):
+033:         mv[i, 0] = 2 * center[0] - xyz[i, 0]
    __pyx_t_8 = 0;
    __pyx_t_9 = -1;
    if (__pyx_t_8 < 0) {
      __pyx_t_8 += __pyx_v_center.shape[0];
      if (unlikely(__pyx_t_8 < 0)) __pyx_t_9 = 0;
    } else if (unlikely(__pyx_t_8 >= __pyx_v_center.shape[0])) __pyx_t_9 = 0;
    if (unlikely(__pyx_t_9 != -1)) {
      __Pyx_RaiseBufferIndexError(__pyx_t_9);
      __PYX_ERR(0, 33, __pyx_L1_error)
    }
    __pyx_t_10 = __pyx_v_i;
    __pyx_t_11 = 0;
    __pyx_t_9 = -1;
    if (__pyx_t_10 < 0) {
      __pyx_t_10 += __pyx_v_xyz.shape[0];
      if (unlikely(__pyx_t_10 < 0)) __pyx_t_9 = 0;
    } else if (unlikely(__pyx_t_10 >= __pyx_v_xyz.shape[0])) __pyx_t_9 = 0;
    if (__pyx_t_11 < 0) {
      __pyx_t_11 += __pyx_v_xyz.shape[1];
      if (unlikely(__pyx_t_11 < 0)) __pyx_t_9 = 1;
    } else if (unlikely(__pyx_t_11 >= __pyx_v_xyz.shape[1])) __pyx_t_9 = 1;
    if (unlikely(__pyx_t_9 != -1)) {
      __Pyx_RaiseBufferIndexError(__pyx_t_9);
      __PYX_ERR(0, 33, __pyx_L1_error)
    }
    __pyx_t_12 = __pyx_v_i;
    __pyx_t_13 = 0;
    __pyx_t_9 = -1;
    if (__pyx_t_12 < 0) {
      __pyx_t_12 += __pyx_v_mv.shape[0];
      if (unlikely(__pyx_t_12 < 0)) __pyx_t_9 = 0;
    } else if (unlikely(__pyx_t_12 >= __pyx_v_mv.shape[0])) __pyx_t_9 = 0;
    if (__pyx_t_13 < 0) {
      __pyx_t_13 += __pyx_v_mv.shape[1];
      if (unlikely(__pyx_t_13 < 0)) __pyx_t_9 = 1;
    } else if (unlikely(__pyx_t_13 >= __pyx_v_mv.shape[1])) __pyx_t_9 = 1;
    if (unlikely(__pyx_t_9 != -1)) {
      __Pyx_RaiseBufferIndexError(__pyx_t_9);
      __PYX_ERR(0, 33, __pyx_L1_error)
    }
    *((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_mv.data + __pyx_t_12 * __pyx_v_mv.strides[0]) )) + __pyx_t_13)) )) = ((2.0 * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_center.data) + __pyx_t_8)) )))) - (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_xyz.data + __pyx_t_10 * __pyx_v_xyz.strides[0]) )) + __pyx_t_11)) ))));
+034:         mv[i, 1] = 2 * center[1] - xyz[i, 1]

How can I get rid of the python interactions?

nos
  • 19,875
  • 27
  • 98
  • 134
  • 1
    You just need to read a bit further in the tutorial: https://cython.readthedocs.io/en/latest/src/userguide/numpy_tutorial.html#tuning-indexing-further - this looks like bounds-checking – DavidW Aug 12 '23 at 06:00
  • Certainly wraparound yes (not bound-checking though both are quite similar). Note that the layout is inefficient. So it will make your Cython code quite slow, and Numpy too. You should transpose the array. Most people use `(N,3)` coordinate preventing SIMD vectorization while `(3,N)` is nearly always (much) faster. – Jérôme Richard Aug 12 '23 at 11:15
  • Thank you guys. Adding the directives indeed remove the python interactions – nos Aug 12 '23 at 13:11

0 Answers0