Cython: unsigned int индексы для массивов numpy дают отличный результат

Question

Я преобразовал в cython функцию python, просто добавив некоторые типы и скомпилировав их. Я получал небольшие числовые различия между результатами функций python и cython. После некоторой работы я обнаружил, что различия связаны с доступом к массиву numpy с использованием unsigned int вместо int.

Я использую неподписанных Int индексов для ускорения доступа в соответствии с: http://docs.cython.org/src/userguide/numpy_tutorial.html#tuning-indexing-further

во всяком случае я думал, что это безвредно для применения неподписанных Ints.

Смотреть этот код:

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

печатает:

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

Почему это случилось !!! это ошибка?

Ok, в соответствии с просьбой здесь является SSCCE с теми же типами и значениями, которые я использовал в своей первоначальной функции

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

печатает

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

Я использую Python 2.7.3 Cython 0,18 и msvc9 экспресс

Answer 1

Я изменил пример в вопросе, чтобы упростить чтение сгенерированного источника C для модуля. Меня интересует только логика, которая создает объекты Python float вместо получения np.float32 объектов из массива response.

Я использую pyximport для компиляции модуля расширения. Он сохраняет сгенерированный файл C в подкаталоге ~/.pyxbld (возможно, %userprofile%\.pyxbld в 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')) 

Вот сгенерированный код C для секции интересов (немного переформатировать, чтобы сделать его более удобным для чтения). Оказывается, когда вы используете переменные индекса C unsigned int, сгенерированный код захватывает данные непосредственно из буфера массива и вызывает PyFloat_FromDouble, который принуждает его к double. С другой стороны, когда вы используете переменные индекса Python int, он принимает общий подход. Он формирует кортеж и вызывает PyObject_GetItem. Этот способ позволяет ndarray правильно оценивать тип 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

Играя с этим на моей машине, я не вижу разницы. Я использую 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 

Первая пара результатов, разница только точность выхода заявления печати. Какую версию Cython вы используете? Индексы крайне маловероятны для ответа, поскольку он просто обращается к фиксированной длине памяти, которую хранит атрибут data массива numpy.