python中浮点精度的最佳实践
Mag*_*eed 3 python floating-point precision double
今晚早些时候，我的一个朋友把这个可爱的问题递给了我。问题说：
在MATLAB中编写程序以检查点是否在三角形内。不要忘记检查该点是否也在边界上。三角点
x=(0,0)，y=(0,1)和z=(1,0)
这个问题并不难解决。想法是找到斜边的方程，并检查该点是否位于三角形的任何边上。但是，检查内部和外部并不那么困难。
我在MATLAB上编写了代码，逻辑似乎很好。但是问题是结果与该逻辑不一致！我开始对我的代码提出质疑，因为我对MATLAB的了解并不熟练。尽管如此，我还是尝试了首选语言Python。
这是我的代码：
def isInsideTriangle(x,y):
if x == 0 or y == 0 or y == 1-x:
print('on the border of the triangle')
elif x > 1 or y > 1 or x < 0 or y < 0 or y > 1-x:
print('outside of the triangle')
print(1-x) # check the value
else:
# verbose these values to double check
print(1-x)
print(y)
print(type(y))
print(type(1-x))
print(y==(1-x))
print('inside of the triangle')
isInsideTriangle(0.2,0.8)
尝试使用这两个值时，控制台上的结果应为on the border。但是，程序说的是inside！我试图之间切换x，并y即isInsideTriangle(0.8,0.2)但是程序输出的预期结果这一次。
这使我认识到，逻辑与浮点精度无关。我在MATLAB上将变量的大小增加到64位精度，并且程序运行正常。
我现在的问题
作为Python专家，如何避免Python中的此类问题的最佳编程实践是什么？我们如何才能特别避免在生产环境中出现此类烦人的问题？
首先，您的逻辑不正确。考虑的情况下x=0.9，y=0.9。显然这是在三角形之外，但不满足任何条件x > 1 or y > 1 or x < 0 or y < 0。
第二，任何涉及相等性比较的浮点算法（例如测试点是否在形状的边界上）都可能会受到精度问题的影响。修改逻辑以测试点是否在边界的小范围内可能会更好。
我建议不要将Decimal类用于本机不是十进制数字的任何内容，例如货币。对小数执行除基本算术之外的任何操作（例如math.sqrt），无论如何都会内部将其转换为浮点数。
前言
您的问题是“ 点是否在规则的多边形内？ ” 的特殊化，因为规则的意思是不像我们在GIS系统中经常看到的那样自相交或多多边形。顺便说一下，因为您要的是三角形，所以我假设多边形是凸的。
有趣的是，交叉产品是解决它的关键。在2D平面中处理矢量时，叉积与该平面正交。要提取的有用信息是：它指向上还是向下？
浮点算术错误将发生，并且当叉积接近零但不等于零时，它将变得很关键，然后它将具有一个符号而不是为null。
要检查您的点是否在多边形内，它仅归结为检查边缘和点之间的所有叉积是否具有相同的符号，例如：sign(h x w) = -1。
同样，检查多边形是否为凸形将要检查连续边的所有叉积具有相同的符号，例如：sign(u x v) = -1。
实作
让我们建立一个小类来检查点是否在规则凸多边形的内部（在边缘或外部）：
import numpy as np
class cpoly:
def __init__(self, points=[[0,0], [0,1], [1,0]], assert_convexity=True):
"""
Initialize 2D Polygon with a sequence of 2D points
"""
self._points = np.array(points)
assert self.p.shape[0] >= 3
assert self.p.shape[1] == 2
assert self.is_convex or not(assert_convexity)
@property
def n(self):
return self.p.shape[0]
@property
def p(self):
return self._points
@property
def is_convex(self):
"""
Check convexity of the polygon (operational for a non intersecting polygon)
"""
return self.contains()
def contains(self, p=None, debug=False, atol=2e-16):
"""
Check if a 2D convex polygon contains a point (also used to assess convexity)
Returns:
-1: Point is oustide the polygon
0: Point is close to polygon edge (epsilon ball)
+1: Point is inside the polygon
"""
s = None
c = False
n = self.n
for k in range(n):
# Vector Differences:
d1 = self.p[(k+1)%n,:] - self.p[k%n,:]
if p:
d2 = p - self.p[k%n,:]
else:
d2 = self.p[(k+2)%n,:] - self.p[(k+1)%n,:]
# Cross Product:
z = np.cross(d1, d2)
if np.allclose(z, 0, atol=atol):
s_ = 0
c = True
else:
s_ = np.sign(z)
# Debug Helper:
if debug:
print("k = %d, d1 = %s, d2 = %s, z = %.32f, s = %d" % (k, d1, d2, z, s_))
# Check if cross product sign change (excluded null, when point is colinear with the segment)
if s and (s_ != s) and not(s_ == 0):
# Nota: Integer are exact if float representable, therefore comparizons are correct
return -1
s = s_
if c:
return 0
else:
return 1
def plot(self, axe=None):
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
if not(axe):
fig, axe = plt.subplots()
axe.plot(self.p[:,0], self.p[:,1], 'x', markersize=10, label='Points $p_i$')
axe.add_patch(Polygon(self.p, alpha=0.4, label='Area'))
axe.set_xlabel("$x$")
axe.set_ylabel("$y$")
axe.set_title("Polygon")
axe.set_aspect('equal')
axe.legend(bbox_to_anchor=(1,1), loc='upper left')
axe.grid()
return axe.get_figure(), axe
该类使用2D点列表进行初始化（默认情况下是您的）。
p = cpoly()
在我的设置中，浮点精度约为：
e = np.finfo(np.double).eps # 2.220446049250313e-16
我们创建用于测试目的的试验数据集：
p = cpoly()
r = [
[0,0], [0,1], [1,0], # Polygon vertices
[0,0.5], [-e,0.6], [e,0.4], [0.1, 0.1], [1,1],
[0.5+e,0.5], [0.3-e,0.7], [0.7+e/10,0.3],
[0, 1.2], [1.2, 0.], # Those points make your logic fails
[0.2,0.8], [0.1,0.9],
[0.8+10*e,0.2], [0.9+10*e,0.1]
]
如果我们调整您的功能以具有兼容的输出：
def isInsideTriangle(x,y):
if x == 0 or y == 0 or y == 1-x:
return 0
elif x > 1 or y > 1 or x < 0 or y < 0 or y > 1-x:
return -1
else:
return 1
然后，我们检查试验点，以查看我们两个函数是否运行良好：
_, axe = p.plot()
cols = {-1: 'red', 0:'orange', 1:'green'}
for x in r:
q1 = p.contains(x)
q2 = isInsideTriangle(*x)
print(q1==q2)
axe.plot(*x, 'o', markersize=4, color=cols[q1])
使用此设置，所有点均已正确分类。但是您可以看到您的算法存在缺陷。主要是以下几行：
if x == 0 or y == 0 or y == 1-x:
拒绝[0, 1.2]和失败[1.2, 0.]。
靠近边缘的点和浮动算术误差
即使对于恰好在边缘的[0.2,0.8]点（如浮点错误）也会导致点错误分类。以下几点将使叉积不完全等于零。在下面查看详细信息：
p.contains([0.2,0.8], debug=True) # True
# k = 0, d1 = [0 1], d2 = [0.2 0.8], z = -0.20000000000000001110223024625157, s = -1
# k = 1, d1 = [ 1 -1], d2 = [ 0.2 -0.2], z = 0.00000000000000005551115123125783, s = 0
# k = 2, d1 = [-1 0], d2 = [-0.8 0.8], z = -0.80000000000000004440892098500626, s = -1
这就是为什么我们必须添加一个半径球atol来检查在给定公差下实际上是否为零：
if np.allclose(z, 0, atol=atol):
s_ = 0
# ...
else:
s_ = np.sign(z)
这意味着我们必须接受将足够靠近边缘的点（从两侧）视为包含在多边形中。这是Float算术固有的，最好的办法是atol为您的应用程序调整到可接受的值。或者，您可以找到另一个不受此问题困扰的逻辑或数据模型。
缩放到精确比例
如果我们缩放到精确比例以查看边缘附近发生了什么，我们得到：
n = 40
lim = np.array([0.5,0.5]) + [-n*e/2, +n*e/2]
x = np.linspace(lim[0], lim[1], 30)
X, Y = np.meshgrid(x, x)
x, y = X.reshape(-1), Y.reshape(-1)
_, axe = p.plot()
axe.set_xlim(lim)
axe.set_ylim(lim)
for r in zip(x,y):
q = p.contains(r)
axe.plot(*r, 'o', color=cols[q], markersize=2)
我们看到一些非常靠近边缘但在多边形内部或外部的点被分类为“在边缘上”。这是由于ε球标准。您还可以观察到点之间的间距不相等（无论我是否使用过linspace），因为不可能表示10为的整数次幂2。
结论
上面的解决方案是您的问题的概括，在中执行O(n)。它似乎过大了，但是它是通用的（适用于任何常规多边形）并且是综合的（它依赖于众所周知的几何概念）。
实际上，该算法只是在通过路径时检查该点是否位于多边形所有边的同一侧。如果是这样，那么它将得出该点在多边形内部的结论！
上面的解决方案由于依赖于浮点运算（请参见point 10）而不受浮点运算错误的影响。幸运的是，使用epsilon球测试，我们可以减轻它。
参考
如果您想更深入地了解有限精度算术，我建议您阅读一本出色的书：《数值算法的准确性和稳定性》，J。Higham。
奖金
以下是所有答案与试验数据集的比较：
我们可以对此软检查所强调的不同类型的“错误”给出一些背景信息：
- 点
11和由于设计缺陷而12被错误分类@MagedSaeed； @MahmoudElshahat在它的逻辑中使用直浮子平等，这就是为什么它返回等价类的点（例如点：不同的结果0到2，多边形分）;@SamMason使用最简单的逻辑和epsilon球测试。它似乎具有最大的错误率，因为它无法区分内部和边界（我们可以忽略所有确切答案是0且函数返回的点1）。IMO，这是在if-then中运行的最佳算法O(1)。另外，它可以很容易地进行更新以考虑3种状态的逻辑。- Point
10是故意设计来挑战边界的逻辑的，该点显然距离多边形较远，且与机器精度有关。这是浮点算术错误变得显着且逻辑更加模糊的地方：可接受且易于处理的距离边缘有多远？