Passing list of numpy arrays to C using cython

Question

I have a list list_of_arrays of 3D numpy arrays that I want to pass to a C function with the template

int my_func_c(double **data, int **shape, int n_arrays)

such that

data[i]  : pointer to the numpy array values in list_of_arrays[i]
shape[i] : pointer to the shape of the array in list_of_arrays[i] e.g. [2,3,4]

How can I call my_func_c using a cython interface function?

My first idea was to do something like below (which works) but I feel there is a better way just using numpy arrays without mallocing and freeing.

# my_func_c.pyx

import numpy as np
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, free

cdef extern from "my_func.c":
    double my_func_c(double **data, int **shape, int n_arrays)

def my_func(list list_of_arrays):
    cdef int n_arrays  = len(list_of_arrays)
    cdef double **data = <double **> malloc(n_arrays*sizeof(double *))
    cdef int **shape   = <int **> malloc(n_arrays*sizeof(int *))
    cdef double x;

    cdef np.ndarray[double, ndim=3, mode="c"] temp

    for i in range(n_arrays):
        temp = list_of_arrays[i]
        data[i]  = &temp[0,0,0]
        shape[i] = <int *> malloc(3*sizeof(int))
        for j in range(3):
            shape[i][j] = list_of_arrays[i].shape[j]

    x = my_func_c(data, shape, n_arrays)

    # Free memory
    for i in range(n_arrays):
        free(shape[i])
    free(data)
    free(shape)

    return x

N.B.

To see a working example we can use a very simple function calculating the product of all the arrays in our list.

# my_func.c

double my_func_c(double **data, int **shape, int n_arrays) {
    int array_idx, i0, i1, i2;

    double prod = 1.0;

    // Loop over all arrays
    for (array_idx=0; array_idx<n_arrays; array_idx++) {
        for (i0=0; i0<shape[array_idx][0]; i0++) {
            for (i1=0; i1<shape[array_idx][1]; i1++) {
                for (i2=0; i2<shape[array_idx][2]; i2++) {
                    prod = prod*data[array_idx][i0*shape[array_idx][1]*shape[array_idx][2] + i1*shape[array_idx][2] + i2];
                }
            }
        }
    }

    return prod;
}

Create the setup.py file,

# setup.py

from distutils.core import setup
from Cython.Build import cythonize
import numpy as np

setup(
    name='my_func',
    ext_modules = cythonize("my_func_c.pyx"),
    include_dirs=[np.get_include()]
    )

Compile

python3 setup.py build_ext --inplace

Finally we can run a simple test

# test.py

import numpy as np
from my_func_c import my_func

a = [1+np.random.rand(3,1,2), 1+np.random.rand(4,5,2), 1+np.random.rand(1,2,3)]

print('Numpy product: {}'.format(np.prod([i.prod() for i in a])))
print('my_func product: {}'.format(my_func(a)))

using

python3 test.py

It seems the question is not clear enough but I'm not sure exactly how to edit it. The output in this case is essentially what we want to see. The cython function gives the same result as numpy (with some small rounding error). Its stochastic so changes each time. What I wanted was a way to deal with lists of numpy arrays in cython without having to resort to mallocing and freeing memory. — rwolst
– rwolst, Commented Sep 15, 2017 at 14:26
@rwolst You could pass a pointer to numpy's inbuilt shape arrays rather than copying them element by element. Other than that I think your approach is pretty good. — DavidW
– DavidW, Commented Sep 15, 2017 at 17:49

DavidW · Accepted Answer · 2017-09-17 21:04:27Z

One alternative would be to let numpy manage your memory for you. You can do this by using numpy arrays of np.uintp which is an unsigned int with the same size as any pointer.

Unfortunately, this does require some type-casting (between "pointer sized int" and pointers) which is a good way of hiding logic errors, so I'm not 100% happy with it.

