如何使用scipy.weave更改pythons代码？(如何更快地执行代码？)

Question

在Python中:我有2个数组:lon,lat(超过500 000个点).

     lon = numpy.array()
     lat = numpy.array()

     flag = True
     while flag:
        lon1 = lon[:-1]
        lon2 = lon[1:]
        lat1 = lat[:-1]
        lat2 = lat[1:]

        '''distance'''
        x = (lon2 - lon1)
        y = (lat2 - lat1)
        d = numpy.sqrt(x * x + y * y)

        min = numpy.min(d)
        if min < 0.015:
            j = numpy.where(d == min)[0][0]

            lon[j] = (lon[j] + lon[j + 1]) / 2
            lat[j] = (lat[j] + lat[j + 1]) / 2

            lon = numpy.delete(lon, j + 1)
            lat = numpy.delete(lat, j + 1)

        else:
            flag = False

这段代码工作得很慢!

1. 请提示,如何更快地执行代码？

2. 我知道有scipy.weave.请提示,如何使用scipy.weave更改代码？

   import scipy
   from scipy import weave
   
   codeC = """
   ???
   """
   scipy.weave.inline(codeC, ['lon', 'lat'], compiler='gcc')
   

谢谢.

更新:

    lon1 = lon[:-1]
    lon2 = lon[1:]
    lat1 = lat[:-1]
    lat2 = lat[1:]

    x = (lon2 - lon1)
    y = (lat2 - lat1)
    d = x * x + y * y

    while True and d.size > 1:

        j = np.argmin(d)
        if d[j] > min_radius:
            break

        lon[j] = (lon[j] + lon[j + 1]) / 2
        lat[j] = (lat[j] + lat[j + 1]) / 2

        '''lon = np.delete(lon, j + 1)
           lat = np.delete(lat, j + 1)
           d = np.delete(d, j)'''
        if j == (d.size - 1):
            lon = lon[:j + 1]
            lat = lat[:j + 1]
            d = d[:j]
        else:
            lon = np.hstack((lon[:j + 1], lon[j + 2:]))
            lat = np.hstack((lat[:j + 1], lat[j + 2:]))
            d = np.hstack((d[:j], d[j + 1:]))


        if j > 0:
            x = lon[j] - lon[j - 1]
            y = lat[j] - lat[j - 1]
            d[j - 1] = x * x + y * y
        if j < (d.size - 1):
            x = lon[j + 1] - lon[j]
            y = lat[j + 1] - lat[j]
            d[j] = x * x + y * y

@ilmiacs,@ Jaime,@ Luke,谢谢!它工作得更快!

但总的来说,代码在500点的数组中运行很长时间...... :(((

Answer 1

您正在重复循环中每次迭代的整个昂贵的计算.一次迭代只删除一个节点,然后在删除一个元素时复制整个数组并重做距离计算.然后你搜索下一个节点.

因此,您的代码至少运行时不是最理想的O(N^2).但是,对于你想要达到的目标, O(N log N)或者可能O(N)就足够了.

如果重新排列算法,一个没有numpy的简单python模块应该足够快.食谱:

1. 避免重新计算整个距离数组.删除一个节点时,只需重新计算两个值.
2. 避免复制.它也很贵.

希望这可以帮助.

编辑:

不重新计算距离数组的示例实现将是:

lon = numpy.array()
lat = numpy.array()

lon1 = lon[:-1]
lon2 = lon[1:]
lat1 = lat[:-1]
lat2 = lat[1:]

'''distance'''
x = (lon2 - lon1)
y = (lat2 - lat1)
d = x * x + y * y

while True:
    min = numpy.min(d)
    if min < 0.015 * 0.015:
        j = numpy.where(d == min)[0][0]
    else:
        break

    lon[j] = (lon[j] + lon[j + 1]) / 2
    lat[j] = (lat[j] + lat[j + 1]) / 2

    lon = numpy.delete(lon, j + 1)  # <-- you spend most of the time here
    lat = numpy.delete(lat, j + 1)  # <-- and here

    d = numpy.delete(d, j)      # <-- and here

    x = lon[j]-lon[j-1]
    y = lat[j]-lat[j-1]
    d[j-1] = x * x + y * y
    x = lon[j+1]-lon[j]
    y = lat[j+1]-lat[j]
    d[j] = x * x + y * y

虽然这会给你一个提升,但你仍然是O(N^2),因为广泛的复制.为避免这种情况,您需要将数组存储在双向链表中.然后删除一个元素O(1),而不是O(N).使用链表,您的最差操作将是找到最小值.你可以通过引入一个排序版本来改进d,你也可以跟踪其中的位置d.您需要跟踪更改.但最后你会得到一个算法O(N log N).

编辑2:

正如@Jaime所指出的那样,比较平方距离也是不必要的.我已经纠正了上面的代码.

编辑3:

在周末,我正在考虑更多关于算法的问题,结果证明这是一个非常有趣的问题,也很有趣.我想分享一下我的想法.

让它工作O(N log N)并不像我想象的那么简单.我建议保留一份清单sorted_d并进行管理.每次删除链接时,我们都必须更新两个条目sorted_d.这意味着,我们必须找到两个相应距离的新位置sorted_d.在正常列表中,找到位置可以通过二分法来完成O(log N).但是,据我所知,二分法在链表中是不可行的O(N).这就是问题.如果有人能证明我错了并且想出一种方法来在一个有序的链表中插入一个元素,我会很高兴O(log N).我个人没有看到它.

但是,有一个替代链表sorted_d.所有人都需要一个B树或Patricia树,它与浮动作为键.在这样的树中插入可以在(非常快)'O(log N)`中完成.但是,我不知道这种树的任何标准实现.标准实现适用于字符串或整数.因此,如果我们有一个从浮点数到整数的有序和无碰撞映射,我们就可以实现它.这种映射通常不存在,但毫无疑问可以为任何特定数据集构建.并且有问题的数据集是特别的:确实存在上限和精度.

希望这可以帮助.如果有进一步的兴趣,我也可以提出我对这个问题的最佳数据结构的看法.