我有两组坐标,pos1.txt
和pos2.txt
(附件)。点的ID对应于行号。我想找到一组索引应用于第一组,以尽可能接近第二组,同时保持第一个位置。scipy.optimize
似乎在这方面做得很好,下面是我目前使用的代码。
然而,我想“放松”优化问题,因为它对我的目的来说太慢了(理想情况下,我会使用500,000个点),而且我真的不知道如何进行:
1.到目前为止,最重要的匹配部分是z和r。我并不真的关心[x, y]
中的位置是否有很大的不同,只要径向距离和z相同。
z/zmax>90%
的位置是完全可以互换的。如果这个区域中的一个点改变了它的径向位置很多或者它的z(同时保持在〉90%),这是可以的。
我附上了ID和误差的当前图,我们可以看到,在整个几何结构中,径向/轴向距离的差异最多可达7%,这高于我的目标值,但速度相当慢。
有没有人能告诉我该怎么做,以及如何减少时间?
下面是代码:
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(20,10))
# First set of points with random order
list1 = np.genfromtxt('pos1.txt', delimiter=',')
ax = fig.add_subplot(231, projection="3d")
s = ax.scatter3D(list1[:,0], list1[:,1], list1[:,2], c=np.arange(list1.shape[0]), cmap='jet', s=5)
plt.colorbar(s)
plt.title('ID1')
# Second set of points, sorted and different positions
list2 = np.genfromtxt('pos2.txt', delimiter=',')
ax = fig.add_subplot(232, projection="3d")
s = ax.scatter3D(list2[:,0], list2[:,1], list2[:,2], c=np.arange(list2.shape[0]), cmap='jet', s=5)
plt.colorbar(s)
plt.title('ID2')
# Calculate third set of points
# With scipy.optimize which indices to apply on first set to have lowest distance to second set
cost = np.linalg.norm(list2[:, np.newaxis, :] - list1, axis=2)
_, indices = scipy.optimize.linear_sum_assignment(cost)
list3 = list1[indices]
ax = fig.add_subplot(233, projection="3d")
list2 = np.genfromtxt('pos2.txt', delimiter=',')
s = ax.scatter3D(list3[:,0], list3[:,1], list3[:,2], c=np.arange(list3.shape[0]), cmap='jet', s=5)
plt.colorbar(s)
plt.title('ID3')
# Calculate error made on distance, radial distance, axial distance
error_dist = np.linalg.norm(list2-list3, axis=1)/np.max(cost)*100
ax = fig.add_subplot(234, projection="3d")
list2 = np.genfromtxt('pos2.txt', delimiter=',')
s = ax.scatter3D(list3[:,0], list3[:,1], list3[:,2], c=error_dist, cmap='jet', s=5)
plt.colorbar(s)
plt.title('Error on distance (%)')
r_dist2 = np.linalg.norm(list2[:, :2], axis=1)
r_dist3 = np.linalg.norm(list3[:, :2], axis=1)
error_r_dist = (r_dist3-r_dist2)/np.max(r_dist2)*100
ax = fig.add_subplot(235, projection="3d")
list2 = np.genfromtxt('pos2.txt', delimiter=',')
s = ax.scatter3D(list3[:,0], list3[:,1], list3[:,2], c=error_r_dist, cmap='jet', s=5)
plt.colorbar(s)
plt.title('Error on radial distance (%)')
error_z_dist = (list3[:,2]-list2[:,2])/(np.max(list3[:,2])-np.min(list3[:,2]))*100
ax = fig.add_subplot(236, projection="3d")
list2 = np.genfromtxt('pos2.txt', delimiter=',')
s = ax.scatter3D(list3[:,0], list3[:,1], list3[:,2], c=error_z_dist, cmap='jet', s=5)
plt.colorbar(s)
plt.title('Error on axial distance (%)')
plt.tight_layout()
plt.savefig('match_pos.png', dpi=600)
plt.show()
我在这个网站上上传了职位文件:
- 示例https://zerobin.net/?1be6d4810e148e20#qflDwnHeACHuWzZKrtFdTxukdsvNSrWK+2cJ9l3XWLQ=
- 示例https://zerobin.net/?ecc23876cd7190e7#9e2AKTWp/mSXjlqEZavwEeslQ5oipN5gM+v3Xp4U1ws=
谢谢你!
1条答案
按热度按时间sdnqo3pr1#
我偶然发现了this post,它详细介绍了如何使用ase.geometry来“更快”地计算距离。然而,下面的代码用了3分钟14秒才完成,而用我的~ 15,000个点来完成numpy.linalg.norm方法用了47秒。我的实现有什么问题吗?