def my_func(list list_of_arrays):
    cdef int n_arrays  = len(list_of_arrays)
    cdef np.uintp_t[::1] data = np.array((n_arrays,),dtype=np.uintp)
    cdef np.uintp_t[::1] shape = np.array((n_arrays,),dtype=np.uintp)
    cdef double x;

    cdef np.ndarray[double, ndim=3, mode="c"] temp

    for i in range(n_arrays):
        temp = list_of_arrays[i]
        data[i]  = <np.uintp_t>&temp[0,0,0]
        shape[i] = <np.uintp_t>&(temp.shape[0])

    x = my_func_c(<double**>(&data[0]), <np.intp_t**>&shape[0], n_arrays)

(I should point out that I've only confirmed that it compiles and not tested it further, but the basic idea should be OK)

The way you've done it is probably a pretty sensible way. One slight simplification to your original code that should work

shape[i] = <np.uintp_t>&(temp.shape[0])

instead of malloc and copy. I'd also recommend putting the frees in a finally block to ensure they get run.

Edit: @ead has helpfully pointed out that the numpy shape is stored as as np.intp_t - i.e. an signed integer big enough to fit a pointer in, which is mostly 64bit - while int is usually 32 bit. Therefore, to pass the shape without copying you'd need to change your C api. Casting help makes that mistake harder to spot ("a good way of hiding logic errors")

@ead's suggestion of using C++ vectors is probably a reasonable alternative too - it has some of the advantages of my "numpy array" solution with less messy casting. I'll let him/her post it though...
In the original solution one also should check, that malloc was successful, i.e. the result isn't 0. So if c++ isn't an option, my trade-off would be casts (your solution) over "many things I should remember and do right".
actually, a point against casting would be: shape-entries are 8byte-integers, but after the cast they will be reinterpreted as 4byte-integers (at least on the most usual systems).
@ead Which is probably why OP was copying them I guess. My issue with casting is mostly that it makes it easy to make that kind of mistake. Thanks for pointing that out
So you proved yourself right:) I think a safer way to do it would be to reinterpret the reserved memory first: cdef np.uintp_t[::1] shape_mem = np.array((n_arrays,),dtype=np.uintp); cdef int** shape=<int**>&shape_mem[0] and then enjoy the advantages of the type safety: now shape[i] = &temp.shape[0] would be a compile error if temp.shape[0] wasn't an int.

ead · Accepted Answer · 2017-09-18 04:31:01Z

I think this is a good pattern to consume C-functionality from C++-code, and it can be also used here and would have two advantages:

Memory management is taken care of.
Thanks to templates, no casting needed, so we still have the safety-net of c's type safety.

To solve your problems you could use std::vector:

import numpy as np
cimport numpy as np
from libcpp.vector cimport vector

cdef extern from "my_func.c":
    double my_func_c(double **data, int **shape, int n_arrays)

def my_func(list list_of_arrays):
    cdef int n_arrays  = len(list_of_arrays)
    cdef vector[double *] data
    cdef vector [vector[int]] shape_mem # for storing casted shapes
    cdef vector[int *] shape  #pointers to stored shapes
    cdef double x
    cdef np.ndarray[double, ndim=3, mode="c"] temp

    shape_mem.resize(n_arrays)  
    for i in range(n_arrays):
        print "i:", i
        temp = list_of_arrays[i]
        data.push_back(&temp[0,0,0])
        for j in range(3):
            shape_mem[i].push_back(temp.shape[j])
        shape.push_back(shape_mem[i].data())

    x = my_func_c(data.data(), shape.data(), n_arrays)

    return x

Also your setup would need a modification:

# setup.py    
from distutils.core import setup, Extension
from Cython.Build import cythonize
import numpy as np

setup(ext_modules=cythonize(Extension(
            name='my_func_c',
            language='c++',
            extra_compile_args=['-std=c++11'],
            sources = ["my_func_c.pyx", "my_func.c"],
            include_dirs=[np.get_include()]
    )))

I prefer to use std::vector.data() over &data[0] because the second would mean undefined behavior for empty data, and that is the reason we need std=c++11 flag.

But in the end, it is for you to decide, which trade-off to make: the additional complexity of C++ (it has it own pitfalls) vs. handmade memory management vs. letting go of type safety for a short moment.

Collectives™ on Stack Overflow

Passing list of numpy arrays to C using cython

2 Answers 2

5 Comments

Comments

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

5 Comments

Comments

Linked

Related