cython function output slightly different from python function output
cython, numpy, python
Solution
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 uses`PyObject_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)
Problem
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?!