cython 函数输出与 python 函数输出略有不同答案

【问题标题】：cython function output slightly different from python function outputcython 函数输出与 python 函数输出略有不同
【发布时间】：2013-03-13 03:05:50
【问题描述】：

我通过向一些变量添加类型将 python 函数转换为 cython 等价物。但是，cython 函数产生的输出与原始 python 函数略有不同。

我在这篇文章中了解了造成这种差异的一些原因 Cython: unsigned int indices for numpy arrays gives different result 但即使我在这篇文章中学到了什么，我仍然无法让 cython 函数产生与 python 相同的结果。

所以我整理了 4 个函数来说明我尝试过的内容。有人可以帮助揭示为什么每个函数我得到的结果略有不同吗？以及如何获得一个与function1返回相同精确值的cython函数？我在下面做了一些cmets：

%%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

结果：

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 代表我在原始 python 函数中所做的操作。 tmp1 是结果。

function2 是我的第一个 cython 版本，它产生的结果略有不同。显然，如果响应数组使用类型化变量 unsigned int 或 int 进行索引，则即使数组的类型是 np.float32_t，结果也会被强制为 double（使用 PyFloat_FromDouble）。但是，如果数组使用 python int 进行索引，则使用函数 PyObject_GetItem，我得到 np.float32_t，这就是函数 1 中发生的情况。所以 function1 中的表达式是使用 np.float32_t 操作数计算的，而 function2 中的表达式是使用双精度计算的。我得到的打印结果与 function1 中的略有不同。

function3 是我第二次尝试获得与 function1 相同的输出。在这里，我使用 unsigned int 索引来访问数组响应，但结果留在 np.float32_t 中间变量上，然后我在计算中使用它们。我得到的结果略有不同。显然打印语句将使用 PyFloat_FromDouble 所以它不能打印 np.float32_t。

然后我尝试更改 python 函数以匹配 cython 函数。 function4 尝试通过将每个表达式中的至少一个操作数转换为浮点数来实现这一点，因此其余操作数也被强制转换为 python 浮点数，这是 cython 中的双精度数，并且表达式是用双精度数计算的，如函数2。函数内部的print和function2一模一样，只是返回值略有不同？！

【问题讨论】：

最多相差10^-9（十亿分之一）。为什么您会惊讶于不同的实现会有如此大的差异？（在什么应用程序中会导致问题？）
通过 hex()，print ", ".join([x.hex() for x in result4]) 打印浮点值。
函数 1 使用 Python float (IEEE double) 值作为 tmp1、tmp2 和 tmp3。函数 3 将它们显式声明为 np.float32_t（IEEE 单一）。他们怎么可能返回同样的东西？让中间类型匹配，然后使用完全不同的最终类型就达不到目的了。
同时，您到底想完成什么？在几乎所有情况下，尽早而不是稍后切换到float 至少会给您带来同样好的结果，如果不是更好的话。除非您有特定的理由相信 function2 在实际值方面是“错误的”或“更糟的”（正如大卫罗宾逊指出的那样，这可能还不够准确，甚至没有意义），什么是问题？
好的。好吧，目的是我首先开发 python 函数，当它们按照我的意愿运行时，我为它创建一个 cython 函数。所以我希望 cython 函数给我与 python 相同的结果。这些小的差异会随着我之后执行的其他操作累积起来，从而为我的系统产生不同的基准测试结果

标签： python numpy cython

【解决方案1】：

如果您使用单精度浮点数，其精度只有 7.225 位十进制数字，我不认为从强制转换到加倍的小差异很重要。

为了澄清您对function2 的描述，如果您使用对象进行索引，Cython 使用PyObject_GetItem 来获取np.float32 标量对象（不是np.float32_t，它只是C float 的typedef ）。如果您改为直接索引到缓冲区，并且 Cython 需要一个对象，它会调用 PyFloat_FromDouble。它需要对象来分配tmp1、tmp2 和tmp3，因为它们没有被键入。

另一方面，在function3 中，您键入了tmp 变量，但它仍然需要创建float 对象来打印并返回结果。如果你改用 NumPy ndarray（见下文），你就不会有这个问题：

在function1中，顺便说一下，当你除以2时，你将结果提升为np.float64。例如：

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

对比

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

即使您确保在def 和cpdef 函数中的所有操作都是float32，最终结果在编译的扩展模块中仍然可能在两者之间有所不同。在以下示例中，我检查了 function1 中的中间结果是否都是 np.float32 对象。在function2 生成的C 中，我检查了没有强制转换为double（或等效的typedef）。然而，这两个函数仍然产生略有不同的结果。我可能不得不深入研究编译的程序集以找出原因，但也许我忽略了一些简单的事情。

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

比较：

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

【讨论】：

好的，谢谢您的回答。实际上，在将括号放在正确的位置之后（tmp2 是分子，tmp3 是分母，所以我的 tmp1 缺少一对括号），对于函数 1 和函数 2 得到相同的确切结果。我看到有很多细节在 cython 文档中并不明显。 type(np.float32(1) / 2) == 有点小门，一个操作数是 np.float32 对象，另一个是 int 对象，结果两者都不是！！
所以，现在我可以匹配结果，但这意味着我必须牢记这一点来编写我的 python 和 cython，而不是只编写一个 python 版本然后获得一个等效的 cython 版本。我想在纯python模式下我也会有这个细微的差异问题？

【解决方案2】：

我们比较一下：

function1 一直保持 float32_t。
function2 在索引时转换为float，使用float 执行中间步骤，然后转换回float32_t 以获得最终结果。
function3 转换为 float，但随后立即返回到 float32_t，并以此执行中间步骤。
function4 转换为float，执行中间步骤，然后将最终结果返回为float。

至于为什么function4 打印与function2 相同的内容，但返回的内容不同：如果您查看类型，很简单。这些值显然足够接近，以至于它们以相同的方式发生在 print 上，但与 repr 以同样的方式不够接近。这并不奇怪，因为它们不是同一类型。

【讨论】：