以最快的方式处理大型数据集的 Hausdorff 距离
Hausdorff distance for large dataset in a fastest way
我的数据集中的行数是 500000+。我需要每个 id 与其他人之间的 Hausdorff 距离。并对整个数据集重复它
我有一个庞大的数据集。这是小部分：
df =
id_easy ordinal latitude longitude epoch day_of_week
0 aaa 1.0 22.0701 2.6685 01-01-11 07:45 Friday
1 aaa 2.0 22.0716 2.6695 01-01-11 07:45 Friday
2 aaa 3.0 22.0722 2.6696 01-01-11 07:46 Friday
3 bbb 1.0 22.1166 2.6898 01-01-11 07:58 Friday
4 bbb 2.0 22.1162 2.6951 01-01-11 07:59 Friday
5 ccc 1.0 22.1166 2.6898 01-01-11 07:58 Friday
6 ccc 2.0 22.1162 2.6951 01-01-11 07:59 Friday
我要计算Haudorff Distance:
import pandas as pd
import numpy as np
from scipy.spatial.distance import directed_hausdorff
from scipy.spatial.distance import pdist, squareform
u = np.array([(2.6685,22.0701),(2.6695,22.0716),(2.6696,22.0722)]) # coordinates of `id_easy` of taxi `aaa`
v = np.array([(2.6898,22.1166),(2.6951,22.1162)]) # coordinates of `id_easy` of taxi `bbb`
directed_hausdorff(u, v)[0]
输出为0.05114626086039758
现在我想为整个数据集计算这个距离。对于所有 id_easy。期望的输出是对角线上带有 0 的矩阵（因为 aaa 和 aaa 之间的距离是 0）：
aaa bbb ccc
aaa 0 0.05114 ...
bbb ... 0
ccc 0
尝试使用 scipy
中的 computed distance
from scipy.spatial.distance import cdist
hausdorff_distance = cdist(df[['latitude', 'longitude']], df[['latitude', 'longitude']], lambda u, v: directed_hausdorff(u, v)[0])
hausdorff_distance_df = pd.DataFrame(hausdorff_distance)
请注意，无论您最终使用哪种方法 - 计算都将花费大量时间，这只是由于数据量巨大。问问自己，你是否真的需要每一对距离。
实际上，这类问题是通过将配对数量限制在可管理的数量来解决的。例如，将您的数据框分割成更小的集合，每个集合仅限于一个地理区域，然后在该地理区域内找到一对距离。
超市使用上述方法来确定新店的位置。他们并没有计算他们拥有的每家商店与竞争对手拥有的每家商店之间的一对距离。首先他们限制了面积，总共只有5-10家店，然后他们才开始计算距离。
你说的是计算 500000^2+ 距离。如果每秒计算 1000 个这样的距离，则需要 7.93 年 才能完成矩阵。我不确定 Hausdorff distance 是否对称，但即使是，也只能为您节省两倍（3.96 年）。
该矩阵还将占用大约 1 TB 的内存。
我建议只在需要时才计算它，或者如果您确实需要整个矩阵，则需要并行计算。从好的方面来说，这个问题很容易解决。比如四核，可以这样拆分问题（伪代码）：
n = len(u)
m = len(v)
A = hausdorff_distance_matrix(u[:n], v[:m])
B = hausdorff_distance_matrix(u[:n], v[m:])
C = hausdorff_distance_matrix(u[n:], v[:m])
D = hausdorff_distance_matrix(u[n:], v[m:])
results = [[A, B],
[C, D]]
其中 hausdorff_distance_matrix(u, v) returns u 和 v 之间的所有距离组合。不过，您可能需要将它分成四个以上的部分。
应用是什么？您可以只根据需要分段计算这些吗？
首先我定义了一个提供一些示例数据的方法。如果您在问题中提供类似的内容，那将会容易得多。在大多数与性能相关的问题中，需要实际问题的大小才能找到最佳解决方案。
在下面的回答中，我将假设 id_easy 的平均大小为 17，并且有 30000 个不同的 ID，这导致数据集大小为 510_000。
创建示例数据
import numpy as np
import numba as nb
N_ids=30_000
av_id_size=17
#create_data (pre sorting according to id assumed)
lat_lon=np.random.rand(N_ids*av_id_size,2)
#create_ids (sorted array with ids)
ids=np.empty(N_ids*av_id_size,dtype=np.int64)
ind=0
for i in range(N_ids):
for j in range(av_id_size):
ids[i*av_id_size+j]=ind
ind+=1
豪斯多夫函数
以下函数是 scipy-source 的略微修改版本。
做了以下修改：
- 对于非常小的输入数组，我注释掉了改组部分（启用 对较大的数组进行改组并在您的真实数据上尝试最好的方法
- 至少在 Windows 上，Anaconda scipy 函数看起来有一些性能问题（比 Linux 慢得多），基于 LLVM 的 Numba 看起来是一致的
- 删除了 Hausdorff 对的索引
为 (N,2) 情况展开的距离循环
#Modified Code from Scipy-source
#https://github.com/scipy/scipy/blob/master/scipy/spatial/_hausdorff.pyx
#Copyright (C) Tyler Reddy, Richard Gowers, and Max Linke, 2016
#Copyright © 2001, 2002 Enthought, Inc.
#All rights reserved.
#Copyright © 2003-2013 SciPy Developers.
#All rights reserved.
#Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
#Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
#Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following
#disclaimer in the documentation and/or other materials provided with the distribution.
#Neither the name of Enthought nor the names of the SciPy Developers may be used to endorse or promote products derived
#from this software without specific prior written permission.
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
#BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
#IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
#OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
#OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@nb.njit()
def directed_hausdorff_nb(ar1, ar2):
N1 = ar1.shape[0]
N2 = ar2.shape[0]
data_dims = ar1.shape[1]
# Shuffling for very small arrays disbabled
# Enable it for larger arrays
#resort1 = np.arange(N1)
#resort2 = np.arange(N2)
#np.random.shuffle(resort1)
#np.random.shuffle(resort2)
#ar1 = ar1[resort1]
#ar2 = ar2[resort2]
cmax = 0
for i in range(N1):
no_break_occurred = True
cmin = np.inf
for j in range(N2):
# faster performance with square of distance
# avoid sqrt until very end
# Simplificaten (loop unrolling) for (n,2) arrays
d = (ar1[i, 0] - ar2[j, 0])**2+(ar1[i, 1] - ar2[j, 1])**2
if d < cmax: # break out of `for j` loop
no_break_occurred = False
break
if d < cmin: # always true on first iteration of for-j loop
cmin = d
# always true on first iteration of for-j loop, after that only
# if d >= cmax
if cmin != np.inf and cmin > cmax and no_break_occurred == True:
cmax = cmin
return np.sqrt(cmax)
计算子集上的 Hausdorff 距离
@nb.njit(parallel=True)
def get_distance_mat(def_slice,lat_lon):
Num_ids=def_slice.shape[0]-1
out=np.empty((Num_ids,Num_ids),dtype=np.float64)
for i in nb.prange(Num_ids):
ar1=lat_lon[def_slice[i:i+1],:]
for j in range(i,Num_ids):
ar2=lat_lon[def_slice[j:j+1],:]
dist=directed_hausdorff_nb(ar1, ar2)
out[i,j]=dist
out[j,i]=dist
return out
示例和时间
#def_slice defines the start and end of the slices
_,def_slice=np.unique(ids,return_index=True)
def_slice=np.append(def_slice,ids.shape[0])
%timeit res_1=get_distance_mat(def_slice,lat_lon)
#1min 2s ± 301 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
这取决于您的确切用例：
- 如果您需要预先计算值，因为您需要快速获得大致答案，则可以 'bin' 您的答案。在最粗略的情况下，假设您有 1.2343 * 10^7 和 9.2342 * 10^18。您显着降低分辨率并直接查找。也就是说，您在预先计算的查找 table 中查找“1E7 和 1E19”。您得到的正确答案精确到 10 的幂。
- 如果您想表征您的全球答案分布，请选择统计上足够数量的随机对并改用它。
- 如果您想绘制答案图表，运行 连续传递更精确的答案。例如，从 50x50 点等间距制作图表。计算完成后，如果用户还在，则开始填充 250x250 分辨率，依此类推。有一些可爱的技巧可以决定在哪里集中更多的计算能力，只需跟随用户的鼠标即可。
- 如果您可能想要某些部分的答案，但直到稍后的计算才知道是哪一个，那么让评估变得惰性。也就是说，尝试
__getitem__() 实际计算出的值。
- 如果您认为自己现在确实需要所有答案，请多想想。这么多活动部件没有问题。你应该在已经存在的数据中找到一些 exploitable 约束，否则数据已经不堪重负了。
继续思考！继续黑客攻击！做笔记。
我的数据集中的行数是 500000+。我需要每个 id 与其他人之间的 Hausdorff 距离。并对整个数据集重复它
我有一个庞大的数据集。这是小部分：
df =
id_easy ordinal latitude longitude epoch day_of_week
0 aaa 1.0 22.0701 2.6685 01-01-11 07:45 Friday
1 aaa 2.0 22.0716 2.6695 01-01-11 07:45 Friday
2 aaa 3.0 22.0722 2.6696 01-01-11 07:46 Friday
3 bbb 1.0 22.1166 2.6898 01-01-11 07:58 Friday
4 bbb 2.0 22.1162 2.6951 01-01-11 07:59 Friday
5 ccc 1.0 22.1166 2.6898 01-01-11 07:58 Friday
6 ccc 2.0 22.1162 2.6951 01-01-11 07:59 Friday
我要计算Haudorff Distance:
import pandas as pd
import numpy as np
from scipy.spatial.distance import directed_hausdorff
from scipy.spatial.distance import pdist, squareform
u = np.array([(2.6685,22.0701),(2.6695,22.0716),(2.6696,22.0722)]) # coordinates of `id_easy` of taxi `aaa`
v = np.array([(2.6898,22.1166),(2.6951,22.1162)]) # coordinates of `id_easy` of taxi `bbb`
directed_hausdorff(u, v)[0]
输出为0.05114626086039758
现在我想为整个数据集计算这个距离。对于所有 id_easy。期望的输出是对角线上带有 0 的矩阵（因为 aaa 和 aaa 之间的距离是 0）：
aaa bbb ccc
aaa 0 0.05114 ...
bbb ... 0
ccc 0
尝试使用 scipy
中的 computed distancefrom scipy.spatial.distance import cdist
hausdorff_distance = cdist(df[['latitude', 'longitude']], df[['latitude', 'longitude']], lambda u, v: directed_hausdorff(u, v)[0])
hausdorff_distance_df = pd.DataFrame(hausdorff_distance)
请注意，无论您最终使用哪种方法 - 计算都将花费大量时间，这只是由于数据量巨大。问问自己，你是否真的需要每一对距离。
实际上，这类问题是通过将配对数量限制在可管理的数量来解决的。例如，将您的数据框分割成更小的集合，每个集合仅限于一个地理区域，然后在该地理区域内找到一对距离。
超市使用上述方法来确定新店的位置。他们并没有计算他们拥有的每家商店与竞争对手拥有的每家商店之间的一对距离。首先他们限制了面积，总共只有5-10家店，然后他们才开始计算距离。
你说的是计算 500000^2+ 距离。如果每秒计算 1000 个这样的距离，则需要 7.93 年 才能完成矩阵。我不确定 Hausdorff distance 是否对称，但即使是，也只能为您节省两倍（3.96 年）。
该矩阵还将占用大约 1 TB 的内存。
我建议只在需要时才计算它，或者如果您确实需要整个矩阵，则需要并行计算。从好的方面来说，这个问题很容易解决。比如四核，可以这样拆分问题（伪代码）：
n = len(u)
m = len(v)
A = hausdorff_distance_matrix(u[:n], v[:m])
B = hausdorff_distance_matrix(u[:n], v[m:])
C = hausdorff_distance_matrix(u[n:], v[:m])
D = hausdorff_distance_matrix(u[n:], v[m:])
results = [[A, B],
[C, D]]
其中 hausdorff_distance_matrix(u, v) returns u 和 v 之间的所有距离组合。不过，您可能需要将它分成四个以上的部分。
应用是什么？您可以只根据需要分段计算这些吗？
首先我定义了一个提供一些示例数据的方法。如果您在问题中提供类似的内容，那将会容易得多。在大多数与性能相关的问题中，需要实际问题的大小才能找到最佳解决方案。
在下面的回答中，我将假设 id_easy 的平均大小为 17，并且有 30000 个不同的 ID，这导致数据集大小为 510_000。
创建示例数据
import numpy as np
import numba as nb
N_ids=30_000
av_id_size=17
#create_data (pre sorting according to id assumed)
lat_lon=np.random.rand(N_ids*av_id_size,2)
#create_ids (sorted array with ids)
ids=np.empty(N_ids*av_id_size,dtype=np.int64)
ind=0
for i in range(N_ids):
for j in range(av_id_size):
ids[i*av_id_size+j]=ind
ind+=1
豪斯多夫函数
以下函数是 scipy-source 的略微修改版本。 做了以下修改：
- 对于非常小的输入数组，我注释掉了改组部分（启用 对较大的数组进行改组并在您的真实数据上尝试最好的方法
- 至少在 Windows 上，Anaconda scipy 函数看起来有一些性能问题（比 Linux 慢得多），基于 LLVM 的 Numba 看起来是一致的
- 删除了 Hausdorff 对的索引
为 (N,2) 情况展开的距离循环
#Modified Code from Scipy-source #https://github.com/scipy/scipy/blob/master/scipy/spatial/_hausdorff.pyx #Copyright (C) Tyler Reddy, Richard Gowers, and Max Linke, 2016 #Copyright © 2001, 2002 Enthought, Inc. #All rights reserved. #Copyright © 2003-2013 SciPy Developers. #All rights reserved. #Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: #Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. #Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following #disclaimer in the documentation and/or other materials provided with the distribution. #Neither the name of Enthought nor the names of the SciPy Developers may be used to endorse or promote products derived #from this software without specific prior written permission. #THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, #BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. #IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, #OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; #OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT #(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @nb.njit() def directed_hausdorff_nb(ar1, ar2): N1 = ar1.shape[0] N2 = ar2.shape[0] data_dims = ar1.shape[1] # Shuffling for very small arrays disbabled # Enable it for larger arrays #resort1 = np.arange(N1) #resort2 = np.arange(N2) #np.random.shuffle(resort1) #np.random.shuffle(resort2) #ar1 = ar1[resort1] #ar2 = ar2[resort2] cmax = 0 for i in range(N1): no_break_occurred = True cmin = np.inf for j in range(N2): # faster performance with square of distance # avoid sqrt until very end # Simplificaten (loop unrolling) for (n,2) arrays d = (ar1[i, 0] - ar2[j, 0])**2+(ar1[i, 1] - ar2[j, 1])**2 if d < cmax: # break out of `for j` loop no_break_occurred = False break if d < cmin: # always true on first iteration of for-j loop cmin = d # always true on first iteration of for-j loop, after that only # if d >= cmax if cmin != np.inf and cmin > cmax and no_break_occurred == True: cmax = cmin return np.sqrt(cmax)
计算子集上的 Hausdorff 距离
@nb.njit(parallel=True)
def get_distance_mat(def_slice,lat_lon):
Num_ids=def_slice.shape[0]-1
out=np.empty((Num_ids,Num_ids),dtype=np.float64)
for i in nb.prange(Num_ids):
ar1=lat_lon[def_slice[i:i+1],:]
for j in range(i,Num_ids):
ar2=lat_lon[def_slice[j:j+1],:]
dist=directed_hausdorff_nb(ar1, ar2)
out[i,j]=dist
out[j,i]=dist
return out
示例和时间
#def_slice defines the start and end of the slices
_,def_slice=np.unique(ids,return_index=True)
def_slice=np.append(def_slice,ids.shape[0])
%timeit res_1=get_distance_mat(def_slice,lat_lon)
#1min 2s ± 301 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
这取决于您的确切用例：
- 如果您需要预先计算值，因为您需要快速获得大致答案，则可以 'bin' 您的答案。在最粗略的情况下，假设您有 1.2343 * 10^7 和 9.2342 * 10^18。您显着降低分辨率并直接查找。也就是说，您在预先计算的查找 table 中查找“1E7 和 1E19”。您得到的正确答案精确到 10 的幂。
- 如果您想表征您的全球答案分布，请选择统计上足够数量的随机对并改用它。
- 如果您想绘制答案图表，运行 连续传递更精确的答案。例如，从 50x50 点等间距制作图表。计算完成后，如果用户还在，则开始填充 250x250 分辨率，依此类推。有一些可爱的技巧可以决定在哪里集中更多的计算能力，只需跟随用户的鼠标即可。
- 如果您可能想要某些部分的答案，但直到稍后的计算才知道是哪一个，那么让评估变得惰性。也就是说，尝试
__getitem__()实际计算出的值。 - 如果您认为自己现在确实需要所有答案，请多想想。这么多活动部件没有问题。你应该在已经存在的数据中找到一些 exploitable 约束，否则数据已经不堪重负了。
继续思考！继续黑客攻击！做笔记。