Cython: gli indici int unsigned per gli array numpy danno un risultato diverso

Question

Ho convertito in cython una funzione python aggiungendo solo alcuni tipi e compilandolo. Stavo ottenendo piccole differenze numeriche tra i risultati delle funzioni python e cython. Dopo un po 'di lavoro ho scoperto che le differenze derivavano dall'accedere ad una matrice numpy usando int unsigned invece di int.

ero utilizzando indici int senza segno per accelerare l'accesso secondo: http://docs.cython.org/src/userguide/numpy_tutorial.html#tuning-indexing-further

comunque ho pensato che fosse innocuo per utilizzare unsigned int.

Vedi questo codice:

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])) 

stampe:?

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

perché questo accade !!! E 'un errore?

Ok, come richiesto: ecco uno SSCCE con gli stessi tipi e valori che ho usato nella mia funzione originaria

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])) 

stampe

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

io uso Python 2.7.3 Cython 0,18 e msvc9 express

Answer 1

Ho modificato l'esempio nella domanda per semplificare la lettura della sorgente C generata per il modulo. Mi interessa solo vedere la logica che crea oggetti Python float invece di ottenere oggetti np.float32 dall'array response.

Sto usando pyximport per compilare il modulo di estensione. Salva il file C generato in una sottodirectory di ~/.pyxbld (probabilmente %userprofile%\.pyxbld su 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')) 

Ecco il codice C generato per la sezione di interesse (un po 'riformattato per rendere più facile la lettura). Si scopre che quando si utilizzano le variabili di indice C unsigned int, il codice generato acquisisce i dati direttamente dal buffer di matrice e chiama PyFloat_FromDouble, che lo costringe a double. D'altra parte, quando si usano le variabili di indice Python int, prende l'approccio generico. Forma una tupla e chiama PyObject_GetItem. In questo modo è possibile che ndarray rispetti correttamente il dtype np.float32.

#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

Giocare con questo sulla mia macchina, non vedo la differenza. Sto usando il notebook con la magia ipython Cython:

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 

La prima coppia di risultati, la differenza è solo la precisione di uscita del comunicato stampa. Quale versione di Cython stai usando? È estremamente improbabile che gli indici effettuino la risposta poiché sta semplicemente accedendo a una lunghezza fissa di memoria che l'attributo dei dati della serie numpy sta memorizzando.