Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to declare an ndarray in cython with a general floating point type

What is the best way to declare a numpy array in cython if It should be able to handle both float and double?

I guess it won't be possible with a memory view since there the datatype is crucial, but for an ndarray is there any way to give it a general float type which would still benifit from the swiftness of cython?

so this is what I would usualy do:

def F( np.ndarray A):
    A += 10

I've seen that there is also:

def F( np.ndarray[np.float32_t, ndim=2] A):
    A += 10

but that again will give a bit size for the type. I've also thought along the lines of creating a memory view inside the function depending on the bit size (32 or 64).

Any thought are highly appreciated


Thank you so much for the tip on the floating type. I've tried it like this

import numpy as np
cimport numpy as np
import cython
cimport cython
from libc.math cimport sqrt, abs
from cython cimport floating

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def Rot_Matrix(np.ndarray[floating, ndim=3] Fit_X,
               np.ndarray[floating, ndim=3] Ref_X,
               weight = None):
    cdef:
        unsigned int t, T = Fit_X.shape[0]
        unsigned int n, N = Fit_X.shape[1]
        np.ndarray[floating, ndim=3] Rot = np.empty((T,3,3))

    return Rot

when I now call the function with two arrays of np.float32 I get the error

ValueError: Buffer dtype mismatch, expected 'float' but got 'double'

If I do not use the the type definition in the brakes for Rot so it reads np.ndarray[floating, ndim=3] Rot = np.empty((T,3,3)) then I get ndarray back and it works fine. Would you happen to have a pointer for me what I'm doing wrong?

like image 215
Magellan88 Avatar asked Apr 28 '14 21:04

Magellan88


1 Answers

Well this is actually really easy with fused types support:

This goes inside your code.

from cython cimport floating

def cythonfloating_memoryview(floating[:, :] array):
    cdef int i, j

    for i in range(array.shape[0]):
        for j in range(array.shape[1]):
            array[i, j] += 10

Of course, there are loads of ways of doing this:

Name this fuzed.pyx. There's no need to compile or run cython on it; it's handled by pyximport. Don't use pyximport for production code, though, as you should typically only ship .c files.

from cython cimport floating
from numpy import float32_t, float64_t, ndarray

ctypedef fused myfloating:
    float32_t
    float64_t

def cythonfloating_memoryview(floating[:, :] array):
    # ...

def cythonfloating_buffer(ndarray[floating, ndim=2] array):
    # ...

def myfloating_memoryview(myfloating[:, :] array):
    # ...

def myfloating_buffer(ndarray[myfloating, ndim=2] array):
    # ...

and here's a little test script:

Name this test.py and run it as a normal Python script:

import pyximport
pyximport.install()

import fuzed
import numpy

functions = [
    fuzed.cythonfloating_memoryview,
    fuzed.cythonfloating_buffer,
    fuzed.myfloating_memoryview,
    fuzed.myfloating_buffer,
]

for function in functions:
    floats_32 = numpy.zeros((100, 100), "float32")
    doubles_32 = numpy.zeros((100, 100), "float64")

    function(floats_32)
    function(doubles_32)

    print(repr(floats_32))
    print(repr(doubles_32))

It's worth noting that fused types are specialised at compilation, and are constant for a particular function invocation. The empty Numpy array you make is always of a double type, but you assign it to either a 32-bit float or a 64-bit float. Here's what you should do:

from cython cimport floating
import numpy

def do_some_things(floating[:] input):
    cdef floating[:] output

    if floating is float:
        output = numpy.empty(10, dtype="float32")
    elif floating is double:
        output = numpy.empty(10, dtype="float64")
    else:
        raise ValueError("Unknown floating type.")

    return output

and some tests to prove the point:

import pyximport
pyximport.install()
#>>> (None, None)

import floatingtest
import numpy

floatingtest.do_some_things(numpy.empty(10, dtype="float32"))
#>>> <MemoryView of 'ndarray' at 0x7f0a514b3d88>
floatingtest.do_some_things(numpy.empty(10, dtype="float32")).itemsize
#>>> 4

floatingtest.do_some_things(numpy.empty(10, dtype="float64"))
#>>> <MemoryView of 'ndarray' at 0x7f0a514b3d88>
floatingtest.do_some_things(numpy.empty(10, dtype="float64")).itemsize
#>>> 8
like image 83
Veedrac Avatar answered Sep 29 '22 11:09

Veedrac