cython function output slightly different from python function output

Question

I have converted a python function into a cython equivalent by adding types to some variables. However, the cython function to produces slightly different output than the original python function.

I've learnt some of the reasons for this difference in this post Cython: unsigned int indices for numpy arrays gives different result But even with what I've learnt in this post I still can't get the cython function to produce the same results as the python one.

So I have put together 4 functions illustrating what I have tried. Could someone help to unveil why I get slightly different results for each function? and how to get a cython function that returns the same exact values as function1? I make some comments bellow:

%%cython
import numpy as np
cimport numpy as np    

def function1(response, max_loc):    
    x, y = int(max_loc[0]), int(max_loc[1])

    tmp1 = (response[y,x+1] - response[y,x-1]) / 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))
    tmp2 = (response[y,x+1] - response[y,x-1])
    tmp3 = 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))

    print tmp1, tmp2, tmp3        
    return tmp1, tmp2, tmp3

cpdef function2(np.ndarray[np.float32_t, ndim=2] response, np.ndarray[np.float64_t, ndim=1] max_loc):
    cdef unsigned int x, y 
    x, y = int(max_loc[0]), int(max_loc[1])

    tmp1 = (response[y,x+1] - response[y,x-1]) / 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))        
    tmp2 = (response[y,x+1] - response[y,x-1])
    tmp3 = 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))     

    print tmp1, tmp2, tmp3        
    return tmp1, tmp2, tmp3


cpdef function3(np.ndarray[np.float32_t, ndim=2] response, np.ndarray[np.float64_t, ndim=1] max_loc):     
    cdef unsigned int x, y 
    x, y = int(max_loc[0]), int(max_loc[1])

    cdef np.float32_t tmp1, tmp2, tmp3
    cdef np.float32_t r1 =response[y,x+1]
    cdef np.float32_t r2 =response[y,x-1]
    cdef np.float32_t r3 =response[y,x]
    cdef np.float32_t r4 =response[y,x-1]
    cdef np.float32_t r5 =response[y,x+1]    

    tmp1 = (r1 - r2) / 2*(r3 - min(r4, r5))  
    tmp2 = (r1 - r2)
    tmp3 = 2*(r3 - min(r4, r5))

    print tmp1, tmp2, tmp3        
    return tmp1, tmp2, tmp3

def function4(response, max_loc):     
    x, y = int(max_loc[0]), int(max_loc[1])

    tmp1 = (float(response[y,x+1]) - response[y,x-1]) / 2*(float(response[y,x]) - min(response[y,x-1], response[y,x+1]))
    tmp2 = (float(response[y,x+1]) - response[y,x-1])
    tmp3 = 2*(float(response[y,x]) - min(response[y,x-1], response[y,x+1]))

    print tmp1, tmp2, tmp3        
    return tmp1, tmp2, tmp3

max_loc = np.asarray([ 15., 25.], dtype=np.float64) 
response = np.zeros((49,49), dtype=np.float32)     
x, y = int(max_loc[0]), int(max_loc[1])

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

result1 = function1(response, max_loc)
result2 = function2(response, max_loc)
result3 = function3(response, max_loc)
result4 = function4(response, max_loc)
print result1
print result2
print result3
print result4

and the results:

0.0821185777156 0.314914 1.04306030273
0.082118573023 0.314914017916 1.04306024313
0.0821185708046 0.314914017916 1.04306030273
0.082118573023 0.314914017916 1.04306024313
(0.082118577715618812, 0.31491402, 1.043060302734375)
(0.08211857302303427, 0.3149140179157257, 1.0430602431297302)
(0.08211857080459595, 0.3149140179157257, 1.043060302734375)
(0.082118573023034269, 0.31491401791572571, 1.0430602431297302)

function1 represents the operations I did in my original python function. tmp1 is the result.

function2 is my first cython version, which produces slightly different results. Apparently, if the response array is indexed with a typed variable, unsigned int or int, the result is coerced to a double (using PyFloat_FromDouble) even if the type of the array is np.float32_t. But if the array is indexed with a python int the function PyObject_GetItem is used instead and I get the np.float32_t, which is what happens in function1. So the expressions in function1 are calculated using np.float32_t operands, whereas the expressions in function2 are calculated using doubles. I get slightly different print out than in function1.

