[Cython] Passing numpy arrays to C using Cython

[Dag Sverre Seljebotn]

First: This question (how to safely pass numpy arrays to C or Fortran 
code using Cython) actually comes up a lot. I'm wondering about perhaps 
making it a feature in Cython itself, borrowing a few syntax things from 
Fortran. (That would clean up a lot of code I'm using myself with Cython 
these days, so it could get realized in some months).

Basically, perhaps one should allow declaring functions like this:

cdef extern from "myheader.h":
   int myfunc(double[0:nrows, 0:ncols] data, int nrows, int ncols)

(here data is really "double* data" in the C header) and then call it 
like this:

cdef np.ndarray[double, ndim=2] arr = x
myfunc(arr)

nrows and ncols would be inferred from the shape of "arr", and Cython 
would transparently make sure contiguous memory is passed, making a copy 
if necesarry. (One could leave out the auto-nrows-ncols-bit for now though.)

This would be backed by the buffer PEP and not be NumPy-specific.

The thing lacking from the syntax above is how to specify whether data 
only needs to be passed in, or whether it is written to. "const 
double[...]" perhaps...

Thoughts?

Then to your question:

Ilmar Wilbers wrote:
> Dag Sverre Seljebotn wrote:
>> Hi Ilmar,
>>
>> Ilmar Wilbers wrote:
>>   
>>> Hi list,
>>>
>>> While trying to wrap an external C library, I have stumbled upon some
>>> question. A lot of them were answered by going through the mailing list
>>> archive, however a few remain.
>>>
>>> I have the following function with arguments:
>>> def tridag(np.ndarray[DTYPE_t, ndim=1, negative_indices=False] a)
>>>
>>> I want to do some assertions on the array a:
>>>     assert a.dtype == np.float32
>>>     assert a.flags['C_CONTIGUOUS']
>>> but keep getting an error  "Attempting to index non-array type 'int'".
>>> Checking the type of a.flags indicates that it is in int. What am I
>>> missing?
>>>     
>> You have encountered one of the less polished aspects of Cython I'm
>> afraid. The flags field is typed as int to give efficient access for
>> or-style flags checking, which precludes the object behaviour you use.
>>
>> To make a long story short, here's your current options:
>> - assert (<object>a).flags[....]
>> - The assertion for datatype is already made if DTYPE_t is float32_t. Use
>> mode="c" as a flag to ndarray to do the other assertion automatically (on
>> entering the function).
>>
>> This is definitely something that should be more done something about
>> (i.e. this should Just Work), I'll try to file an enhancement ticket for
>> it.
>>
>>   
> OK, thank you, Dag Sverre.
> 
> I realize that the mode is not really important for 1D arrays, but was 
> merely working my way up to matrices.

Mode is important for 1D arrays in this setting -- if you have strided 
memory, you cannot pass it on to C. I.e.

x = arange(10)
x = x[::2]

Here no memory is changed on the second line, only a stride is set 
skipping every other element.

> This requires the explicit use of a = zeros(n, numpy.float32) in the 
> calling Python code, which I was hoping to avoid.

You can also use arr.astype(), which will make a copy without filling 
with zeros. If that isn't flexible enough (cannot set C or Fortran 
mode?), one can also use np.empty, which only allocates memory and 
doesn't fill with zeros either.

> My question is, can I use an alternative to
> cdef float * p_a = <np.float32_t*>a.data
> 
> in case the arrays have dtype float64, for instance, for simply writing
> cdef float * p_a = <np.float64_t*>a.data
> 
> doesn't help. Aka: is it possible to not make assumptions on the array 
> dtype? Using np.float without 32 or 64 works fine, but the data sent to 
> the C function is wrong. I guess float in C has the same precision as 
> float32 in numpy and double in C has the same precision as float64 in numpy?

If you simply write "np.ndarray[double] arr", you WILL get a runtime 
exception if arr is assigned something that does not contain a C double. 
  So this is already handled automatically and you do not need to code 
the check yourself -- simply do

cdef np.ndarray[double] a = x
cdef double * p_a = <double*>a.data

Yes, usually double == float64.


-- 
Dag Sverre