Removing neighbors in a point cloud I have written a program to optimize a point cloud in dependency of their distances to each other. The code works very well for smaller number of points. For 1700 points it takes ca. 6 minutes. But I have 300000 points in the point cloud. The code is still running after almost 30 hours. 

I am sure there is a pythonic way to optimze the code. But I do not know how can I reduce the calculation time and make the performance better. I have read alot about multithreding, chunk size etc. But I think there is no need for this size of data. I think, I am using the memory in a very bad way. Am I right? Here is the code:

import pickle
import math
import itertools
import gc

def dist(p0, p1):
    return math.sqrt((p0[0] - p1[0])**2 + (p0[1] - p1[1])**2)

def join_pair(points, r):
    for p, q in itertools.combinations(points, 2):
        if dist(p, q) &lt; r:
            points.remove(q)
            return True
    return False


if __name__ == "__main__":
    filename = 'points.txt'
    filename2 = 'points_modified.txt'
    mynumbers = []
    with open(filename) as f:
        for line in f:
            mynumbers.append([float(n) for n in line.strip().split(' ')])
        mynumbers = sorted(mynumbers, key=lambda tup: (tup[1], tup[0]))
    gc.collect()    
    while join_pair(mynumbers, 3.6):
        pass

    with open (filename2,"w")as fp2:
        for j in range(len(mynumbers)):
            fp2.write(' '.join(map(repr, mynumbers[j]))+'\n')


and the input data is like the following:

661234.58 5763766.03 72.63
661254.81 5763765.08 75.04
661244.86 5763764.14 74.99
661234.90 5763763.21 74.94
661225.13 5763762.29 74.89
1. Introduction

It is not a good plan to leave code running for 30 hours. Until you have analyzed the code, and have figured out how its runtime varies according to the size of the input, then you have no idea how long it will take. Maybe it will take 31 hours, but maybe it will take a million hours.

In this situation you need to (i) analyze the algorithm to get an expression for its runtime complexity (that is, how its runtime depends on the size of the input); (ii) measure the performance of the code on small inputs; (iii) use the results of (i) and (ii) to extrapolate to the runtime for the full data set. 

2. Review


There are no docstrings. What are these functions supposed to do?
dist can be simplified using math.hypot:

def dist(p, q):
    "Return the Euclidean distance between points p and q."
    return math.hypot(p[0] - q[0], p[1] - q[1])

The name join_pairs is misleading because it does not join anything. I am sure you could pick a better name.
It would be useful (for testing and performance measurement) to have a function that calls join_pairs in a loop until it returns False. I am going to use this function:

def sparse_subset1(points, r):
    """Return a maximal list of elements of points such that no pairs of
    points in the result have distance less than r.

    """
    points = points[:]
    while join_pairs(points, r):
        pass
    return points


The reason that I am taking a copy of points (rather than modifying the original) is so that I can call this function for testing or performance measurement without destroying the input data. See below for some examples.


3. Performance


Here's a test case of 1,000 random points:

import random
POINTS = [[random.random() for _ in range(2)] for _ in range(1000)]


The code in the post takes about 11 seconds to reduce this to a sparse subset of 286 points at least 0.04 apart:

&gt;&gt;&gt; from timeit import timeit
&gt;&gt;&gt; timeit(lambda:sparse_subset1(POINTS, 0.04), number=1)
10.972541099996306


Why does it take so long? The problem is that join_pairs removes just one point each time it is called. This means that you have to call it in a loop, but each time you call it it starts over comparing the same pairs of points that were already compared in the previous iteration. As a result of this algorithm the runtime is \$O(n^3)\$: that is, the runtime grows like the cube of the length of the input.

We can see this cubic behaviour by measuring the time \$t\$ taken by the function on a range of sizes \$n\$ and finding a cubic equation \$t = an^3\$ of best-fit:

&gt;&gt;&gt; import numpy as np
&gt;&gt;&gt; n = np.arange(100, 1001, 100)
&gt;&gt;&gt; t = [timeit(lambda:sparse_subset1(POINTS[:i], 0.04), number=1) for i in n]
&gt;&gt;&gt; a, = np.linalg.lstsq(n[:, np.newaxis] ** 3, t)[0]


We can plot this on a log-log graph using matplotlib:

&gt;&gt;&gt; import matplotlib.pyplot as plt
&gt;&gt;&gt; plt.xlabel("n")
&gt;&gt;&gt; plt.ylabel("t (seconds)")
&gt;&gt;&gt; plt.loglog(n, t, 'r+', label="data")
&gt;&gt;&gt; plt.loglog(n, a * t**3, label="least squares fit to $t = an^3$")
&gt;&gt;&gt; plt.legend()
&gt;&gt;&gt; plt.grid(alpha=0.25)
&gt;&gt;&gt; plt.show()




