More efficient way to get diagonal length per coordinate

user_7885426 0
Viewed 162

I have an array of x and y values (coordinates) representing matches and for each of these x,y's I want to know the length of the diagonal it is part of. For example let's take these coordinates

Data description

coords = np.asarray([[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]])
# [[0 0]
#  [0 7]
#  [1 1]
#  [1 6]
#  [2 2]
#  [2 5]
#  [3 3]
#  [3 4]
#  [4 4]]

We can transform it to a matrix but this is too inefficient in my case with enermous tables (for instance, scipy todia() will throw an inefficient warning; see below). Anyway let's make the matrix to make the problem more clear:

[[1 0 0 0 0 0 0 1]
 [0 1 0 0 0 0 1 0]
 [0 0 1 0 0 1 0 0]
 [0 0 0 1 1 0 0 0]
 [0 0 0 0 1 0 0 0]]

Goal
Looking at the table above we see two diagonals (or one diagonal and one antidiagonal). For each position of the diagonal I want to know the length of the diagonal it is part of, so a table like this:

# x, y, diag length
[[0 0 5]
 [1 1 5]
 [2 2 5]
 [3 3 5]
 [4 4 5]
 [3 4 4]
 [2 5 4]
 [1 6 4]
 [0 7 4]]

Inefficient solution
I figured I could represent this data in a sparse scipy matrix while this gives the desired result transforming the sparse matrix to a diagonal coordinate matrix is already inefficient for 100 diagonals let alone for the thousands I have.

from scipy.sparse import dia_matrix, coo_matrix
coords = np.asarray([[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]])

# Create the scipy coord matrix
x = coords[:,0]
y = coords[:,1]
tot_elem = coords.shape[0]*2
data = np.repeat(1, len(x))
co_mat = coo_matrix( (data, (x, y)), shape=(max(x)+1, max(y)+1))

# Get the diagonal matrix
dia_mat = dia_matrix(co_mat).tocoo()
diag_coords = np.column_stack((dia_mat.row, dia_mat.col))

# Get the consecutive values to put them to lengths
difs = np.diff(diag_coords[:, 1])
cuts = [0] + list(np.where(difs != 1)[0] + 1) + [diag_coords.shape[0]]
sizes = np.diff(cuts)
sizes = np.repeat(sizes, sizes)

# Combine with the original coords
dia_sizes = np.column_stack((dia_mat.row, dia_mat.col, sizes))
print(dia_sizes)

*Just realized a coordinate can be part of both a diagonal and antidiagonal, in this case I can report both or only report the length of the longest diagonal - which my solution does not take care of :(

EDIT: More efficient solution

Looking at the todia() code here I noticed they use a smart trick to see if points are on a diagonal, namely x-y should be the same for points on the same diagonal. However, this is not true for the anti-diagonal. So I assume the opposite, x + y does give us poinst on the same antidiagonal. Using this I came up with the code which already is much faster than using scipy.

import numpy as np

coords = np.asarray([[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]])
x = coords[:,0]
y = coords[:,1]

# Get the diagonal (inspired by scripy todia code)
ks1 = y - x

# Unlike scipy, I think we can do the same by summing to get the anti-diagonal
ks2 = y + x

# Sort these to get the groups in the same diagonal
idx = np.argsort(ks1)
anti_idx = np.argsort(ks2)

def get_dia_len(arr,ori):
    sizes = np.diff([0] + list(np.where(np.diff(arr)!= ori)[0] + 1) + [arr.shape[0]])
    size_arr = np.repeat(sizes, sizes)
    return size_arr

# Get the diagonal lengths, i.e. cut at changing values and get the gaps between them
norm_sizes = get_dia_len(x[idx],1)
anti_sizes = get_dia_len(y[anti_idx],-1)

# Gather this in a table
norm = np.column_stack([x[idx], y[idx], norm_sizes])
anti = np.column_stack([x[anti_idx], y[anti_idx], anti_sizes])
dia_coord = np.concatenate((norm, anti))

# We only have a diagonal when we have >1 value
dia_coord = dia_coord[dia_coord[:, -1] > 1]
print(dia_coord)

Have been bending my head around this for a while and curious to see if someone has a smart way to solve this :)

diagonal scipy matrix numpy python
1 Answers

One approach could be to loop through the coordinates and construct 45 degree lines through each point (assuming that's what "diagonal" means), and then remove from the coords list any points that lie on this line -

This function computes points on the 45 degree line of the fixed point and returns only those points which are on the coords list

coords = [[0,0], [0,7], [1,1], [1,6], [2,2], [2,5], [3,3],[3,4], [4,4]]
coords = [tuple(_) for _ in coords]

def get_y(x, fixed_point, allowed_slopes=(1, -1), coords=coords.copy()):
    coords = [tuple(_) for _ in coords]
    x_fixed, y_fixed = fixed_point
    possible_y = [y_fixed + slope*(x - x_fixed) for slope in allowed_slopes]
    possible_coords = [(x, y) for y in possible_y]
    available_coords = list(set(possible_coords) & set(coords))
    return available_coords
print(get_y(1, (0,0)))
#[(1, 1)]
print(get_y(6, (0,0)))
#[] because (6, 6) is not on coords

And then we can loop through coords while removing all points that are on the same line. Using the list.pop ensures that we don't have to unnecessarily compute diagonals multiple times for the same group of points

idx = 0
grouped_points = list()
while coords:
    group = list()
    fixed_point = coords.pop()
    print(f'fixed_point is now {fixed_point}')
    group.append(fixed_point)
    print(f'group is now {group}')
    available_x = set([x for (x, y) in coords])
    print(f'available_x is now {available_x}')
    for x in available_x:
        pt, *_ = get_y(x, fixed_point)
        print(f'pt is now {pt}')
        if pt and pt in coords:
            group.append(pt)
            coords.remove(pt)
        print(f'coords is now {coords}')
        print(f'group is now {group}')
    print(idx, group, sep='\t')
    grouped_points.append(group)
    idx += 1

And then append the lengths to the output to get the desired result

grouped_points = [(*pt, len(group)) for group in grouped_points for pt in group]
print(*grouped_points, sep='\n')
#(4, 4, 5)
#(0, 0, 5)
#(1, 1, 5)
#(2, 2, 5)
#(3, 3, 5)
#(3, 4, 4)
#(0, 7, 4)
#(1, 6, 4)
#(2, 5, 4)

Timing this using timeit shows that this solution is about 10X faster for this set of coords

user_42488420
Related