Cython: unsigned int indices dla numpy array daje inny wynik

Question

Przekonwertowałem do cytonu funkcję Pythona po prostu dodając niektóre typy i kompilując je. Otrzymywałem małe różnice liczbowe między wynikami funkcji python i cython. Po pracy odkryłem, że różnice wynikają z dostępu do tablicy numpy przy użyciu unsigned int zamiast int.

używałem unsigned int indeksów w celu przyspieszenia dostępu według: http://docs.cython.org/src/userguide/numpy_tutorial.html#tuning-indexing-further

i tak myślałem, że to nieszkodliwe używać niepodpisanych ints.

Zobacz ten kod:

cpdef function(np.ndarray[np.float32_t, ndim=2] response, max_loc): 
    cdef unsigned int x, y 
    x, y = int(max_loc[0]), int(max_loc[1]) 
    x2, y2 = int(max_loc[0]), int(max_loc[1]) 
    print response[y,x], type(response[y,x]), response.dtype 
    print response[y2,x2], type(response[y2,x2]), response.dtype 
    print 2*(response[y,x] - min(response[y,x-1], response[y,x+1])) 
    print 2*(response[y2,x2] - min(response[y2,x2-1], response[y2,x2+1])) 

drukuje:

0.959878861904 <type 'float'> float32 
0.959879 <type 'numpy.float32'> float32 
1.04306024313 
1.04306030273 

Dlaczego tak się dzieje !!! czy to błąd?

Ok, wnioskowane Oto SSCCE z tych samych typów i wartości, które użyłem w mojej pierwotnej funkcji

cpdef function(): 
    cdef unsigned int x, y 
    max_loc2 = np.asarray([ 15., 25.], dtype=float) 
    cdef np.ndarray[np.float32_t, ndim=2] response2 = np.zeros((49,49), dtype=np.float32)  
    x, y = int(max_loc2[0]), int(max_loc2[1]) 
    x2, y2 = int(max_loc2[0]), int(max_loc2[1]) 

    response2[y,x] = 0.959878861904 
    response2[y,x-1] = 0.438348740339 
    response2[y,x+1] = 0.753262758255 


    print response2[y,x], type(response2[y,x]), response2.dtype 
    print response2[y2,x2], type(response2[y2,x2]), response2.dtype 
    print 2*(response2[y,x] - min(response2[y,x-1], response2[y,x+1])) 
    print 2*(response2[y2,x2] - min(response2[y2,x2-1], response2[y2,x2+1])) 

drukuje

0.959878861904 <type 'float'> float32 
0.959879 <type 'numpy.float32'> float32 
1.04306024313 
1.04306030273 

używam Python 2.7.3 Cython 0,18 i msvc9 express

Answer 1

Zmodyfikowałem przykład w pytaniu, aby ułatwić odczytanie wygenerowanego źródła C dla modułu. Interesuje mnie tylko logika, która tworzy obiekty Python float zamiast obiektów np.float32 z tablicy .

Używam pyximport do kompilacji modułu rozszerzenia. Zapisuje wygenerowany plik C w podkatalogu ~/.pyxbld (prawdopodobnie %userprofile%\.pyxbld w systemie Windows).

import numpy as np 
import pyximport 
pyximport.install(setup_args={'include_dirs': [np.get_include()]}) 

open('_tmp.pyx', 'w').write(''' 
cimport numpy as np 
cpdef function(np.ndarray[np.float32_t, ndim=2] response, max_loc): 
    cdef unsigned int p_one, q_one 
    p_one = int(max_loc[0]) 
    q_one = int(max_loc[1]) 
    p_two = int(max_loc[0]) 
    q_two = int(max_loc[1]) 
    r_one = response[q_one, p_one] 
    r_two = response[q_two, p_two] 
''') 

import _tmp 
assert(hasattr(_tmp, 'function')) 

Oto wygenerowany kod C dla interesującego nas odcinka (nieco sformatowany, aby ułatwić czytanie). Okazuje się, że gdy używasz zmiennych indeksowych C unsigned int, wygenerowany kod pobiera dane bezpośrednio z bufora tablicowego i wywołuje PyFloat_FromDouble, co wymusza na double. Z drugiej strony, gdy używasz zmiennych indeksowych Python int, przyjmuje to podejście ogólne. Tworzy krotkę i wywołuje PyObject_GetItem. W ten sposób pozwala ndarray poprawnie honorować typ dtype.

#define __Pyx_BufPtrStrided2d(type, buf, i0, s0, i1, s1) \ 
    (type)((char*)buf + i0 * s0 + i1 * s1) 

    /* "_tmp.pyx":9 
*  p_two = int(max_loc[0]) 
*  q_two = int(max_loc[1]) 
*  r_one = response[q_one, p_one]    # <<<<<<<<<<<<<< 
*  r_two = response[q_two, p_two] 
*/ 
    __pyx_t_3 = __pyx_v_q_one; 
    __pyx_t_4 = __pyx_v_p_one; 
    __pyx_t_5 = -1; 

    if (unlikely(__pyx_t_3 >= (size_t)__pyx_bshape_0_response)) 
    __pyx_t_5 = 0; 
    if (unlikely(__pyx_t_4 >= (size_t)__pyx_bshape_1_response)) 
    __pyx_t_5 = 1; 

    if (unlikely(__pyx_t_5 != -1)) { 
    __Pyx_RaiseBufferIndexError(__pyx_t_5); 
    { 
     __pyx_filename = __pyx_f[0]; 
     __pyx_lineno = 9; 
     __pyx_clineno = __LINE__; 
     goto __pyx_L1_error; 
    } 
    } 

    __pyx_t_1 = PyFloat_FromDouble((
    *__Pyx_BufPtrStrided2d(
     __pyx_t_5numpy_float32_t *, 
     __pyx_bstruct_response.buf, 
     __pyx_t_3, __pyx_bstride_0_response, 
     __pyx_t_4, __pyx_bstride_1_response))); 

    if (unlikely(!__pyx_t_1)) { 
    __pyx_filename = __pyx_f[0]; 
    __pyx_lineno = 9; 
    __pyx_clineno = __LINE__; 
    goto __pyx_L1_error; 
    } 

    __Pyx_GOTREF(__pyx_t_1); 
    __pyx_v_r_one = __pyx_t_1; 
    __pyx_t_1 = 0; 

    /* "_tmp.pyx":10 
*  q_two = int(max_loc[1]) 
*  r_one = response[q_one, p_one] 
*  r_two = response[q_two, p_two]    # <<<<<<<<<<<<<< 
*/ 
    __pyx_t_1 = PyTuple_New(2); 

    if (unlikely(!__pyx_t_1)) { 
    __pyx_filename = __pyx_f[0]; 
    __pyx_lineno = 10; 
    __pyx_clineno = __LINE__; 
    goto __pyx_L1_error; 
    } 

    __Pyx_GOTREF(((PyObject *)__pyx_t_1)); 
    __Pyx_INCREF(__pyx_v_q_two); 
    PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_q_two); 
    __Pyx_GIVEREF(__pyx_v_q_two); 
    __Pyx_INCREF(__pyx_v_p_two); 
    PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_v_p_two); 
    __Pyx_GIVEREF(__pyx_v_p_two); 

    __pyx_t_2 = PyObject_GetItem(
    ((PyObject *)__pyx_v_response), 
    ((PyObject *)__pyx_t_1)); 

    if (!__pyx_t_2) { 
    __pyx_filename = __pyx_f[0]; 
    __pyx_lineno = 10; 
    __pyx_clineno = __LINE__; 
    goto __pyx_L1_error; 
    } 

    __Pyx_GOTREF(__pyx_t_2); 
    __Pyx_DECREF(((PyObject *)__pyx_t_1)); 
    __pyx_t_1 = 0; 
    __pyx_v_r_two = __pyx_t_2; 
    __pyx_t_2 = 0; 

Answer 2

Zabawa z tym na moim komputerze, nie widzę różnicy. Używam notebook ipython z Cython magii:

In [1]: 

%load_ext cythonmagic 

In [12]: 

%%cython 

import numpy as np 
cimport numpy as np 

cpdef function(np.ndarray[np.float32_t, ndim=2] response, max_loc): 
    cdef unsigned int x, y 
    x, y = int(max_loc[0]), int(max_loc[1]) 
    x2, y2 = int(max_loc[0]), int(max_loc[1]) 
    #return 2*(response[y,x] - min(response[y,x-1], response[y,x+1])), 2*(response[y2,x2] - min(response[y2,x2-1], response[y2,x2+1])) 
    print response[y,x], type(response[y,x]), response.dtype 
    print response[y2,x2], type(response[y2,x2]), response.dtype 
    print 2*(response[y,x] - min(response[y,x-1], response[y,x+1])) 
    print 2*(response[y2,x2] - min(response[y2,x2-1], response[y2,x2+1])) 

In [13]: 

a = np.random.normal(size=(10,10)).astype(np.float32) 
m = [3,2] 
function(a,m) 

0.586090564728 <type 'float'> float32 
0.586091 <type 'numpy.float32'> float32 
4.39655685425 
4.39655685425 

Pierwsza para wyników, różnica jest tylko precyzja Wyjście instrukcji print. Jakiej wersji Cython używasz? Indeksy są bardzo mało prawdopodobne, aby wpłynęły na odpowiedź, ponieważ ma ona tylko dostęp do stałej długości pamięci, która jest przechowywana w atrybucie danych tablicy numpy.