function3 is my second cython attempt trying to get the same output as function1. Here I use unsigned int indices to access the array response but the results are left on np.float32_t intermediate variables that then I use in the calculation. I get slightly different result. Apparently the print statement will use PyFloat_FromDouble so it won't be able to print the np.float32_t.

Then I tried to change the python function to match the cython one. function4 tries to achieve this by converting to float at least one operand in each expression so the rest of operands are coerced to python float too, which is a double in cython, and the expression is calculated with doubles, as in function2. The prints inside the function are exactly the same as function2, but the returned values are slightly different?!

Answer 1

If you're using single-precision floats, which only have 7.225 decimal digits of precision, I wouldn't expect a small variance from coercion to double to matter much.

To clarify your description of function2, if you index with an object, Cython usesPyObject_GetItem to get an np.float32 scalar object (not an np.float32_t, which is just a typedef for a C float). If you instead index directly into the buffer, and Cython needs an object, it calls PyFloat_FromDouble. It needs objects to assign tmp1, tmp2, and tmp3, since they aren't typed.

In function3, on the other hand, you typed the tmp variables, but it still needs to create float objects to print and return the results. If you instead use a NumPy ndarray (see below), you won't have that problem:

In function1, by the way, you promote the result to np.float64 when you divide by 2. For example:

>>> type(np.float32(1) / 2)
<type 'numpy.float64'>

vs.

>>> type(np.float32(1) / np.float32(2))
<type 'numpy.float32'>

Even if you ensure that all operations are float32 in both the def and cpdef functions, the end result might still vary between the two in the compiled extension module. In the following example I checked that the intermediate results in function1 are all np.float32 objects. In the generated C of function2 I checked that there was no cast to double (or equivalent typedef). Yet these two functions still produce slightly different results. I'd probably have to dive into the compiled assembly to figure out why, but maybe I'm overlooking something simple.

def function1(response, max_loc):    
    tmp = np.zeros(3, dtype=np.float32)
    x, y = int(max_loc[0]), int(max_loc[1])
    tmp[0] = (((response[y,x+1] - response[y,x-1]) / np.float32(2)) *
             (response[y,x] - min(response[y,x-1], response[y,x+1])))
    tmp[1] = response[y,x+1] - response[y,x-1]
    tmp[2] = 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))

    print tmp[0], tmp[1], tmp[2]
    return tmp

cpdef function2(np.ndarray[np.float32_t, ndim=2] response, max_loc):
    cdef np.ndarray[np.float32_t, ndim=1] tmp = np.zeros(3, dtype=np.float32)
    cdef unsigned int x, y
    x, y = int(max_loc[0]), int(max_loc[1])
    tmp[0] = (((response[y,x+1] - response[y,x-1]) / <np.float32_t>2) *
             (response[y,x] - min(response[y,x-1], response[y,x+1])))
    tmp[1] = response[y,x+1] - response[y,x-1]
    tmp[2] = 2*(response[y,x] - min(response[y,x-1], response[y,x+1]))

    print tmp[int(0)], tmp[int(1)], tmp[int(2)]
    return tmp

Comparison:

>>> function1(response, max_loc)
0.0821186 0.314914 1.04306
array([ 0.08211858,  0.31491402,  1.0430603 ], dtype=float32)

>>> function2(response, max_loc)
0.0821186 0.314914 1.04306
array([ 0.08211857,  0.31491402,  1.0430603 ], dtype=float32)

Answer 2

Let's compare:

• function1 stays float32_t all the way through.
• function2 converts to float upon indexing, does the intermediate steps with float, then converts back to float32_t for the final results.
• function3 converts to float, but then immediately back to float32_t, doing the intermediate steps with that.
• function4 converts to float, does the intermediate steps, then returns the final results as float.

As for why function4 prints the same thing as function2, but returns something different: If you look at the types, it's simple. The values are apparently close enough that they happen to print the same way, but not close enough to repr the same way. Which isn't surprising given that they're not the same type.