Skip to content

Instantly share code, notes, and snippets.

@jjerphan
Last active June 29, 2021 12:09
Show Gist options
  • Save jjerphan/4906c0a0ae35c76e83b409741fa8c13a to your computer and use it in GitHub Desktop.
Save jjerphan/4906c0a0ae35c76e83b409741fa8c13a to your computer and use it in GitHub Desktop.
Python interaction when using memory views as attributes on Cython classes
%%cython --annotate
#cython: boundscheck=False
#cython: wraparound=False
#cython: cdivision=True
## Adapted from sklearn.neighbors.Mahalanobis
# https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/neighbors/_dist_metrics.pyx#L669
import numpy as np
cimport numpy as np
np.import_array()
ctypedef np.float64_t DTYPE_t
ctypedef np.intp_t ITYPE_t
from libc.math cimport sqrt
DTYPE = np.float64
ITYPE = np.intp
cdef class SnippedMahalanobisDistance():
cdef:
DTYPE_t[:, ::1] mat
DTYPE_t[::1] vec
DTYPE_t p
ITYPE_t size
def __cinit__(self):
self.p = 2
self.vec = np.zeros(1, dtype=DTYPE, order='c')
self.mat = np.zeros((1, 1), dtype=DTYPE, order='c')
self.size = 1
def __init__(self, V=None, VI=None):
if VI is None:
if V is None:
raise ValueError("Must provide either V or VI "
"for Mahalanobis distance")
VI = np.linalg.inv(V)
if VI.ndim != 2 or VI.shape[0] != VI.shape[1]:
raise ValueError("V/VI must be square")
self.mat = np.asarray(VI, dtype=float, order='C')
self.size = self.mat.shape[0]
# we need vec as a work buffer
self.vec = np.zeros(self.size, dtype=DTYPE)
cdef inline DTYPE_t rdist(self,
const DTYPE_t* x1,
const DTYPE_t* x2,
ITYPE_t size) nogil except -1:
cdef DTYPE_t tmp, d = 0
cdef np.intp_t i, j
# compute (x1 - x2).T * VI * (x1 - x2)
# There is interaction with CPython from here...
for i in range(size):
self.vec[i] = x1[i] - x2[i]
for i in range(size):
tmp = 0
for j in range(size):
tmp += self.mat[i, j] * self.vec[j]
d += tmp * self.vec[i]
# ... to there
return d
cdef inline DTYPE_t snipped_mahalanobis(self,
DTYPE_t[:, ::1] mat,
DTYPE_t[::1] vec,
const DTYPE_t* x1,
const DTYPE_t* x2,
ITYPE_t size) nogil except -1:
# A variant which does not have any interaction
# with CPython
cdef DTYPE_t tmp, d = 0
cdef np.intp_t i, j
# compute (x1 - x2).T * VI * (x1 - x2)
for i in range(size):
vec[i] = x1[i] - x2[i]
for i in range(size):
tmp = 0
for j in range(size):
tmp += mat[i, j] * vec[j]
d += tmp * vec[i]
return d
@jjerphan
Copy link
Author

jjerphan commented Jun 29, 2021

Generated code for SnippedMahalanobisDistance.rdist:

+47:     cdef inline DTYPE_t rdist(self,
 48:                               const DTYPE_t* x1,
 49:                               const DTYPE_t* x2,
 50:                               ITYPE_t size) nogil except -1:
+51:         cdef DTYPE_t tmp, d = 0
 52:         cdef np.intp_t i, j
 53: 
 54:         # compute (x1 - x2).T * VI * (x1 - x2)
 55: 
 56:         # There is interaction with CPython from here...
+57:         for i in range(size):
+58:             self.vec[i] = x1[i] - x2[i]
    if (unlikely(!__pyx_v_self->vec.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 58, __pyx_L1_error)}
    __pyx_t_4 = __pyx_v_i;
    *((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) ( /* dim=0 */ ((char *) (((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) __pyx_v_self->vec.data) + __pyx_t_4)) )) = ((__pyx_v_x1[__pyx_v_i]) - (__pyx_v_x2[__pyx_v_i]));
  }
 59: 
+60:         for i in range(size):
+61:             tmp = 0
+62:             for j in range(size):
+63:                 tmp += self.mat[i, j] * self.vec[j]
      if (unlikely(!__pyx_v_self->mat.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 63, __pyx_L1_error)}
      __pyx_t_7 = __pyx_v_i;
      __pyx_t_8 = __pyx_v_j;
      if (unlikely(!__pyx_v_self->vec.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 63, __pyx_L1_error)}
      __pyx_t_9 = __pyx_v_j;
      __pyx_v_tmp = (__pyx_v_tmp + ((*((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) ( /* dim=1 */ ((char *) (((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) ( /* dim=0 */ (__pyx_v_self->mat.data + __pyx_t_7 * __pyx_v_self->mat.strides[0]) )) + __pyx_t_8)) ))) * (*((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) ( /* dim=0 */ ((char *) (((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) __pyx_v_self->vec.data) + __pyx_t_9)) )))));
    }
+64:             d += tmp * self.vec[i]
    if (unlikely(!__pyx_v_self->vec.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 64, __pyx_L1_error)}
    __pyx_t_4 = __pyx_v_i;
    __pyx_v_d = (__pyx_v_d + (__pyx_v_tmp * (*((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) ( /* dim=0 */ ((char *) (((__pyx_t_54_cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9_DTYPE_t *) __pyx_v_self->vec.data) + __pyx_t_4)) )))));
  }
 65: 
 66:         # ... to there
+67:         return d

Screenshot 2021-06-29 at 11-45-14 Cython _cython_magic_a64c371eaaa29cbf2ef07aea5d70c4ef7f1ec7b9 pyx

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment