加速python代码我可以向量化double for循环吗？

Question

我是python新手。我正在使用dbscan代码对一些改变。现在代码运行良好，但速度非常慢。所以我发现我必须从我的代码。这里是代码的一部分：

class Point:
    def __init__(self, x = 0, y = 0, visited = False, isnoise = False):
        self.x = x  
        self.y = y  
        self.visited = False  
        self.isnoise = False  

    def show(self):  
        return self.x, self.y 

    def dist(self, p1, p2):  
        #Calculate the great circle distance between two points on the earth (specified in decimal degrees)return distance between two point  
        # convert decimal degrees to radians 
        dlat = radians(p2.x-p1.x)
        dlon = radians(p2.y-p1.y)
        a = sin(dlat/2) * sin(dlat/2) + cos(radians(p1.x))* cos(radians(p2.x)) * sin(dlon/2) * sin(dlon/2)
        c = 2 * atan2(sqrt(a), sqrt(1-a))
        d = 6371 * c
        return d 

    def distanceQuery(self,neighbor_pts):
        dista=[]
        for i in range(len(neighbor_pts)):
          for j in range(i+1,len(neighbor_pts)):
            z=self.dist(neighbor_pts[i],neighbor_pts[j])
            dista.append(z)
        return max(dista)

distanceQuery函数正在使用double for循环。我有什么办法可以把这个去掉吗？我能把这个双循环矢量化吗？由于这是群集代码，因此需要附加一些步骤。我已经读过，numpy数组的工作方式与python list的工作方式不同。附加numpy数组效率低下。在

编辑：

所以可以向量化。但这是代码的另一部分，在这里，我检查了某些条件，然后进行了附加。在

^{pr2}$

现在如果我也向量化相邻点。我要解决附加的问题？所以每个点都将附加到neighbour_points中，然后它将生成一个distanceQuery。这个过程也是迭代的一部分。所以这里也有两个循环，我只想确保在numpy数组中追加不会效率低下

Answer 1

import numpy as np

def dist(p1, p2): 
    # Initially, p1.shape() == (n, 2) and p2.shape() == (m, 2)
    # Now, p1.shape() == (1, n, 2) and p2.shape() == (m, 1, 2)
    p1 = p1[np.newaxis, :, :]
    p2 = p2[:, np.newaxis, :]

    # get all the vectory things
    from numpy import sin, cos, radians, sqrt, arctan2 as atan2 

    # do the same math as before, but use `p[..., 0]` instead of `p.x` etc
    dlat = radians(p2[..., 0] - p1[..., 0])
    dlon = radians(p2[..., 1] - p1[..., 1])
    a = sin(dlat/2) * sin(dlat/2) + cos(p1[..., 0])*cos(p2[..., 0]) * sin(dlon/2) * sin(dlon/2)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = 6371 * c
    return d 

def distanceQuery(neighbor_pts):
    return np.max(dist(neighbor_pts, neighbor_pts))

例如：

^{pr2}$

时间安排：

def dist_slow(p1, p2):
    """your function, adjusted to take an array instead of a `Point`"""
    from math import radians, cos, sqrt, atan2
    # compute the distance for all possible pairs
    dlat = radians(p2[0]-p1[0])
    dlon = radians(p2[1]-p1[1])

    a = sin(dlat/2) * sin(dlat/2) + cos(radians(p1[0]))*cos(radians(p2[0])) * sin(dlon/2) * sin(dlon/2)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = 6371 * c
    return d

def query_iter(p):
    return max(dist_slow(p1, p2) for p1, p2 in itertools.combinations(p, 2))

def query_orig(p):
    dista=[]
    for i in range(len(p)):
      for j in range(i + 1, len(p)):
        z = dist_slow(p[i], p[j])
        dista.append(z)
    return max(dista)

def query_mine(p):
    return dist(p, p).max()

然后：

>>> points = np.random.rand(1000, 2)
>>> timeit query_orig(points)
1 loops, best of 3: 7.94 s per loop
>>> timeit query_iter(points)
1 loops, best of 3: 7.35 s per loop
>>> timeit query_mine(points)
10 loops, best of 3: 150 ms per loop

Answer 2

你可以用numpy ufunc以“向量”的形式做任何事情：

from numpy import radians, sin, cos, sqrt, arctan2
from numpy import random

def max_dist(p1x,p1y,p2x,p2y):
    # give them "orthogonal" shape
    p1x = p1x.reshape(p1x.size,1)
    p1y = p1y.reshape(p1y.size,1)
    p2x = p2x.reshape(1,p2x.size)
    p2y = p2y.reshape(1,p2y.size)

    # compute the distance for all possible pairs
    dlat = radians(p2x-p1x)
    dlon = radians(p2y-p1y)

    a = sin(dlat/2) * sin(dlat/2) + cos(radians(p1x))*cos(radians(p2x)) * sin(dlon/2) * sin(dlon/2)
    c = 2 * arctan2(sqrt(a), sqrt(1-a))
    d = 6371 * c

    return d.max()


if __name__=='__main__':
    # generate random samples
    N = 1000
    p1x,p1y,p2x,p2y = random.rand(4,N)

    print 'max_dist=',max_dist(p1x,p1y,p2x,p2y)