The output of the cython function is slightly different from the output of the python function

Question

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

I found out some reasons for this difference in this Cython post : unsigned int indices for numpy arrays give different results But even with what I learned in this post, I still can't get the cython function to produce the same results as python.

So, I put together 4 functions illustrating what I tried. Can someone help to reveal why I get slightly different results for each function? and how to get a cython function that returns the exact exact values ​​as function1? I make a few comments below:

%%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, yx, 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, yx, 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 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 that I performed in my original python function. The result is tmp1.

function2 is my first version of cython, which gives slightly different results. Apparently, if the response array is indexed with a typed variable, unsigned int or int, the result is forced to double (using PyFloat_FromDouble), even if the array type is np.float32_t. But if the array is indexed using python int, the PyObject_GetItem function is used instead, and I get np.float32_t, which happens in function1. Thus, the expressions in function 1 are evaluated using np.float32_t operands, while the expressions in function 2 are calculated using doubles. I get a slightly different fingerprint than in function1.

function3 is my second cython attempt, trying to get the same result as function1. Here I use unsigned int indexes to access the response of the array, but the results remain on the np.float32_t intermediate variables, which are then used in the calculation. I get a slightly different result. Obviously, the print statement will use PyFloat_FromDouble, so it will not be able to print the np.float32_t file.

Then I tried to change the python function to match cython. function4 tries to achieve this by converting the float to at least one operand in each expression so that the rest of the operands are forcibly applied to the pyat float, which is double in cython, and the expression is evaluated using doubles, as in function2. The fingerprints inside the function are exactly the same as function2, but the return values ​​are slightly different ?!

Answer 1

If you use single-point floats that have only 7.225 decimal digits of precision, I would not expect a slight deviation from coercion to double.

To clarify your description of function2 , if you are indexing an object, Cython uses PyObject_GetItem to get the scalar object np.float32 (and not np.float32_t , which is just a type for C float ). If you pointed directly to the buffer, and Cython needs an object, it calls PyFloat_FromDouble . He needs objects to assign tmp1 , tmp2 and tmp3 , since they are not printed.

In function3 , on the other hand, you entered tmp variables, but you still need to create float objects to print and return the results. If you use NumPy ndarray (see below), you will not have this problem:

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

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

against.

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

Even if you guarantee that all float32 operations in the def and cpdef , the end result may still change between the two in the compiled extension module. In the following example, I checked that the intermediate results in function1 are np.float32 objects. In the generated C function2 I checked that there was no translation on double (or equivalent typedef). However, these two functions still give slightly different results. I probably would have to dive into the assembly assembly to understand why, but maybe I will miss 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 passes float32_t to the end.
• function2 converts to float when indexing whether it performs intermediate steps with float , and then returns to float32_t for final results.
• function3 converts to float , and then immediately returns to float32_t , performing intermediate steps with this.
• function4 converts to float whether it performs the intermediate steps, and then returns the final results as a float .

As for why function4 prints the same thing as function2 , but returns something else: if you look at types, it's simple. The values ​​are apparently close enough to match print equally, but not so close to repr . Which is not surprising, given that they are not of the same type.