2
\$\begingroup\$

This is a follow up to a previously asked Finding the closest pair of points divide-and-conquer question . The original code changed significantly based on answers from @AJNeufeld therefore I created a new question per codereview guidelines for additional help.

After incorporating suggestions from answers to the original question linked above, the run-time improved by 30% which is significant. Nevertheless, the code is still slow and I think there is room for further improvements. My hunch is that find_min_distance_in_rec is slowing the code down. Below is the current code. Again, the code has been stress tested so I am confident that it's correct, however, it's slow.

#Uses python3
import math
import statistics as stats

# helper functions:
def two_point_distance(p0,p1):
    # returns distance between two (x,y) pairs
    return math.sqrt( ((p0[0]-p1[0])*(p0[0]-p1[0])) + 
                     ((p0[1] - p1[1])*(p0[1] - p1[1])) )

def combine_xy(x_arr,y_arr):
    # combine x_arr and y_arr to combined list of (x,y) tuples 
    return list(zip(x_arr,y_arr))

def find_closest_distance_brute(xy_arr):
    # brute force approach to find closest distance 
    dmin = math.inf
    for i, pnt_i in enumerate(xy_arr[:-1]):      
        dis_storage_min = min( two_point_distance(pnt_i, pnt_j) for pnt_j in xy_arr[i+1:])      
        if dis_storage_min < dmin:
            dmin = dis_storage_min  
    return dmin

def calc_median_x(xy_arr):
    # return median of x values in list of (x,y) points
    return stats.median( val[0] for val in xy_arr )

def filter_set(xy_arr_y_sorted, median, distance):
# filter initial set such than |x-median|<= distance
    return [ val for val in xy_arr_y_sorted if abs(val[0] - median) <= distance ]

def x_sort(xy_arr):
    # sort array according to x value
    return sorted(xy_arr, key=lambda val: val[0])

def y_sort(xy_arr):
    # sort array according to y value
    return sorted(xy_arr, key=lambda val: val[1])


def split_array(arr_x_sorted, arr_y_sorted,median):
    # split array of size n to two arrays of n/2
    # input is the same array twice, one sorted wrt x, the other wrt y
    leq_arr_x_sorted = [ val for val in arr_x_sorted if val[0] < median ]
    geq_arr_x_sorted = [ val for val in arr_x_sorted if val[0] > median ]
    eq_arr_x        = [ val for val in arr_x_sorted if val[0] == median ]

    n = len(eq_arr_x)//2
    leq_arr_x_sorted = leq_arr_x_sorted + eq_arr_x[:n]
    geq_arr_x_sorted = eq_arr_x[n:] + geq_arr_x_sorted

    leq_arr_y_sorted = [ val for val in arr_y_sorted if val[0] < median ]
    geq_arr_y_sorted = [ val for val in arr_y_sorted if val[0] > median ]
    eq_arr_y        = [ val for val in arr_y_sorted if val[0] == median ]

    n = len(eq_arr_y)//2
    leq_arr_y_sorted = leq_arr_y_sorted + eq_arr_y[:n]
    geq_arr_y_sorted = eq_arr_y[n:] + geq_arr_y_sorted

    return leq_arr_x_sorted, leq_arr_y_sorted, geq_arr_x_sorted, geq_arr_y_sorted

def find_min_distance_in_rec(xy_arr_y_sorted,dmin):
    # takes in array sorted in y, and minimum distance of n/2 halves
    # for each point it computes distance to 7 subsequent points
    # output min distance encountered

    dmin_rec = dmin

    if len(xy_arr_y_sorted) == 1:
        return math.inf

    if len(xy_arr_y_sorted) > 7:            
        for i, pnt_i in enumerate(xy_arr_y_sorted[:-7]):
            dis_storage_min = min(two_point_distance(pnt_i, pnt_j) 
                                 for pnt_j in xy_arr_y_sorted[i+1:i+1+7])
            if dis_storage_min < dmin_rec:
                dmin_rec = dis_storage_min

        dis_storage_min = find_closest_distance_brute(xy_arr_y_sorted[-7:])
        if dis_storage_min < dmin_rec:
            dmin_rec = dis_storage_min
    else:
        for k, pnt_k in enumerate(xy_arr_y_sorted[:-1]):      
            dis_storage_min = min( two_point_distance(pnt_k, pnt_l) 
                                for pnt_l in xy_arr_y_sorted[k+1:])      
            if dis_storage_min < dmin_rec:
                dmin_rec = dis_storage_min 

    return dmin_rec             

def find_closest_distance_recur(xy_arr_x_sorted, xy_arr_y_sorted):
    # recursive function to find closest distance between points
    if len(xy_arr_x_sorted) <=3 :
        return find_closest_distance_brute(xy_arr_x_sorted)

    median = calc_median_x(xy_arr_x_sorted)
    leq_arr_x_sorted, leq_arr_y_sorted , grt_arr_x_sorted, grt_arr_y_sorted = split_array(xy_arr_x_sorted, xy_arr_y_sorted, median)

    distance_left = find_closest_distance_recur(leq_arr_x_sorted, leq_arr_y_sorted)
    distance_right = find_closest_distance_recur(grt_arr_x_sorted, grt_arr_y_sorted)
    distance_min = min(distance_left, distance_right)

    filt_out = filter_set(xy_arr_y_sorted, median, distance_min)
    distance_filt = find_min_distance_in_rec(filt_out, distance_min)

    return min(distance_min, distance_filt)

