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)]
pts = [ [x[i],y[i]] for i in range(N)]

P0 = pts[np.random.randint(0,N)]
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()

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:

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:

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
  • \$\begingroup\$ Ow, thanks so much. I'm learning a lot with your comments. I didn't know about argmax. And, yes, much faster your code. I think it will be useful when I try with n=5000 and K=1000. \$\endgroup\$ – Sigur Nov 4 '17 at 19:27
  • \$\begingroup\$ I observed that you computed distances using ((p0 - points)**2).sum(axis=1) instead of using scipy.spatial.distance as I did. I suppose that using some predefined function could be faster than doing all math. \$\endgroup\$ – Sigur Nov 4 '17 at 19:41
  • 1
    \$\begingroup\$ @Sigur Btw, for N, K = 5000, 1000, your code takes 64s and mine 0.12s :) \$\endgroup\$ – Graipher Nov 4 '17 at 19:44
  • 1
    \$\begingroup\$ @ImportanceOfBeingErnest As noted in the answer, this code, just as the OP's code,does not produce the set of points which are farthest apart, just a set of points fairly far apart (but it saves a lot of time by doing this). It does this by choosing a random first point and then choosing the point farthest away from the current point. Just like with the traveling salesman problem (where a solution is to always choose to travel to the closest city next), this does not produce the optimal solution but it does in general produce a fairly good solution, without having to try all permutations. \$\endgroup\$ – Graipher May 22 '18 at 14:20
  • 1
    \$\begingroup\$ @ImportanceOfBeingErnest: I added a paragraph near the beginning pointing this out. Hope it is clearer now. \$\endgroup\$ – Graipher May 22 '18 at 19:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.