尝试使用 Numpy 检查圆内的点

Question

我是 numpy 的新手，我的任务是处理以下情况：我需要创建两个随机整数（0 + 1 之间）的 numpy 数组。一个 numpy 数组表示 x 坐标，另一个表示 y 坐标。然后我需要使用 squareroot(x^2 + y^2)

我目前只是尝试对我的数组进行平方并将它们加在一起。可能是一个非常简单的任务是给我无穷无尽的麻烦。

import matplotlib.pyplot as plt
import numpy as np

plots = 100

dataOne = np.random.random(size = plots)
dataTwo = np.random.random(size = plots)

circle = plt.Circle((0,0), 1, alpha = 0.1)
plt.gca().add_patch(circle)
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.show()

squareDataOne = dataOne ** 2
squareDataTwo = dataTwo ** 2

if squareDataOne + squareDataTwo < 1:
    print("Now check square-root!")

我不断收到一条错误消息：包含多个元素的数组的真值不明确。使用 a.any() 或 a.all()。谁能解释为什么 Python/Numpy 不喜欢这个？我被告知尝试使用布尔表达式对数组进行切片。任何人都可以提供有关将其合并到我给定代码中的最佳方法的建议吗？ 对新手的任何建议或提示都表示赞赏。

Answer 1

squareDataOne 看起来像这样：

[7.43871942e-02 2.73007883e-01 5.23115388e-03 6.57541340e-01
 3.08779564e-01 1.24098667e-02 5.08258990e-01 6.52590269e-01
 8.90656103e-02 3.76389212e-02 2.12513661e-01 2.79683875e-01
 7.76233370e-01 6.48353342e-02 8.01663208e-01 8.69331480e-01
 4.34903542e]

squareData2 看起来很相似。 if 语句中的表达式：

squareDataOne + squareDataTwo < 1

产生这个数组：

[ True False  True  True  True  True  True  True False False  True False
  True False False  True  True  True  True  True False  True False False
  True  True False  True  True  True  True  True  True  True  True  True
  True  True False False  True False]

因此，您的 if 语句需要 True 或 False 值，并且正在获取此数组。错误消息告诉您 Python 不知道如何将此数组转换为单个 True 或 False 值。

我不太了解您的代码逻辑，不知道您需要做什么来解决此问题。显然你有很多数据，但你希望决定一个二元事件；如果你应该打印“现在检查平方根！”或不。我不知道你应该怎么做。

Answer 2

然后我需要检查这些点是否落在一个圆圈内 半径一使用 squareroot(x^2 + y^2)

你可以使用array filtering

pt_norm = (squareDataOne + squareDataTwo)
r_inside_circle = np.sqrt(pt_norm[pt_norm < 1])

这将为您提供落在r_inside_circle 中的圆内的所有点的半径。随着您增加 plots 的值，您会看到 (4.0*len(r_inside_circle))/len(dataOne) 将接近 PI。

Answer 3

我想你想替换这个：

if squareDataOne + squareDataTwo < 1:
    print("Now check square-root!")

类似：

# Calculate radii
radii = (squareDataOne + squareDataTwo)**0.5

# Create a boolean array (True/False values) for filtering
incircle = radii < 1

# Print the number of points inside the circle (Python 3.6+ for the f-string)
print(f"There are {np.sum(incircle)} points inside the circle")

# Plot the points within the circle
plt.scatter(dataOne[incircle], dataTwo[incircle])

dataOne[incircle] 和dataTwo[incircle] 将仅从每个数组中提取incircle 的值为True 的元素（即[x, y] 坐标对）。

Answer 4

# Change number of random points to 5.
plots = 5
...
inside_circle = squareDataOne + squareDataTwo < 1
print(inside_circle)

output:
[False  True  True False False]

它表明inside_circle是一个dtype为“bool”的numpy数组，表示每个点是否在圆内，而不是单个“bool”变量，这就是pythonif语句抛出错误的原因。您可以将inside_circle中的所有“bool”求和以获得圆内的点数（当对bool数组求和时，numpy将True视为1，False视为0）。

plots = 100
...
num_inside = np.sum(inside_circle)
# Check if all points are inside.
if num_inside == plots:
    print('All inside.')

# If what you really want to achieve is to estimate the area of one quarter of a circle
# using Monte Carlo Sampling.
print(num_inside / plots)    # Your estimate.
print(np.pi / 4)             # Actual area.

output:
0.79
0.7853981633974483

估计很接近，对吧？

Answer 5

以下是使用半径**2 和使用边界框进行预过滤的示例。这些优化对大型集合很有帮助。

使用边界框和 Radius**2 进行预过滤的速度是原来的 10 倍。此代码在 0.128 秒内执行：

import time
import numpy as np

plots = 10000000
radius = 1
center = np.array([0,0])

all_points = np.random.randint(low = -5, high = 5, size = plots*2).reshape(-1, 2)
print('All Points: ' + str(len(all_points)))

start_time = time.time()

# Limit expensive circle search to only points within box bounding circle
bounding_box_lowerleft = center - np.array([radius, radius])
bounding_box_upperright = center + np.array([radius, radius])
in_box_points = all_points[(all_points[:, 0] < bounding_box_upperright[0]) & (all_points[:,0] > bounding_box_lowerleft[0]) & (all_points[:, 1] < bounding_box_upperright[1]) & (all_points[:,1] > bounding_box_lowerleft[1])]
    
print('In Box: ' + str(len(in_box_points)))

# Use squared to avoid cost of sqrt
radius_squared = radius**2
in_circle_points = in_box_points[((in_box_points[:, 0] - center[0])**2 + (in_box_points[:, 1] - center[1])**2) <= radius_squared ]

print('In Circle: ' + str(len(in_circle_points)))
print('Elapsed time: ' + str(time.time() - start_time) + ' seconds')

Radius**2 快 4 倍。此代码在 0.349 秒内执行：

import time
import numpy as np

plots = 10000000
radius = 1
center = np.array([0,0])

all_points = np.random.randint(low = -5, high = 5, size = plots*2).reshape(-1, 2)
print('All Points: ' + str(len(all_points)))

start_time = time.time()

# Use squared to avoid cost of sqrt
radius_squared = radius**2
in_circle_points = all_points[((all_points[:, 0] - center[0])**2 + (all_points[:, 1] - center[1])**2) <= radius_squared ]

print('In Circle: ' + str(len(in_circle_points)))
print('Elapsed time: ' + str(time.time() - start_time) + ' seconds')

没有任何优化。此代码在 1.185 秒内执行：

import time
import numpy as np

plots = 10000000
radius = 1
center = np.array([0,0])

all_points = np.random.randint(low = -5, high = 5, size = plots*2).reshape(-1, 2)
print('All Points: ' + str(len(all_points)))

start_time = time.time()

in_circle_points = all_points[((all_points[:, 0] - center[0])**2 + (all_points[:, 1] - center[1])**2)**0.5 <= radius ]

print('In Circle: ' + str(len(in_circle_points)))
print('Elapsed time: ' + str(time.time() - start_time) + ' seconds')