def find_closest_point(x_arr, y_arr):
    # input is x,y points in two arrays, all x's in x_arr, all y's in y_arr
    xy_arr = combine_xy(x_arr,y_arr)
    xy_arr_x_sorted = x_sort(xy_arr)
    xy_arr_y_sored = y_sort(xy_arr)

    min_distance = find_closest_distance_recur(xy_arr_x_sorted, xy_arr_y_sored)

    return min_distance
\$\endgroup\$
2
\$\begingroup\$

The original code changed significantly based on answers

I think that the plural is inappropriate there: the changes seem to be solely based on AJNeufeld's answer. Not a single one of the issues I pointed out has been addressed, not even the bug:

  1. My understanding is that the end result should be four arrays such that leq_arr_x_sorted and leq_arr_y_sorted contain the same points sorted on different projections. But the dmy_x / dmy_y separation of points whose x coordinate is median doesn't seem to guarantee that those points will be sent to the same half in the x-projection and the y-projection. In order to guarantee that an additional pre-condition is needed: most obviously that in arr_x_sorted ties are broken by y; but that pre-condition is not produced by x_sort.

This bug can be demonstrated easily by adding

    print(leq_arr_x_sorted, leq_arr_y_sorted, geq_arr_x_sorted, geq_arr_y_sorted)

before the return of split_array and calling with the simple test case find_closest_point([1, 2, 2, 3], [4, 3, 2, 1]).

\$\endgroup\$
  • \$\begingroup\$ I edited the question to give credit to AJNeufeld, thanks for the suggestion. Regarding your points for improvement, I commented on your answer under original question \$\endgroup\$ – mghah Jul 18 '18 at 20:32
2
\$\begingroup\$

I guess I didn't make this clear in my answer on the original question.

The minimum of the square-root of the square distances is the square-root of the minimum of the square distances. You can save time by calculating the square-root only when needed.

Replace this (and all calls to it):

def two_point_distance(p0,p1):
    # returns distance between two (x,y) pairs
    return math.sqrt( ((p0[0]-p1[0])*(p0[0]-p1[0])) + 
                     ((p0[1] - p1[1])*(p0[1] - p1[1])) )

with this:

def square_distance(p0, p1):
    dx = p0[0] - p1[0]
    dy = p0[1] - p1[1]
    return dx * dx  +  dy * dy

For example, find_closest_distance_brute() becomes:

def find_closest_distance_brute(xy_arr):
    # brute force approach to find closest distance
    dist_sqr_min = math.inf
    for i, pnt_i in enumerate(xy_arr[:-1]):      
        d_sqr_min = min( square_distance(pnt_i, pnt_j) for pnt_j in xy_arr[i+1:])      
        if d_sqr_min < dist_sqr_min:
            dist_sqr_min = d_sqr_min

    return math.sqrt(dist_sqr_min)   # Only calculate square-root of final value

Or more pythonically:

def find_closest_distance_brute(xy_arr):
    # brute force approach to find closest distance
    dist_sqr_min = min((square_distance(pnt_i, pnt_j) for i, pnt_i in enumerate(xy_arr[:-1])
                                                      for pnt_j in xy_arr[i+1:]),
                       default=math.inf)
    return math.sqrt(dist_sqr_min)   # Only calculate square-root of final value

Regarding find_min_distance_in_rect():

def find_min_distance_in_rec(xy_arr_y_sorted,dmin):

    dmin_rec = dmin

    if len(xy_arr_y_sorted) == 1:
        return math.inf

What happens if len(xy_arr_y_sorted) == 0? A better test would be <= 1.

if len(xy_arr_y_sorted) > 7:            
    # ... complicated code involving lots of 7's.
else:
    for k, pnt_k in enumerate(xy_arr_y_sorted[:-1]):      
        dis_storage_min = min( two_point_distance(pnt_k, pnt_l) 
                            for pnt_l in xy_arr_y_sorted[k+1:])      
        if dis_storage_min < dmin_rec:
            dmin_rec = dis_storage_min 

This else: code looks exactly like find_closest_distance_brute() code. In fact, you can replace it with a call to that function.

if len(xy_arr_y_sorted) > 7:            
    # ... complicated code involving lots of 7's.
else:
    dmin_rec = find_closest_distance_brute(xy_arr_y_sorted)      

I don't like # ... complicated code involving lots of 7's.. I'll rework it in a later edit.


See Peter Taylor's answer ... especially the part about median = calc_median_x(xy_arr_x_sorted), for an O(n) to O(1) improvement.


As pointed out a second time by Peter Taylor, your code has a bug; xy_arr_x_sorted and xy_arr_y_sorted are supposed to have the same content, but due to the bug, may not. He has suggested a change to fix the bug. I have a different approach: eliminate xy_arr_x_sorted entirely.

In find_closest_distance_recur(), you have:

    if len(xy_arr_x_sorted) <=3 :
        return find_closest_distance_brute(xy_arr_x_sorted)

Since the arrays are supposed to contain the same content, you could replace this with:

    if len(xy_arr_y_sorted) <= 3:
        return find_closest_distance_brute(xy_arr_y_sorted)

Once you do that, you no longer have any useful references to the [1] members of the contents of xy_arr_x_sorted. Instead, you could replace it with a sorted list of just the x-coordinates.

Once that is done, the partitioning code should be revisited. There is no difference between duplicate median values, so a simple slice operation can divide the x-coordinate list in two with zero copying required.

\$\endgroup\$

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.