So extrapolating from this analysis and measurement, we estimate that if \$n=300000\$ then the time taken will be roughly:

&gt;&gt;&gt; a * 300000**3
339550658.803593


That's about ten years! Looking at the graph, it seems that the slope of the data is not quite as steep as the line of best-fit, so maybe the true exponent is a bit less than 3 and the extrapolated time somewhat less than ten years. But as a rough analysis this tells us that there's no hope for it finish in a few more hours, and so we have to find a better approach.
It would be better to do the work in a single pass over the list of points, like this:

def sparse_subset2(points, r):
    """Return a maximal list of elements of points such that no pairs of
    points in the result have distance less than r.

    """
    result = []
    for p in points:
        if all(dist(p, q) &gt;= r for q in result):
            result.append(p)
    return result


This builds up the modified point cloud in the list result and returns it (instead of modifying the input points). The code considers each point in the input, and adds it to the result so long as it is sufficiently far from all the points that have been previously added to the result.

To check that this is correct, we can compare the output against the original version of the code:

&gt;&gt;&gt; sparse_subset1(POINTS, 0.04) == sparse_subset2(POINTS, 0.04)
True


This revised code takes \$O(nm)\$, where \$n\$ is the number of points in the input and \$m\$ is the number of points in the output. In the worst case \$m = n\$ and the code has to compare every pair of points. On the 1000-points test case, this is hundreds of times faster than the code in the post:

&gt;&gt;&gt; timeit(lambda:sparse_subset2(POINTS, 0.04), number=1)
0.04830290700192563

To further improve the performance, we need a spatial index, that is, a data structure supporting efficient nearest-neighbour queries. For this problem, I'm going to use the Rtree package from the Python Package Index, which implements the R-tree data structure. An R-tree maintains a collection of rectangles (or hyper-rectangles if the search space has more than two dimensions), and can efficiently find the rectangles in the collection that intersect with a query rectangle.

import rtree

def sparse_subset3(points, r):
    """Return a maximal list of elements of points such that no pairs of
    points in the result have distance less than r.

    """
    result = []
    index = rtree.index.Index()
    for i, p in enumerate(points):
        px, py = p
        nearby = index.intersection((px - r, py - r, px + r, py + r))
        if all(dist(p, points[j]) &gt;= r for j in nearby):
            result.append(p)
            index.insert(i, (px, py, px, py))
    return result


When this code appends a point to the result list, it also adds a corresponding rectangle to the R-tree. This rectangle has zero width and height, so it represents a single point. When considering a point p for inclusion in the result set, the code queries a rectangle centred on p with width and height 2 * r. All points within a distance r of p are inside this rectangle, but there might also be some points in the corners of the rectangle that are more distant than r. That's why the code has to further check each point in the query rectangle.

Again, we need to check that this is correct:

&gt;&gt;&gt; sparse_subset1(POINTS, 0.04) == sparse_subset3(POINTS, 0.04)
True


An R-tree takes \$O(\log m)\$ to find an intersecting rectangle in an index of \$m\$ rectangles, so we expect the runtime of sparse_subset3 to be \$O(n\log m)\$, where \$n\$ is the number of points in the input and \$m\$ is the number of points in the output.

Let's make a much bigger test case, with 100,000 points:

POINTS = [[random.random() for _ in range(2)] for _ in range(10**5)]


The performance of sparse_subset2 and sparse_subset3 depends on \$m\$, the number of points in the result. We expect that when r is big, \$m\$ will be small and so there won't be much difference. Here we use r=0.05 and get a subset of 287 points, and the two algorithms have almost identical runtime:

&gt;&gt;&gt; timeit(lambda:sparse_subset2(POINTS, 0.05), number=1)
3.5980678769992664
&gt;&gt;&gt; timeit(lambda:sparse_subset3(POINTS, 0.05), number=1)
3.598255231976509


But if we reduce r to 0.01, getting a subset of 5994 points, then the R-tree shines:

&gt;&gt;&gt; timeit(lambda:sparse_subset2(POINTS, 0.01), number=1)
87.03980600900832
&gt;&gt;&gt; timeit(lambda:sparse_subset3(POINTS, 0.01), number=1)
4.6150341169850435


(I didn't try sparse_subset1 on this test case as it will take many hours.)

The worst case for sparse_subset3 is when r is zero and so all the points in the input end up in the result:

&gt;&gt;&gt; timeit(lambda:sparse_subset3(POINTS, 0), number=1)
10.394442071992671


Extrapolating from this, you should be able to run your 300,000-point case in under a minute, which is a big saving on ten years.