rankordinal.py

def rankordinal(points):
    n, m = points.shape
    P = np.zeros((n, m), dtype=int)
    R = np.zeros((n, m), dtype=int)
    
    P[:, 0] = np.lexsort(np.flip(points.T, 0))
    R[P[:, 0], 0] = np.linspace(0, n-1, n)

    for i in range(1, m):
        P[:, i] = sorted(P[:, i-1], key=lambda x: points[x, i])
        R[P[:, i], i] = np.linspace(0, n-1, n)
        
    ranks = np.zeros(n, dtype=int)
    for x in P[:, 0]:
        p, r = max((R[x, o], o) for o in range(m))
        for y in P[p+1:, r]:
            if ranks[x] != ranks[y]:
                continue
                
            if np.all(R[x, :] < R[y, :]):
                ranks[y] += 1
    
    return ranks