Farthest point algorithm in Python I am learning Python and I tried to do a program to find K points from a set of N points, farthest as much as possible (I'm not sure if this is the case).

I'd like to know if my code is not redundant, or if the code can be improved.

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import euclidean

def dist_ponto_cj(ponto,lista):
    return [ euclidean(ponto,lista[j]) for j in range(len(lista)) ]

def ponto_mais_longe(lista_ds):
    ds_max = max(lista_ds)
    idx = lista_ds.index(ds_max)
    return pts[idx]

N = 80
K = 40
farthest_pts = [0]*K
print( 'N = %d, K = %d' % (N, K))

#x=[ np.random.randint(1,N) for p in range(N)]
#y=[ np.random.randint(1,N) for p in range(N)]
x=np.random.random_sample((N,))
y=np.random.random_sample((N,))
pts = [ [x[i],y[i]] for i in range(N)]

P0 = pts[np.random.randint(0,N)]
farthest_pts[0]=P0
ds0 = dist_ponto_cj(P0,pts)

ds_tmp = ds0
#print ds_tmp
for i in range(1,K):
    #PML = 
    farthest_pts[i] = ponto_mais_longe(ds_tmp)
    ds_tmp2 = dist_ponto_cj(farthest_pts[i],pts)
    ds_tmp = [ min(ds_tmp[j],ds_tmp2[j]) for j in range(len(ds_tmp))]
    print ('P[%d]: %s' % (i,farthest_pts[i]))

#print farthest_pts

xf=[ farthest_pts[j][0] for j in range(len(farthest_pts))]
yf=[ farthest_pts[j][1] for j in range(len(farthest_pts))]

fig, ax = plt.subplots()
plt.grid(False)
plt.scatter(x,y,c='k',s=4)
plt.scatter(xf,yf,c='r',s=4)
plt.show()
First, a style comment. On the internet, especially in programming, and in particular on this website, English is the lingua franca. So you should avoid mixing other languages and English. This way your code is the most transferable, re-usable and readable.



Second, a comment on the algorithm itself. Your algorithm (and this includes any changes I make to it down below) does not actually find the set of points furthest apart from each other (i.e. the optimal solution). What it does is generate a solution where points tend to be far apart from each other. This is similar to the common strategy for the traveling salesman problem, where you choose to always travel to the closest (unvisited) city next.

This algorithm has the advantage that it does not need to try all combinations, usually quickly leads to a good enough solution, and is very easy to implement. It has the disadvantage that it is on average about a quarter less far apart than the optimal solution and might even return the worst possible solution for some cases.

You might want to look at the other heuristic solutions in that link for different algorithms.



Finally, let's see if the code you have can be improved.

Since you are already using numpy, you should take more advantage of it. It's power lies in using its internal functions, which are implemented and executed in C, independent of the Python interpreter.

As a start, I put your code doing the actual calculations into a function, this way it is re-usable and testable:

@timeit
def op(pts, N, K):
    farthest_pts = [0] * K

    P0 = pts[np.random.randint(0, N)]
    farthest_pts[0] = P0
    ds0 = dist_ponto_cj(P0, pts)

    ds_tmp = ds0
    for i in range(1, K):
        farthest_pts[i] = ponto_mais_longe(ds_tmp)
        ds_tmp2 = dist_ponto_cj(farthest_pts[i], pts)
        ds_tmp = [min(ds_tmp[j], ds_tmp2[j]) for j in range(len(ds_tmp))]
        # print ('P[%d]: %s' % (i, farthest_pts[i]))
    return farthest_pts


if __name__ == "__main__":
    N, K = 80, 40
    pts = np.random.random_sample((N, 2))
    farthest_pts = op(pts, N, K)
    ...


I also used pts = np.random.random_sample((N, 2)) to directly calculate x and y. This means that to plot the points afterwards, you need to use array indexing: plt.scatter(pts[:, 0], pts[:, 1], c='k', s=4).

The @timeit is a decorator*, that prints out the time spent in that particular function whenever it is run. Your code takes about 0.05 seconds on my machine. This is our baseline.

Finally, I put the calling code into a if __name__ == "__main__": guard to allow importing parts of this script from other scripts.



The first thing I would change in your code is the calculation of distances. Here we can use the fact that numpy can operate on the whole array in parallel and just write:

def calc_distances(p0, points):
    return ((p0 - points)**2).sum(axis=1)


Next, here is a way to implement your algorithm using more numpy functions:

@timeit
def graipher(pts, K):
    farthest_pts = np.zeros((K, 2))
    farthest_pts[0] = pts[np.random.randint(len(pts))]
    distances = calc_distances(farthest_pts[0], pts)
    for i in range(1, K):
        farthest_pts[i] = pts[np.argmax(distances)]
        distances = np.minimum(distances, calc_distances(farthest_pts[i], pts))
    return farthest_pts


I also start with pre-assigning an empty array, but using numpy.zeros. I then choose a random point and choose the next furthest point the same way you do.I use numpy.argmax to basically do what your ponto_mais_longe function does, namely return the index of the maximal value. numpy.minimum returns the minimal value for each element in the two given sequences.

This implementation takes about 0.0004 seconds on my machine (so almost 100 times faster).

I don't think this can be sped-up further, because each iteration of the for loop depends on the previous iteration. The only way to speed it up further is to use a different algorithm.



Side note: My timeit decorator looks like this:

from time import clock

def timeit(func):
    def wrapper(*args, **kwargs):
        starting_time = clock()
        result = func(*args, **kwargs)
        ending_time = clock()
        print('Duration: {}'.format(ending_time - starting_time))
        return result
    return wrapper